A New Precursor Integral Method for Solving Space-Dependent Kinetic Equations in Neutronic and Thermal-Hydraulic Coupling System

The accurate prediction of the neutronic and thermal-hydraulic coupling system transient behavior is important in nuclear reactor safety analysis, where a large-scale nonlinear coupling system with strong stiﬀness should be solved eﬃciently. In order to reduce the stiﬀness and huge computational cost in the coupling system, the high-performance numerical techniques for solving delayed neutron precursor equation are a key issue. In this work, a new precursor integral method with an exponential approximation is proposed and compared with widely used Taylor approximation-based precursor integral methods. The truncation errors of exponential approximation and Taylor approximation are analyzed and compared. Moreover, a time control technique is put forward which is based on ﬂux exponential approximation. The procedure is tested in a 2D neutron kinetic benchmark and a simpliﬁed high-temperature gas-cooled reactor-pebble bed module (HTR-PM) multiphysics problem utilizing the eﬃcient Jacobian-free Newton–Krylov method. Results show that selecting appropriate ﬂux approximation in the precursor integral method can improve the eﬃciency and precision compared with the traditional method. The computation time is reduced to one-ninth in the HTR-PM model under the same accuracy when applying the exponential integral method with the time adaptive technique.


Introduction
e reactor is a system involving many physical phenomena, such as neutron kinetics and thermal-hydraulics. Specifically, the neutron flux distribution directly influences the heat release and affects the thermal-hydraulics in the reactor [1]. Hence, the accurate and efficient solution of the neutron field is the first but crucial step for the successful simulation of the whole coupling system and has attracted persistent attention in the field of nuclear reactor safety analysis [2]. e variation of neutron flux and delayed neutron precursor always has a smaller time scale than thermal-hydraulic variables but may lead to a significant change of heat release in the core. For the delayed neutron precursor, six additional equations are usually required to describe its transient behavior, leading to the significantly enhanced scale of the whole neutronic/ thermal-hydraulic coupled equation system. Since the delayed neutron precursor equations are only simple ordinary differential equations without the space operator, it could be analytically solved the delayed neutron precursor concentration by integrating their equations over time. As a result, the amounts of precursor concentration variables are eliminated, and the number of unknowns is significantly reduced. e analytical precursor concentrations consist of the integral of neutron flux over time and then eliminate precursor variables in the neutron flux equation. erefore, the integral of neutron flux overtime should be accurately calculated.
e linear flux hypothesis based on Taylor expansion approximation is widely used to calculate the integral of neutron flux over time [3][4][5][6]. A new exponential form flux is proposed in this work to calculate the integral of neutron flux over time, inspired by the ideas of the stiffness confinement method (SCM) [7] and frequency transform method [8]. e truncation errors of the proposed exponential approximation and Taylor approximation are analyzed and compared. Moreover, since the coefficient of the exponential function reflects the time scale of neutron kinetics, a time step control method is also derived to reduce calculation cost and ensure precision. In order to pursue the high computational performance, an advanced Jacobian-free Newton-Krylov (JFNK) method [9] is employed to solve space-dependent neutron kinetic diffusion equations and the neutronic/thermal-hydraulic coupled system. e remaining content is organized as follows. In Section 2, the outline of the numerical method for the space-dependent kinetic equations and precursor integral method is introduced. Exponential flux approximation as well as Taylor approximation methods are also derived in this section. Numerical results and discussions with regard to the precursor integral method are given in Section 3, and the simulation results of the exponential flux approximation method with time control are appended. Finally, conclusion remarks are provided in Section 4.

Mathematical Background
e mathematical background is presented in this section. In detail, the space-dependent neutron kinetic equations, including a neutron diffusion equation and delayed neutron equations, are described in Section 2.1. e precursor integral method with different approximated expansion methods is presented and discussed in Section 2.2. Expansion functions are employed in this section, including the Taylor series and a new exponential function. Error analysis of both expansion functions is discussed in Section 2.3. An adaptive time step method is proposed through the exponential function. In Section 2.4, an efficient fully implicit algorithm, JFNK, is introduced.

Space-Dependent Neutron Kinetic Equations.
e timedependent Boltzmann transport equations and time-dependent precursor equations are used to describe the neutron kinetic problems. e equations are as follows: zt � ∇ · D g (r, t)∇ϕ g (r, t) − Σ t,g (r, t)ϕ g (r, t) vΣ f,g′ (r, t)ϕ g , (r, t) + zC m (r, t) zt � β m G g′�1 vΣ f,g′ (r, t)ϕ g′ (r, t) − λ m C m (r, t) . (2) Equation (1) is discretized in space by using the conventional finite difference method. Compared with other truncation errors, the spatial discretization error is treated as high-order term and can be ignored.
To discretize the time derivatives in equations (1) and (2), we invoke the simple first-order approximation as follows: In the numerical simulation, the time intervals t ∈ [0, T] are divided into several intervals from t 0 � 0, and the Δt n is the time step length for step n. e selection of time step length depends on the accuracy of the time discretization scheme and the dynamic time scale of physical quantities. A suitable time step should be chosen to make a balance between accuracy and efficiency. For the explicit scheme, the time step always suffers from stability limitation. As a result, the time step is usually quite small to satisfy the stability requirement, leading to a relatively large computational cost. erefore, the fully implicit scheme is used in this work due to its advantage of the unconditional stability [10]. e number of delayed neutron variables is slightly more than the number of multigroup neutron fluxes, and occupies a large part in the solution of the whole system. As shown in equation (2), the precursor concentration is only related to the local neutron fission rate. erefore, the solution system contains a large number of simple ordinary differential equations. Numerous research studies discussed the numerical solution of equations (1) and (2). e other main problem of these analysis consists in the socalled stiffness of the system [11]. According to Chao's work [7], there is an abrupt turn in neutron flux, reflecting the stiffness characteristic due to the orders of magnitude difference between the time coefficients in equation (4). Stiff differential equations frequently arise in physical situation characterized by the existence of greatly differing time constants [12].

Precursor Integral Method.
As mentioned above, the stiffness and delayed neutron variables make the neutron kinetic problem difficult to solve. e precursor integral method was born to solve the problem. Precursor concentrations are not sensitive to the stiffness of the problem [10].
is phenomenon indicates the possibility of proposing a method to avoid the difficulty of stiffness in evaluation of C (n+1) m and confine it to that of ϕ (n+1) ( r → , g). Integrating equation (2) and eliminating variables C (n+1) m in time-space neutron kinetic equations can effectively reduce stiffness in equations. e precursor concentrations for current time (n + 1) are easily derived from the integral: 2 Science and Technology of Nuclear Installations An analytical neutron flux distribution ϕ( r → , g, t) is applied to solve the integral, and C (n+1) m can be solved analytically by integral equation (4). Since precursor concentrations are not sensitive to the stiffness [13], C (n+1) m is acceptable during simulation. It is noted that λ m C m is actually a nonlinear term according to equation (4). e key to this method is precisely calculating integral in equation (4). In general, the numerical approximation of the integral term can be realized by using expansion functions. e precursor integral method has two benefits. One is that the discretized time step can be enlarged because of the decoupling of stiffness due to the elimination of precursor variables. Although the stiffness is still present, the precursor concentration variables are eliminated which is beneficial for efficient solving. Moreover, the solving equations become nonlinear equations because of the existence of nonlinear terms in equation (4).

Taylor Approximation of Neutron Flux.
e unknown functions are commonly approximated by Taylor series. Applying expansion to neutron flux at initial of time step t n , Considering the zero-order and first-order expansions of neutron flux, equation (4) can be rewritten as Solving the integral is difficult because of the time derivate term. When applying implicit scheme, the term can be approximated by forward difference: Integrate equation (6) in time t n ≤ t ≤ t (n+1) : Precursor concentrations can be analytically acquired through neutron flux expansion. If zero-order neutron flux Taylor expansion is used in equation (5), we can obtain 2.2.2. Exponential Approximation of Neutron Flux. As mentioned above, flux approximation technique is a key to the precursor integral method. e same as the precursor integral method, precursor concentration is analytical acquired utilizing trial functions. Inspired by the idea of the SCM method, the neutron kinetic frequency is defined: Without considering the feedback effect, equations (1) and (2) are first-order linear differential equation with constant coefficients. e time characteristics of neutron flux are usually in exponential form. Since the neutron flux changes from the beginning of the time step ϕ n , it is assumed that its solution at current time step is Science and Technology of Nuclear Installations e coefficient preceding the exponential functions is the weak time-dependent part of neutron flux. Supposing this part is flat flux distribution in the time step, we have e frequency w g ( r → , t − t n ) used in equation (11) still has time dependence, which can be difficult to handle in the implicit method. erefore, an approximation is used, which is the average value in time step: Substitute equation (11) into equation (4), and integral on time t n ≤ t ≤ t (n+1) : By substituting the formula into the transient neutron diffusion equation, the neutron space-time kinetic problem can be completely solved. Because the average frequency w g ( r → ) contains variables at the current time ϕ (n+1) , it should be emphasized that each iteration average frequency w g ( r → ) must be updated.

Truncation Error Analysis.
In this section, the precision of flux approximation and time discretization are analyzed within the context of a neutron kinetic problem. e partial differential equations (PDEs) are defined in the above section, and here, we discuss the origin of truncation error. Equation (1) can be rewritten as Two types of time discretization are discussed in this work: A first-order accurate method (backward Euler): A second-order accurate method (Crank-Nicolson): A truncated Taylor series will help to reveal the appearance of time discretization errors, for example, We use the simplified notation _ ϕ ≡ (dϕ/dt).It is understood that any time integration method can fail to maintain its designed or asymptotic, convergence rate if the time step is too large. e truncation error in the expansion would become dominant and could not be omitted when time is roughly discretized. e truncation errors for different time discretization schemes are derived.
Backward Euler: ere is a nonlinear term in integral precursor equation: As mentioned earlier, the nonlinear terms are linearized by different neutron flux approximation methods. Hence, the truncation error occurs and depends on the level of linearization. When adopting the backward Euler method, the truncation error is It can be seen from equation (22) when applying flux approximation that an additional truncation error is involved.
e higher-order flux approximations permit accurate precursor concentration as well as smaller truncation error. As demonstrated above, the error is introduced by the deviation between neutron flux approximation ϕ( r → , t ′ ) and actual value ϕ( r → , t ′ ) in the integrals. e truncation error of and when using zero-order Taylor approximation, the truncation error is Here, Δt, t n stands for length of time step and last time point, respectively. Utilizing the integral characteristics of exponential function, the truncation error introduced by zero-order Taylor neutron flux approximation is After similar derivation, the truncation error of the firstorder Taylor approximation can be obtained: Replaced by exponential expansion, the actual neutron flux is expanded by Taylor series: Substituting ω by equation (13) gives Bring equations (28) and (29) into equation (31) and the truncation error can be analyzed according to above results: Two kinds of truncation errors are discussed in this section, time discretization error and flux approximation error, receptively. Each truncation error is a function of time step Δt, which means that the selection of time step will affect errors in varying degrees. For the flux approximation method, exponential and first-order Taylor expansion will provide high-order precision. erefore, the focus of our work is to analyze and compare the accuracy and efficiency of two approximation methods.

Time Step Control Method in Exponential Approximation.
Time step control method is widely used in neutron program as described in SPANDEX [14]. It is almost impossible to predict the truncation error accurately in neutron kinetic calculations, especially for the implicitly discretized diffusion equation. Within the tolerance of introduced truncation error, we proposed an adaptive time step method to reduce computational cost as much as possible. e adaptive time step method assumes an exponential time dependence in neutron flux, which is identical with exponential approximation in the precursor integral method. Applying backward Euler time discretization and fully implicit formulation in neutron kinetic equations, Equation assumes an exponential time dependence e ωt in neutron flux. is is feasible since the neutron flux can be expressed as a series of exponential functions in the point reactor neutron dynamic method. After a sufficient time, the solution will be dominated by asymptotic period T ≡ s −1 0 (s 0 is the least negative root). At that time, ω is equal to s −1 0 . Two forms of neutron flux expression can be obtained by exponential expansion and variable derivation: e truncation error ∈ of neutron flux is e value of relative power is important in validation of transient neutron kinetic cases. e estimated time step is derived in which relative power error ∈ power is taken into account: Science and Technology of Nuclear Installations For the selection of ω( r → , g), the average frequency in Section 2.2.2 is a good estimation which describes the dynamic changes in neutron flux. In order to make the estimation more conservative, it is often helpful to multiply an empirical coefficient η em in equation (35).
Equation (34) analytically constructs a roughly estimation of truncation error. Even though the truncation error is obtained through numerous approximations and based on several assumptions, it can reflect the order of magnitude about relative power error. Equation (35) gives the length of time steps given an expected truncation error. Different from traditional time control strategy which utilizes dynamic time scale and empirical functions, this new adaptive time step method could select time step on the basis of expectant precision. Meanwhile, this method could work with exponential expansion in the precursor integral method.

JFNK Method.
Implicit formulations are well-suited to maintain the asymptotic balance in complex multiple timescale problems, because all relevant terms in the partial differential equation (PDE) system can be approximated at the same time level [8]. However, implicit algorithms usually are related to nonlinear PDEs and require nonlinear solver. e residual equations F → ( x → ) � 0 obtained from discretization of the nonlinear system are used in the Jacobianfree Newton-Krylov method. e nonlinear terms are solved using a Newton-Modified equation, which requires inverting the Jacobian system Here, ∈ a and ∈ r are an absolute and a relative tolerance, respectively. e inversion is performed iteratively using Krylov methods [15] (typically, GMRES), which could be facilitated by a Jacobian-free implementation and preconditioner. e Jacobian-free scheme uses the finite difference approach to approximate the matrix-vector production without building the Jacobian matrix explicitly. e preconditioner is based on reformulating the Jacobian system as JP −1 Pδ x → � F → (right preconditioning), where P −1 approximates the inversed Jacobian matrix, and it should be computational inexpensive.

Implementation.
In order to verify the efficient of the precursor integral method and compare the different flux approximation methods in integral, especially the exponential form, we have developed the space-dependent neutron kinetic solver based on diffusion theory and the multigroup approximation. In the kinetic solver, the finite volume method is used for spatial discretization. Note that implicit scheme is applied in our work and can achieve stability with the expense of increased computational cost. e main sources of two truncation errors, i.e., truncation error in time discretization and delayed neutron integral technique, are analyzed in space-dependent neutronic kinetic calculation. In the previous section, we demonstrated that flux approximations may have a significant impact on the accuracy of a numerical scheme. e precursor integral method is tested with the above flux approximation method and the results are compared with those obtained from other methods. At the same time, the importance of two kinds of truncation errors is compared. A time control method is implemented in HTR-PM transient, and the result of CPU time is shown compared with the fixed time step method. e problems include step reactivity insertion and ramp input. In the following section, the examples are described separately.
For solving partial differential equations with nonlinear term, the preconditioned JFNK method is the core of the algorithm. An efficient physic-based preconditioning method is implemented in this solver, and we use incomplete LU factorization to obtain the inverse of preconditioning matrix.

TWIGL Benchmark.
TWIGL seed-blanket reactor kinetic benchmark is implemented as a verification calculation to confirm the performance of the precursor integral method for the space-dependent kinetic equations. Two reactivity insertion cases are included in the TWIGL benchmark problem [16]. e original TWIGL bench problem consists of a two-dimensional core, as shown in Figure 1, and material properties are shown in Table 1.
In addition, the spatial mesh widths for x and y directions are Δx � 1 cm and Δy � 1 cm in the configurations. e performance of each method is evaluated by relative power.

Reactivity Insert Cases.
Accurate calculate neutron relative power is essential in neutron kinetic problems, which needs to properly simulate the changes in neutron flux and precursor concentration. e relative power is given by Reactivity introduction is one of the most classical neutron dynamic problems. e TWIGL benchmark provides two ways of reactivity introduction. e first case simulates a step reactivity insertion in a thermal reactor without extern neutron source. In this case, a,2 (t) � a,2 (0) − 0.0035(t � 0.0s). e other reactivity input model which has linear reactivity function of time (the ramp input) is considered. is ramp reactivity takes the form a,2 (t) � a,2 (0) × (1 − 0.11667t)(t ≤ 0.2s). Numerical simulations were performed for three precursor treatments: independent variable, Taylor expansion approximation, and exponential expansion approximation. Backward Euler time discretization is applied in these case and time steps are selected as 10 ms and 20 ms. Five reference points are selected in the following tables.
From Tables 2 and 3, the results indicate the following: (1) Compared with the method which treats precursor concentrations as independent variables, the neutron precursor integral method is superior in terms of the accuracy, but an appropriate neutron flux approximation should be employed in analytical precursor concentration expression. (2) e results of ramp reactivity are not as accurate as those of step reactivity with the same time step length. is suggests that the length of time step should be changed for considering different ways of reactivity insertion.
(3) Different flux approximations in the precursor integral method introduce truncation errors of different sizes. Higher-order flux approximations decrease the error of the results by a magnitude. (4) e exponential approximation gives the most accurate result among precursor integral methods in step reactivity insertion. As mentioned above, firstorder Taylor and exponential approximations generate second-order time precision, and the result of these methods shows good accuracy.
e above simulations of reactivity introduction problems test relative power changes in a period with the same time step. is suggests that a shorter time step size is necessary for accurate simulation. e influence of time step selection on simulation accuracy is also a problem that needs to be paid attention. Considering different time step selections among simulations, a time period of 0 ∼ 50 ms in these two reactivity input cases is chosen. Because of the lack of reference solution, the result of the independent precursor variable method of minimum time step is chosen as base solution. e numerical results are shown in Tables 4 and 5.
It can be seen that all methods produce larger truncation errors as time step increases from Tables 4 and 5, which is due to the fact that the expression with truncation error contains the time step term Δt. A time step length greater than 10 ms simulates relative power error close to 10%. Considering the efficiency and accuracy of the calculation, the time step used for subsequent calculation is 10ms.
Previous simulations have compared several different precursor integral methods using different flux approximate methods and different length of time steps. is part focuses on the influence of time discretization on the calculation accuracy, and thus, the Crank-Nicolson method and backward Euler method are applied in simulations. e    Figure 2. It can be observed that the numerical result from the low-order method (backward Euler) provides large truncation error compared with the high-order method (Crank-Nicolson). Also from Figure 2, it can be observed that the error of the low-order time discretization method is relatively large in the time of reactivity introduction, due to the fact that backward Euler scheme cannot accurately represent the violent change of neutron flux. However, when neutron reaction rate is stable, the relative error of the method will be less. It is also important to point out that the higher-order time discretization scheme has the most significant effect on the reduction of relative error.
A computational efficiency study was also performed in these two reactivity insertion cases. e computation time which treats the precursor as independent variables is much larger than the precursor integral method. is is due to an overlarge solution vector dimension used in the independent variable method causes high computation cost, and this will result in too much inner iterations. Figure 3 shows that the decrease in time step causes exponential CPU time rise. Crank-Nicolson scheme does not bring too much    computational burden compared with backward Euler scheme, and this is mainly because the difference of two methods is only whether they inherit the previous time step's information. e precursor integral method with exponential flux approximation shows robustness and a relatively high computational efficiency for the numerical simulations of TWIGL neutron kinetic problems.

HTR-PM Model.
In the previous section, the precursor integral method is proved to be of higher accuracy, less computing time, and memory storage than any other methods. In addition, the first-order Taylor approximation and exponential approximation methods could still maintain a high accuracy even at a larger time step. e selection of time step depends on the dynamical time scale and characteristic of the nonlinear system. Some special numerical characteristics could occur when coupling neutron with thermal-hydraulic [1]. e high-temperature gas-cooled reactor-pebble bed module (HTR-PM), which is under construction in Shidao Bay, Shandong Province, consists of two nuclear modules, and one shared steam turbine [17]. Each module contains its own reactor core, helical coiled one-through steam generator, helium turbine, and vertical pipes. e flowing helium is heated from 250 ∘ C to 750 ∘ C in the reactor region, where many complex phenomena related to the neutron physics and thermal-hydraulics exist [18]. erefore, the present modeling mainly focuses on these issues, especially the solution of neutron and physical quantities affecting the precursors significantly. In terms of the algorithm, the JFNK method is employed to solve this multiphysics nonlinear problem. ere are many elements including pebble bed, reflector, carbon bricks, and helium flow region in the HTM reactor model. e artificial homogenization sections generated by the dedicated nuclear physics software V.S.O.P [19] are used to characterize different material composition. e neutron energy is divided into four energy groups. e reactor core is modeled as a two-dimensional region in r-z coordinate, in which the left boundary is treated as a symmetric boundary, while the remained boundaries are considered as vacuum boundaries. In the solution, the cross-sections are updated according to the temperature field from thermal-hydraulic calculation [20]: where T s is solid (pebble bed area and reflector) temperature, B 1 , B 2 , B 3 , B 4 , and B 5 are constant coefficients. Four-group diffusion equations are used to solve neutron flux. Six sets of delayed neutron precursor are considered in neutron kinetic simulation. Decay constant and delayed neutron portion are shown in Table 6.

3.3.2.
ermal-Hydraulic Model. In the thermal-hydraulic model, the HTR-PM core is divided into 7 regions, such as pebble bed, gas plenums in the top, and bottom reflectors. ese regions are characterized by different materials. e properties of these materials and some other constitutive relations for thermal-hydraulic calculation are obtained from KTA safety standards [21] and the simplified version of empirical correlations in TINTE [22]. As to the numerical technique, the finite difference method is employed to discretize the governing equations with regard to the thermal-hydraulic quantities and neutron physics quantities on the same grid.
where T s , T f , and p are solid temperature, fluid temperature, and pressure, respectively, μ r and μ z are mass fluxes, λ s , λ f , α, C p , and ε are solid effective thermal conductivity coefficient, fluid effective thermal conductivity coefficient, convective heat transfer coefficient, fluid heat capacity, and porosity, and W, ρ, and g are resistance factor, density, and acceleration of gravity, respectively. e coolant channel model in HTR-PM is simplified into one dimension because the diameter of the tube is rather small. General design parameters are shown in Table 7.

Numerical Result.
As stated previously, the solution of thermal-hydraulic field and neutron physics field in the HTR-PM reactor model is performed on the same spatial mesh shown in Figure 4.
is nonuniform mesh is constructed following the high-temperature gas-cooled reactor safety analysis program TINTE, which is verified in the design and analysis of HTR-PM.
A transient case is investigated in this paper. By amplifying the fission cross-section at the steady state, a positive reactivity ρ � 0.31595 is introduced into the HTR-PM reactor core model. e time duration of this transient is 0.8 s. Compared with the simulation results at the steady state shown in Figure 5, it is observed that this transient leads to a significant change to the solid temperature, fluid temperature field, outlet temperature of the core, and power in the core, as shown in the figure. Figure 6.
In this transient, the total power of the reactor has increased to 1.6 times within 0.8 s. e power approaches an almost steady value at the end of this transient, indicating the reactor does not reach the prompt criticality. ree methods, namely, zero-order Taylor expansion, first-order Taylor expansion, and exponential expansion are employed to approximate the flux in the precursor integral method, where two time steps, i.e., Δt � 100 ms and 10 ms, under a backward Euler time discretization scheme are adopted to solve this transient problem. e comparisons of    Table 8. In these two tables, the reference result is derived from a case which treats precursors as independent variables and uses a rather small time step Δt � 0.1 ms. From Table 8, it is observed that all methods lead to a well-agreed result to the reference solution, which is due to these three approximations of λ m C m are all accurate. However, for a larger time step, the result derived from zero-order Taylor expansion deviates a lot from the reference solution because of excessive truncation error, which is shown in Table 8. In addition, it is also found that the results from two higher-order approximations of neutron flux, viz, first-order Taylor expansion and exponential expansion are almost identical. erefore, it is concluded that these two methods are more stable and qualitatively correct than the zero-order Taylor approximation.
To evaluate the effect of time discretization on the solution, another scheme, i.e., Crank-Nicolson, is also employed. e relative power error of exponential expansion and first-order Taylor expansion using different time steps (Δt � 0.001 s, 0.01 s, 0.05 s, and 0.1 s) under both time discretization schemes is shown in Figure 7. In the figure, it is indicated that backward Euler (BE) and Crank-Nicolson (CN) schemes are both accurate at a small time step, but CN is more accurate than BE at a larger time step due to a second-order truncation error in time discretization. Even in the largest time step (Δt � 0.1 s), relative power error is restricted within 1% using the CN scheme. Moreover, it is also shown that the BE scheme introduces a large error at the beginning of the transient, but decreases rapidly with time proceeding due to the stable asymptotic period of reactor. Compared to the BE scheme, the CN scheme leads to a rather smaller relative power error at any time due to a higher precision in time discretization. In addition, no obvious difference is observed between exponential approximation and first-order Taylor approximation using the  same second-order time discretization, which means that time discretization dominates the truncation error.

Time
Step Control Method. e time step control technique stated in Section 2.4 is verified in this section. A reactivity insertion model identical to that in HTR-PM is used. In this transient, the time step is determined by the criterion shown as equation (35). Figure 8 shows the variation of estimated time step with time of different time step. It is of interest that the exponential approximation of flux is adopted here since the time adapt method is derived from this approximation.
In Figure 8, it is noted that the estimated time step increases continuously as time proceeds, which means that the neutron fluxes change dramatically at the beginning of the transient and tend to be stable gradually. To maintain the accuracy of the results, a small time step is required initially, while a large time step is also feasible after almost 0.6 s because of a relatively stable exponential variation of neutron flux. In addition, it is also observed that the estimated time step corresponding to Δt � 50 ms is the smallest compared to Δt � 1 ms and 10 ms, which is due to a ω increase of truncation error under large time step. erefore, the smaller the Δt adapt , the smaller the estimated time step. Figure 9 shows the comparison of relative power error of the time adapt method and fixed time steps, i.e., Δt � 1 ms, 10 ms, 50 ms, and 100 ms. Both CPU time and number of time step are included in this figure. For the time adapt method, the initial time step is 1ms to resolve fast change of neutron flux and other time steps are determined by equation (5). Figure 9 indicates that the time adapt method can realize a relatively low and stable error. In addition, it is found a continuous drop of error at fixed time step, which means that the fixed time step does not respond to the change of neutron flux well. After almost 0.6 s, the time step is larger because the neutron flux changes slowly then and small time step is no longer required to reduce calculation time.
e results indicate that the time adapt method is a powerful tool to balance the computational cost and precision in neutron kinetic calculation.   e initial time step of the time adapt method is 1 ms, and equation (35) is used to determine the next time step. Figure 9 indicates that the time adapt method can maintain errors at a relatively stable value. Fixed time step methods result in a continuous error drop which reflects in this figure.
is also means that the fixed time step method has not responded to the change of neutron flux. e result is that the time adapt method initially takes smaller time steps in order to resolve fast change in neutron flux. After approximately 0.6 second into the transient, the time adapt method takes longer time steps because the neutron flux is now changing more slowly and the very small time step is no longer required in order to reduce calculation time. e time adapt method used here can balance computational cost and precision in neutron kinetic calculation.

Conclusion
e work presents an exponential flux approximation technique in the precursor integral method, which is shown to be advantageous for solving space-dependent neutron kinetic equations. Because the neutron flux is approximated by exponential function, the coefficients in the function reflect the reactor asymptotic period. is means that the time step control algorithm can be based on the coefficients, and an adaptive time step method is proposed which is instructed by an predicted error. Truncation error analysis among different flux approximations used in the precursor integral method is performed in this work. e method of treating the precursor as independent variables is appended as a contrast. Nonlinear term appeared, and efficient nonlinear solution method, JFNK, is utilized.
ere are two sources leading to numerical errors in transient calculation, viz, time discretization and delayed neutron integral technique. Analytical insight on the effects of the time discretization and neutron flux approximation method is obtained. e numerical verification of the method above is carried out in TWIGL neutronic benchmark. Besides, more realistic applications are extended. HTR-PM multiphysics problem is presented with a thermal-hydraulic and a fourgroup neutron kinetic model in 2D cylindrical coordinate. e main remarks can be drawn as follows: (i) e traditional method of treating the delayed neutron precursor as independent variables causes expansive calculation cost compared with the precursor integral method. e CPU time can be reduced by more than half in using precursor integral methods. (ii) First-order Taylor approximation and exponential approximation for neutron flux can provide sufficient accuracy of calculation and the stability of numerical results in neutron kinetic problems. e comprehensive performance of exponential approximation is better. (iii) Time discretization dominates the truncation error of transient calculation. e time discretization of second-order accuracy is suggested in the calculation of neutron dynamics. (iv) e time adapt method proposed in this work is proven to be efficient in HTR-PM transient calculation. is method is based on exponential flux approximation and has the features of small error and low computational cost for solving neutron kinetic equations.
In the future work, more cases will be tested in neutronic and thermal-hydraulic coupling problem, in which precursor integral method with exponential approximation is used for solving neutron kinetic problem.