Abstract
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 stiffness should be solved efficiently. In order to reduce the stiffness 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 flux exponential approximation. The procedure is tested in a 2D neutron kinetic benchmark and a simplified high-temperature gas-cooled reactor-pebble bed module (HTR-PM) multiphysics problem utilizing the efficient Jacobian-free Newton–Krylov method. Results show that selecting appropriate flux approximation in the precursor integral method can improve the efficiency 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.
1. Introduction
The 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]. The 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. The analytical precursor concentrations consist of the integral of neutron flux over time and then eliminate precursor variables in the neutron flux equation. Therefore, the integral of neutron flux overtime should be accurately calculated.
The linear flux hypothesis based on Taylor expansion approximation is widely used to calculate the integral of neutron flux over time [3–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]. The 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.
The 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.
2. Mathematical Background
The 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. The 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.
2.1. Space-Dependent Neutron Kinetic Equations
The time-dependent Boltzmann transport equations and time-dependent precursor equations are used to describe the neutron kinetic problems. The equations are as follows:
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 are divided into several intervals from , and the is the time step length for step . The 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. Therefore, the fully implicit scheme is used in this work due to its advantage of the unconditional stability [10].
The 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. Therefore, the solution system contains a large number of simple ordinary differential equations. Numerous research studies discussed the numerical solution of equations (1) and (2). The other main problem of these analysis consists in the so-called 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].
2.2. Precursor Integral Method
As mentioned above, the stiffness and delayed neutron variables make the neutron kinetic problem difficult to solve. The precursor integral method was born to solve the problem. Precursor concentrations are not sensitive to the stiffness of the problem [10]. This phenomenon indicates the possibility of proposing a method to avoid the difficulty of stiffness in evaluation of and confine it to that of . Integrating equation (2) and eliminating variables in time-space neutron kinetic equations can effectively reduce stiffness in equations. The precursor concentrations for current time are easily derived from the integral:
An analytical neutron flux distribution is applied to solve the integral, and can be solved analytically by integral equation (4). Since precursor concentrations are not sensitive to the stiffness [13], is acceptable during simulation. It is noted that is actually a nonlinear term according to equation (4). The 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.
The 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).
2.2.1. Taylor Approximation of Neutron Flux
The unknown functions are commonly approximated by Taylor series. Applying expansion to neutron flux at initial of time step ,
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 :
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. The 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. The time characteristics of neutron flux are usually in exponential form. Since the neutron flux changes from the beginning of the time step , it is assumed that its solution at current time step is
The 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
The frequency used in equation (11) still has time dependence, which can be difficult to handle in the implicit method. Therefore, an approximation is used, which is the average value in time step:
Substitute equation (11) into equation (4), and integral on time :
By substituting the formula into the transient neutron diffusion equation, the neutron space-time kinetic problem can be completely solved. Because the average frequency contains variables at the current time , it should be emphasized that each iteration average frequency must be updated.
2.3. Truncation Error Analysis
In this section, the precision of flux approximation and time discretization are analyzed within the context of a neutron kinetic problem. The 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 .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. The truncation error in the expansion would become dominant and could not be omitted when time is roughly discretized.
The truncation errors for different time discretization schemes are derived. Backward Euler: Crank–Nicolson:
There 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. The 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 and actual value in the integrals. The truncation error of isand when using zero-order Taylor approximation, the truncation error is
Here, 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 first-order 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 , 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. Therefore, the focus of our work is to analyze and compare the accuracy and efficiency of two approximation methods.
2.4. 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.
The 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 in neutron flux. This 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 ( is the least negative root). At that time, is equal to . Two forms of neutron flux expression can be obtained by exponential expansion and variable derivation:
The truncation error of neutron flux is
The value of relative power is important in validation of transient neutron kinetic cases. The estimated time step is derived in which relative power error is taken into account:
For the selection of , 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 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.
2.5. JFNK Method
Implicit formulations are well-suited to maintain the asymptotic balance in complex multiple time-scale 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. The residual equations obtained from discretization of the nonlinear system are used in the Jacobian-free Newton–Krylov method. The nonlinear terms are solved using a Newton–Modified equation, which requires inverting the Jacobian system , where is so-called Jacobian matrix. The updated scheme is , and the iteration proceeds until . Here, and are an absolute and a relative tolerance, respectively. The inversion is performed iteratively using Krylov methods [15] (typically, GMRES), which could be facilitated by a Jacobian-free implementation and preconditioner. The Jacobian-free scheme uses the finite difference approach to approximate the matrix-vector production without building the Jacobian matrix explicitly. The preconditioner is based on reformulating the Jacobian system as (right preconditioning), where approximates the inversed Jacobian matrix, and it should be computational inexpensive.
3. Analysis and Computational Results
3.1. 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. The 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. The 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. The 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.
3.2. 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]. The 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 directions are in the configurations. The performance of each method is evaluated by relative power.
3.2.1. 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. The relative power is given by
Reactivity introduction is one of the most classical neutron dynamic problems. The TWIGL benchmark provides two ways of reactivity introduction. The first case simulates a step reactivity insertion in a thermal reactor without extern neutron source. In this case, . The other reactivity input model which has linear reactivity function of time (the ramp input) is considered. This ramp reactivity takes the form . 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)The results of ramp reactivity are not as accurate as those of step reactivity with the same time step length. This 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)The exponential approximation gives the most accurate result among precursor integral methods in step reactivity insertion. As mentioned above, first-order Taylor and exponential approximations generate second-order time precision, and the result of these methods shows good accuracy.
The above simulations of reactivity introduction problems test relative power changes in a period with the same time step. This suggests that a shorter time step size is necessary for accurate simulation. The 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 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. The 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 . A time step length greater than 10 ms simulates relative power error close to . Considering the efficiency and accuracy of the calculation, the time step used for subsequent calculation is .
Previous simulations have compared several different precursor integral methods using different flux approximate methods and different length of time steps. This 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. The relative power errors of two reactivity insertion cases are investigated with .
The numerical results of the simulations are shown in 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)
(b)
A computational efficiency study was also performed in these two reactivity insertion cases. The computation time which treats the precursor as independent variables is much larger than the precursor integral method. This 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. The precursor integral method with exponential flux approximation shows robustness and a relatively high computational efficiency for the numerical simulations of TWIGL neutron kinetic problems.
(a)
(b)
3.3. 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. The 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].
The 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. The flowing helium is heated from to in the reactor region, where many complex phenomena related to the neutron physics and thermal-hydraulics exist [18]. Therefore, 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.
3.3.1. Reactor Physics Model
There are many elements including pebble bed, reflector, carbon bricks, and helium flow region in the HTM reactor model. The artificial homogenization sections generated by the dedicated nuclear physics software V.S.O.P [19] are used to characterize different material composition. The neutron energy is divided into four energy groups. The 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 is solid (pebble bed area and reflector) temperature, 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. Thermal-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. These regions are characterized by different materials. The 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. Solid temperature equation: Fluid temperature equation: Mass flux equation:where are solid temperature, fluid temperature, and pressure, respectively, are mass fluxes, are solid effective thermal conductivity coefficient, fluid effective thermal conductivity coefficient, convective heat transfer coefficient, fluid heat capacity, and porosity, and are resistance factor, density, and acceleration of gravity, respectively.
The 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.
3.3.3. 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. This 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 is introduced into the HTR-PM reactor core model. The time duration of this transient is . 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 . The power approaches an almost steady value at the end of this transient, indicating the reactor does not reach the prompt criticality.
Three 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., , under a backward Euler time discretization scheme are adopted to solve this transient problem. The comparisons of the results are shown in 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 . 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 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. Therefore, 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. The relative power error of exponential expansion and first-order Taylor expansion using different time steps () 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 (), relative power error is restricted within 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.
(a)
(b)
3.3.4. Time Step Control Method
The 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 is the smallest compared to , which is due to a increase of truncation error under large time step. Therefore, the smaller the , 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., . Both CPU time and number of time step are included in this figure. For the time adapt method, the initial time step is 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 , the time step is larger because the neutron flux changes slowly then and small time step is no longer required to reduce calculation time. The results indicate that the time adapt method is a powerful tool to balance the computational cost and precision in neutron kinetic calculation.
(a)
(b)
The relative errors of power compared with the fixed time step and the time adapt method are shown in Figure 9. CPU times and number of time step in each method are also provided in the figure. The 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. This also means that the fixed time step method has not responded to the change of neutron flux. The 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. The time adapt method used here can balance computational cost and precision in neutron kinetic calculation.
4. Conclusion
The 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. This 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. The method of treating the precursor as independent variables is appended as a contrast. Nonlinear term appeared, and efficient nonlinear solution method, JFNK, is utilized. There 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. The 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 four-group neutron kinetic model in 2D cylindrical coordinate. The main remarks can be drawn as follows:(i)The traditional method of treating the delayed neutron precursor as independent variables causes expansive calculation cost compared with the precursor integral method. The 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. The comprehensive performance of exponential approximation is better.(iii)Time discretization dominates the truncation error of transient calculation. The time discretization of second-order accuracy is suggested in the calculation of neutron dynamics.(iv)The time adapt method proposed in this work is proven to be efficient in HTR-PM transient calculation. This 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.
Notations
| : | Groupwise fraction |
| : | Fission spectra |
| : | Decay constants of delayed neutron groups |
| : | Truncation error |
| : | Residual function |
| : | Jacobian matrix |
| : | Neutron velocity |
| : | Frequency coefficient |
| : | Prompt neutron yield |
| : | Scalar neutron flux |
| : | Fission cross-section |
| : | Fission cross-section |
| : | Dynamical time scale |
| : | Delayed neutron precursor concentration |
| : | Diffusion coefficient |
| : | Number of delayed neutron groups |
| : | Number of neutron energy groups. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the Chinese National Natural Science Foundation Project (11505102 and 11375099), Chinese National S&T Major Project (2018ZX06 902013), and IAEA CRP (I31020).