The aim of this work is to carry out the numerical simulation of fission product (FP) behavior inside the reactor building under loss of coolant accident (LOCA) using MATLAB. For this purpose, a kinetic model has been developed and implemented in MATLAB to study the behavior of in-containment FPs during postulated LOCA for typical 1000 MW pressurized water reactor (PWR). A continuous release of the FPs from the reactor pressure vessel (RPV) has been implemented with coolant retention. The in-containment FP behavior is influenced by containment atmosphere and containment safety systems. The sensitivity analysis and removal rate of airborne isotopes with the containment spray system have been studied for various spray activation time, spray failure time, droplet size and spray pH value. The droplet size and pH value of the spray system effectively remove the airborne isotopes. The alkaline (sodium thiosulfate, Na2S2O3) spray solution and spray with pH 9.5 have similar scrubbing properties for iodine. However, the removal rate from the containment spray system has been found an approximately inverse square of droplet diameter (1/d2).
- fission products
- sodium thiosulfate
- droplet diameter
The nuclear reactor systems are sufficiently complex that there could be the possibility of an accident followed by the release of fission product (FPs). Such a release could require multiple failures of safety systems and barriers. In the case of a break in the hot/cold leg in a pressurized water reactor (PWR), coolant and energies are first released from the reactor coolant system to the containment through the break. The FP also released along with the coolant. This type of accident usually occurs in the high-pressure cold leg. The worst condition of such an uncontrolled break is the guillotine type of break. In such type of accident, the envelope of primary systems is breached . If such an accident is not controlled by safety systems, then such accidents may transform into the severe accident.
In severe accidents, FP is released during the progression of accidents . Owing to the strong influence of thermal hydraulics on FP release and transportation, FP release and transport mechanism is very complicated and complex. The FP behavior inside the containment is the fundamental of the source term. The source term results are the outputs of level 2 PSA , which are necessary for radiological assessments and consequences. The dominant FPs that constitute in hazardous effects can be categorized as noble gasses (Xe, Kr), volatile (I, Cs, Te), semi-volatile (Ru, Ag, Ba, Sr., Tc, Rh) and nonvolatile (Nb, Zr, Y, Pd, La, Mo, Tc, Nd, Ce) ([4, 5]). The aerosols are 129Te, 127Te, 105Rh, 103Ru, 105Ru, 137Cs, 138Cs, 89Sr, 90Sr and 140Ba. These isotopes release in the particulate form, and going through agglomeration and nucleation process, they form aerosols . However, iodine may transform into volatile species and possess a complex chemistry . The common organic form of iodine is available in chemical forms as CH3I, CsI and HI . The behavior of FP is highly influenced by the in-containment atmosphere, heat loads, containment pressure and steam generation rate. The containment is installed with the spray system and cooling fans to prevent the early over-pressurization due to the heat load. The containment spray system is significant in enhancing the early depletion of radionuclides during early in-vessel release phase from the containment atmosphere. The spray system is automatically activated, as an emergency designed device to prevent containment integrity .
The FP release from a nuclear power plant (NPP) is known as a key factor affecting both the design of safety equipment and safety evaluation, including safety and risk assessment . Experimental research on FP release behavior was conducted by many investigators [11, 12, 13]. Many experiments had played a significant role in understanding the behavior of aerosols, FPs, iodine chemistry, and transportation under accident situations [13, 14, 15, 16, 17]. The Phébus-FP project  was the most impressive program initiated to study the behavior of FP. The main objectives of this project were (1) to minimize the uncertainty in source term evaluation, (2) to study the FPs, structural and control rod material release transportation and deposition from the degraded core through coolant, and (3) and behavior of FP inside the containment building [19, 20]. Meanwhile, several analytical and computational codes were developed. ASTEC is one of the most popular codes used to study the behavior of FPs in severe accident conditions . MELCOR along with MACCS can be used to assess FP release and assessment of radiological consequence . MAAP is the most popular tool to calculate severe accident source term, and its quick calculation is its prime character. Therefore, MAAP code is widely used in the nuclear industry .
Moreover, the numerical simulation of FP activity has been carried out by several researchers.  have developed an analytical model FIPRAP “FP Release Analysis Program” for the numerical simulation of FPs released from the fuel. The FIPRAP code can estimate the volatile FPs released from the nuclear fuel under changing irradiation conditions with the incorporation of all physical phenomena and fulfill the requirements of fuel designing, performance, degradation and source term estimation codes. Lewis et al.  have presented a review of FPs release modeling in support of fuel failure monitoring analysis for the characterization and allocation of defected fuel. A generalized model for FP transport in the fuel-to-sheath gap was given by . Koo et al.  have proposed a model describing pallet oxidation and bubble formation at grain boundaries, their interlinkages, and release into exposed surfaces. Avanov  described a good model for the description release of FPs from the porous ceramic fuel, its leakage from cladding and mixing with the primary coolant. Tucker and white  have proposed an analytical model for the estimation of FPs from ceramic UO2 fuel. In this model, the PF leakage probabilities from the fuel interior through grain are figured out. These probabilities strongly depend on the interconnectivity of pores in the ceramic fuel. Awan et al.  have also carried out the numerical simulation of FP activity in the reactor primary coolant. The proposed developed model is hybrid and analyzes the static and dynamic FP activity in the primary coolant of the reactor .
The goal of this chapter is to carry out the numerical simulation of FP behavior inside the reactor containment building under LOCA using MATLAB. The calculation process of iodine and other FPs is shown in Figure 1. A semi-kinetic model has been developed and implemented in MATLAB to carry out the sensitivity analysis of FPs during postulated LOCA for typical 1000 MW PWR. The kinetic model is presented in section II, which contains the deterministic as well as the kinetic approach. The deterministic computational methodology and computational steps flow chart are described in Section III. Next, the flow chart of model and implementation of model in MATLAB are described in Section IV. The examples and outcomes of the simulation results are presented in Section V. Finally, Section VI is the conclusion.
2. In-containment fission product release model
Figure 1 shows the process of release of FPs from fuel to cladding, cladding to coolant and then to the containment. In this work, a 1000-MW pressurized water reactor (PWR) has been considered with the design specification as shown in Table 1. The PWR system along with the containment system is shown in Figure 2. We have developed a real-time kinetic model to simulate the FP behavior inside the containment. The analytical model is a set of coupled ordinary differential equations (ODEs). The FP activity inside the reactor containment building and on the surfaces and walls of the containment is governed by the following sets of ODEs [8, 32, 33].
where i indicates the isotope, whereas V and S indicate the volumetric and surface activities of ith isotope. The puff release of FP is mv (t) = fx × ff × fp × fc × Ac/V g.m−3. The values of various parameters used in these simulations are listed in Table 2.
|Average fuel enrichment wt%||2.4%|
|Specific power (MWth/kg U)||33.3|
|Power density (MWth/m3)||66.6|
|System pressure (MPa)||15.166|
|System pressure (MPa)||14.96|
|Coolant flow (kg/s)||17387.7|
|Core height (m)||12.41|
|Core active region (m)||3.65|
|Core diameter (m2)||3.81|
|Control rod assemblies||69|
|Fuel rod outer diameter (cm)||1.092|
|Rod pitch (cm)||1.443|
|Fuel assembly matrix||15 × 15|
|Coolant inlet temperature (k)||564.81|
|Coolant outlet temperature (k)||592.98|
|Control rod material||Ag (80%)-In (15%)-Cd (5%)|
|Containment free volume||V (m3)||57,600|
|Containment free surface||S (m2)||~34,374|
|Leakage rate||Lr (m3/s)||14.15|
|Core damage fraction||fc (%)||35%|
|Fuel release fraction||ff (%)||9.0 × 10−1|
|Water release fraction||fp (%)||3.00 × 10−1–1.0 × 102|
|Recirculation rate||Rres (m3/s)||1–5|
|Recirculation filtration efficiency||ηres (%)||10–90%|
|Exhaust filter efficiency||ηex (%)||90–98|
|Fraction of immediately released radioisotopes||fx (%)||2.0 × 10−1|
|Mixing rate||wx (s−1)||0.1–1.0|
|Spray flow rate||F (m3/s)||0.1–1.0|
|Droplet size||d (micron)||100–1000|
|Deposition velocity (ud)||(m/s)||5.5 × 10−4|
|Resuspension rate||s−1||≤2.3 × 10−6|
|Partition coefficient||H||200 (for pH 5.0), 5000 (for pH 9.5) and 10,000 (Na2S2O3)|
|Molar weights of solvent||Ml (g/mole)||18.01528|
|Temperature||T (K)||80 + 273.15|
|Molecular volume of I2||υ (cm3/g)||71.5|
|Spray flow rate||F (m3/sec)||0.35|
|Degree of solvent||x||2.6 for H2O|
2.1. Kinetic source of fission product
The (1−fx) exp.(−wxt) is the airborne FP activity released along with the coolant with mixing rate wx. Where K is the normalization constant and expressed as follows. The overall radioactive mass inventory, including kinetic and static parts, is depicted in Eq. (6).
2.2. Fission product removal with spray
The removal of iodine and aerosols from the containment with the spray system can be expressed as depicted in Eqs. (7) and (8), where mri and mra are the removal rates of iodine and aerosols, respectively.
3. Deterministic computational methodology
Several steps are involved in the simulation of FP behavior inside the reactor building starting from the generation of FP in fuel along with the fuel burn-up. Leakage of FP into the coolant and then from the coolant to containment along with the leakage of coolant. The computational steps are listed in Figure 3. A two-stage methodology has been adopted: (1) evaluation of activity in the core just before the accident and (2) kinetic quantification of airborne activity under confined conditions. The core activity has been evaluated at for one complete fuel cycle to get maximum core activity. The behavior of airborne FP activity has been quantified for loss of the coolant accident (LOCA) under NUREG-1465  and regulatory guide 1.183  assumptions. The developed model uses subroutine functions containing coupled ODEs and Runge–Kutta (RK) method. The ODEs (Eqs. (1)–(12)) are implemented in MATLAB. The system of ODEs (Eqs. (1), (3), (7), (8)) is coupled and solved numerically using the Runge–Kutta (RK) method in this program.
The RK numerical provides efficient time-domain solution, yielding static as well as dynamic values of FPAs corresponding to about 84 different dominant FPs. The computational cycle starts with the initialization of the variables with t = 0. In the time loop, the values of FPAs inside the containment building are calculated using RK scheme for each next time step. The program allows performing these calculations for spray system operation.
4. Implementation in MATLAB
The above equations can be implemented in MATLAB. The flow chart of the MATLAB program is shown in Figure 4. In the first step, the physical constant and parameters are defined, and the time array and droplet size are determined by the user.
% MATLAB Program for In-containment Fission product program by Khurram Mehboob
% Date : 08-07-2017
clear; clc; clear all;
Global Hi Lr V S vd dec r Rr neu EI h Klcm Kgcm d Ea fr H y00 Q y t I Ac D Core_I
Cont_A QQ f x fc B wx YY Sorc wx1
tn = input('Enter end time = tn = '); h = input('Enter stepsize = h = '): t = (0:h:tn); % time array
for d1=100: 100: 1000; % particle diameter (microns)
d = d1*1e-4; % particle diameter (cm)
k=d1/100; % Droplet control Factors for printing
fx = 0.20; % activity immediately available in the containment air
fc = 0.35; % core damage fraction.
H =10000; % partition coefficient for iodine
Rr = 4.719; % Recirculation flow rate
Lr = 14.15; % leakage rate
wx = 0.01; % mixing rate
In the second step, the fixed variables are loaded from an input text file. The input text file contains the output data from the ORIGEN2.2 code that contains data for 84 different FPs.
V = input2(1,1); % free volume of the containment
S = input2(2,1); % free surface of the containment
No_iso = input2(5,1); % no of isotopes
Hi = input2(6,1); % height of spray system 40.0 m
fr = input2(10,1); % Recirculation filtration efficiency
wx1 = wx/10.0; % mixing rate of fission products in coolant
K = (wx*wx1)/(wx-wx1); % normalization constant
for i = 1:No_iso; % loop to read isotope data
Ac = input(i+10,4); % fission product activity in the core.
for j = 6:7
Ac = Ac*input(i+10,j); % Ac*ff*fw
dec = input2(i+10,2); % decay constant
vd = input2(i+10,08); % deposition velocity
r = input2(i+10,09); % resuspension rate
neu = input2(i+10,10); % heap filter efficiency
In the third step, the fixed variables for Eqs. (9)–(12) are solved for droplet by using the data listed in Table 3. The values of parameters (T, x, Ml, ul Sc, Re) are the constant variable for the static containment atmosphere.
T = 80+273.15; % Temperature (K)
x = 2.6; % H2O
Ml = 18.01528 ; % molar weight of solvent (g/mole )
v = 71.5; % molar volume of diffusing substance (I2)
Sc = 1.742; % Schmidt number.
Re = 1.95; % Reynold number.
c1 = T-281.615; c2 = (T-281.615)^2; c3 = 8078.4+ c2; c4 = c3^0.5; C = 2.1482*(c1+c4)-120;
ul = (100/C)*2.41908832931* 14.8816390000001; %conversion of Centipoise to g/cm.s
Next, the fixed variables are used to determine the dynamics variables (DL, Dmix, KG, KL). The Eqs. (9)–(12) are solved using the static data calculated above. Also, the kinetic release of FP is Q(t) that is determined as the function of time.
DI2= (((7.4)*(10^-8))*((x*Ml)^0.5)*T)/(ul*(v^0.6)) ; %diffusion coefficient if I2 (cm2/s)
DIMix = 0.00035*0.258064; %diffusion coeffi of I2 in steam ( cm2/s)
Kgcm = (DIMix/d)*(2.0+0.6*(Re^(1/2))*(Sc^(1/3))); % liquid phase trans coefficient ( cm/s )
Klcm = ((2/3)*(3.1416^2)*(DI2))/(d); % gas phase trans coefficient (cm/s)
vt = ((Re*2387.4)/(d1))* 0.508; % terminal velocity (cm/s)
te = (40*100)/(2*vt); % exposure time (s)
A = Kg*(te); BB = d*(H+(Kg/Kl));
EI = (1 - exp(-6*(A/BB))); % FP removal rate (s-1)
Q = ((1-fx)*B*Sorc*(exp(-wx1*t)- exp(-wx*t))); % kinetic source
Next, the initial conditions for airborne and surface activities are implemented, for example, mv (t) = fx × ff × fp × fc × Ac/V g.m−3 and ms(t) = 0.0. The kinetic and static parameter values are implemented in coupled equation written as subfunction diffeq, and Rang–Kutta fourth-order method is implemented by calling odeRK4 subroutine function.
y00 = fx*fc*Ac/V; y0 = [y00, 0, 0];
[t,y] = odeRK4(@diffeq, tn, h, y0);
YY(:,i)= y(:,1) ; % YY=y(length(y),:);
I(i)= input2(i+10,1); % Atoms number
D(i)= input2(i+10,2); % decay constant
Core_I(i)= input2(i+10,4); % Activity in the core
Cont_A(i)= y00; % immediate released activity in containment
function dy = diffeq(t,y)
global Lr V S vd dec r F EI H Rr fr fx Ac fc wx B wx1 Q Sorc
if t <=700
F = 0.0; % input2(8,1); % spray flow rate (0.1-2.0 m3/s)(950 m3/h=0.264m3/s)
Q = ((1-fx)*B*Sorc*(exp(-wx1*t)- exp(-wx*t)));
dy(2)= vd*y(1)- r*y(2);
dy(3)= Q - ((EI*H*F)/V)*y(3);
The Range–Kutta fourth-order method is implemented by calling the odeRK4 subroutine function. The function is capable of solving N number of coupled ODEs at separate time steps. The implementation if RK4 method is shown below.
% implementation of Range kutta 4thorder method.
function [t,y] = odeRK4(diffeq,tn ,h, y0)
t = (0:h:tn); % time vector with spacing h
nt = length(t); % no. of elements in vector t
neq = length(y0);
y = zeros(nt, neq); % prelocation of y for speed
y(1,:) = y0(:)';
h2=h/2; h3=h/3; h6=h/6;
k1 = zeros(neq,1); k2= k1; k3 = k1; k4= k1; ytemp = k1;
told = t(j-1); yold = y(j-1,:)';
k1= feval(diffeq,told, yold);
for n= 1:neq
ytemp(n) = yold(n) + h2*k1(n);
k2= feval(diffeq,told+ h2, ytemp);
for n= 1:neq
ytemp(n) = yold(n) + h2*k2(n);
k3= feval(diffeq,told +h2, ytemp);
for n= 1:neq
ytemp(n) = yold(n) + h*k3(n);
k4= feval(diffeq,told+h, ytemp);
for n= 1:neq
5. Some results from numerical simulation
5.1. Volumetric fission product inventory
The core inventory for typical 1000 MW PWR has been evaluated by ORIGEN 2.2 code which is used by our model as a subroutine. A 35% core damage has been considered and 20% (fx) as the puff release. While the rest of radioactive mass release along with coolant with mixing rate wx = 0.01 s−1. The FP release inside the containment building as a function of time is depicted in Figure 5. The volumetric radioactive mass found to increase during first 300 s then starts decreasing with a constant rate. The cesium found to be dominant with 100 times higher than the other radioactive released masses. The Krypton gas is found to be 15% higher in magnitude with Xenon gas. However, the other isotopes show the similar behavior but with less in magnitude.
5.2. Puff release (fx) effect on in-containment FPA
The LOCA is due to the uncontrolled leakage from coolant piping (hot leg or cold leg). The coolant burst release generates the immediate escape of radioiodine into containment with the rapture. The volumetric activity of 137Xe inside the containment has been simulated for various values of instantaneous burst release (i.e., fx = 10–70%) of total activity inside the core. The simulation results are depicted in Figure 6. The results indicated that with the higher percentage of instantaneous release (fx = 50–70%), the activity in containment slightly increased and then decreased linearly.
However, a less fraction of burst release (fx = 10–30% of total activity), the activity inside the containment first increases and then starts decreasing after approaching to the maximum value. As the value of fx decreases, the peak shifts toward higher timescale. The peak becomes more prominent with small values of fx. This happens due to competition between fx term in the initial condition and (1-fx)(exp -wx) term in source term (Eq. (4)). The behavior of 137Xe for various values of instantaneous release (fx) with mixing rate wx = 0.01/s explains the clearer picture of airborne Xenon (Figure 6).
5.3. Mixing rate (wx) effect on in-containment FPA
The delayed core released fraction (1-fx)(exp -wxt) for various values of mixing rate wx which contributes to airborne volumetric activity through coolant has been simulated. The volumetric activity of 131I inside the containment for various values of mixing rate wx from 0.005 to 1.0 s−1 has been numerically simulated. It has been observed that a slight change in the value of wx results in astonishing variation in airborne activity. The results are shown in Figure 7. For wx = 1.0 s−1, highest magnitude has been observed with a shift of peak toward lower timescale value. However, the 131I has been observed to reduce almost linearly with mixing rate 1.0 s−1 (Figure 7). The in-containment volumetric concentration of 131I increases and then starts decreasing gradually for mixing rate value from wx = 1.0 to 0.03 s−1; however, it remains constant for mixing rate 0.005 s−1. This is because of trivial mixing of iodine in the coolant. A higher magnitude of 131I activity has been observed with higher mixing rate.
5.4. Fission product activity with spray system
The primary purpose of the spray system is to mitigate the FP exposure to the environment and to maintain the containment integrity. In this work, we have studied the effect of the spray system in mitigating the radioactive masses (gaseous and particles) released during in-vessel release phase. During loss of coolant accident, the temperature and pressure inside the containment start raising. It reached to 80° pressure and reached to 7.533 psi within few minutes. The simulation has been carried out by assuming the containment temperature at 80°C and pressure at 7.533 psi and the spray with pH 5.0 and 9.5 and the alkaline spray. The spray system is found to have minimum effect on noble gasses and reduces iodine and other radioactive particles effectively. The spray system is started at 100, 500, 1000, and 1500 s after the release time. The effect of spray system activation time on noble gasses and iodine is shown in Figure 8.
The iodine concentration first increases exponentially in containment and immediately starts reducing with the activation of the spray system. This is because of computation between the continuous source of radioactive iodine (Pi) coming from reactor pressure vessel (RPV) and removal of iodine with the containment spray system (Figure 8). The volumetric iodine starts reducing exponents after the activation of the spray system. The droplet size has been assumed 800 microns for the simulation. However, the noble gasses are observed to be unaffected with a spray system. The response of airborne iodine to the failure of the spray system is depicted in Figure 9.
Figure 9 indicates that the premature failure of the system (t = 100 s) does not affect the airborne iodine concentration. The slight decrease in airborne iodine concentration is seen if the spray system is failed or malfunctions at 500 s. However, the spiracles are operated for a longer period, for example, 1000–1500 s during the in-vessel phase. The airborne concentration of iodine is reduced significantly. The failure of the system at 1000 and 15,000 s caused the regaining of airborne iodine (Figure 9).
5.5. Droplet diameter and pH effect on in-containment FPA
The droplet collection efficiency of spray also depends on the containment atmospheric temperature, pressure and spiracle pH value. The affect of spray water pH value and an alkaline spray has been simulated. The results are depicted in Figure 10. The results showed that the higher pH spray solution (pH 9.5) and alkaline solution (Na2S2O3) have similar removal characteristics for airborne iodine. The iodine removal rates of Boric (pH 5.0), NaOH (pH 9.5) and alkaline solution with different droplet sizes are shown in Figure 11.
The removal rate has been seen to decrease exponentially with the increase in droplet size for an alkaline solution (Figure 11). However, for a spray solution with pH value 5.0, the removal rate decreases in a linear manner. The in-containment volumetric mass under the atmospheric conditions with a temperature of 80°C and 7.533 Psi. Assuming 35% core damage and 20% burst release. The rest of mass is assumed to release along with the coolant with mixing rate wx = 0.01 s−1. The mass concentration of tellurium is simulated for droplet sizes (100–1000 microns). The containment spray system is activated with the initiation of an accident with a constant flow rate of 0.2 m3/s.
Simulation results showed that the droplet size is quite effective to reduce the airborne FPs. It has been observed that the concentration of airborne tellurium decreases with a decrease in droplet size (Figure 12). The peak concentration in Tellurium mass reaches to a maximum concentration at a longer time with the higher droplet diameter. The magnitude of maximum concentration has been found the approximately inverse square of droplet diameter (1/d2). The containment spray system removal rate for iodine versus droplet diameter is depicted in Figure 11. The maximum removal rate has been found 452 s−1 with alkaline solution spray with a droplet size of 100 micrometers. The removal rate is found to decrease exponentially as the droplet diameter increases. 44.7 s−1 removal rate has been seen for 1000-micron diameter droplet size for pH 9.5 and alkaline spray solution with spray flow rate 0.35 m3/s. The gas phase and liquid phase coefficients play a vital role in absorption efficiency. Both gas- and liquid-phase mass transfer coefficients (KG and KL) decrease drastically with an increase in droplet size (Figure 13). However, the gas- and liquid-phase mass transfer coefficients (KG and KL) are also related to the inverse square of droplet diameter (1/d2).
This chapter has presented the numerical simulation of FP activity inside the reactor containment building under LOCA. The numerical simulation of in-containment FPs against the mixing rate, puff release, droplet diameter, spray pH value and spray performance has been simulated. The results indicate that the mixing rate of FPs in coolant significantly affects airborne FP activity inside the containment. The higher pH spray solution (9.5) and spray with sodium thiosulfate (Na2S2O3) have observed similar scrubbing properties. The droplet size is significantly important for removal of FP. There is a higher tendency of FP to interact with airborne particles (Figure 11) with, due to their higher values of liquid- and gas-phase mass transfer coefficients (KL and KG) (Figure 13). Therefore, the acceptance criteria of droplet size have been suggested between 600 and 800 microns, with pH value higher than 7.0 which delivers higher removal rate. The earlier the containment spray system has operated, the airborne concentration will be minimum (Figure 8). However, the delay operation caused the higher airborne concentration of radioactive mass. It has also been observed if the containment spray system is failed during in-vessel phase a regain will be caused in radioactive mass (Figure 9). Based on our work, we are suggesting 600–8500 microns mean droplet diameter containment spray system should be used to get maximum radiation hazard safety. Moreover, from our results, we can conclude that the spray system should be operated within 500 s after the accident and should be operated more than 3000 s (whole in vessel phase). The uncertainties in simulated results depend on generally available data in the literature.
This article was funded by the deanship of scientific research (DSR) at King Abdul Aziz University, Jeddah. The author, therefore, acknowledges with thanks DSR for technical and financial support.