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.

function PWR\_Fission\_Product % 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) %=======Control Variables==================== 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.

load 'input.txt'

%=======Fixed variables==============

Figure 4. Flow diagram of computer program.

No\_iso = input2(5,1); % no of isotopes

V = input2(1,1); % free volume of the containment S = input2(2,1); % free surface of the containment

Numerical Simulation of Fission Product Behavior Inside the Reactor Containment Building Using MATLAB

http://dx.doi.org/10.5772/intechopen.70706

51

Hi = input2(6,1); % height of spray system 40.0 m fr = input2(10,1); % Recirculation filtration efficiency Numerical Simulation of Fission Product Behavior Inside the Reactor Containment Building Using MATLAB http://dx.doi.org/10.5772/intechopen.70706 51

Figure 4. Flow diagram of computer program.

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 [8] and regulatory guide 1.183 [32] 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

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

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

% MATLAB Program for In-containment Fission product program by Khurram Mehboob

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

tn = input('Enter end time = tn = '); h = input('Enter stepsize = h = '): t = (0:h:tn); % time array

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.

program allows performing these calculations for spray system operation.

defined, and the time array and droplet size are determined by the user.

%================================================%

%================================================

for d1=100: 100: 1000; % particle diameter (microns) %=======Control Variables==================== d = d1\*1e-4; % particle diameter (cm)

fc = 0.35; % core damage fraction.

Lr = 14.15; % leakage rate wx = 0.01; % mixing rate

H =10000; % partition coefficient for iodine Rr = 4.719; % Recirculation flow rate

k=d1/100; % Droplet control Factors for printing

fx = 0.20; % activity immediately available in the containment air

4. Implementation in MATLAB

50 Numerical Simulations in Engineering and Science

Cont\_A QQ f x fc B wx YY Sorc wx1

function PWR\_Fission\_Product

% Date : 08-07-2017

clear; clc; clear all;

program.



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.

%Initial Condition

end

if t <=700

else F= 0.35;

end

end

y00 = fx\*fc\*Ac/V; y0 = [y00, 0, 0];

[t,y] = odeRK4(@diffeq, tn, h, y0);

function dy = diffeq(t,y)

dy(2)= vd\*y(1)- r\*y(2); dy(3)= Q - ((EI\*H\*F)/V)\*y(3);

neq = length(y0);

h2=h/2; h3=h/3; h6=h/6;

y(1,:) = y0(:)';

for j=2:nt

%================================================

%================================================

%================================================

Cont\_A(i)= y00; % immediate released activity in containment

condition for containment spray is controlled in subroutine function here.

F = 0.0; % input2(8,1); % spray flow rate (0.1-2.0 m3/s)(950 m3/h=0.264m3/s)

global Lr V S vd dec r F EI H Rr fr fx Ac fc wx B wx1 Q Sorc

The subroutine function containing Eqs. (1), (3), (4) and (7) or (8) is depicted below. The

Numerical Simulation of Fission Product Behavior Inside the Reactor Containment Building Using MATLAB

http://dx.doi.org/10.5772/intechopen.70706

53

dy(1)= -dec\*y(1)-(Lr/V)\*y(1)-vd\*(S/V)\*y(1)+r\*(S/V)\*y(2)-fr\*(Rr/V)\*y(1)-((F\*H\*EI)/V)\*y(3)+((1-fx)

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.

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

Q = ((1-fx)\*B\*Sorc\*(exp(-wx1\*t)- exp(-wx\*t)));

The implementation if RK4 method is shown below.

% implementation of Range kutta 4thorder method.

t = (0:h:tn); % time vector with spacing h nt = length(t); % no. of elements in vector t

y = zeros(nt, neq); % prelocation of y for speed

k1 = zeros(neq,1); k2= k1; k3 = k1; k4= k1; ytemp = k1;

\*B\*Ac\*(fc/V)\*((exp(-wx1\*t))-exp(-wx\*t)));

function [t,y] = odeRK4(diffeq,tn ,h, y0)


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.


Next, the initial conditions for airborne and surface activities are implemented, for example, mv (t) = fx � ff � fp � fc � Ac/V g.m�<sup>3</sup> and ms(t) = 0.0. The kinetic and static parameter values are implemented in coupled equation written as subfunction diffeq, and Rang–Kutta fourthorder method is implemented by calling odeRK4 subroutine function.

```
%Initial Condition
```
wx1 = wx/10.0; % mixing rate of fission products in coolant

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

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

Next, the initial conditions for airborne and surface activities are implemented, for example, mv (t) = fx � ff � fp � fc � Ac/V g.m�<sup>3</sup> and ms(t) = 0.0. The kinetic and static parameter values are implemented in coupled equation written as subfunction diffeq, and Rang–Kutta fourth-

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)

EI = (1 - exp(-6\*(A/BB))); % FP removal rate (s-1)

Q = ((1-fx)\*B\*Sorc\*(exp(-wx1\*t)- exp(-wx\*t))); % kinetic source

order method is implemented by calling odeRK4 subroutine function.

Ac = input(i+10,4); % fission product activity in the core.

K = (wx\*wx1)/(wx-wx1); % normalization constant

for i = 1:No\_iso; % loop to read isotope data

%===================================%

52 Numerical Simulations in Engineering and Science

Ac = Ac\*input(i+10,j); % Ac\*ff\*fw

static containment atmosphere.

x = 2.6; % H2O

A = Kg\*(te); BB = d\*(H+(Kg/Kl));

Sorc= Ac\*fc/V;

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

T = 80+273.15; % Temperature (K)

Sc = 1.742; % Schmidt number. Re = 1.95; % Reynold number.

FP is Q(t) that is determined as the function of time.

Ml = 18.01528 ; % molar weight of solvent (g/mole ) v = 71.5; % molar volume of diffusing substance (I2)

for j = 6:7

end

%================================================ 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 end

The subroutine function containing Eqs. (1), (3), (4) and (7) or (8) is depicted below. The condition for containment spray is controlled in subroutine function here.

```
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)
 else
 F= 0.35;
end
Q = ((1-fx)*B*Sorc*(exp(-wx1*t)- exp(-wx*t)));
dy(1)= -dec*y(1)-(Lr/V)*y(1)-vd*(S/V)*y(1)+r*(S/V)*y(2)-fr*(Rr/V)*y(1)-((F*H*EI)/V)*y(3)+((1-fx)
*B*Ac*(fc/V)*((exp(-wx1*t))-exp(-wx*t)));
dy(2)= vd*y(1)- r*y(2);
dy(3)= Q - ((EI*H*F)/V)*y(3);
end
```
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;
for j=2:nt
```

```
told = t(j-1); yold = y(j-1,:)';
  k1= feval(diffeq,told, yold);
  for n= 1:neq
     ytemp(n) = yold(n) + h2*k1(n);
  end
  k2= feval(diffeq,told+ h2, ytemp);
  for n= 1:neq
     ytemp(n) = yold(n) + h2*k2(n);
  end
  k3= feval(diffeq,told +h2, ytemp);
  for n= 1:neq
     ytemp(n) = yold(n) + h*k3(n);
  end
  k4= feval(diffeq,told+h, ytemp);
  for n= 1:neq
     y(j,n)= yold(n)+h6*(k1(n)+k4(n))+h3*(k2(n)+k3(n));
  end
 end
end
```