Exact solution of the one-dimensional neutron diffusion kinetic equation with one delayed precursor concentration in Cartesian geometry

1.   Introduction

The numerical methods ^[1] and the analytical approximations ^[2,3,4] for the solution of the kinetic neutron diffusion equations with delayed neutron precursor concentrations in the nuclear reactor have been of interest in both nuclear physics and reactor design for many safety considerations. In the last few years special attention to analytical solutions of linear and nonlinear problems by the scientific community. In addition to the mathematical elegance, the analytical solution and the explicit dependence on the parameters in their analytical expressions possess an adequate way to generate benchmark solutions to validate computational code results. Moreover, the analytical solution in some sense eliminates mathematical effort of the error evaluation required by various numerical techniques. In addition to computing effort required in numerical techniques. Obviously, The External and scattering slowing down terms are neglected due to the assumption of one energy group, but the fission and the absorption term are included in Eq (1). Moreover, we clarify the following:

Analytical solutions of the coupled three-dimensional multi-group time-dependent neutron transport equation with thermal-hydraulics feedback including generalized neutron diffusion source terms, scattering source term, absorption source term, external source term, and fission source are prohibitively difficult. So approximate methods are employed. This article is concerned with the most common approximation to the time-dependent transport equation i.e., the time dependent neutron diffusion equation.

Moreover, some simplifications to the physical problem are essential to be made to get the required analytical solution, namely, one dimensional, one delayed neutron precursor concentration with no source term as in the current article. Otherwise, numerical approximate methods are to be employed to solve more generalized and realistic reactor dynamics equations with thermal hydraulics feedback although it needs too much calculations effort to achieve the required accuracy ^[5,6]. Also, many attempts for illustrating the benchmark results for computational codes validation for the analytical solutions neutron diffusion kinetic equations can be found in reference ^[7].

In this work, keeping us in the tracking of searching analytical solution, we apply the Laplace transform technique to solve the one-dimensional neutron diffusion kinetics equation in Cartesian coordinates. The Laplace transform technique is a well-established methodology to solve analytically linear differential equations for a broad class of problems in the field of physics and engineering ^[8,9,10]. As a matter of fact, analytical solutions assure that no approximation is done along the solution derivation.

The main idea of this approach relies on the transformation of partial differential equations of the original physical problem into ordinary differential equation which can be solved exactly then by applying inverse Laplace transform we can get the solution of the original physical problem. In this work, keeping us in the line of searching analytical solutions we present an exact closed form solution for the neutron diffusion kinetic equation in Cartesian geometry. In this article we specialize the application, without losing generality, to the one-dimensional, monoenergetic diffusion kinetic equation with one delayed neutron precursor. The solution of the current article affirms that the methodology can be considered as cornerstone and straightforward manner for more realistic physical problems. This is means that later we can consider a six delayed neutron precursor and several energy groups as well as dealing with multidimensional problems. Some numerical results have been given to complete and validate the solution obtained.

2.   The physical problem

Assuming monoenergetic neutrons and only one delayed neutron concentration we can write the one-dimensional neutron diffusion equation in cartesian geometry is given by

For $t > 0$ and $0 < x < L$, subjected to the vacuum boundary condition $\phi (0, t) = \phi (L, t) = 0$ and the initial conditions:

Here $\phi (x, t)$ denotes the neutron flux, ${\phi }_{o}$ is the neutron flux at the initial time $\left(t = 0\right)$, $C\left(x, t\right)$ is the delayed neutron concentration and $, V$, $D, {{{\Sigma } }}_{a}\beta , \nu , {{{\Sigma } }}_{f}$ and $\lambda$ are the standard neutron parameters. Normally, the flux and delayed neutron precursor concentrations in space-time dynamics problems calculate their initial state values by solving the corresponding neutron diffusion steady-state model. Although in the current calculations we assume that they are simply given by the definitions initial conditions appeared in Eq (3). This is obvious because we are focused on getting a new explicit mathematical solution to the physical problem.

3.   The exact solution

For convenience let us rewrite Eqs (1) and (2) as follow:

Where for simplicity we put

The neutron flux Eq (4) can be solved exactly by the help of Laplace transform as follow. Taking Laplace transform to both sides of the partial differential Equation (4) with help of Eq (6), we get

Rearranging and making use of initial condition (3),

Again, using Laplace transform for the Precursor concentration Eq (5) with help of Eq (6) yields

Rearranging and making use of initial condition (3)

Where

Now, substituting Eq (10) into Eq (8) and rearranging we get,

Which equivalent to the known differential equation of the form,

Where

The general solution of Eq (13) is as follow

Substituting the boundary conditions appearing in Eq (3) after subjected to Laplace transformation we get the following forms for the unknown functions ${h}_{1}\left(s\right), {h}_{2}\left(s\right)$

Now, substituting unknown functions ${h}_{1}\left(s\right), \text{}{h}_{2}\left(s\right)$ from Eq (17) and, $g\left(s\right), \omega \left(s\right)$ from Eqs (14), (15) into Eq (16) yields

Now to make a possible apply for the inverse Laplace transform to Eq (18) we have to find first the poles and residues of Eq (18) as follow:

3.1.   Poles calculations of $\widetilde {\phi }(x, s)$

From Eq (18) we first have to find poles by setting,

Using definition of $\omega \left(s\right)$ from Eq (14) to find the first order poles of the Eq (19) yields ${\omega }^{2}\left(s\right) = 0$, consequently $\omega \left({s}_{1}\right) = \omega \left({s}_{2}\right) = 0$,

where

Also setting,

Consequently,

Substituting $\omega \left(s\right)$ from Eq (22) into Eq (14) and solving for s yields,

Where

3.2.   Residues calculation of $\widetilde {\phi }(x, s)$

Noting that $\omega \left({s}_{1}\right) = \omega \left({s}_{2}\right) = 0$ it can be easily checked that,

Similarly,

Now for $s = {s}_{3}ands = {s}_{4}$ $\mathrm{l}$ et us reform $\widetilde {\phi }(x, s)$ by substituting from Eqs (14) and (22) into (18) to get the following

Using the last form of $\widetilde {\phi }(x, s)$ appearing in Eq (27) and use of definition of $\omega \left(s\right)$ from Eq (22) we get

Similarly,

Where

and ${s}_{3} = {s}_{3}\left(n\right), \text{}{s}_{4} = {s}_{4}\left(n\right)$.

Applying the residue theorem method and using Eqs (25), (26), (28) and (29) we finally get,

Where

Now using the convolution theorem to find inverse Laplace transform of Eq (10) we finally get:

Where, ${s}_{3}$ and ${s}_{4}$ as defined in Eq (23).

Concluding, we must observe that the solution of the problem (1) is well determined by the Eqs (31) and (34) where all coefficients and functions appeared in these Equations are well defined in the previous discussion.

4.   Numerical results

In order to demonstrate the feasibility of methodology to handle the diffusion kinetic equation, we apply the proposed method to solve the problem with the following parameters ^[2]: $D$ = 0.96343, $V$ = 1.103497×10⁷, ${{{\Sigma } }}_{a}$ = 1.58430×10⁻², $\nu {\mathit{\Sigma } }_{f}$ = 3.33029×10⁻², L = 22.9, $\beta$ = 0.0045 and, $\lambda = 0.08$. We begin presenting the convergence of the results encountered for the kinetic equation to the proposed exact solution. In Table 1 we display the results attained by this methodology for the neutron flux. We vary the number of terms of the proposed Laplace transform method. In order to compare the number of terms needed to achieve the desired precision of the solution in literature, Table 1 presents the values of the neutron flux by increasing the number N of terms of the expanded series in both spatial and time variables at position $x = 11.45cm, t = 1sec$. From a simple inspection of the displayed results, we promptly figure out a coincidence of eight significant digits between the values for the neutron flux. Therefore, we may say that the fictitious problem solved converges to the exact solution of the diffusion kinetic equation.

From the results presented in Table 1, one observes that with 200 terms one attains the same accuracy. A comment is in order here, it seems that the current calculations is more efficient than reference ^[3] since there are only 200 terms were required to get results stable up to the eight digits. We show the numerical convergence of the results obtained for the neutron flux increasing the number of terms in the series solution up to 200 terms. Given a closer look to Table 1, we readily realize that we reach an accuracy of eight significant digits, when we make the summation of two hundred terms in the series solution. This discussion gives us confidence to affirm that the results appearing in Table 1 for the neutron flux has an accuracy of eight significant digits.

In Figures 1 and 2, we show the numerical results by plotting both neutron flux and precursor concentration flux as function of the spatial variable and time along the domain. Given a closer look to these figures, as expected, we promptly observe the asymptotic behavior of the neutron flux as time goes to infinity, we mean the time depending on solution tends to the stationary solution with the time increasing. It should be noted here that Figures 1 and 2 are actually plotted generated in the domain [0.0005, 30] for t, i.e., the behaviors of the flux and the precursor at the initial time t = 0 are not included, however, MATHEMATICA generated these figures as starting with 0, may be because MATHEMATICA considered the value 0.0005 as a zero, for the fact that 0.0005 is very close to 0. These are the reasons that flux appears as a parabolic shape in Figure 1 because the flux is, naturally, not a constant at our starting point t = 0.0005 (while the flux is only constant at t = 0 which is not actually displayed in Figure 1). It should also note that when considering the domain of t starting with 0, MATHEMATICA generated the curve of the flux as a sharp decreasing function just after t = 0 which is not physically acceptable.

Another important notice here is that at t = 0, the expressions 31 and 34 for the flux and precursor, respectively, become infinite series of functions of the spatial variable x, however, these infinite series can be easily summed to the corresponding constant initial conditions in Eq (3) (proofs can be easily accomplished using properties of Fourier series). But on the other hand, when calculating such infinite series (at t = 0) of functions of the spatial variable x using MATHEMATICA one can easily observe that the results don't give the right-hand side of the initial conditions in Eq (3), (really, some complex terms appeared in the results by MATHEMATICA). This was really an obstacle. This an additional reason behind our choice to plot the figures by avoiding the behaviors at the initial time.

Although considering the initial time as t = 0 was really an obstacle for MATHEMATICA to generate the right curves, the change of the starting point for the time to be 0.0005 instead of 0 overcome this difficulty. Accordingly, the curves of the flux and the precursor in Figures 1 and 2 appeared as acceptable from the physical point of view. The author would like to reveal that there are many similar problems in physics have the same difficulties when considering boundary value problems (BVPs) consisting of PDEs having the same types of the present ICs and BCs ^[11].

5.   Conclusions

In this work, we determined the one-dimensional, monoenergetic diffusion kinetic solutions with one delayed neutron precursor concentration in Cartesian geometry using the technique of Laplace Transform. Finally motivated by the promising results attained by the methodology, in a forthcoming paper we extend the present case to work out the diffusion kinetic equation in analytical fashion.

Besides the elegance and the analytical feature of the proposed exact analytical solution, we bolster this affirmative recalling that this approach is quite general in the sense that it means this technique can be applied in a straightforward manner to more realistic physical problems, particularly the ones considering six delayed neutron group and also multigroup model with up two hundred groups of energy or to a multi-regions slab problem, each one with its distinct and specific physical parameters.

Regarding the aptness of the current technique to generate benchmark solutions for computational codes validation. In fact, we claim that we can get exact results with a prescribed accuracy by increasing the number of terms summation in the series solution, paying attention to the expected observed oscillatory behavior of the convergence. This will open pathways to solve the global calculus of criticality of a nuclear reactor nucleus, applying continuity conditions for the flux and density of neutrons.

Conflict of interest

The author declares that there are no conflicts of interest regarding the publication of this paper.