Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2009 / Article

Research Article | Open Access

Volume 2009 |Article ID 249162 |

F. Faghihi, M. Saidi Nezhad, "A MATLAB-Linked Solver to Find Fuel Depletion in a PWR, a Suggested VVER-1000 Type", Mathematical Problems in Engineering, vol. 2009, Article ID 249162, 18 pages, 2009.

A MATLAB-Linked Solver to Find Fuel Depletion in a PWR, a Suggested VVER-1000 Type

Academic Editor: Shijun Liao
Received11 Sep 2009
Accepted30 Nov 2009
Published15 Mar 2010


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 [13], 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𝑑𝑥1(𝑡)𝑑𝑡=𝑥1=𝑝11(𝑡)𝑥1++𝑝1𝑛(𝑡)𝑥𝑛,𝑑𝑥2(𝑡)𝑑𝑡=𝑥2=𝑝21(𝑡)𝑥1++𝑝2𝑛(𝑡)𝑥𝑛𝑑𝑥𝑛(𝑡)𝑑𝑡=𝑥𝑛=𝑝𝑛1(𝑡)𝑥1++𝑝𝑛𝑛(𝑡)𝑥𝑛,(1.1) where initial values of the dependent variables are:𝑥𝑖𝑡0=𝛼𝑖,𝛼𝑖=const.,(𝑖=1,2,,𝑛).(1.2) 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 𝑥1(𝑡) as an example, user should write “1" and then a new window will be opened and 𝑥1(𝑡) 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 𝑘=0,1,2,,𝑁), and therefore we can assume average of 𝑝𝑖𝑗s in each step-size:𝑝𝑖𝑗𝑡𝑘=𝑝𝑖𝑗,where𝑝𝑖𝑗=𝑝𝑖𝑗𝑡𝑘+𝑝𝑖𝑗𝑡𝑘+12.(1.3) 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.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 𝑇2 molecules. Thus the formed (𝑡𝜇) reacts with 𝐷2, 𝐷𝑇, or 𝑇2 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 𝜔0𝑠. 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 𝜔𝑠=(1𝑅)𝜔0𝑠.

Now, consider a homogeneous media in which the 𝜇CF is carried out [79]. 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 (𝜇CF) are given by𝑑𝑁𝜇(𝑡)𝑑𝑡=𝜆0𝑁𝜇(𝑡)𝜆𝑎𝑐𝑑𝜙𝑚𝑁𝜇(𝑡)𝜆𝑎𝑐𝑡𝜙𝑚𝑁𝜇(𝑡)+1𝜔𝑠𝜆𝜇𝑑𝑡𝑁𝑑𝑡𝜇(𝑡),𝑑𝑁𝜇𝑑(𝑡)𝜆𝑑𝑡=0+1𝜏𝑁𝜇𝑑(𝑡)+𝜆𝑎𝑐𝑑𝜙𝑚𝑁𝜇(𝑡)𝜆𝑑𝑡𝑐𝑡𝜙𝑚𝑁𝜇𝑑(𝑡),𝑑𝑁𝜇𝑡(𝑡)𝜆𝑑𝑡=0+1𝜏𝑁𝜇𝑡(𝑡)+𝜆𝑎𝑐𝑡𝜙𝑚𝑁𝜇(𝑡)+𝜆𝑑𝑡𝑐𝑡𝜙𝑚𝑁𝜇𝑑(𝑡)𝜆𝑑𝑡𝜇𝑐𝑑𝜙𝑚𝑁𝜇𝑡(𝑡),𝑑𝑁𝑑𝑡𝜇(𝑡)𝑑𝑡=𝜆𝑑𝑡𝜇𝑐𝑑𝜙𝑚𝑁𝜇𝑡(𝑡)𝜆𝑓𝑑𝑡𝜇𝑁𝑑𝑡𝜇(𝑡)𝜆0𝑁𝑑𝑡𝜇(𝑡),𝑑𝑁𝑛(𝑡)𝑑𝑡=𝜆𝑓𝑑𝑡𝜇𝑁𝑑𝑡𝜇(𝑡),(2.1) 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 𝜆0 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).

Process Parameter Value

Relative concentration of deuterium 𝑐 𝑑 0.5
Relative concentration of tritium 𝑐 𝑡 0.5
Muon decay constant 𝜇 0 0 . 4 5 5 × 1 0 6 s 1
Muonic-atom formation rate 𝜆 𝑎 1 × 1 0 1 0 s 1
Media ion density 𝜙 𝑚 Liquid hydrogen density
Sticking coefficient 𝜔 𝑠 0.008
Fusion rate of 𝑑 𝑡 𝜇 molecule 𝜆 𝑓 𝜇 𝑑 𝑡 1 . 1 × 1 0 1 2 s 1
Formation rate of 𝑑 𝑡 𝜇 molecule 𝜆 𝜇 𝑑 𝑡 2 . 8 × 1 0 8 s 1
Time of leakage 𝜏 7 . 0 8 × 1 0 7 s
Muon exchange rate between 𝑑 𝜇 and 𝑡 𝜇 atoms 𝜆 𝑑 𝑡 1 × 1 0 1 0 s 1

We have solved these coupled dynamics equations in time range of [0 to 2.2×106sec], the muon life-time, using our MATLAB-linked solver with the following initial conditions:𝑁𝜇(𝑡=0)=1,𝑁𝑘(𝑡=0)=0,𝑘=𝜇𝑑,𝜇𝑡,𝑑𝑡𝜇,𝑛.(2.2) According to the coupled dynamical equation (2.1) and also Figure 1, the calculated neutrons, 𝑁𝑛(𝑡=2.2×106s) for each one inserted muon, corresponds to the muon cycling rate of the explained cold fusion (𝜇CF), 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, 𝑁𝑛(𝑡=2.2𝜇s)110 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 110×17.6=1.94GeV. For producing one muon, due to an accelerator, we must expense about 4GeV 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 UO2 pellets, and uranium consist of 235U and 238U isotopes (neglecting 234U 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 U -235 atom. The process at which U -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 238U concentration is represented by 𝑁28, and 239Pu by 𝑁49. The other plutonium isotopes such as 240Pu, 241Pu, and 242Puare denoted by 𝑁40, 𝑁41, 𝑁42, respectively. All interested materials and isotopes are balanced as follows.235U depletion = Absorption of thermal neutrons in the 235U cause fission.238U depletion = Absorption of thermal neutrons in 238U and absorption of resonance neutron in the 238U to produce 239Pu, and absorption of fast neutrons in 238U to cause fission.239Pu production = Absorption of thermal neutrons in the 238U+ Absorption of resonance neutrons from 235U fission in 238U+Absorption of resonance neutrons from 239Pu fission in 238U+ Absorption of resonance neutrons from 241Pu fission in 238U-Absorption of thermal neutrons in 239Pu+ Absorption of fast neutrons from 235U,238U, and 239Pu fissions in 238U.240Pu production = Neutron absorption in the 239Pu to produce 240Pu minus neutron absorption in the 240Pu to produce 241Pu.241Pu production = Neutron absorption in the 240Pu to produce 241Pu minus neutron absorption in the 241Pu to produce 242Pu.Fission fragments production = Fission yields of 235U + Fission yields of 238U + Fission yields of 239Pu+ Fission yields of 241Pu.

Therefore, according to our defined parameters, the coupled first-order differential equations which describes plutonium and uranium isotopes concentrations are given as:𝑑𝑁25(𝑡)𝑑𝑡=𝜎𝑎25𝜙𝑁25(𝑡),𝑑𝑁28(𝑡)𝑑𝑡=𝜎𝑎28𝜙𝑁28(𝑡)𝜎𝑎28𝜖𝑃1𝜙𝑁28(𝑡)𝜎𝑎28𝜖𝜙𝑁28(𝑡),𝑑𝑁49(𝑡)𝑑𝑡=𝜎𝑎28𝜙𝑁28(𝑡)+𝜂25𝜖𝑃1(1𝑝)𝜎𝑎25𝜙𝑁25(𝑡)+𝜂49𝜖𝑃1(1𝑝)𝜎𝑎49𝜙𝑁49(𝑡)+𝜂41𝜖𝑃1(1𝑝)𝜎𝑎41𝜙𝑁41(𝑡)𝜎𝑎49𝜙𝑁49+𝛼(𝑡),281+𝛼28𝜖1𝜂28𝜂125𝜎𝑎25𝑁25(𝑡)+𝜂49𝜎𝑎49𝑁49(𝑡)+𝜂41𝜎𝑎41𝑁41(𝑡)𝜙,𝑑𝑁40(𝑡)=𝛼𝑑𝑡49𝜎𝑎491+𝛼49𝜙𝑁49(𝑡)𝜎𝑎40𝜙(𝑡)𝑁40(𝑡),𝑑𝑁41(𝑡)𝑑𝑡=𝜎𝑎40𝜙𝑁40(𝑡)𝜎𝑎41𝜙𝑁41(𝑡),(2.3) where 𝜙(𝑡) is the average thermal neutron flux of the core (we have considered 𝜙(𝑡) is independent of time and equal to a constant), and 𝜎𝑎𝛼𝛽(𝛼=2,4 and 𝛽=0,5,8,9) 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 UO2 fuel and there are, in fact, no plutonium isotopes at starting time. Thus the only ones should be determined are 𝑁28 and 𝑁25 at the time of reactor startup. On the other hand uranium element is composed of two isotopes of U -235 and U-238, in which in a typical PWR type the fuel is enriched to about averaged value of 3.5 percent where the initial values are given in Table 2.

Material Mass, kg Concentration, N( × 1 0 2 4 ), atoms/ c m 3

U O 2 932460.0217391
2 3 5 U 3058.0430.0007700
2 3 8 U 84314.6200.0209645

NuclideSubscript 𝜎 𝑎 (barn) 𝜂 𝛼

2 3 5 U 255561.960.2398
2 3 8 U 28 2.2342 2.3432 0.1907
2 3 9 P u 49 1618.2 1.86 0.5430
2 4 0 P u 40 2616.8
2 4 1 P u 41 1567.3 2.223 0.3765

𝜖 is fast fission factor = 1 . 0 4 7 6 .
𝑝 is resonance scape probability = 0 . 7 7 2 5 .
𝑃 1 is fast to 2 3 8 U resonance probability = 0 . 9 8 8 9 .
𝜂 is neutrons produced per neutrons absorbed.
𝛼 is ratio of capture to fission cross section.

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 3.5×1013neutrons/cm2s.(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 7300 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 U-235, U-238, important plutonium isotopes, fission fragments burnup/or buildup and also End of Cycle (EOC) masses are illustrated in Table 4.

Material BOCEOC

2 3 5 U 30582632
2 3 8 U 8431483604
P u 0210
Fission fragments0926

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 U-235 fission; that is, they have a finite yield from fission. Consider, for example, that U-235 fission rate is 𝑁𝑓𝜎𝑓𝜙(𝑡), in which 𝑁𝑓 is U-235 atom density, 𝜎𝑓 is the effective U-235 microscopic fission cross section, and 𝜙(𝑡) is the time dependent neutron flux within the core. So as a result, this amount of U-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 U-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.

NuclideNuclide indicationHalf-lifeAbsorption cross section (b)Direct yield from U - 2 3 5

1 4 7 N d 𝑁 1 11.1 d00.0236
1 4 7 P m 𝑁 2 2.62 y845.72 to 1 4 8 m P m , 448.23 to 1 4 8 P m 0
1 4 7 S m 𝑁 3 274.2 0
1 4 7 m P m 𝑁 4 42 d (7 % to 1 4 8 P m , 93 % to 1 4 8 S m )31964 0
1 4 8 P m 𝑁 5 5.4 d13858 0
1 4 9 P m 𝑁 6 53.1 h1105.6 0.0113
1 4 9 S m 𝑁 7 73635 0
1 5 0 P m 𝑁 8 2.7 h0 0
1 5 0 S m 𝑁 9 158.38 0
1 5 1 S m 𝑁 1 0 87 y9734.5 0.0044
1 5 2 S m 𝑁 1 1 813.01 0.00281

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 U-235 are written as follows: 14760Nd𝑑𝑁1𝑑𝑡=𝑦1𝑁𝑓𝜎𝑓𝜎𝜙1𝜙+𝜆1𝑁1,(2.4)14761Pm𝑑𝑁2𝑑𝑡=𝑦2𝑁𝑓𝜎𝑓𝜙+𝜆1𝑁1𝜎2𝜙+𝜆2𝑁2,(2.5)14762Sm𝑑𝑁3𝑑𝑡=𝑦3𝑁𝑓𝜎𝑓𝜙+𝜆2𝑁2𝜎3𝜙𝑁3,(2.6)148𝑚61Pm𝑑𝑁4𝑑𝑡=𝜎24𝜙𝑁2𝜎4𝜙+𝜆4𝑁4,(2.7)14861Pm𝑑𝑁5𝑑𝑡=𝜆45𝑁4+𝜎25𝜙𝑁2𝜎5𝜙+𝜆5𝑁5,(2.8)14961Pm𝑑𝑁6𝑑𝑡=𝑦6𝑁𝑓𝜎𝑓𝜙+𝜎4𝜙𝑁4+𝜎5𝜙𝑁5𝜎6𝜙+𝜆6𝑁6,(2.9)14962Sm𝑑𝑁7𝑑𝑡=𝑦7𝑁𝑓𝜎𝑓𝜙+𝜆6𝑁6𝜎7𝜙𝑁7,(2.10)15061Pm𝑑𝑁8𝑑𝑡=𝜎6𝜙𝑁6𝜎8𝜙+𝜆8𝑁8,(2.11) where in (2.7), 𝜎24 is radiative capture cross section for 147Pm(𝑛,𝛾)148𝑚Pm reaction and in a similar manner, in (2.8) 𝜎25 is for 147Pm(𝑛,𝛾)148Pm reaction. Also in (2.8) 𝜆45 is a radioactive decay constant that 148𝑚Pm, as a result of a 𝛽 decay, disintegrates to 148Pm. Moreover, 𝑦𝑖 indicates direct fission yield for 𝑖-th species from U-235 thermal fission, and finally, (2.15) is a rate equation for U-235 atom density in which 𝜎𝑓 is U-235 effective fission cross section. This equation implies that U-235 atom density decreases as an exponential function as time goes on. 𝑁𝑓 is time-dependent U-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.

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 xenon135 isotope. It is the most important neutron absorber (poison) in a typical PWR type such that it is produced directly from U-235 nucleus fission and indirectly from decay of Te-135 chain. Its decay chain is in the formFission13552Te13553I13554Xe13555Cs13556Ba.(2.16) In addition, 13554Xe is a great neutron absorber and under the neutron flux within the reactor will be changed into 13654Xe. 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 U-235 atom. They are as follows:

135Te𝑑𝑇𝑑𝑡=𝑦𝑇𝑁𝑓𝜎𝑓𝜙𝜆𝑇𝑇(2.17)135I𝑑𝐼𝑑𝑡=𝜆𝑇𝑇𝜆𝐼𝐼(2.18)135Xe𝑑𝑋𝑑𝑡=𝑦𝑋𝑁𝑓𝜎𝑓𝜙+𝜆𝐼𝐼𝜎𝑋𝑋𝜙𝜆𝑋𝑋(2.19)135Cs𝑑𝐶𝑑𝑡=𝜆𝑋𝑋𝜎𝐶𝐶𝜙𝜆𝐶𝐶(2.20)135Ba𝑑𝐵𝑑𝑡=𝜆𝐶𝐶(2.21)235U𝑑U𝑑𝑡=𝜎𝑎𝑁𝑓𝜙(2.22) in which 𝑦𝑇 and 𝑦𝑋 stand for fission yields of 135Te and 135Xe, respectively. Also, 𝜆 are associated radioactive decay constants; 𝜎 are associated absorption cross-section, and explicitly 𝜎𝑓 and 𝜎𝑎 are fission and absorption cross-sections for U-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 (50%) and therefore a xenon peak occurs after 11 hours, but xenon steady-state concentration is more than previous period as we expected.

NuclideNuclide indicationHalf lifeRadioactive decay constant s 1 Absorption cross-section (b)Direct yield from U - 2 3 5

1 3 5 T e 𝑇 2 9 s 0 . 0 2 3 9 0 0 . 0 6 0 9
1 3 5 I 𝐼 6 . 7 h 2 . 8 7 × 1 0 5 00
1 3 5 X e 𝑋 9 . 2 h 2 . 0 9 × 1 0 5 2 . 6 4 × 1 0 6 0.0032
1 3 5 C s 𝐶 3 × 1 0 6 yr 7 . 3 × 1 0 1 5 17.20
1 3 5 B a 𝐵 stable000

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.


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


load H



save moham O

The second program is

function coeficent


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


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


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:');


B(1,w)=input('Enter Yi: ');



   % t0 and t1 are the start and end points of time interval


    T0=input('Insert Start-point of the computations: ');


    T1=input('Insert End-point of the computations: ');



    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('**press any key to continue computations**')




  J=input('Which variables you want for plotting? ');



   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('**press any key to continue computations**')



[t,y]=ode23(@coupled,[T0 T1],B)

J=input('Which variables you want for plotting? ');


elseif A==3

  disp('Software will use variable order Adams-Bashforth-Moulton PECE


  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('**press any key to continue computations**')



    [t,y]=ode113(@coupled,[T0 T1],B)

  J=input('which variables you want for plotting? ');



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('**press any key to continue computations**')



          [t,y]=ode15s(@coupled,[T0 T1],B)

    J=input('Which variables you want for plotting? ');


    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('**press any key to continue computations**')



    [t,y]=ode23s(@coupled,[T0 T1],B)

  J=input('Which variables you want for plotting? ');


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('**press any key to continue computations**')



    [t,y]=ode23t(@coupled,[T0 T1],B)

J=input('Which variables you want for plotting? ');


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('**press any key to continue computations**')



          [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


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.


  1. F. Faghihi and K. Hadad, “Numerical solutions of coupled differential equations and initial values using Maple software,” Applied Mathematics and Computation, vol. 155, no. 2, pp. 563–572, 2004. View at: Publisher Site | Google Scholar
  2. F. Faghihi, “Coupled differential equations and initial values using softwares; an applied program,” in Proceedings of the International Conference on Applied Mathematics and Numerical Analysis (ICNAAM '05), Rhodes, Greece, September 2005. View at: Google Scholar
  3. F. Faghihi, “Numerical minimization procedures of the adiabatic approximation wavefunction,” Applied Mathematics and Computation, vol. 162, no. 3, pp. 1013–1021, 2005. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  4. F. Faghihi and M. R. Nematollahi, “Nuclear fuel depletion analysis using MATLAB software,” International Journal of Modern Physics C, vol. 17, no. 6, pp. 805–815, 2006. View at: Google Scholar
  5. MATLAB Software Help, Version, License No. 205273, Release 2004.
  6. K. Nagamine, T. Matsuzaki, K. Ishida, S. N. Nakamura, and N. Kawamura, “Implications of the recent D-T μCF experiments at RIKEN-RAL and near-future directions,” Hyperfine Interactions, vol. 119, no. 1–4, pp. 273–280, 1999. View at: Google Scholar
  7. F. Faghihi and M. R. Eskandari, “Nuclear fusion rate study of A muonic molecule via nuclear threshold resonances,” International Journal of Modern Physics E, vol. 14, no. 8, pp. 1213–1221, 2005. View at: Publisher Site | Google Scholar
  8. M. R. Eskandari and F. Faghihi, “Forced fusion in the excited state of dtμ muonic-molecule and its possible drivers,” International Journal of Modern Physics C, vol. 14, no. 10, pp. 1281–1293, 2003. View at: Publisher Site | Google Scholar
  9. M. R. Eskandari, F. Faghihi, and R. Gheisari, “μCF study of D/T and H/D/T mixtures in homogeneous and inhomogeneous medium, and comparison of their fusion yields,” International Journal of Modern Physics C, vol. 13, no. 5, pp. 689–705, 2002. View at: Publisher Site | Google Scholar
  10. M. Benedict, T. H. Pigford, and H. W. Levi, Nuclear Chemical Engineering, McGraw-Hill, New York, NY, USA, 1981.
  11. M. Venkatesulu and P. D. N. Srinivas, “Solutions of nonstandard initial value problems for a first order ordinary differential equation,” Journal of Applied Mathematics and Simulation, vol. 2, no. 4, pp. 225–237, 1989. View at: Google Scholar | Zentralblatt MATH | MathSciNet
  12. S. Sekar, “Analysis of linear and nonlinear stiff problems using the RK-Butcher algorithm,” Mathematical Problems in Engineering, vol. 2006, Article ID 39246, 15 pages, 2006. View at: Publisher Site | Google Scholar | MathSciNet
  13. I. P. Gavrilyuk, M. Hermann, M. V. Kutniv, and V. L. Makarov, “Difference schemes for nonlinear BVPs using Runge-Kutta IVP-solvers,” Advances in Difference Equations, vol. 2006, Article ID 12167, 29 pages, 2006. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  14. J. R. Lamarsh, Introduction to Nuclear Engineering, Prentice Hall, Upper Saddle River, NJ, USA, 2001.
  15. M. Legua, I. Morales, and L. M. Sánchez Ruiz, “Resolution of first- and second-order linear differential equations with periodic inputs by a computer algebra system,” Mathematical Problems in Engineering, vol. 2008, Article ID 654820, 7 pages, 2008. View at: Publisher Site | Google Scholar
  16. A. J. Arenas, G. González-Parra, L. Jódar, and R.-J. Villanueva, “Piecewise finite series solution of nonlinear initial value differential problem,” Applied Mathematics and Computation, vol. 212, no. 1, pp. 209–215, 2009. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet

Copyright © 2009 F. Faghihi and M. Saidi Nezhad. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

More related articles

 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

Article of the Year Award: Outstanding research contributions of 2020, as selected by our Chief Editors. Read the winning articles.