# Generalization of the Analytical Exponential Model for Homogeneous Reactor Kinetics Equations

Authors: Abdallah A. Nahla; Mohammed F. Al-Ghamdi
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

Mathematical form for two energy groups of three-dimensional homogeneous reactor kinetics equations and average one group of the precursor concentration of delayed neutrons is presented. This mathematical form is called “two energy groups of the point kinetics equations.” We rewrite two energy groups of the point kinetics equations in the matrix form. Generalization of the analytical exponential model (GAEM) is developed for solving two energy groups of the point kinetics equations. The GAEM is based on the eigenvalues and the corresponding eigenvectors of the coefficient matrix. The eigenvalues of the coefficient matrix are calculated numerically using visual FORTRAN code, based on Laguerre’s method, to calculate the roots of an algebraic equation with real coefficients. The eigenvectors of the coefficient matrix are calculated analytically. The results of the GAEM are compared with the traditional methods. These comparisons substantiate the accuracy of the results of the GAEM. In addition, the GAEM is faster than the traditional methods.


---

## Body

## 1.  Introduction

In real nuclear reactor problems that involve safety considerations, the coefficients of the space-time neutron diffusion equations depend upon the neutron flux density, precursor concentration of delayed neutrons groups, space, time, and temperature feedback. The numerical methods and the analytical approximations for solving the space-time neutron diffusion equations in the nuclear reactor have been of interest in nuclear physics and reactor design.Three-dimensional nodal mixed dual method was used for the nuclear reactor kinetics equations with improved quasistatic model and the semiimplicit scheme to solve the precursor equations [1]. A rigorous weight function in the quasistatic method was presented to solve the time-dependent multigroup diffusion equations for subcritical systems with external sources [2]. Parallelised Krylov methods were applied to the improved quasistatic approach in addition to the direct, implicit time difference, approach for solving space-time dependent multigroup neutron diffusion equations [3]. The nodal diffusion method was presented to solve space-time neutron kinetics using the finite element, primal and mixed hybrid nodal methods [4]. A one-step implicit method and a nodal modal method were studied for the solution of the time-dependent neutron diffusion equations which-based on a hexagonal spatial mesh [5–7]. An adiabatic kinetics model was developed for a boiling water reactor (BWR) core simulator AETNA [8]. Mathematical modeling of the space-time kinetics phenomena in advanced heavy water reactor, a 920 MW thermal, vertical pressure tube type thorium-based nuclear reactor, was presented using nodal modal method [9]. Class of Padé and cut-product approximations were applied to solve the two energy groups space-time nuclear reactor kinetics equations in two dimensions [10]. Adaptive Matrix Formation (AMF) method was presented and applied to homogeneous, symmetric heterogeneous, and nonsymmetric heterogeneous reactors in the cases of two and three dimensions [11]. The generalized Runge-Kutta method was developed for solving the multigroup, multidimensional, static, and transient neutron diffusion kinetics equations [12]. Computation accuracy and efficiency of a power series analytic method were presented for the time-space-dependent neutron diffusion equations with adiabatic heat up and doppler feedback [13]. The multigroup integro-differential equations of the neutron diffusion kinetics were presented and solved numerically in multislab geometry with the use of the progressive polynomial approximation [14].The different kinds of Padé approximations were applied to the nonlinear point kinetics equations including temperature feedback via analytical inversion method [15, 16]. The straightforward recurrence relation of a power series was presented to solve the point kinetics equations in the presence of delayed neutrons with Newtonian feedback [17]. The stability of generalized Runge-Kutta method of order four was studied for the stiff system point kinetics equations with temperature feedback reactivity [18, 19]. The analytical exponential model was generalized to solve the point kinetics equations of beryllium and heavy water moderated reactors with step, ramp, and temperature feedback reactivities [20]. Power series solution (PWS) was applied for solving the point kinetics equations with lumped model temperature feedback [21]. Recently, new analytical solution, based on the roots of inhour equation, the eigenvalues of the coefficient matrix, was presented and applied for the point kinetics equations with the constant, ramp, and temperature feedback reactivities [22].In this work, the two energy groups of the point kinetics equations and the average one group of delayed neutron are introduced. The generalization of the analytical exponential model is calculated using the eigenvalues and the eigenvectors of the coefficient matrix. The eigenvalues and the corresponding eigenvectors of the coefficient matrix are calculated. The results of the generalization of the analytical exponential model are compared with the results of the traditional methods.


## 2.  Mathematical Model

The two energy groups neutron diffusion equations governing the dynamic groups diffusion neutron flux and average one group of the precursor concentration of delayed neutrons behavior are written as follows [23–25]:
(2.1)1v1∂∂tΦ1(r,t)=D1∇2Φ1(r,t)-(Σa1+Σs12)Φ1(r,t)+νΣf1(1-β)Φ1(r,t)+νΣf2(1-β)Φ2(r,t)+λC(r,t),1v2∂∂tΦ2(r,t)=D2∇2Φ2(r,t)-Σa2Φ2(r,t)+Σs12Φ1(r,t),∂∂tC(r,t)=β(νΣf1Φ1(r,t)+νΣf2Φ2(r,t))-λC(r,t),
where Φ1(r,t) and Φ2(r,t) are fast and thermal neutron flux, C(r,t) is average one group of the precursor concentration of delayed neutrons, D1 and D2 are fast and thermal diffusion coefficients, Σa1 and Σa2 are fast and thermal absorption cross-sections, Σf1 and Σf2 are fast and thermal fission cross-sections, Σs12 is the scattering cross-section from fast-to-thermal neutron, ν is the neutron fission, v1 and v2 are fast and thermal neutron speed, λ is decay constant of average one group of delayed neutrons, and β is the fraction of delayed neutrons.The neutron flux and the precursor concentration of delayed neutrons can be written as separable functions of space and time as(2.2)Φ1(r,t)=v1n1(t)Ψ̂(r),Φ2(r,t)=v2n2(t)Ψ̂(r),C(r,t)=c(t)Ψ̂(r),
where n1(t) and n2(t) are the fast and thermal neutron density, c(t) is the precursor concentration density of delayed neutrons and Ψ̂(r) is the fundamental function which can be determined from the following:
(2.3)∇2Ψ̂(r)+B2Ψ̂(r)=0,
where B2 is the material buckling.Substituting (2.2) and (2.3) into (2.1) yields
(2.4)dn1(t)dt=-D1B2v1n1(t)-(Σa1+Σs12)v1n1(t)+νΣf1(1-β)v1n1(t)+νΣf2(1-β)v2n2(t)+λc(t),dn2(t)dt=-D2B2v2n2(t)-Σa2v2n2(t)+Σs12v1n1(t),dc(t)dt=β(νΣf1v1n1(t)+νΣf2v2n2(t))-λc(t).
Let us consider that l1=1/v1νΣf1 and l2=1/v2νΣf2 are fast and thermal generation time of neutrons, L12=D1/Σa1 and L22=D2/Σa2 are fast and thermal diffusion length, k1=νΣf1/Σa1[1+L12B2] and k2=νΣf2/Σa2[1+L22B2] are fast and thermal multiplication factor, ρ1=(k1-1)/k1 and ρ2=(k2-1)/k2 are fast and thermal reactivities, κ=v1Σs12, μ1=β/l1 and μ2=β/l2.Equations (2.4) are rewritten in the following form:
(2.5)dn1(t)dt=(ρ1l1-μ1-κ)n1(t)+(1l2-μ2)n2(t)+λc(t),dn2(t)dt=κn1(t)+(ρ2-1l2)n2(t),dc(t)dt=μ1n1(t)+μ2n2(t)-λc(t).
Aboanber [26] presented two-energy group two-point kinetic model of reflected reactor. So, (2.5) can be called “two energy groups of the point kinetics equations”.The initial values of this system,dn1(0)/dt=dn2(0)/dt=dc(0)/dt=0, take the following form
(2.6)n1(0)=[1-ρ2]l2κn2(0),c(0)=μ1λn1(0)+μ2λn2(0),
and the material buckling takes the form:
(2.7)B2=-Σa1+Σs12-νΣf12D1-Σa22D2+(Σa1+Σs12-νΣf12D1+Σa22D2)2-(Σa2(Σa1+Σs12-νΣf1)-Σs12νΣf2D1D2).


## 3.  Generalization of the Analytical Exponential Model

Let us rewrite (2.5) in matrix form as follows
(3.1)ddtΨ(t)=AΨ(t),
where
(3.2)Ψ(t)=[n1(t)n2(t)c(t)],A=[-αζλκ-η0μ1μ2-λ],α=μ1+κ-ρ1l1,ζ=1l2-μ2,η=1-ρ2l2.
Using the exponential integrator [27, 28], we can determine the general solution of (3.1) as follows:
(3.3)Ψ(t)=Bexp⁡(∫t0tAdt),
where B is the integral constant.This solution includes the exponential function of the matrixA which can be calculated using generalization of the analytical exponential model [20, 29] as
(3.4)Ψ(t)=∑k=13Bkexp⁡(ωkt)Uk,
where Bk are the constants, ωk are the eigenvalues of the matrix A and Uk are the corresponding eigenvectors of the matrix A.Using the initial conditionΨ(0)=Ψ0 to determine the constants Bk=VkTΨ0, (3.4) becomes as follows
(3.5)Ψ(t)=∑k=13exp⁡(ωkt)UkVkTΨ0,
where, Vk are the eigenvectors of matrix AT which corresponding to eigenvalues ωk and satisfies the normalization condition UkTVk=1. The eigenvectors Uk and Vk of the coefficient matrix A and AT, respectively, are calculated analytically in the next subsections.Equation (3.5) gives good results only for a constant matrix A. But for variable matrix A, the solution of the equation (3.1) can be takes the form
(3.6)Ψ(tm+1)=∑k=13exp⁡(hωk)UkVkTΨ(tm),
where, the values of the eigenvalues ωk and the corresponding eigenvectors Uk and Vk are calculated at time tm+(h/2), and h=tm+1-tm is the time step interval [tm,tm+1].Equations (3.5) and (3.6) represent the solutions of the two energy groups of point kinetics equation with average one group of delayed neutrons. These solutions represent a generalization of the analytical exponential model for step and variable coefficient matrix A, respectively.


### 3.1.  Eigenvalues of the Matrix A

The eigenvalues of the matrixA are determined from the following algebraic equation:
(3.7)ω3+ω2(α+η+λ)+ω[αη+αλ+ηλ-κζ-λμ1]+αηλ-κζλ-κλμ2-ηλμ1=0.
The equation (3.7) is called “inhour equation.” This equation has a 3rd degree algebraic equation. To determine the roots of (3.7), we use the visual FORTRAN code, based on Laguerre’s numerical method, to calculate the roots of an algebraic equation with real coefficients.


### 3.2.  Eigenvectors of the Matrix A

The eigenvectors of the matrixA are determined analytically as follows:
(3.8)Uk=[1κ(ωk+η)μ1(ωk+η)+κμ2(ωk+λ)(ωk+η)].


### 3.3.  Eigenvectors of the Matrix AT

The eigenvectors of the matrixAT are determined analytically as follows:
(3.9)Vk=11+(κζ(ωk+λ)+κλμ2)/(ωk+λ)(ωk+η)2+(λμ1(ωk+η)+κλμ2)/(ωk+λ)2(ωk+η)×[1ζ(ωk+λ)+λμ2(ωk+λ)(ωk+η)λ(ωk+λ)],
where (1+(κζ(ωk+λ)+κλμ2)/(ωk+λ)(ωk+η)2+(λμ1(ωk+η)+κλμ2)/(ωk+λ)2(ωk+η))-1 is the normalization factor.


## 3.1.  Eigenvalues of the Matrix A

The eigenvalues of the matrixA are determined from the following algebraic equation:
(3.7)ω3+ω2(α+η+λ)+ω[αη+αλ+ηλ-κζ-λμ1]+αηλ-κζλ-κλμ2-ηλμ1=0.
The equation (3.7) is called “inhour equation.” This equation has a 3rd degree algebraic equation. To determine the roots of (3.7), we use the visual FORTRAN code, based on Laguerre’s numerical method, to calculate the roots of an algebraic equation with real coefficients.


## 3.2.  Eigenvectors of the Matrix A

The eigenvectors of the matrixA are determined analytically as follows:
(3.8)Uk=[1κ(ωk+η)μ1(ωk+η)+κμ2(ωk+λ)(ωk+η)].


## 3.3.  Eigenvectors of the Matrix AT

The eigenvectors of the matrixAT are determined analytically as follows:
(3.9)Vk=11+(κζ(ωk+λ)+κλμ2)/(ωk+λ)(ωk+η)2+(λμ1(ωk+η)+κλμ2)/(ωk+λ)2(ωk+η)×[1ζ(ωk+λ)+λμ2(ωk+λ)(ωk+η)λ(ωk+λ)],
where (1+(κζ(ωk+λ)+κλμ2)/(ωk+λ)(ωk+η)2+(λμ1(ωk+η)+κλμ2)/(ωk+λ)2(ωk+η))-1 is the normalization factor.


## 4. Results and Discussions

To check the accuracy of the results of the mathematical model of two energy groups of the homogeneous reactor kinetics equations, it is applied to the three-dimensional homogeneous reactor. This reactor is a bare homogeneous cube of side length 200 cm on each side with two neutron groups and one precursor group. The parameters of the three-dimensional homogeneous reactor are taken from [11] and shown in Table 1. The mathematical model of two energy groups of the point kinetics equations are solved by the generalization of the analytical exponential model (GAEM). Two different types of problems; step and ramp reactivities, are applied.Table 1
Parameters of three-dimensional homogeneous reactor [11].
ParameterGroup 1Group 2Diffusion coefficientD(cm)1.351.08Absorption cross-sectionΣa(cm-1)0.0013820.0054869Neutron fissionν (neutron)2.412.41Fission cross-sectionΣf(cm-1)0.0002420.00408Scattering cross-sectionΣsg,g+1(cm-1)0.00230.0Velocityv(cm/s)3.0×1072.2×105Decay constantλ(s-1)0.08Fraction of delayed neutronβ0.0064Keff0.8952858


### 4.1.  Step Reactivity

The thermal neutron flux at the center point of the two energy groups of the homogeneous reactor kinetics equations with a positive step reactivity,ΔΣa2=-0.0000369, is calculated in Table 2. The results of the GAEM for the two energy groups of the point kinetics equations have an accuracy comparable with the results of the adaptive matrix formation (AMF) [11], 3DKIN [30] and the exact solution for the space-time nuclear reactor kinetics equations.Table 2
Thermal neutron flux at the center point of three-dimensional homogeneous reactor with a positive step reactivity.
Time (s)3DKINAMFGAEMExacth=0.001h=0.0001h=0.001h=0.00010.00.8160.8163560.8163560.8163560.8160.051.1241.1272641.1273201.1273001.1270.101.4061.4075631.4075871.4075381.4070.151.6601.6607981.6607821.6606951.6600.201.8901.8901791.8901121.8899791.8900.302.2892.2883962.2881892.2879492.2880.402.6222.6212012.6207962.6204322.620The fast and thermal neutron flux of the two energy groups of the point kinetics equations with a negative step reactivity,ΔΣa2=+0.0000369, is calculated in Table 3. This table substantiates the accuracy of the results of the GAEM compared with the results of the AMF [11].Table 3
Neutron flux at the center point of three-dimensional homogeneous reactor with a negative step reactivity.
Time (s)AMFGAEMFastThermalFastThermal0.02.2289290.8163562.2289290.8163560.051.6454810.6010481.6455130.6010600.101.4094370.5139691.4095140.5139980.151.3132170.4784791.3132830.4785040.201.2732140.4637311.2732610.4637480.301.2475080.4542681.2475190.4542720.401.2396790.4514001.2396610.451394Figure1. shows the fast and thermal neutron flux at three values for positive step reactivities, ΔΣa2=-0.0000369,-2(0.0000369) and -3(0.0000369). And, the fast and thermal neutron flux at three types of negative step reactivities, ΔΣa2=+0.0000369,+2(0.0000369) and +3(0.0000369), are shown in Figure 2.Figure 1
Fast and thermal neutron flux at a positive step reactivity.Figure 2
Fast and thermal neutron flux at a negative step reactivity.The calculation times, CPU times, of the GAEM to the AMF method are 0.43 s: 0.79 s at the same time step using personal desktop computer (processor: Core 2 Duo, 2.93 GHz). This means that the GAEM for two energy groups of the point kinetics equations is faster than the AMF for the space-time nuclear reactor kinetics equations.


### 4.2.  Ramp Reactivity

To substantiate the accuracy and the efficiency of the generalization of the analytical exponential model for solving two energy groups of the point kinetics equations with linear ramp reactivity, it is applied to three dimensional homogeneous reactor. The fast and thermal neutron flux for a positive ramp reactivity,ΔΣa2=-(0.0000369/0.4)t, are calculated in Table 4. This table is demonstrates the comparison between the results of the GAEM and AMF.Table 4
Neutron flux at the center point of three dimensional homogeneous reactor with a positive ramp reactivity.
Time (s)AMFGAEMFastThermalFastThermal0.02.2289290.8163562.2289290.8163560.052.2767180.8339952.2767120.8339930.102.3990850.8791512.3990480.8791370.152.5800620.9459312.5799720.9458980.202.8166161.0332172.8164561.0331580.303.4841461.2795223.4837731.2793850.404.5298381.6653544.5290851.665078The fast and thermal neutron flux of the two energy groups of the point kinetics equations with a negative ramp reactivity,ΔΣa2=+(0.0000369/0.4)t, is calculated in Table 5.Table 5
Neutron flux at the center point of three dimensional homogeneous reactor with a negative ramp reactivity.
Time (s)AMFGAEMFastThermalFastThermal0.02.2289290.8163562.2289290.8163560.052.1822040.7991092.1821700.7990970.102.0729540.7587952.0728920.7587730.151.9366150.7084881.9365280.7084560.201.7955650.6564461.7954600.6564070.301.5411710.5625891.5410330.5625390.401.3405550.4885821.3403960.488524The fast and thermal neutron flux at different types of positive ramp reactivities,ΔΣa2=-(0.0000369t),-2(0.0000369t) and -3(0.0000369t), is shown in Figure 3. And also, the fast and thermal neutron flux at different types of negative ramp reactivities, ΔΣa2=+(0.0000369t),+2(0.0000369t) and +3(0.0000369t), are shown in Figure 4.Figure 3
Fast and thermal neutron flux at a positive ramp reactivity.Figure 4
Fast and thermal neutron flux at a negative ramp reactivity.For ramp reactivity, The calculation times, CPU times, of the GAEM to the AMF method are 0.45 s: 0.81 s. This means, the GAEM is faster than the AMF.


## 4.1.  Step Reactivity

The thermal neutron flux at the center point of the two energy groups of the homogeneous reactor kinetics equations with a positive step reactivity,ΔΣa2=-0.0000369, is calculated in Table 2. The results of the GAEM for the two energy groups of the point kinetics equations have an accuracy comparable with the results of the adaptive matrix formation (AMF) [11], 3DKIN [30] and the exact solution for the space-time nuclear reactor kinetics equations.Table 2
Thermal neutron flux at the center point of three-dimensional homogeneous reactor with a positive step reactivity.
Time (s)3DKINAMFGAEMExacth=0.001h=0.0001h=0.001h=0.00010.00.8160.8163560.8163560.8163560.8160.051.1241.1272641.1273201.1273001.1270.101.4061.4075631.4075871.4075381.4070.151.6601.6607981.6607821.6606951.6600.201.8901.8901791.8901121.8899791.8900.302.2892.2883962.2881892.2879492.2880.402.6222.6212012.6207962.6204322.620The fast and thermal neutron flux of the two energy groups of the point kinetics equations with a negative step reactivity,ΔΣa2=+0.0000369, is calculated in Table 3. This table substantiates the accuracy of the results of the GAEM compared with the results of the AMF [11].Table 3
Neutron flux at the center point of three-dimensional homogeneous reactor with a negative step reactivity.
Time (s)AMFGAEMFastThermalFastThermal0.02.2289290.8163562.2289290.8163560.051.6454810.6010481.6455130.6010600.101.4094370.5139691.4095140.5139980.151.3132170.4784791.3132830.4785040.201.2732140.4637311.2732610.4637480.301.2475080.4542681.2475190.4542720.401.2396790.4514001.2396610.451394Figure1. shows the fast and thermal neutron flux at three values for positive step reactivities, ΔΣa2=-0.0000369,-2(0.0000369) and -3(0.0000369). And, the fast and thermal neutron flux at three types of negative step reactivities, ΔΣa2=+0.0000369,+2(0.0000369) and +3(0.0000369), are shown in Figure 2.Figure 1
Fast and thermal neutron flux at a positive step reactivity.Figure 2
Fast and thermal neutron flux at a negative step reactivity.The calculation times, CPU times, of the GAEM to the AMF method are 0.43 s: 0.79 s at the same time step using personal desktop computer (processor: Core 2 Duo, 2.93 GHz). This means that the GAEM for two energy groups of the point kinetics equations is faster than the AMF for the space-time nuclear reactor kinetics equations.


## 4.2.  Ramp Reactivity

To substantiate the accuracy and the efficiency of the generalization of the analytical exponential model for solving two energy groups of the point kinetics equations with linear ramp reactivity, it is applied to three dimensional homogeneous reactor. The fast and thermal neutron flux for a positive ramp reactivity,ΔΣa2=-(0.0000369/0.4)t, are calculated in Table 4. This table is demonstrates the comparison between the results of the GAEM and AMF.Table 4
Neutron flux at the center point of three dimensional homogeneous reactor with a positive ramp reactivity.
Time (s)AMFGAEMFastThermalFastThermal0.02.2289290.8163562.2289290.8163560.052.2767180.8339952.2767120.8339930.102.3990850.8791512.3990480.8791370.152.5800620.9459312.5799720.9458980.202.8166161.0332172.8164561.0331580.303.4841461.2795223.4837731.2793850.404.5298381.6653544.5290851.665078The fast and thermal neutron flux of the two energy groups of the point kinetics equations with a negative ramp reactivity,ΔΣa2=+(0.0000369/0.4)t, is calculated in Table 5.Table 5
Neutron flux at the center point of three dimensional homogeneous reactor with a negative ramp reactivity.
Time (s)AMFGAEMFastThermalFastThermal0.02.2289290.8163562.2289290.8163560.052.1822040.7991092.1821700.7990970.102.0729540.7587952.0728920.7587730.151.9366150.7084881.9365280.7084560.201.7955650.6564461.7954600.6564070.301.5411710.5625891.5410330.5625390.401.3405550.4885821.3403960.488524The fast and thermal neutron flux at different types of positive ramp reactivities,ΔΣa2=-(0.0000369t),-2(0.0000369t) and -3(0.0000369t), is shown in Figure 3. And also, the fast and thermal neutron flux at different types of negative ramp reactivities, ΔΣa2=+(0.0000369t),+2(0.0000369t) and +3(0.0000369t), are shown in Figure 4.Figure 3
Fast and thermal neutron flux at a positive ramp reactivity.Figure 4
Fast and thermal neutron flux at a negative ramp reactivity.For ramp reactivity, The calculation times, CPU times, of the GAEM to the AMF method are 0.45 s: 0.81 s. This means, the GAEM is faster than the AMF.


## 5.  Conclusions

Two energy groups of the point kinetics equations are represented a mathematical form for three dimensional homogeneous reactor kinetics equations and average one group of precursor delayed neutrons. The mathematical form was rewritten in the matrix form. Generalization of the analytical exponential model (GAEM) was developed to solve the new mathematical matrix form for three dimensional homogeneous reactor. This GAEM based on the eigenvalues and the corresponding eigenvectors of the coefficient matrix.The results of GAEM were compared with the AMF, 3DKIN, and the exact solution for three dimensional homogeneous reactor. This comparison substantiated the accuracy of the results of the GAEM for different types of reactivities. In addition, the GAEM for the mathematical form, two energy groups of the point kinetics equations, is faster than any other traditional methods. So, this mathematical form may be one of the most important models for the neutron diffusion equations in the nuclear reactors.



---

*Source: 282367-2012-05-16.xml*