Solution of Point Reactor Neutron Kinetics Equations with Temperature Feedback by Singularly Perturbed Method

Abstract

The singularly perturbed method (SPM) is proposed to obtain the analytical solution for the delayed supercritical process of nuclear reactor with temperature feedback and small step reactivity inserted. The relation between the reactivity and time is derived. Also, the neutron density (or power) and the average density of delayed neutron precursors as the function of reactivity are presented. The variations of neutron density (or power) and temperature with time are calculated and plotted and compared with those by accurate solution and other analytical methods. It is shown that the results by the SPM are valid and accurate in the large range and the SPM is simpler than those in the previous literature.

1. Introduction

The analysis of variation of neutron density (or power) and reactivity with time under the different conditions is an important content of nuclear reactor physics or neutron kinetics [1–7]. Some important achievements on the supercritical transient with temperature feedback with big (𝜌₀>𝛽) or small (𝜌₀<𝛽) reactivity inserted have been approached through the effort of many scholars [7–12]. The studies on the delayed supercritical transient with small reactivity inserted and temperature feedback are introduced in the related literature [13–15], in which the explicit function of density (or power) and reactivity with respect to time is derived mainly with decoupling method, power prompt jump approximation, precursor prompt jump approximation, temperature prompt jump approximation [10, 16], and so forth. From the detailed analysis and comparison of the results in the early and recent literature [7, 12, 14], it is found that some results have certain limit and rather big error under the particular conditions. In present work, the variation law of power, reactivity, and precursor density with respect to time at any level of initial power is obtained by the singularly perturbed method (SPM). All the results are compared with those obtained by the numerical solution which tend to the accurate solution under very small time step size [17]. It is proved that the SPM is correct and reliable and is simpler than the analytical methods by the related literature.

2. Theoretical Derivation

The point reactor neutron kinetics equations with one group of delayed neutrons are [3, 4] 𝑑𝑛(𝑡)=𝜌(𝑡)𝑑𝑡−𝛽𝑙𝑛(𝑡)(𝑡),(𝑡)+𝜆𝐶𝑑𝐶=𝛽𝑑𝑡𝑙𝑛(𝑡)−𝜆𝐶(𝑡), (1)(2) where 𝑛 is the average neutron density, 𝑡 is the time, 𝜌 is the reactivity, 𝛽 is the total fraction of the delayed neutron, 𝑙  is the prompt neutron lifetime, 𝜆 is the radioactive decay constant of delayed neutron precursor, and 𝐶 is the average density of delayed neutron precursor. When multiplied with a certain coefficient, 𝑛 represents the power. It is assumed that the reactor has a negative temperature coefficient of reactivity 𝛼 (𝛼>0) when a small step reactivity 𝜌₀(<β) is inserted. Consider the temperature feedback, and the real reactor reactivity is 𝜌=𝜌₀−𝛼𝑇, (3) where 𝑇 is the temperature increment of the reactor, namely, 𝑇=𝑇_𝑠−𝑇₀, where 𝑇_𝑠 and 𝑇₀ are the instantaneous temperature and initial temperature, respectively. After the reactivity 𝜌₀ is inserted into the reactor the adiabatic model is still employed [3, 15]; then we have 𝑑𝑇𝑑𝑡=𝐾_𝑐𝑛(𝑡), (4) where 𝐾_𝑐 is the reciprocal of thermal capacity of reactor.

Combining (3) and (4) results in 𝑑𝜌𝑑𝑡=−𝛼𝐾_𝑐𝑛(𝑡). (5)

Substituting (2) into the derivative of (1) with respect to 𝑡 yields 𝑙𝑑²𝑛𝑑𝑡²=𝜌−𝛽𝑑𝑛𝑑𝑡+𝑛𝑑𝜌𝑑𝑡+𝜆𝛽𝑛−𝑙𝜆²𝐶. (6)

Substituting  𝜆𝐶 obtained from (1) into (6) and simplifying it yields 𝑑²𝑛𝑑𝑡²=𝜌−𝛽𝑙−𝜆𝑑𝑛+𝑛𝑑𝑡𝑙𝑑𝜌+𝑑𝑡𝜆𝜌𝑛𝑙. (7) The transient process is supposed to begin at 𝑡=0 and 𝜌(0)≠0 [15, 17], so the initial conditions can be given as 𝜌(0)=𝜌₀,  𝑛(0)=𝑛₀,  (𝑑𝑛/𝑑𝑡)|_𝑡=0=𝜌₀𝑛₀/𝑙,  (𝑑²𝑛/𝑑𝑡²)|_𝑡=0=((𝜌₀−𝛽)/𝑙)(𝜌₀𝑛₀/𝑙)−(𝛼𝐾_𝑐𝑛²₀/𝑙), where 𝑛₀ is the initial neutron density (or power).

For |(𝜌−𝛽)/𝑙|≫𝜆, 𝛽/𝑙≫𝜆 [3], (1) and (2) are stiff equations. According to the singularly perturbed method [18], solution of 𝑛(𝑡) includes the inner solution 𝑛_𝑓(𝑡) in inner part and the outer solution 𝑛_𝑙(𝑡) in outer part. In this paper both inner and outer solutions are approximated to be zero order: 𝑛(𝑡)=𝑛_𝑓(𝑡)+𝑛_𝑙(𝑡). (8) The initial conditions are 𝑛₀=𝑛_𝑓0+𝑛_𝑙0,𝑛^′(0)=𝑛^′_𝑓(0)+𝑛^′_𝑙(0)=𝜌₀𝑛₀𝑙,𝑛^′′(0)=𝑛_𝑓^′′(0)+𝑛_𝑙^′′(0)=𝜌₀−𝛽𝑙𝜌₀𝑛₀𝑙−𝛼𝐾_𝑐𝑛²₀𝑙, (9)(10)(11) where 𝑛_𝑓0 is the initial value of inner solution and 𝑛_𝑙0 is the initial value of outer solution.

Substituting (8) into (7) and (5), respectively, yields 𝑑²𝑛_𝑓+𝑛_𝑙𝑑𝑡²=𝜌−𝛽𝑙𝑑𝑛−𝜆^𝑓+𝑛_𝑙+𝑛𝑑𝑡^𝑓+𝑛_𝑙𝑙𝑑𝜌+𝑛𝑑𝑡𝜆𝜌_𝑓+𝑛_𝑙𝑙,𝑑𝜌𝑑𝑡=−𝛼𝐾_𝑐𝑛_𝑓(𝑡)+𝑛_𝑙(𝑡). (12)

In the outer part the inner solution attenuates to zero, and 𝑑𝑛_𝑓/𝑑𝑡 and 𝑑²𝑛_𝑓/𝑑𝑡² can be neglected, namely, 𝑑𝑛_𝑓/𝑑𝑡≈0 and 𝑑²𝑛_𝑓/𝑑𝑡²≈0; therefore, in the outer part (12) can be simplified as follows: 𝑑²𝑛_𝑙𝑑𝑡²=𝜌−𝛽𝑙𝑑𝑛_𝑙+𝑛𝑑𝑡𝑙^𝑙𝑑𝜌𝑑𝑡−𝜆𝑑𝑛_𝑙+𝑑𝑡𝜆𝜌𝑛_𝑙𝑙,𝑑𝜌𝑑𝑡=−𝛼𝐾_𝑐𝑛_𝑙. (13)(14)

Because 𝑛_𝑙(𝑡) varies slowly and 𝑙≈10⁻⁴ s, compared to other terms, the term on the left side of (13) can be neglected, and then we have 𝛽−𝜌+𝜆𝑙𝑑𝑛_𝑙=𝑑𝑡𝑑𝜌𝑛𝑑𝑡+𝜆𝜌_𝑙. (15)

Combining (14) and (15) results in 𝑑𝑛_𝑙=𝑑𝜌𝛼𝐾_𝑐𝑛_𝑙−𝜆𝜌𝛼𝐾_𝑐.𝛽−𝜌+𝜆𝑙 (16)

Integrating (16) subjective to the initial conditions 𝜌(0)=𝜌₀ and 𝑛_𝑙(0)=𝑛_𝑙0 yields 𝑛_𝑙=𝜆𝜌²₀−𝜌²2𝛼𝐾_𝑐+𝛽−𝜌+𝜆𝑙𝛽−𝜌₀+𝜆𝑙𝑛𝛽−𝜌+𝜆𝑙^𝑙0. (17)

Substituting (17) into (14) leads to −𝑑𝜌=𝜆𝜌𝑑𝑡²⁰−𝜌²+2𝛼𝐾_𝑐𝛽−𝜌₀+𝜆𝑙𝑛_𝑙02.𝛽−𝜌+𝜆𝑙 (18)

With the initial conditions 𝜌(0)=𝜌₀, the solution of (18) is 𝑡=𝛽+𝜆𝑙𝜆𝜌₁𝜌ln¹−𝜌𝜌₁𝜌+𝜌¹+𝜌₀𝜌₁−𝜌₀+1𝜆𝜌ln²¹−𝜌²₀𝜌²₁−𝜌², (19) where 𝜌₁=𝜌²₀+2𝛼𝐾_𝑐(𝛽−𝜌₀+𝜆𝑙)𝑛_𝑙0/𝜆.

In the inner part 𝑛_𝑙(𝑡) is assumed to be constant 𝑛_𝑙(𝑡)=𝑛_𝑙(0)=𝑛_𝑙0,   and 𝑑𝑛_𝑙(𝑡)/𝑑𝑡≈(𝑑𝑛_𝑙(𝑡)/𝑑𝑡)|_𝑡=0,  𝑑²𝑛_𝑙(𝑡)/𝑑𝑡²≈(𝑑²𝑛_𝑙(𝑡)/𝑑𝑡²)|_𝑡=0,  so (12) can be simplified as follows, respectively: 𝑑²𝑛_𝑓𝑑𝑡²+𝑛_𝑙^′′(0)=𝜌−𝛽𝑙−𝜆𝑑𝑛_𝑓𝑑𝑡+𝑛^′_𝑙(0)+𝑛_𝑓+𝑛_𝑙0𝑙𝑑𝜌,𝑑𝑡+𝜆𝜌𝑑𝜌𝑑𝑡=−𝛼𝐾_𝑐𝑛_𝑓(𝑡)+𝑛_𝑙0. (20)(21)

In addition, the temperature feedback is not fast enough to affect the reactivity, and in the inner part it is assumed that 𝜌(𝑡)=𝜌₀,  𝑑𝜌(𝑡)/𝑑𝑡≈−𝛼𝐾_𝑐𝑛₀.  From (20) we can get 𝑙𝑑²𝑛_𝑓𝑑𝑡²=𝜌₀−𝛽−𝜆𝑙𝑑𝑛_𝑓+𝑑𝑡𝜆𝜌₀−𝛼𝐾_𝑐𝑛₀𝑛_𝑓+𝑛_𝑙0+𝜌₀−𝛽−𝑙𝜆𝑛^′_𝑙(0)−𝑛_𝑙^′′(0)𝑙. (22)

The solution of (22) is 𝑛_𝑓=−𝑛_𝑙0−𝜌₀−𝛽−𝑙𝜆𝑛^′_𝑙(0)−𝑛_𝑙^′′(0)𝑙𝜆𝜌₀−𝛼𝐾_𝑐𝑛₀+𝐶₁1exp×𝜌2𝑙₀+−𝛽−𝑙𝜆𝜌₀−𝛽−𝑙𝜆²+4𝑙𝜆𝜌₀−𝛼𝐾_𝑐𝑛₀𝑡+𝐶₂1exp×𝜌2𝑙₀−−𝛽−𝑙𝜆𝜌₀−𝛽−𝑙𝜆²+4𝑙𝜆𝜌₀−𝛼𝐾_𝑐𝑛₀𝑡. (23)

The fast varying part 𝑛_𝑓(𝑡) is assumed to attenuate to zero in the inner part, so 𝐶₁(𝑡)=0 and −𝑛_𝑙0−𝜌₀−𝛽−𝑙𝜆𝑛^′_𝑙(0)−𝑛_𝑙^′′(0)𝑙𝜆𝜌₀−𝛼𝐾_𝑐𝑛₀=0. (24)

Then we can get the fast varying part 𝑛_𝑓(𝑡) in the inner part as follows: 𝑛_𝑓=𝑛_𝑓01exp×𝜌2𝑙₀−−𝛽−𝑙𝜆𝜌₀−𝛽−𝑙𝜆²+4𝑙𝜆𝜌₀−𝛼𝐾_𝑐𝑛₀𝑡. (25)

From (25) we have 𝑛^′_𝑓0=𝑛_𝑓0𝜌2𝑙⁰−−𝛽−𝑙𝜆𝜌₀−𝛽−𝑙𝜆²+4𝑙𝜆𝜌₀−𝛼𝐾_𝑐𝑛₀. (26)

Combining (9), (11), (13), (24), and (26) results in𝑛_𝑓0=2𝑙𝛼𝐾_𝑐𝑛₀−𝜆𝜌₀𝑛₀/𝜌₀−𝛽−𝑙𝜆²−2𝜌₀𝑛₀/𝜌₀−𝛽−𝑙𝜆2𝑙𝛼𝐾_𝑐𝑛₀−𝜆𝜌₀/𝜌₀−𝛽−𝑙𝜆²−1−1+4𝑙𝜆𝜌₀−𝛼𝐾_𝑐𝑛₀/𝜌₀−𝛽−𝑙𝜆²,𝑛_𝑙0=2𝜌₀𝑛₀/𝜌₀−𝛽−𝑙𝜆−𝑛₀1−1+4𝑙𝜆𝜌₀−𝛼𝐾_𝑐𝑛₀/𝜌₀−𝛽−𝑙𝜆²2𝑙𝛼𝐾_𝑐𝑛₀−𝜆𝜌₀/𝜌₀−𝛽−𝑙𝜆²−1−1+4𝑙𝜆𝜌₀−𝛼𝐾_𝑐𝑛₀/𝜌₀−𝛽−𝑙𝜆². (27)

Substituting (27) into (25) and (17) can get 𝑛_𝑓(𝑡) and 𝑛_𝑙(𝑡); then the neutron density (or power) will be obtained by (8).

Combining (1) and (2) results in 𝑑(𝑛+𝐶)=𝜌𝑑𝑡𝑙𝑛. (28)

Eliminating the time variable in (5) and (28) leads to 𝑑(𝑛+𝐶)𝜌𝑑𝜌=−𝛼𝐾_𝑐𝑙. (29)

Integrating (29) with the initial conditions 𝑛(0)=𝑛₀, 𝐶(0)=𝛽𝑛₀/𝜆𝑙 and 𝜌(0)=𝜌₀ yields 𝐶(𝑡)=𝜌²₀−𝜌²2𝛼𝐾_𝑐𝑙+𝛽1+𝑛𝜆𝑙⁰−𝑛(𝑡). (30)

Equations (8), (17), (25), (27), and (30) are the new analytical expressions derived by this paper.

3. Calculation and Analysis

The PWR with fuel ²³⁵U is taken as an example with parameters 𝛽=0.0065, 𝑙=0.0001 s, 𝜆=0.0774 1/s,  𝐾_𝑐=0.05 K/MW·s, and 𝛼=5×10⁻⁵ 1/K [8, 13]. For the reactor with the initial power 1 MW, while reactivity 𝜌₀=0.5𝛽 and 𝜌₀=0.8333𝛽 is inserted, respectively, the variations of reactivity, temperature, power with time, and power with reactivity are presented in Figures 1, 2, 3, and 4. The curves with smaller change are for 𝜌₀=0.5𝛽 and the curves with larger change are for 𝜌₀=0.8333𝛽. The solid line notes the accurate solution by the best basic function method with very small step size [17]. The results of this paper and the accurate solution as well as the temperature prompt jump (TPJ) method in the literature [10] are almost the same and hard to distinguish in Figures 1–4. The short dashed-dot line and long dashed-dot line represent the results of precursor prompt jump (PrPJ) method in the literature [9] and the small parameter (SmP) method in the literature [15], respectively. The dashed line notes the results of power prompt jump (PPJ) method in the literature [16]. The difference is caused by the approximate treatment to obey the methods of PrPJ, SmP, PPJ, and SmP. Furthermore the variation in the vicinity of prompt supercritical process is also calculated and is shown in Figures 5 and 6. The correct results cannot be obtained by the small parameter method (SmP) in the vicinity of prompt supercritical process and are not shown in Figures 5 and 6.

Figure 1 

Variation of output power with time while inserting step reactivity 𝜌0=0.5𝛽 and 𝜌0=0.8333𝛽.

Figure 2 

Variation of total reactivity with time while inserting step reactivity 𝜌0=0.5𝛽 and 𝜌0=0.8333𝛽.

Figure 3 

Variation of temperature rise of reactor with time while inserting step reactivity 𝜌0=0.5𝛽 and 𝜌0=0.8333𝛽.

Figure 4 

Variation of output power with reactivity while inserting step reactivity 𝜌0=0.5𝛽 and 𝜌0=0.8333𝛽.

Figure 5 

Variation of output power with time while inserting step reactivity 0.995𝛽.

Figure 6 

Variation of output power with time while inserting step reactivity 0.998𝛽.

(1) From Figures 1–6 it can be concluded that very good results cannot be obtained by the precursor prompt jump (PrPJ) method to calculate the delayed supercritical progress with small step reactivity and temperature feedback.

(2) For small step reactivity, the results by the small parameter (SmP) method are close to those by the power prompt jump (PPJ) method and are better than those by the precursor prompt jump (PrPJ) method, but the accuracy of results by the small parameter method decreases with the increase of the reactivity inserted. The power is negative when the small parameter method is used to calculate the transient process in the vicinity of prompt supercritical state. From Figures 1–4, it can be seen that the small parameter method is more suitable for the calculation of reactivity and temperature increase than for that of power.

(3) The results are quite precise using the power prompt jump (PPJ) method for the delayed supercritical process, but the main problem compared to the accurate solution is that some displacement exists along time axis. Furthermore it should be pointed out that each power peak value obtained by the precursor prompt jump (PrPJ) method, power prompt jump (PPJ) method, or small parameter (SmP) method is lower than that obtained by the accurate solution or singularly perturbed method (SPM) see Figure 1.

(4) From Figures 1–4 it can be also found that the temperature prompt jump method (TPJ) and the singularly perturbed method (SPM) in this paper are the two most precise methods for the delayed supercritical process with small step reactivity and temperature feedback. However from Figures 5 and 6 it can be seen that as the reactivity inserted increases to the vicinity of prompt supercritical process, the total discrepancy of power by the TPJ method is larger than that by the SPM or PPJ method, and the irrelevant phenomena that the power jumps at first and then decreases monotonously from the peak will appear in the TPJ method as shown in Figure 6.

4. Conclusions

The analytical expressions of power (or neutron density), reactivity, the precursor power (or density), and temperature increase with respect to time are derived for the delayed supercritical process with small reactivity (𝜌₀<𝛽) and temperature feedback by the singularly perturbed method. Compared with the results by the accurate solution and other methods in the literature, it is shown that the singularly perturbed method (SPM) in this paper is valid and accurate in the large range and is simpler than those in the previous literature. The method in this paper can provide a new theoretical foundation for the analysis of reactor neutron dynamics.

Acknowledgments

This research is supported by the National Natural Science Foundation of China (Project no. 11301540) and the Natural Science Foundation of Naval University of Engineering.