Abstract
Coupled first-order IVPs are frequently used in many parts of engineering and sciences. We present a “solver” including three computer programs which were joint with the MATLAB software to solve and plot solutions of the first-order coupled stiff or nonstiff IVPs. Some applications related to IVPs are given here using our MATLAB-linked solver. Muon catalyzed fusion in a D-T mixture is considered as a first dynamical example of the coupled IVPs. Then, we have focused on the fuel depletion in a suggested PWR including poisons burnups (xenon-135 and samarium-149), plutonium isotopes production, and uranium depletion.
1. Introduction
Coupled first-order IVPs are frequently used in many parts of engineering and sciences [1–3], and we presented a package seems to be useful for researchers to solve IVPs [4]. It is possible to describe many dynamical problems using IVPs; MATLAB is the best software for engineers and applied scientists to solve the problems numerically, specially solving IVPs. In our early studies, we have utilized a numerical “MATLAB-linked solver” to calculate stiff or nonstiff first-order coupled IVPs using MATLAB software [4], and reader can find these programs in the appendix. The well-known numerical methods such as Runge-Kutta, Rosenbrock, Classical method, Taylor series, Adams-Bashforth are used to solve IVPs using our MATLAB-linked solver [4, 5].
The main aim of the present research is to give a MATLAB-linked solver to solve first-order coupled differential equation which is used in many subjects of the nuclear engineering. Therefore, in the present study, some dynamical problems (which mathematically are coupled first-order IVPs) are studied as examples of the present solver ability. First we explain Muon Catalyzed Fusion (CF) and find the fusion cycling rate. Then, we focus on the poisons, including xenon135 and samarium-149, burnups in a suggested 1000 MWe PWR as well as its plutonium isotopes build up. Their solutions are given using our “MATLAB-linked solver.”
Basically, consider a first-order coupled IVPs such as where initial values of the dependent variables are: Here, is used for the independent variable and may refer to time in a dynamical problems, and stands for dependent variables.
Coupled IVPs with constant coefficients. First, we consider the IVPs with constant coefficients, or in other words constant , and we illustrate our package procedures to solve and plot the calculated dependent variables. Three programs were written and connected to the MATLAB, software where these programs will be run in the MATLAB's editorial page, by running DEPLET.m M-file. Some questionnaires should be answered by the user such as the following:(i)Entering the number of differential equations (unknowns).(ii)Inserting initial values of (dependent variables).(iii)Inserting start and end points of the computations, or in another words independent variables interval.(iv)The type of coupled differential equation should be specified. The answer includes “Stiff” or “Nonstiff” cases.(v)The next question is the method in which user wants for executing. Answers includes “ode45 method", “ode23 method", “ode113 method" for the nonstiff case, and also “ode15s method", “ode23s method", “ode23t method", “ode23tb method" for stiff case so that for more information about these MATLAB commands, refer to the MATLAB help [5]. By clicking on each solvers, a short review on the specified numerical method will be given. Finally, user should insert s and execution would be begun. After execution, dependent variables ()s will be computed according to the desired numerical method, and user can plot s in the computational interval. For plotting as an example, user should write “1" and then a new window will be opened and will be plotted versus the computed interval.
Coupled IVPs with variable coefficients. The so-called package (called DEPLET m-file) can be extended to the IVPs with variable coefficients or in another words . In this case, the computational interval should be divided into many step-size intervals so that variations of s are ignored in each step-size (each step-size is equal to , where is a desired grid interval and therefore , where ), and therefore we can assume average of s in each step-size: In this case again we return into the coupled differential equations with constant coefficients which can be solved by the DEPLET.m. Calculated , in each interval, will be used as an initial conditions for the next step, and therefore by combining given solutions the full-solutions would be obtained.
2. Some Dynamical Problems Related to Coupled IVPs
2.1. Muon Catalyzed Fusion System
The basic process of the muon catalyzed fusion in a - mixture as depicted in the upper part of Figure 1, can be summarized as follows [6]. After high-energy injection and then stopping and decreasing its energy in a - mixture, either () or a () atom is formed, with a probability more or less, proportional to the relative concentration of () and (). Because of the difference between () and () in the binding energies of their atomic states, in () undergoes a transfer reaction to tritium yielding () during a collision with the surrounding tritium in either - or molecules. Thus the formed () reacts with , , or to form a muonic-molecule at a rate of followed by a fusion reaction occurring from a molecular state of the (). The fusion takes place and a 14-Mev neutron and a 3.6-Mev -particle are emitted. After the fusion reaction inside the () molecule, most of the are liberated to participate in a second CF cycle. There is however some small fraction of the which are captured by the recoiling positively charged . The probability of forming an ()+ ion is called the initial sticking probability . Once the () is formed, the can be stripped from the () ion where it is stuck and liberated again. This process is called regeneration, with a corresponding fraction . Thus, in the form of either a nonstuck or one regenerated from () can participate in a second CF cycle, leading to an effective sticking parameter .
Now, consider a homogeneous media in which the CF is carried out [7–9]. The ion density of the media , in another words tritium and deuterium concentrations, atomic and molecular formation rates (, and ), and fusion decay rates () are known as the constants due to a fixed temperature of the media and therefore are assumed to be independent of time. Therefore, according to the physical model and also Figure 1, the first-order linear coupled dynamical equations for the Muon Catalyzed Fusion system () are given by where and are the relative concentration of deuterium and tritium, respectively. The muonic-atoms formation rates are given by (), and the muonic-molecule formation rate of is given by . The fusion rate is shown by . The possible leakage rate of muonic-atoms is proportional to and , and also is the muon decay constant. In (2.1), is the dimensionless ion density of the media and is proportional to liquid hydrogen density (close to 0.07). As said before, is the muon effective sticking coefficient. Table 1 gives values of the constants for solving (2.1).
We have solved these coupled dynamics equations in time range of [0 to ], the muon life-time, using our MATLAB-linked solver with the following initial conditions: According to the coupled dynamical equation (2.1) and also Figure 1, the calculated neutrons, for each one inserted muon, corresponds to the muon cycling rate of the explained cold fusion (), or which are proportional to the number of fusion (i.e., each neutron corresponds to one fusion: ).
At the end of running, neutron concentration is plotted in our calculated time interval and is given in Figure 2. According to Figure 2, is the muon cycling rate of the mentioned Muon Catalyzed Fusion system. Each fusion gives 17.6Me 5 energy so that total obtained energy is about . For producing one muon, due to an accelerator, we must expense about energy so that the above cold fusion is not commercial. Increasing temperature as well as more ion density concentration together with decreasing muon sticking coefficient has been some ideas for obtaining commercial cold fusion which is under research. Another comments are taken into account which are beyond the scope of present research.
2.2. Plutonium Build Up in a Nuclear Pressurized Water Reactor
Consider a PWR which has been operating in a suggested time interval. In a PWR reactors, nuclear fuel is O pellets, and uranium consist of and isotopes (neglecting isotope). The most important isotopes of interest that have been produced as a result of uranium fuel depletion, are two isotopes of plutonium element, that is, Pu-239 and Pu-241. Because they can also be employed as fuel, like as -235 atom. The process at which -238 fresh fuel is converted into plutonium isotopes is shown in Figure 3. The appropriate magnitude of each neutron capture cross-section for corresponded nuclear (n,) reaction is given on each arrow in terms of barns [10]. According to the foregoing chain, a set of coupled first-order ordinary differential equation can be established to give the time-dependent concentration of some of interesting isotopes. This is done via the conservation of mass principle, that is, production rate minus consumption rate equals the net rate of change of isotope concentration. The concentration is represented by , and by . The other plutonium isotopes such as , , and are denoted by , , , respectively. All interested materials and isotopes are balanced as follows. depletion = Absorption of thermal neutrons in the cause fission. depletion = Absorption of thermal neutrons in and absorption of resonance neutron in the to produce , and absorption of fast neutrons in to cause fission. production = Absorption of thermal neutrons in the + Absorption of resonance neutrons from fission in +Absorption of resonance neutrons from fission in + Absorption of resonance neutrons from fission in -Absorption of thermal neutrons in + Absorption of fast neutrons from ,, and fissions in . production = Neutron absorption in the to produce minus neutron absorption in the to produce . production = Neutron absorption in the to produce minus neutron absorption in the to produce .Fission fragments production = Fission yields of + Fission yields of + Fission yields of + Fission yields of .
Therefore, according to our defined parameters, the coupled first-order differential equations which describes plutonium and uranium isotopes concentrations are given as: where is the average thermal neutron flux of the core (we have considered is independent of time and equal to a constant), and and are the nuclear thermal microscopic absorption cross section which refer to the desired isotopes. Also, other parameters are described in Table 3. As it is seen, the set of (2.3) are IVP, so it is apparent that we must know the values of atom densities at the initiation of fuel irradiation. But, it was stated that the core is initially loaded with fresh fuel and there are, in fact, no plutonium isotopes at starting time. Thus the only ones should be determined are and at the time of reactor startup. On the other hand uranium element is composed of two isotopes of - and -, in which in a typical PWR type the fuel is enriched to about averaged value of percent where the initial values are given in Table 2.
To solve the set of above IVP, (2.3), we make some simplifying assumptions in the first iteration such as the following. (i)Effective cross sections remain constant throughout the core and during fuel lifetime.(ii)Average neutron flux within the core is constant and is considered to be equal to .(iii)The time duration at which reactor fuel has to be replaced with the fresh fuel, due to neutronic and/or thermal hydraulic reasons, is about hours.
Using data given in Table 3, the set of (2.3) are solved using our presented MATLAB-linked solver. The solver was run and gave our desired results. Beginning of cycle (BOC) masses of -235, -238, important plutonium isotopes, fission fragments burnup/or buildup and also End of Cycle (EOC) masses are illustrated in Table 4.
2.3. Samarium-149 Build Up in a Nuclear Pressurized Water Reactor
The fission fragments are highly radioactive which undergo and emissions. Some of the fission fragments are highly neutron absorber materials and strongly affect neutronic balance within the core as if they act as a neutron poison. They tend to capture a neutron and form a nucleus which contains a neutron more. So as will be seen, as time goes on, fission fragments would be converted to some other atoms and it is necessary to make an estimation of about their atom density (number of atoms per unit volume within the fuel), with respect to irradiation time. According to the foregoing discussion, it is expected to have a completely different fuel at the end of fuel life with that originally loaded within the core. In most cases of interest, such as study of fission products poisoning, involved isotopes form a radioactive and neutron reacting chain in which its members are linked together via decay and reactions. Also some members of the chain are produced directly from -235 fission; that is, they have a finite yield from fission. Consider, for example, that -235 fission rate is , in which is -235 atom density, is the effective -235 microscopic fission cross section, and is the time dependent neutron flux within the core. So as a result, this amount of -235 atoms are undergoing fission per unit volume of the fuel per unit time. We define here fission yield of -th species, , as the ratio of -th atoms produced to -235 atoms undergoing fission. Consequently, constant formation rate of the -th nuclide per unit volume could be written as .
As shown in Figure 4, some members of the chain will have two different probable modes of disappearance, depending on whether the decay or neutron capture is more probable, they tend to make two completely different nuclei. This state of affair is taken into account in writing the rate equations for some nuclei. Using the so-called chain, we can develop appropriate ”rate equation” for individual nuclide per unit volume. Before this, we show some characteristics of the involved isotopes of the Sm-149 chain in Table 5.
According to Figure 4, a set of 12 coupled ordinary first-order differential equations that describe the rate of change of each of the 11 nuclei in the Sm-149 fission-product chain as well as -235 are written as follows: where in (2.7), is radiative capture cross section for reaction and in a similar manner, in (2.8) is for reaction. Also in (2.8) is a radioactive decay constant that , as a result of a decay, disintegrates to . Moreover, indicates direct fission yield for -th species from -235 thermal fission, and finally, (2.15) is a rate equation for -235 atom density in which is -235 effective fission cross section. This equation implies that -235 atom density decreases as an exponential function as time goes on. is time-dependent -235 atom density.
The set of equations of (2.4) through (2.15) should be solved simultaneously to give desired result and when the matrix of coefficient is established and further investigated, it turns out that this set of equations is a nonstiff one and here, is then solved using the Runge-Kuta method [11]. Using our MATLAB-linked solver as well as data given in Table 5, our calculated results are shown in Figure 5.
(a)
(b)
(c)
2.4. Xenon135 Build Up in a Nuclear Pressurized Water Reactor
Another poison of our interest, as the greatest fission product in a nuclear reactor, is xenon isotope. It is the most important neutron absorber (poison) in a typical PWR type such that it is produced directly from -235 nucleus fission and indirectly from decay of Te-135 chain. Its decay chain is in the form In addition, is a great neutron absorber and under the neutron flux within the reactor will be changed into . Similar to the previous samarium poisoning subsection and according to the so-called neutron absorbing and radioactive decay chain, we develop a set of six coupled first-order differential equations that describe the rate of change of each of the five nuclei in the fission-product chain as well as again -235 atom. They are as follows:
in which and stand for fission yields of and , respectively. Also, are associated radioactive decay constants; are associated absorption cross-section, and explicitly and are fission and absorption cross-sections for -235 nucleus, respectively. Constant values that were appeared in the set of equations of (2.17) through (2.22) are given in Table 6. Equations (2.17) through (2.22) are again coupled IVPs and are stiff case [12, 13]. These equations are solved simultaneously to give desired results; which we have focused on the Xe-135 concentration in the case of reactor power variation and results are given in Figure 6. In the first period, reactor operates at full power and xenon concentration increases toward to a constant value after about 40 hours. In the second period, reactor is shut down and therefore, xenon peaks after about 11 hours and then decreases. In the third period, reactor operates at full power and a same manner as like as period 1 for xenon behavior are obtained. In the fourth period, reactor operates at half of nominal power () and therefore a xenon peak occurs after 11 hours, but xenon steady-state concentration is more than previous period as we expected.
3. Conclusion
We have presented a computer package to solve first-order IVPs with constant and variable coefficients using MATLAB software, in which the solution of a given stiff or nonstiff coupled differential equations with known initial values were found and plotted. In the present paper, some well-known nuclear engineering dynamical problems, related to the IVPs, were given. A major application of IVPs to a real problem is the fuel depletion in a suggested PWR, where it is computed by the present MATLAB-linked “solver”. We used matrices form such as with known initial values in each case. But we have focused on the constant matrix, where its elements are multiplication of neutronic flux and material cross sections. Our results are good compared with the well-known texts [10, 14]. Obviously, our approach should be extended to a variable or coupled IVPs with variable coefficients for more accuracy, for instance, cross section is not fixed during fuel depletion [15, 16].
Our aim here is to bring out a MATLAB-linked solver for researchers to solve coupled IVPs numerically where it is appeared frequently in many cases of nuclear engineering problems. Reader may refer to the appendix to find our written MATLAB-linked solver program.
Appendix
Three programs which are named COUPLED, COEFFICIENT, and DEPLET must be written as an M-file and then saved in the work directory of the MATLAB software. The first program is
function dy = coupled(t,y) format('long','e') global Di
dy = zeros(Di,1); % a column vector
%disired variable
load moham if O==1 run coefficient
end
load H
dy=H*y
O=O+1;
save moham O
The second program is
function coeficent
disp(’************************************************************’)
disp(’Coupled Differential Equations is computed in the form of
Dy=H*y.’) disp(’ You should ENTER the H(n,m) array.’)
disp(’*********************************’) load Di
for n=1:Di
for m=1:Di
disp(’In the following, you can find desired (n,m) to RUN:’)
disp([n m])
H(n,m)=input(’Please ENTER value of H(n,m) for the above given
(n,m):’);
end end save H H
The 3rd program is:
This program compute and plot set of Coupled Differential
Equations and Inintial Values(IVP) Using MATLAB commands. To
start computation one must enetr number of unknowns and equations,
constants and choose desired numerical method. function deplet
disp('*********************************************')
Di=input('**Please ENTER the Number of Differential
Equations(Unknowns): ') save Di Di O=1; save moham O
%–––––––––––––––––––––––-
% xi's are initial conditions for unknowns.
B=zeros(1,Di); for w=1:Di;
disp('Insert initial values of Yi where i is:');
disp([w])
B(1,w)=input('Enter Yi: ');
end
%––––––––––––––––––––––
% t0 and t1 are the start and end points of time interval
disp('******************************************')
T0=input('Insert Start-point of the computations: ');
disp('****************************************')
T1=input('Insert End-point of the computations: ');
disp('************************************')
%––––––––––––––––––––––––
P=menu('What is the type of your coupled differential equation?
','NonStiff Equations','Stiff Equation');
if P==1;
A=menu('Which method you want for executing?',
'ode45 method','ode23 method',' ode113 method');
if A==1;
disp('This methos is Based on an explicit Runge-Kutta (4,5)
formula, the Dormand-Prince pair...')
disp('It is a one-step solver - in computing, it needs only ')
disp('the solution at the immediately preceding time point,.')
disp('In general, ode45 is the best function to
apply as "first try" for most problems.')
disp('###############################')
disp('**press any key to continue computations**')
disp('###############################')
pause
[t,y]=ode45(@coupled,[T0:1:T1],B);
J=input('Which variables you want for plotting? ');
plot(t,y(:,J))
%–––––––––––––––––––––––––––
elseif A==2
disp('This method is Based on an explicit Runge-Kutta (2,3) pair ')
disp('of Bogacki and Shampine. It may be more efficient')
disp('than ode45 at crude tolerances and in the ')
disp('presence of mild stiffness.')
disp('Like ode45, ode23 is a one-step solver.')
disp('###############################')
disp('**press any key to continue computations**')
disp('###############################')
pause
[t,y]=ode23(@coupled,[T0 T1],B)
J=input('Which variables you want for plotting? ');
plot(t,y(:,J))
elseif A==3
disp('Software will use variable order Adams-Bashforth-Moulton PECE
solver.')
disp('It may be more efficient than ode45 at stringent')
disp('tolerances and when the ODE function is particularly ')
disp('expensive to evaluate. ode113 is a multistep')
disp('solver - it normally needs the solutions')
disp('at several preceding time points ')
disp('###############################')
disp('**press any key to continue computations**')
disp('###############################')
pause
[t,y]=ode113(@coupled,[T0 T1],B)
J=input('which variables you want for plotting? ');
plot(t,y(:,J))
end
elseif P==2
G=menu('Which method you want for executing',
'ode15s method','ode23s','ode23t','ode23tb')
if G==1
disp('Software will use Variable-order solver based on the ')
disp('numerical differentiation formulas (NDFs).')
disp('Optionally it uses the backward differentiation formulas')
disp('BDFs, (also known as Gear method).')
disp('Like ode113, ode15s is a multistep solver.')
disp('If you suspect that a problem is stiff ')
disp('or if ode45 failed or was very inefficient, try ode15s')
disp('###############################')
disp('**press any key to continue computations**')
disp('###############################')
pause
[t,y]=ode15s(@coupled,[T0 T1],B)
J=input('Which variables you want for plotting? ');
plot(t,y(:,J))
elseif G==2
disp('This method is Based on a modified Rosenbrock formula of
order 2.')
disp('Because it is a one-step solver, ')
disp('it may be more efficient than ode15s ')
disp('at crude tolerances. It can solve some')
disp('kinds of stiff problems for which ode15s is not effective.')
disp('###############################')
disp('**press any key to continue computations**')
disp('###############################')
pause
[t,y]=ode23s(@coupled,[T0 T1],B)
J=input('Which variables you want for plotting? ');
plot(t,y(:,J))
elseif G==3
disp('software wil use an implementation of the trapezoidal rule ')
disp('using a "free" interpolant.')
disp('Use this solver if the problem')
disp('is only moderately stiff and you')
disp('need a solution without numerical damping.')
disp('###############################')
disp('**press any key to continue computations**')
disp('###############################')
pause
[t,y]=ode23t(@coupled,[T0 T1],B)
J=input('Which variables you want for plotting? ');
plot(t,y(:,J))
elseif G==4
disp('software will use an implementation of TR-BDF2,')
disp('an implicit Runge-Kutta formula with ')
disp('a first stage that is a trapezoidal ')
disp('rule step and a second stage that is a')
disp('backward differentiation formula of ')
disp('order 2. Like ode23s, this solver may ')
disp('be more efficient than ode15s at crude tolerances.')
disp('###############################')
disp('**press any key to continue computations**')
disp('###############################')
pause
[t,y]=ode23tb(@coupled,[T0 T1],B)
J=input('Which variables you want for plotting? ');
plot(t,y(:,J)) end disp('If you want to RUN this code again, you
must rewrite (re-Enter) options.') disp('TO start again, RUN
deplet.m') end
Acknowledgments
This work is supported under academic Grant no. 88-GR-ENG-6. The corresponding author wishes to acknowledge the Research Council of the Shiraz University for their financial support.