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

Authors: F. Faghihi; M. Saidi Nezhad
Journal: Mathematical Problems in Engineering
 (2009)

---

## 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 aD-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.


---

## Body

## 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(1.1)dx1(t)dt=x1′=p11(t)x1+⋯+p1n(t)xn,dx2(t)dt=x2′=p21(t)x1+⋯+p2n(t)xn…⋮dxn(t)dt=xn′=pn1(t)x1+⋯+pnn(t)xn,
where initial values of the dependent variables are:(1.2)xi(t0)=αi,αi=const.,(i=1,2,…,n).
Here, t is used for the independent variable and may refer to time in a dynamical problems, and xi(t) stands for dependent variables.Coupled IVPs with constant coefficients. First, we consider the IVPs with constant coefficients, or in other words constant pij, 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 ofxi (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 pijs and execution would be begun. After execution, dependent variables (xi(t))s will be computed according to the desired numerical method, and user can plot xi(t)s in the computational interval. For plotting x1(t) as an example, user should write “1" and then a new window will be opened and x1(t) 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 pij(t). In this case, the computational interval [ab] should be divided into many step-size intervals so that variations of pijs are ignored in each step-size (each step-size is equal to h=(b-a)/N, where  N is a desired grid interval and therefore tk=a+hk, where k=0,1,2,…,N), and therefore we can assume average of pijs in each step-size:(1.3)pij(tk)=p̅ij,wherep̅ij=pij(tk)+pij(tk+1)2.
In this case again we return into the coupled differential equations with constant coefficients which can be solved by the DEPLET.m. Calculated x(tk), 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 aD-T 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 D-T mixture, either (dμ) or a (tμ) atom is formed, with a probability more or less, proportional to the relative concentration of D (cd) and T (ct). Because of the difference between (dμ) and (tμ) in the binding energies of their atomic states, μ- in (dμ) undergoes a transfer reaction to tritium yielding (tμ) during a collision with the surrounding tritium in either D-T or T2 molecules. Thus the formed (tμ) reacts with D2, DT, or T2 to form a muonic-molecule dtμ at a rate of λdtμ followed by a fusion reaction occurring from a molecular state of the (dtμ). The fusion takes place and a 14-Mev neutron and a 3.6-Mev α-particle are emitted. After the fusion reaction inside the (dtμ) 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 ωs0. 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 R. 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 ωs=(1-R)ωs0.Figure 1
Reaction cycle of the Muon Catalyzed Fusion,μCF.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 (λa, and λdtμ), and fusion decay rates (λdtμf) 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(2.1)dNμ(t)dt=-λ0Nμ(t)-λacdϕmNμ(t)-λactϕmNμ(t)+(1-ωs)λμdtNdtμ(t),dNμd(t)dt=-(λ0+1τ)Nμd(t)+λacdϕmNμ(t)-λdtctϕmNμd(t),dNμt(t)dt=-(λ0+1τ)Nμt(t)+λactϕmNμ(t)+λdtctϕmNμd(t)-λdtμcdϕmNμt(t),dNdtμ(t)dt=λdtμcdϕmNμt(t)-λdtμfNdtμ(t)-λ0Ndtμ(t),dNn(t)dt=λdtμfNdtμ(t),
where cp and ct are the relative concentration of deuterium and tritium, respectively. The muonic-atoms formation rates are given by λi (i=μd,μt), and the muonic-molecule formation rate of dtμ is given by λdtμ. The dtμ fusion rate is shown by λdtμf. The possible leakage rate of muonic-atoms is proportional to Nμd/τ and Nμt/τ, and also λ0 is the muon decay constant. In (2.1), ϕm is the dimensionless ion density of the media and is proportional to liquid hydrogen density (close to 0.07). As said before, ωs is the muon effective sticking coefficient. Table 1 gives values of the constants for solving (2.1).Table 1
Constant values for solving (2.1) in 500°C temperature and media with one liquid hydrogen density. Also, it is supposed that there are equal concentrations of deuterium and tritium in the media.
ProcessParameterValueRelative concentration of deuteriumcd0.5Relative concentration of tritiumct0.5Muon decay constantμ00.455×106s-1Muonic-atom formation rateλa1×1010s-1Media ion densityϕmLiquid hydrogen densitySticking coefficientωs0.008Fusion rate ofdtμ moleculeλμdtf1.1×1012s-1Formation rate ofdtμ moleculeλμdt2.8×108s-1Time of leakageτ7.08×10-7sMuon exchange rate betweendμ and tμ atomsλdt1×1010s-1We have solved these coupled dynamics equations in time range of [0 to2.2×10-6sec], the muon life-time, using our MATLAB-linked solver with the following initial conditions:(2.2)Nμ(t=0)=1,Nk(t=0)=0,k=μd,μt,dtμ,n.
According to the coupled dynamical equation (2.1) and also Figure 1, the calculated neutrons, Nn(t=2.2×10-6s) 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: d+t⇀α+n).At the end of running, neutron concentration is plotted in our calculated time interval and is given in Figure2. According to Figure 2, Nn(t=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 ωs 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.Figure 2
:dNn(t)/dt is plotted versus time using MATLAB-linked solver. In our real physical conditions and our suggested model here Nn(t=2.2μs)≃110, is the muon cycling rate of the mentioned μCF.


### 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 isUO2 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 N28, and 239Pu by N49. The other plutonium isotopes such as 240Pu, 241Pu, and 242Puare denoted by N40, N41, N42, 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 of235U + Fission yields of 238U + Fission yields of 239Pu+ Fission yields of 241Pu.Figure 3
U-238 is converted to plutonium isotopes through the above chains.Therefore, according to our defined parameters, the coupled first-order differential equations which describes plutonium and uranium isotopes concentrations are given as:(2.3)dN25(t)dt=-σa25ϕN25(t),dN28(t)dt=-σa28ϕN28(t)-σa28ϵP1ϕN28(t)-σa28ϵϕN28(t),dN49(t)dt=σa28ϕN28(t)+η25ϵP1(1-p)σa25ϕN25(t)+η49ϵP1(1-p)σa49ϕN49(t)+η41ϵP1(1-p)σa41ϕN41(t)-σa49ϕN49(t),+α281+α28ϵ-1η28-1[η25σa25N25(t)+η49σa49N49(t)+η41σa41N41(t)]ϕ,dN40(t)dt=α49σa491+α49ϕN49(t)-σa40ϕ(t)N40(t),dN41(t)dt=σa40ϕN40(t)-σa41ϕN41(t),
where ϕ(t) is the average thermal neutron flux of the core (we have considered ϕ(t) is independent of time and equal to a constant), and σaαβ(α=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 N28 and N25 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.Table 2
Masses and atom densities of each fuel isotopes at the time of reactor startup.
MaterialMass, kgConcentration,N(×1024), atoms/cm3UO2932460.0217391235U3058.0430.0007700238U84314.6200.0209645Table 3
Effective properties of nuclide for thermal neutron in 1000 MWe PWR (3000 MW thermal power).
NuclideSubscriptσa(barn)ηα235U255561.960.2398238U282.23422.34320.1907239Pu491618.21.860.5430240Pu402616.8241Pu411567.32.2230.3765ϵ is fast fission factor=1.0476.p is resonance scape probability=0.7725.P1 is fast to 238U resonance probability=0.9889.η 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 to3.5×1013neutrons/cm2·s.(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 about7300 hours.Using data given in Table3, 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.Table 4
Fuel composition in a PWR after 1 year irradiation. All values are in terms of kg.
MaterialBOCEOC235U30582632238U8431483604Pu0210Fission 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 (n,γ) 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 Nfσfϕ(t), in which Nf is U-235 atom density, σf is the effective U-235 microscopic fission cross section, and ϕ(t) 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 i-th species, yi, as the ratio of i-th atoms produced to U-235 atoms undergoing fission. Consequently, constant formation rate of the i-th nuclide per unit volume could be written as Pi=yiNfσfϕ(t).As shown in Figure4, 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.Table 5
Nuclear properties forSm decay chain.
NuclideNuclide indicationHalf-lifeAbsorption cross section (b)Direct yield fromU-235147NdN111.1 d00.0236147PmN22.62 y845.72 to148mPm, 448.23 to 148Pm0147SmN3∞274.20147mPmN442 d (7% to 148Pm, 93% to 148Sm )319640148PmN55.4 d138580149PmN653.1 h1105.60.0113149SmN7∞736350150PmN82.7 h00150SmN9∞158.380151SmN1087 y9734.50.0044152SmN11∞813.010.00281Figure 4
The fission-product chain leading to Sm-149.According to Figure4, 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:
(2.4)60147Nd:dN1dt=y1Nfσfϕ-(σ1ϕ+λ1)N1,(2.5)61147Pm:dN2dt=y2Nfσfϕ+λ1N1-(σ2ϕ+λ2)N2,(2.6)62147Sm:dN3dt=y3Nfσfϕ+λ2N2-(σ3ϕ)N3,(2.7)61148mPm:dN4dt=σ24ϕN2-(σ4ϕ+λ4)N4,(2.8)61148Pm:dN5dt=λ45N4+σ25ϕN2-(σ5ϕ+λ5)N5,(2.9)61149Pm:dN6dt=y6Nfσfϕ+σ4ϕN4+σ5ϕN5-(σ6ϕ+λ6)N6,(2.10)62149Sm:dN7dt=y7Nfσfϕ+λ6N6-σ7ϕN7,(2.11)61150Pm:dN8dt=σ6ϕN6-(σ8ϕ+λ8)N8,(2.12)62150Sm:dN9dt=σ7ϕN7+λ8N8-σ9ϕN9,(2.13)62151Sm:dN10dt=y10Nfσfϕ+σ9ϕN9-(σ10ϕ+λ10)N10,(2.14)62152Sm:dN11dt=y11Nfσfϕ+σ10ϕN10-σ11ϕN11,(2.15)92235U:dNfdt=-σfϕNf,
where in (2.7), σ24 is radiative capture cross section for 147Pm(n,γ)148mPm reaction and in a similar manner, in (2.8) σ25 is for 147Pm(n,γ)148Pm reaction. Also in (2.8) λ45 is a radioactive decay constant that 148mPm, as a result of a β decay, disintegrates to 148Pm. Moreover, yi indicates direct fission yield for i-th species from U-235 thermal fission, and finally, (2.15) is a rate equation for U-235 atom density in which σf is U-235 effective fission cross section. This equation implies that U-235 atom density decreases as an exponential function as time goes on. Nf 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.(a) Uranium-235 depletion. (b) Samarium-149 atom density variation as a function of irradiation time. (c) Atom density variation of other important members of the samarium chain.
(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 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 form(2.16)Fission⇀52135Te⇀53135I⇀54135Xe⇀55135Cs⇀56135Ba.
In addition, 54135Xe is a great neutron absorber and under the neutron flux within the reactor will be changed into 54136Xe. 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:(2.17)Te135:dTdt=yTNfσfϕ-λTT(2.18)I135:dIdt=λTT-λII(2.19)Xe135:dXdt=yXNfσfϕ+λII-σXXϕ-λXX(2.20)Cs135:dCdt=λXX-σCCϕ-λCC(2.21)Ba135:dBdt=λCC(2.22)U235:dUdt=-σaNfϕ
in which yT and yX stand for fission yields of 135Te and 135Xe, respectively. Also, λ are associated radioactive decay constants; σ are associated absorption cross-section, and explicitly σf and σa 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.Table 6
Nuclear properties of fission products of mass135.
NuclideNuclide indicationHalf lifeRadioactive decay constants-1Absorption cross-section (b)Direct yield fromU-235135TeT29 s0.023900.0609135II6.7 h2.87×10-500135XeX9.2 h2.09×10-52.64×1060.0032135CsC3×106 yr7.3×10-1517.20135BaBstable000Figure 6
Xenon135 atom density as well as power variations in the four-time interval is given.


## 2.1. Muon Catalyzed Fusion System

The basic process of the muon catalyzed fusion in aD-T 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 D-T mixture, either (dμ) or a (tμ) atom is formed, with a probability more or less, proportional to the relative concentration of D (cd) and T (ct). Because of the difference between (dμ) and (tμ) in the binding energies of their atomic states, μ- in (dμ) undergoes a transfer reaction to tritium yielding (tμ) during a collision with the surrounding tritium in either D-T or T2 molecules. Thus the formed (tμ) reacts with D2, DT, or T2 to form a muonic-molecule dtμ at a rate of λdtμ followed by a fusion reaction occurring from a molecular state of the (dtμ). The fusion takes place and a 14-Mev neutron and a 3.6-Mev α-particle are emitted. After the fusion reaction inside the (dtμ) 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 ωs0. 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 R. 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 ωs=(1-R)ωs0.Figure 1
Reaction cycle of the Muon Catalyzed Fusion,μCF.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 (λa, and λdtμ), and fusion decay rates (λdtμf) 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(2.1)dNμ(t)dt=-λ0Nμ(t)-λacdϕmNμ(t)-λactϕmNμ(t)+(1-ωs)λμdtNdtμ(t),dNμd(t)dt=-(λ0+1τ)Nμd(t)+λacdϕmNμ(t)-λdtctϕmNμd(t),dNμt(t)dt=-(λ0+1τ)Nμt(t)+λactϕmNμ(t)+λdtctϕmNμd(t)-λdtμcdϕmNμt(t),dNdtμ(t)dt=λdtμcdϕmNμt(t)-λdtμfNdtμ(t)-λ0Ndtμ(t),dNn(t)dt=λdtμfNdtμ(t),
where cp and ct are the relative concentration of deuterium and tritium, respectively. The muonic-atoms formation rates are given by λi (i=μd,μt), and the muonic-molecule formation rate of dtμ is given by λdtμ. The dtμ fusion rate is shown by λdtμf. The possible leakage rate of muonic-atoms is proportional to Nμd/τ and Nμt/τ, and also λ0 is the muon decay constant. In (2.1), ϕm is the dimensionless ion density of the media and is proportional to liquid hydrogen density (close to 0.07). As said before, ωs is the muon effective sticking coefficient. Table 1 gives values of the constants for solving (2.1).Table 1
Constant values for solving (2.1) in 500°C temperature and media with one liquid hydrogen density. Also, it is supposed that there are equal concentrations of deuterium and tritium in the media.
ProcessParameterValueRelative concentration of deuteriumcd0.5Relative concentration of tritiumct0.5Muon decay constantμ00.455×106s-1Muonic-atom formation rateλa1×1010s-1Media ion densityϕmLiquid hydrogen densitySticking coefficientωs0.008Fusion rate ofdtμ moleculeλμdtf1.1×1012s-1Formation rate ofdtμ moleculeλμdt2.8×108s-1Time of leakageτ7.08×10-7sMuon exchange rate betweendμ and tμ atomsλdt1×1010s-1We have solved these coupled dynamics equations in time range of [0 to2.2×10-6sec], the muon life-time, using our MATLAB-linked solver with the following initial conditions:(2.2)Nμ(t=0)=1,Nk(t=0)=0,k=μd,μt,dtμ,n.
According to the coupled dynamical equation (2.1) and also Figure 1, the calculated neutrons, Nn(t=2.2×10-6s) 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: d+t⇀α+n).At the end of running, neutron concentration is plotted in our calculated time interval and is given in Figure2. According to Figure 2, Nn(t=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 ωs 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.Figure 2
:dNn(t)/dt is plotted versus time using MATLAB-linked solver. In our real physical conditions and our suggested model here Nn(t=2.2μs)≃110, is the muon cycling rate of the mentioned μCF.


## 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 isUO2 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 N28, and 239Pu by N49. The other plutonium isotopes such as 240Pu, 241Pu, and 242Puare denoted by N40, N41, N42, 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 of235U + Fission yields of 238U + Fission yields of 239Pu+ Fission yields of 241Pu.Figure 3
U-238 is converted to plutonium isotopes through the above chains.Therefore, according to our defined parameters, the coupled first-order differential equations which describes plutonium and uranium isotopes concentrations are given as:(2.3)dN25(t)dt=-σa25ϕN25(t),dN28(t)dt=-σa28ϕN28(t)-σa28ϵP1ϕN28(t)-σa28ϵϕN28(t),dN49(t)dt=σa28ϕN28(t)+η25ϵP1(1-p)σa25ϕN25(t)+η49ϵP1(1-p)σa49ϕN49(t)+η41ϵP1(1-p)σa41ϕN41(t)-σa49ϕN49(t),+α281+α28ϵ-1η28-1[η25σa25N25(t)+η49σa49N49(t)+η41σa41N41(t)]ϕ,dN40(t)dt=α49σa491+α49ϕN49(t)-σa40ϕ(t)N40(t),dN41(t)dt=σa40ϕN40(t)-σa41ϕN41(t),
where ϕ(t) is the average thermal neutron flux of the core (we have considered ϕ(t) is independent of time and equal to a constant), and σaαβ(α=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 N28 and N25 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.Table 2
Masses and atom densities of each fuel isotopes at the time of reactor startup.
MaterialMass, kgConcentration,N(×1024), atoms/cm3UO2932460.0217391235U3058.0430.0007700238U84314.6200.0209645Table 3
Effective properties of nuclide for thermal neutron in 1000 MWe PWR (3000 MW thermal power).
NuclideSubscriptσa(barn)ηα235U255561.960.2398238U282.23422.34320.1907239Pu491618.21.860.5430240Pu402616.8241Pu411567.32.2230.3765ϵ is fast fission factor=1.0476.p is resonance scape probability=0.7725.P1 is fast to 238U resonance probability=0.9889.η 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 to3.5×1013neutrons/cm2·s.(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 about7300 hours.Using data given in Table3, 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.Table 4
Fuel composition in a PWR after 1 year irradiation. All values are in terms of kg.
MaterialBOCEOC235U30582632238U8431483604Pu0210Fission 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 (n,γ) 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 Nfσfϕ(t), in which Nf is U-235 atom density, σf is the effective U-235 microscopic fission cross section, and ϕ(t) 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 i-th species, yi, as the ratio of i-th atoms produced to U-235 atoms undergoing fission. Consequently, constant formation rate of the i-th nuclide per unit volume could be written as Pi=yiNfσfϕ(t).As shown in Figure4, 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.Table 5
Nuclear properties forSm decay chain.
NuclideNuclide indicationHalf-lifeAbsorption cross section (b)Direct yield fromU-235147NdN111.1 d00.0236147PmN22.62 y845.72 to148mPm, 448.23 to 148Pm0147SmN3∞274.20147mPmN442 d (7% to 148Pm, 93% to 148Sm )319640148PmN55.4 d138580149PmN653.1 h1105.60.0113149SmN7∞736350150PmN82.7 h00150SmN9∞158.380151SmN1087 y9734.50.0044152SmN11∞813.010.00281Figure 4
The fission-product chain leading to Sm-149.According to Figure4, 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:
(2.4)60147Nd:dN1dt=y1Nfσfϕ-(σ1ϕ+λ1)N1,(2.5)61147Pm:dN2dt=y2Nfσfϕ+λ1N1-(σ2ϕ+λ2)N2,(2.6)62147Sm:dN3dt=y3Nfσfϕ+λ2N2-(σ3ϕ)N3,(2.7)61148mPm:dN4dt=σ24ϕN2-(σ4ϕ+λ4)N4,(2.8)61148Pm:dN5dt=λ45N4+σ25ϕN2-(σ5ϕ+λ5)N5,(2.9)61149Pm:dN6dt=y6Nfσfϕ+σ4ϕN4+σ5ϕN5-(σ6ϕ+λ6)N6,(2.10)62149Sm:dN7dt=y7Nfσfϕ+λ6N6-σ7ϕN7,(2.11)61150Pm:dN8dt=σ6ϕN6-(σ8ϕ+λ8)N8,(2.12)62150Sm:dN9dt=σ7ϕN7+λ8N8-σ9ϕN9,(2.13)62151Sm:dN10dt=y10Nfσfϕ+σ9ϕN9-(σ10ϕ+λ10)N10,(2.14)62152Sm:dN11dt=y11Nfσfϕ+σ10ϕN10-σ11ϕN11,(2.15)92235U:dNfdt=-σfϕNf,
where in (2.7), σ24 is radiative capture cross section for 147Pm(n,γ)148mPm reaction and in a similar manner, in (2.8) σ25 is for 147Pm(n,γ)148Pm reaction. Also in (2.8) λ45 is a radioactive decay constant that 148mPm, as a result of a β decay, disintegrates to 148Pm. Moreover, yi indicates direct fission yield for i-th species from U-235 thermal fission, and finally, (2.15) is a rate equation for U-235 atom density in which σf is U-235 effective fission cross section. This equation implies that U-235 atom density decreases as an exponential function as time goes on. Nf 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.(a) Uranium-235 depletion. (b) Samarium-149 atom density variation as a function of irradiation time. (c) Atom density variation of other important members of the samarium chain.
(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 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 form(2.16)Fission⇀52135Te⇀53135I⇀54135Xe⇀55135Cs⇀56135Ba.
In addition, 54135Xe is a great neutron absorber and under the neutron flux within the reactor will be changed into 54136Xe. 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:(2.17)Te135:dTdt=yTNfσfϕ-λTT(2.18)I135:dIdt=λTT-λII(2.19)Xe135:dXdt=yXNfσfϕ+λII-σXXϕ-λXX(2.20)Cs135:dCdt=λXX-σCCϕ-λCC(2.21)Ba135:dBdt=λCC(2.22)U235:dUdt=-σaNfϕ
in which yT and yX stand for fission yields of 135Te and 135Xe, respectively. Also, λ are associated radioactive decay constants; σ are associated absorption cross-section, and explicitly σf and σa 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.Table 6
Nuclear properties of fission products of mass135.
NuclideNuclide indicationHalf lifeRadioactive decay constants-1Absorption cross-section (b)Direct yield fromU-235135TeT29 s0.023900.0609135II6.7 h2.87×10-500135XeX9.2 h2.09×10-52.64×1060.0032135CsC3×106 yr7.3×10-1517.20135BaBstable000Figure 6
Xenon135 atom density as well as power variations in the four-time interval is given.


## 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(d/dt)[X]j=[P]ij[X]j with known initial values in each case. But we have focused on the constant [P]ij 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 [P]ij 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.



---

*Source: 249162-2010-03-15.xml*