On Singular Perturbation of Neutron Point Kinetics in the Dynamic Model of a PWR Nuclear Power Plant

: This short communication makes use of the principle of singular perturbation to approximate the ordinary differential equation (ODE) of prompt neutron (in the point kinetics model) as an algebraic equation. This approximation is shown to yield a large gain in computational efﬁciency without compromising any signiﬁcant accuracy in the numerical simulation of primary coolant system dynamics in a PWR nuclear power plant. The approximate (i.e., singularly perturbed) model has been validated with a numerical solution of the original set of neutron point-kinetic and thermal–hydraulic equations. Both models use variable-step Runge–Kutta numerical integration.


Introduction
The neutron point-kinetics equation with six (or less) delayed groups has been traditionally used for modeling the dynamics of neutron power in nuclear power plants (e.g., pressurized water reactor (PWR)) [1,2]. The neutron point kinetic model is a well-known and experimentally validated model, which is extensively used in both academia and industry. In the integrated (e.g., lumped-parameter) representation of PWR plants, the point-kinetics model is coupled with thermal-hydraulic and fuel heat transfer models to simulate the steady-state and transient behaviors of the reactor and the primary coolant system. This integrated system model, which is represented by (nonlinear) ordinary differential equations (ODEs), is solved numerically because there is no analytical solution, in general. Such numerical simulations are time-consuming especially if the set of ODEs is stiff. Therefore, a challenge in the transient analysis of nuclear plants for design of (real-time) monitoring and (active) control systems is to construct a dynamic model that would be computationally efficient and yet serve the purpose at hand. Since the thermal-hydraulic and fuel heat transfer models act as low-pass filters that attenuate the high-frequency transients of prompt neutrons, the accuracy of predicted transients are not significantly affected if the transients due to prompt neutron are eliminated, provided that net deviations in the total reactivity are small (e.g., in the order of those that may occur during reactor-following/turbine-leading operations of PWR plants [3,4]). A fortiori, when the turbine load changes, due to the thermal capacitance of the pressurized water held in the primary system, which acts a low-pass filter for the secondary system, the states of the primary coolant presents slow dynamics ensuring that there is sufficient time for reactor control.
Over the last several decades, researchers have applied the concept of singular perturbation for numerical analysis of neutron-point-kinetics and thermal-hydraulic dynamics in nuclear power plants (e.g., [1,[5][6][7][8]). However, in the current state-of-the-art (e.g., [9][10][11][12][13]), different types of nuclear reactors are still being modeled by including the dynamics of prompt neutrons along with the dynamics of six (or less) delayed groups and temperature reactivity feedback from fuel and coolant. This short communication addresses the above issue by making use of the principle of singular perturbation [14,15] to approximate the differential equation of prompt neutron as an algebraic equation, which yields a large gain in computational efficiency without compromising any significant accuracy of the simulation results. The main objective here is to demonstrate to practicing engineers, with a rigorous (and yet simple) example, that this approximation of the dynamics of prompt neutrons in the point kinetics model as an algebraic equation indeed significantly improves the speed of dynamic simulation without any noticeable loss of accuracy. This information is very important for control system design in PWR and other types of commercial nuclear power plants.

Lumped Parameter Model for Simulation
The equations of point kinetics including the dynamics of prompt neutrons and N c delayed groups (1 ≤ N c ≤ 6) are modeled as: where the terms in Equations (1), (2) and others are defined in the list of abbreviations at the end of this short communication.
At a steady state (ss), having the total reactivity ρ ss = 0, dn(t) dt | ss = 0, and dc i (t) dt | ss = 0, it follows from Equation (2) that: For convenience of analysis, the relative neutron density n r (t) n(t) The heat transfer from fuel to coolant and the net heat removal from the coolant are modeled as: The state equations for the lumped fuel temperature T f and lumped coolant temperature T l are obtained as: where the average reactor power at time t is obtained in terms of its initial value as: P a (t) = P a (0)n r (t) (10) and the total reactivity due to control rod and temperature feedback from the coolant and fuel is: Under a steady-state normal operation, the reactivity is balanced to be zero by the loaded fuel, the control rod position, and the coolant Boron concentration. The coolant and fuel reactivity feedbacks are affected by deviations from the respective temperature reference points. In this short communication, the reference temperatures are set as the initial fuel and average coolant temperatures, respectively; and the control rod reactivity is set as zero before any reactivity insertion. Table 1 lists numerical values of the pertinent parameters and initial conditions, which will allow reproduction of the simulation results presented in this short communication.  Figure 1 shows the dynamics of the reactor parameters after three consecutive steps of control rod reactivity insertion and withdrawal by numerically solving Equations (4) to (11) with Runge-Kutta-Fehlberg method [16], which is an adaptive stepsize version of the standard (fixed stepsize) Runge-Kutta method; the variable stepsize reduces the local truncation error to the order of five. It is noted that the reactivity insertion rate is 3.07 cents per second, which is ∼0.0002 reactivity per second. The simulation code, written in Python 3, takes ∼7 s to run on a PC with 4th generation 2.2 GHz Intel Core i7 CPU and 16 GB 1600 MHz DDR3 memory; there are 19,245 steps with the tolerance value of 1 × 10 −6 for the estimated local truncation error. Profiles of relative neutron density n r , delayed neutron precursor concentration c r,i , i = 1, · · · , 6, relative average fuel temperature T f ,r and relative outlet coolant temperature T l,r with integration step size of 1 ms and 0.9 million steps.

Enhacement of Numerical Efficiency by Singularly Perturbed Point Kinetics
This section shows how the dynamics of prompt neutrons can be approximated with an algebraic equation, generated by singular perturbation [14,15], and eventually, the step size of integration is dramatically increased to enhance computational efficiency. However, before embarking on construction of singularly perturbed neutron point kinetics, it is necessary to provide a succinct mathematical background, which is presented in the next subsection.

Background of Singular Perturbation
While standard perturbation methods are generally applied to differential equations that are smoothly dependent on a small parameter ε, singular perturbations [14,15] are characterized by discontinuous dependence of system properties on ε, as seen below in the system state equations of a full-order model: where the functions ξ and η in the initial conditions smoothly depend on the parameter ε and the initial time t 0 ∈ (0, t 1 ]. If ε is sufficiently small (i.e., 0 < ε 1), then the time scales of Equations (12) and (13) are treated as slow (t) and fast (τ), respectively, because ε d dt (13), then the transients of z last for a very short time as they have very high frequency contents due to the fact that dz dt = g/ε. Having ε = 0 in Equation (13) makes the transients of z instantaneous whenever g = 0.
Setting the parameter ε = 0, which is called singular perturbation, causes an abrupt and significant change (e.g., reduction in the dimension of the state space) in dynamical behavior of the full-order system, because the fast-time system (13) degenerates to an algebraic (or transcedental) equation, which is assumed to have at least one isolated real root [14,15].
The reduced-order system is obtained by setting ε = 0 as: wherex(t) is the solution of the following equation; andz represents the quasi-steady-state of z when x =x. The above action may give rise to a discontinuity problem that can be circumvented by analysis in multiple time scales, which is the essence of singular perturbation.

Remark 2.
The solutions,x(t) andz(t) of the reduced-order system (i.e., ε = 0) in Equations (16) and (15)), respectively, are expected to be different from x(t) and z(t) in the full-order system in Equations (12) and (13)), respectively. While the unperturbed variable z may start from an arbitrary initial value at time t 0 , the quasi-steady-statez is not allowed to have an arbitrary initial value at any time t. Therefore,z(t) cannot be a uniform approximation of z(t, ε); however, if the relation (x(t, ε) − x(t)) ∼ O(ε) is guaranteed to hold uniformly on an interval [t 0 , t 1 ], then it follows that (z(t, ε) −z(t)) ∼ O(ε) over the same interval, because where a function θ : R → R belongs to the class O(ε) if lim ε→0 |θ(ε)| |ε| exists and is equal to a non-negative real number.
The concept of the so-called "boundary layer" [14,15] is now introduced to investigate the fast transients in terms of exponential stability of the equilibrium point in the state space. With a change of variable for convenience: y z − h(t, x), which shifts the quasi-steady-state of z fromz to the origin, the full-order model in Equations (12) and (13) are rewritten in the (x, y) space as: state of Equation (18) is y = 0, which is now substituted into Equation (17) to yield the reduced-order model in Equation (16). By setting ε = 0 in the fast time τ t−t 0 ε scale, Equation (18) reduces to: which is called the "boundary layer" model for singular perturbation.
During the "boundary-layer" interval [t 0 , , then the function z should asymptotically approachz provided that the following conditions, laid out in Tikhonev Theorem [14,15], are satisfied. 1.
The functions f and g in Equations (12) and (13), respectively, and their first partial derivatives with respect to (x, z, ε) and the first partial derivative of g with respect to t are continuous.

3.
The function h(t, x) in Equation (15) and the Jacobian [∂g(t, x, z, 0)/∂z] have continuous first partial derivatives with respect to their arguments.

4.
The reduced-order system in Equation (16) has a unique solutionx(t) for t ∈ [t 0 , t 1 ] within a compact subset of the solution space.

5.
The origin in the state space of Equation (19) is an exponentially stable equilibrium of the boundary-layer system.

Singularly Perturbed Neutron Point Kinetics
This subsection undertakes the task of singular perturbation of Equation (4) and examines the consequences of this approximation. To reveal the (fast-transient) time constant of neutron kinetics, a time-dependent small parameter ε is defined as: so that Equation (4) can be rewritten as: With a step reactivity insertion ρ r (0 + ) = 2 × 10 −4 by the control rod, and after a relatively long time, T (∼300 s), during which temperature reactivity feedback takes place, the reactor regain steady-state, ρ(t) = 0 ∀t > T. The values of ε(t) before the transients, at the beginning of the transients, and long after the transients are calculated using Equation (20) as: 0.0154, 0.0159, and 0.0154 respectively, as shown in Figure 2. This implies that the meanε of the variable ε(t) is always a small number and the ratio of standard deviation to mean for ε(t) is much smaller than 1, which implies that ε(t) ≈ε ∀t ∈ [t 0 , t 1 ] during the transients. The small value of ε(t) also indicates that the neutron dynamics has a small time scale. Singular perturbation (i.e., setting ε(t) = 0) approximates the ODE in Equation (21) by an algebraic equation as: (22) is a consequence of filtering out fast transients from the original dynamical system. Thus, the boundary layer (i.e., details of the fast time-scale behavior) are eliminated as seen in Figure 3. It is also seen that, beyond the boundary layer, the results of approximation are consistent with the original solution.

Remark 3. The approximation in Equation
The boundary layer model can now be set in terms of the analytical solution under the assumption that the parameters of slow-dynamics system are constant. Solving Equation (21), the following closed-form solution is obtained: Remark 4. The closed-form solution in Equation (23) is interpreted as adding a fast decay term to the algebraic Equation (22), which makes up for the missing transients. Figure 3 compares the results generated from Equations (4), (22) and (23). It is noted that Equation (23) only works under the step reactivity insertion case. An example of oscillating reactivity insertion, where the decay term may not die out asymptotically, is presented below.

Example: Sinusoidally Oscillating Reactivity Insertion
Let us consider a situation in which the position of the reactor control rod is subjected to small fluctuations around its mean value due to variations in the electrical power generation in the reactor-following/turbine-leading mode of plant operation. This situation is simulated by vibratory motion of the reactor control rod that causes approximately sinusoidal reactivity insertion and withdrawal at 1 Hz frequency with the magnitude of 2 × 10 −4 . To this end, Equations (4) and (22), which represent the original version and singularly perturbed version of prompt neutron model, respectively, are used in the simulation. The computation time and the number of integration steps for the simulation using Equation (4) are 412 milliseconds and 1661, respectively. In contrast, the simulation using Equation (22) takes computation time of ∼50 milliseconds and the number of integration steps is 197. The simulation results are shown in Figure 4, and the relative error n r −ñ r n r is plotted in Figure 5. It is calculated that the absolute value of the relative error is bounded by ∼0.31%, which is considered to be very accurate from the perspectives of numerical simulation. The number of integration steps in the simulation using the algebraic relation in Equation (22) is approximately one order of magnitude less than that using the ODE in Equation (21) or Equation (4).  for approximation of relative neutron density n r (see Equation (4) and Equation (22)).

Discussion and Conclusions
This short communication has demonstrated how the fast dynamics of prompt neutron kinetics can be singularly perturbed to yield an algebraic equation with one order of reduction in the simulation time and without any noticeable loss of accuracy (e.g., peak error in the order of ∼0.31% during transients) and then the error of approximation rapidly approaches zero. Stability of transients is guaranteed if net deviations in the total reactivity are small in the sense that the conditions of Tikhonev Theorem [14,15] are satisfied (see the end of Section 3.1).
Reducing the order of the primary coolant system model by singular perturbation of the reactor prompt kinetics equation is indeed sufficient to investigate its impact on the dynamics of the entire PWR plant. The rationale is that the primary coolant system in a PWR is thermally coupled with the secondary coolant system through the steam generator that, having a reasonably large thermal capacitance, would act as a low-pass filter. Since the primary coolant system practically filters fast transients of prompt neutrons, responses of the secondary coolant system will be even less sensitive to singular perturbation.
Author Contributions: conceptualization, A.R.; methodology, A.R., X.C.; software, X.C.; writing-original draft preparation, A.R., X.C.; supervision, A.R.; funding acquisition, A.R. All authors have read and agreed to the published version of the manuscript. These authors contributed equally to this work.