Frontiers | The Advanced Multilevel Predictor-Corrector Quasi-Static Method for Pin-Resolved Neutron Kinetics Simulation (2024)

Introduction

Since the high-performance computing clusters have significant advances recently, the state-of-the-art computer simulation for nuclear reactors is three-dimensional (3D) whole-core time-dependent modeling with high-fidelity pin-resolved features. The nuclear industry relied on the simulation technique to understand many complicated processes and possibly decrease safety conservatism for design accidents, thus increasing nuclear power’s overall costs (Shen et al., 2019). Meanwhile, a significant challenge of the dramatically computational cost has happened to the direct simulation utilizing conventional 3D neutron transport techniques, the 3D complete nuclear reactor core. The real total number of numerical unknowns for a typical reactor core is much too large, approaching 10¹⁵ for steady-state simulation but significantly more for time-dependent simulation (Collins et al., 2016). As a replacement for the direct 3D reactor simulation, a viable solution uses the two-dimensional (2D)/one-dimensional (1D) method. The 2D/1D approach, which is commonly used in high-fidelity codes to solve the forward transport equation, employs two-dimensional heterogeneous transport computations in the radial direction and a lower-order transport calculation in the 1D axial direction, such as the CRX code (Cho et al., 2002), DeCART (Joo et al., 2004), nTRACER (Jung et al., 2013), MPACT (Collins et al., 2016; Kochunas et al., 2017; Shen et al., 2019), NECP-X (Liu et al., 2018), Tiger-3D (Wu, 2014), PANX (Zhang, et al., 2017a; Zhang, et al., 2017b) and PROTEUS-MOC (Zhang et al., 2019). The 2D/1D scheme has been successful in actual reactor applications in those high-fidelity codes, where the Method of Characteristics (MOC) was most often used for solving the 2D radial problem, and for the 1D axial computation, a variety of techniques are employed.

When utilizing huge time steps to reduce the number of transport options, it is challenging to maintain accuracy (Zhu et al., 2016), considering the time-dependent transient analysis’s high computational intensiveness. Among the most frequently utilized and preferred techniques to efficiently solve the time-dependent Boltzmann equation is the quasi-static method (Henry, 1958; Henry and Curlee, 1958), which was then modified as the Improved Quasi-Static (IQS) method (Ott and Meneley, 1969). The basic idea of the quasi-static method is to assume that the neutron flux can be factorized into the amplitude and shape function because the flux amplitude varies considerably faster than the flux shape. Thus, the shape function in the IQS technique may be solved using a modified time-dependent Boltzmann equation with a specified amplitude function. And for the amplitude function, the Exact Point-Kinetics Equations (EPKEs) are usually used, which are obtained by combining the time-dependent Boltzmann equation and the known shape function. The shape and amplitude functions would then be solved using a shape constraint function until the iteration converged (Zhu et al., 2016).

A fine-mesh/coarse-mesh based IQS method was provided to extend the nodal/EPKE based IQS method by introducing a coarse-mesh-wise amplitude function to replace the whole-core amplitude function, which was named Multigrid Amplitude Function (MAF) and was later implemented in a transport-based transient solver (Ban et al., 2012; Shaner et al., 2013; Tsujita et al., 2013; Tsujita et al., 2020). In addition, a factorization technique known as the Predictor-Corrector Quasi-static Method (PCQM) has lately gained popularity in addition to the IQS/MAF method (Kao and Henry, 1989). Rather than solving the shape and amplitude functions iteratively, the PCQM directly calculates the neutron flux in the predictor step and corrects the flux using the amplitude function evaluated in the corrector step, resulting in improved accuracy and computational efficiency the traditional IQS method (Caron et al., 2015). To overcome the unacceptable computing burden of the Multi-Group (MG) Coarse Mesh Finite Difference (CMFD) calculation with the fine time steps in PCQM, the MPACT team at the University of Michigan proposed a new multilevel transient solver named the Transient MultiLevel (TML) method, where the first level couples the transient solver and the MG CMFD solver, and the second level involves the coupling between the MG CMFD solver and the EPKE solver to capture the flux variation in the fine time range (Zhu et al., 2016).

Even so, the computational burden of the MG CMFD is still large and achieves almost 70% of the transient transport burden even with the TML method in the MPACT code (Zhu, 2016). To decrease the cost of the MG CMFD, a two-level generalized equivalence theory-based CMFD (gCMFD) acceleration system was developed, in which an analogous One-Group (1G) CMFD is created to efficiently accelerate the MG CMFD solutions and the 3D whole-core transport computation (Hao et al., 2018; Kang et al., 2020). Then, Shen et al. from the University of Michigan implemented the 1G CMFD method in the MPACT code as a new scheme of TML named TML-4. The TML-4 scheme reduced the total run time of the original TML scheme by at least 16%, even 47% for certain large-scale, full-core problems (Shen et al., 2021). However, the TML-4 remains the original TML’s time step structure, in which the transport level still has a significant computational expense.

The Advanced Multilevel Predictor-Corrector Quasi-static Method (AML-PCQM) scheme is proposed in this work to achieve further efficiency performance with the same accuracy.

The AML-PCQM scheme involves the coupling among four TFSP solvers, including the transient transport solver, the MG CMFD solver, the 1G CMFD solver, and the EPKE solver. For the time step structure, the AML-PCQM scheme makes 1G CMFD level replace the MG CMFD level totally at the second finest time step and expands the time steps of both transport and MG CMFD levels. With the new time step structure, the computational burdens of transport and MG CMFD can be effectively reduced, while the 1G CMFD can help maintain the overall accuracy. The TML and AML-PCQM schemes are implemented in the high-fidelity neutron transport code HNET to verify the functional performance. The HNET code is developed in C language with Message Passing Interface (MPI) parallel strategy, whose transient solver is based on a 2D/1D transport solver accelerated by the two-level time-dependent gCMFD (generalized equivalence theory-based coarse mesh finite difference) technique (Kang et al., 2020). Comparison and analysis of TML and AML-PCQM are provided using two heterogeneous single assembly problems with control rod withdrawal events.

The paper is structured as follows. Transient Methodology provides a comprehensive description of the transient solver at all four levels, i.e., the transport, the MG CMFD, the 1G CMFD, and EPKE transient formulations. The implementation of the AML-PCQM for all levels in HNET is described in Iteration Strategy of the AML-PCQM Scheme. Numerical Results demonstrates the comparison and the evaluation of the numerical results for the heterogeneous single assembly problems with control rod withdrawal events. Finally, Conclusion covers current performance and problems and a work plan for the future.

Transient Methodology

This section shows the specifics of the equations for the four levels of the transient solution. The 3D transient transport equations are presented first, followed by the 3D MG CMFD transient equations, the 3D 1G CMFD transient equations, and the EPKEs.

3D Transport Transient Fixed Source Equation

The 3D transient transport solution begins with the 3D time-dependent multi-group neutron transport equation with isotropic scattering approximation in the Cartesian coordinate system, as shown below.

As well as the delayed neutron precursors are determined as follows,

where φ is the angular flux and C_k is the delayed neutron precursor density for the delayed group k. S_F is the total fission neutron source, and its value is modified by the calculated eigenvalue in the steady-state calculation, and S_s,g is the scatter neutron source for the angle $\overrightarrow{\Omega }$ and the group g. χ_p and χ_d are the prompt and delayed neutron spectrums, respectively. β_k is the delayed neutron fraction and λ_k is decay constant of the k group delayed neutron precursor.

The time-dependent angular flow may be represented using the exponential transformation as,

where the ${\alpha }^n$ can be obtained by the power p of time t_n-1 and t_n-2 as,

Then the left term of Eq. 1A can be transformed into the following form,

Then,

For a given Δt_n at time step n, Eq. 3D may be further discretized for time using the first order Backward Euler technique (implicit method), as shown below,

where

Regarding the precursor equation,

Assuming that the fission source is linear over a time step,

Then, integrating precursor equations over a time step (t_n-1, t_n),

where

If the precursor equations are plugged into the Boltzmann transient equation with a first-order approximation of the fission source in one step, the transient fixed source equation is obtained,

where

In this case, an isotropic approximation of the time-dependent angular flux in the source terms is used for different time steps, which has been shown to be acceptable by some state-of-the-art time-dependent transport solvers (Zhu et al., 2016), and then part of the source terms can be simplified further as,

where

Substituting Eq. 9A into Eq. 8A,the functional transient fixed source equation is therefore as follows,

Eq. 10 may therefore be solved using any conventional steady-state neutron transport solver. A 2D/1D methodology is utilized in this study, in which the 2D radial MOC (method of characteristics) and 1D axial NEM (nodal expansion method) are linked with transverse leakage terms. These techniques’ concepts are described in the referenced article (Kang et al., 2020).

3D MG CMFD Transient Equation

The time-dependent diffusion equation with precursor equations for the CMFD formulation is provided in Eq. 11.

With similar approximations in the transient transport equations, the CMFD TFSP equation can be derived by incorporating the transport TFSP equation (Eq. 10) over 4π in the coarse mesh as follows,

where, c is the index of the center node and the neighbor node corresponds to the index q. n is the index of time step and the l indicates the iteration number. The ${\overline{S}}_{dc}^{n-1}$and ${\overline{S}}_{dt}^{n-1}$are delay neutron sources from the last time step and will be constant in a certain time step. The $\nu {\overline{\Sigma }}_f^n$, ${\overline{\Sigma }}_{g,g^{\prime }}^n$ and ${\overline{\chi }}^n$are the average fission, scattering cross-sections, and fission spectrum with standard definitions in nuclear reactor physics, while the ${\overline{\Sigma }}_t^n$ represents the average value in a node of the altered total cross-section. The $k_{eff}^{ss}$ is the eigenvalue of the steady-state transport calculation.

According to the CMFD methodology, the interface current from node c to node q in Eq. 11 can be quantified as follows,

where the “D-tilde” terms on the interface from node c to node q and the one from node q to node c are defined respectively as,

An albedo boundary condition is utilized for the boundary node treatment, and the “D-tilde” is provided as,

where h_c,q is the thickness of node c in the direction from node c to node q. The D_g is the standard diffusion coefficient. The albedo α has different values for each different boundary condition, e.g., 0.5 for the vacuum boundary condition and 0 for the reflective boundary condition. Additionally, $f_{c,q}^{dis}$and$f_{q,c}^{dif}$ are referred to the Nodal Discontinuity Factor (NDF) and Modified Diffusion coefficient Factor (MDF), respectively, quantified using the surface current, surface flux, and average flux information obtained from the radial MOC and axial NEM calculations. The variables’ specifics may be found in the referenced article (Xu et al., 2012; Hao et al., 2018).

Because the “D-tilde” is used to compel the interface to be the same as produced by higher-order techniques in the 2D/1D solution, the CMFD’s node average fluxes and interface currents will be equal with the transport solutions after global convergence is reached. Therefore, an analogous 3D MG CMFD TFSP linear system is built using the new definition of “D-tilde” and hom*ogenized XSs given by the radial planer MOC and the axial two-node NEM in the transport TFSP solutions, as follows,

where W, E, S, N, T and B represent the node’s west, east, north, south, top, and bottom surfaces. In this situation, the fixed source solver in this work may solve the conventional steady-state fixed source transport equation with an extra transient source component for each time step. The CMFD TFSP iteration must continue to use the eigenvalue from the steady-state computation with no updates.

3D 1G CMFD Transient Equation

Because the original multi-group CMFD has a high condition number, and the condition number becomes much larger if the Wielandt shift is applied to speed the inverse power iteration, the one-level MG CMFD linear system may converge considerably more slowly (Hao et al., 2018). However, since the 1G CMFD linear system is considerably less expensive to solve, it is advantageous to use the fission source from the 1G CMFD and to update the MG nodal scalar flux to decrease the computing load of the MG CMFD calculation (Kang et al., 2020).

The cross-sections, flux, and current information from the MG CMFD TFSP linear system create the 1G CMFD TFSP linear system. Except for the unique handling of “D-tilde,” all other terms in the 1G CMFD TFSP are derived by compressing the corresponding terms in the MG CMFD TFSP over all energy groups. After calculating the coefficients of the 1G CMFD TFSP, the 1G nodal scalar flux is computed by solving

where

The discontinuity factor and diffusion coefficient correction factor in 1G CMFD TFSP are calculated using the interface current from MG CMFD, which varies from the MG CMFD TFSP. Once the MG current at an interface is known, the 1G discontinuity factor and diffusion coefficient correction factor may be calculated, as described in the referenced article (Hao et al., 2018). The 1G CMFD TFSP linear system may therefore be built from the MG CMFD TFSP.

In addition, to solve the MG and 1G CMFD TFSP linear systems, a novel, efficient parallel RSILU preconditioned GMRES (Xu et al., 2019) solver has been developed. RSILU preconditioned GMRES may achieve excellent parallelization efficiency without multi-color ordering and has substantial benefits in reducing iterations and computational cost in parallel computing to find answers to the CMFD fixed source issue.

Exact Point Kinetics Equation

The EPKE can be obtained by integrating the 3D time-dependent diffusion equation with the adjoint MG CMFD fluxes as the weighting function, which is written as

and

where p is the core amplitude function, and ζ_k is the adjoint flux weighted precursor concentration for the delayed group k. The reactivity, delayed neutron fractions, neutron generation time, and delayed neutron constants all have the regular definition as follows,

where the cumulative spectrum of all fission neutrons is obtained by steady-state transport as following

the factor F(t) is defined as

and the total effective delayed neutron fraction is the summation of all delayed neutron groups as follows,

Then, Eq. 17 is solved using precisely the exact discretization for the transient transport equation in 3D Transport Transient Fixed Source Equation.

Iteration Strategy of the AML-PCQM Scheme

The AML-PCQM is based on the standard PCQM, while the essential idea is inspired by the TML method in MPACT. The time step setting for all levels is illustrated in Figure 1. Through this transient solver, the combination of three-level spatial grids and two-level energy grids can maintain a consistent accuracy and minimize the overall computational burden for the transient simulation.

FIGURE 1

In the AML-PCQM scheme, there are four solutions with three coupling levels, in which a new iteration strategy is applied to capture the neutron fluxes varies on each specific level. The overall flowchart is given in Figure 2; the four perpendicular blocks present the iterations strategy of the ATML algorithm. Each vertical block presents one level of the ATML solution. The most left vertical blocks represent the transient transport iteration scheme, in which the predictor angular flux shape on the sub-pin level is assumed to be accurate. Then, as shown in the second left vertical blocks, the scalar flux on the pin mesh is obtained by solving MG CMFD steps and correcting the transport solution’s pin-wise amplitude function.

FIGURE 2

In the same way, the energy shape function by MG CMFD on the pin mesh is presumed accurate, and the 1G scalar flux corrects the whole-energy-space amplitude from 1G CMFD steps. Finally, the core-wise shape function is predicted by 1G CMFD, and the core-wise amplitude is corrected by the finest EPKE steps in the most right vertical blocks. The details of the coupling scheme are introduced as follows.

The first coupling level is between the transport TFSP solution and the MG CMFD TFSP solution. According to the PCQM theory, the angular flux can be factored into an amplitude and shape function in a coarse mesh i.

where P is the amplitude function, which is flat in the coarse mesh, and ψ represents the angular flux shape function on the fine meshes.

Since the amplitude and shape function definition, an artificial constraint for the shape function is necessary. Here the integral of the spatial shape function in the single coarse mesh is required to be unity as,

After integrating the Eq. 19 over angle and space in the \coarse mesh and using the Eq. 20 as the constraint, the spatial amplitude function in the coarse mesh happens to be identical to the corresponding CMFD scalar flux. Thus, the angular flux can be corrected with the MG CMFD TFSP solution as in Eq. 21.

In Eq. 21, the term ${\phi }_g^{Predictor}(\overrightarrow{r},\overrightarrow{\Omega },t)$ is the predictor angular flux by the transport TFSP solution on the transport step $\Delta t_n$, the term ${\overline{\varphi }}_{i,g}^{Predictor}(\overrightarrow{r},t)$ is the predictor coarse-mesh scalar flux hom*ogenized from the transport TFSP solution, and the term${\overline{\varphi }}_{i,g}^{Corrector}(\overrightarrow{r},t)$is the corrector coarse-mesh scalar flux by the MG CMFD TFSP solution. The MG CMFD level is solved using the finer time step $\Delta {t^{\prime }}_n$with the linearly interpolated multi-group coefficients during $\Delta t_n$.

The second coupling level is for the MG CMFD and 1G CMFD. The multi-group coarse-mesh scalar flux is factored into an amplitude and shape function in a certain energy interval, but for the 1G CMFD the certain energy interval means the whole energy space. Then the factoring equation becomes as follows,

where ${\psi }^{\prime }$ represents the coarse-mesh flux shape in the energy space.

The normalization condition is

Then the energy amplitude function can be found identical to the corresponding 1G CMFD scalar flux. Similar to Eq. 21, the multi-group scalar flux can be corrected with the 1G CMFD TFSP solution as in Eq. 24. The 1G CMFD level is solved using the finer time step $\Delta {t^{″}}_n$with the linearly interpolated one-group coefficients during $\Delta {t^{\prime }}_n$.

The last level is the coupling of the 1G CMFD and EPKEs, and the factorization is still necessary, as shown in Eq. 25.

where $\overline{\varphi }$, p and ${\psi }^{″}$represent the 1G CMFD scalar flux, the core-wise amplitude, and the spatially dependent shape function, respectively.

The constraint here is obtained using the multi-group coarse-mesh scalar flux since the EPKE is integrated from the MG CMFD, as shown in Eq. 26.

where$<>$is the integration operator overall spatial regions and energy groups and C is a constant as in standard PCQM.

As a result, the 1G CMFD scalar flux can be corrected by Eq. 27. And the coefficients in the EPKE are linearly interpolated during the step $\Delta {t^{″}}_n$.

Numerical Results

Numerical results are presented in this section. Two heterogeneous single assembly problems are used to verify the performance of the AML-PCQM: 1) a UO₂ assembly problem based on C5G7-TD to check the ability of 1G CMFD acceleration for iterations and overall scheme, 2) a 51-group single assembly problem to analyze the performance of AML-PCQM for the control rod withdrawal event. The HNET code simulates all cases with the 2.30GHz Intel Xeon E5-2699 v3 CPU. The 2D/1D hybrid scheme is performed for transient transport in which the radial 2D calculation uses the planer MOC solver, and the axial 1D equation is solved with the NEM.

Problem 1: A UO₂ Assembly Problem Based on C5G7-TD

Problem 1 is based on a UO₂ assembly from the C5G7-TD benchmark’s TD4 exercise (Hou et al., 2017). The macro cross-sections are shown in a 7-group structure, whereas the kinetics parameters are presented in an 8-delayed group neutron representation. A ray spacing of 0.05cm with 48 azimuthal angles and a Tabuchi-Yamamoto polar quadrature (Yamamoto et al., 2007) using 3 polar angles per half-space are utilized for the planer MOC ray-tracing module. All active pin cells are divided into 32 flat source regions consisting of 3 fuel rings and 1 moderator rings with 8 azimuthal divisions, and the reflector cells use the 6 × 6 type grid. For the 3D configuration, the 2D geometry mesh is extruded with eight 10-cm-thick layers, eight 5-cm-thick layers, and two 20-cm-thick axial layers for both the top and bottom reflector regions. The vacuum boundary condition is applied to the core’s axial boundaries, while the radial boundaries are all reflective. The configuration is provided in Figure 3.

FIGURE 3

Six cases of Problem 1 are simulated for the preliminary verification of the AML-PCQM scheme. The transient event in Problem 1 is a prompt withdrawal of 24 control rods at the beginning of the 0.12s transient. The assembly is partly rodded for the initial condition, where the control rods stay at 20cm of insertion. Cases 1.1 and 1.2 are with the pure transient transport scheme with the time step of 0.2ms. All other cases are with a 5ms time step for the transient transport level. The case 1.3 is provided as the original PCQM case with transport and EPK method. Cases 1.4 and 1.5 use the TML method with a recommended ratio for each PCQM level in MPACT (Zhu et al., 2016). Cases 1.6 and 1.7 present the AML-PCQM scheme with four levels of different time steps, where the 5ms transport time step and the 0.2ms EPKE time step are identical to other mentioned cases, but the time step for MG CMFD and 1G CMFD are different.

Cases 1.1 and 1.4 aim to reference the pure transient transport and the TML scheme from MPACT. Then, the case 1.2 and 1.5 are intended to present the 1G CMFD acceleration for the MG CMFD iterations in case 1.1 and 1.4. Further, the case 1.6 and 1.7 are run to show the preliminary performance of the AML-PCQM scheme. Table 1 provides the details of these cases, where the “N_MG”, “N_1G”, and “N_PK” represent the number of MG CMFD steps per transport step, the number of 1G CMFD steps per MG CMFD step, and the number of EPKE steps per MG or 1G CMFD step, respectively.

TABLE 1

The fractional core power history is shown in Figure 4 in terms of accuracy. The relative errors of the fractional core power in cases 1.2 to 1.7 are also provided in Figure 4, where case 1.1 results are performed as the reference. Also, the reactivity history is presented in Figure 5. Finally, the Root Mean Squared Error (RMSE) of the fractional core power history, the solver run time, and the iteration numbers for each case are summarized in Table 2, where the index “MG” and “1G” represent MG CMFD and 1G CMFD respectively.

FIGURE 4

FIGURE 5

TABLE 2

For the accuracy of the fractional core power, the relative errors are all less than 0.2%, except the standard PCQM case 1.3. The reason is that the reactivity insertion caused by the prompt rod withdrawal does not change as a ramping line during the first several steps, which makes the linear interpolation of the EPK parameters lead to higher relative errors. However, the multilevel quasi-static cases present better reactivity results because the middle-level CMFD solvers capture the flux change in the pin meshes during the middle time steps.

As indicated, reference case 1.1 requires a considerable computational time of 535s for all solvers, which is about 6 times larger than other PCQM cases because of the large number of time steps for pure transport. On the other hand, in case 1.2, the MG CMFD solver’s iterations and run time effectively decrease due to the 1G CMFD acceleration for the MG CMFD solution, which has been discussed in the referenced paper (Kang et al., 2020). First, the standard PCQM results are presented as case 1.3 for the comparison to the multilevel quasi-static cases. Then, cases 1.4 and 1.5 are compared to check the ability of the 1G CMFD iteration accelerator for the TML scheme, where the CMFD solvers’ rum time and the total solvers’ time are reduced by about 87.9 and 64.7%, respectively. Finally, in the case 1.6 and 1.7, several 1G CMFD steps replace specific MG CMFD steps, which both achieve a partial reduction of the MG CMFD’s expense.

The multilevel quasi-static cases show the reasonably good performance of the reference case, which indicates that the TML and AML-PCQM applied in HNET perform superiorly in capturing the amplitude and shape functions’ evolution. Although the AML-PCQM has shown its advantage in decreasing the MG CMFD run time, it is notable that the total solver run time seems to be affected reversely between the TML cases and AML-PCQM cases. The main reason is that the energy group number of C5G7 cross-sections is not many enough for the 1G CMFD level to access an obvious advantage on efficiency. For the same reason, the MG CMFD solver’s run time becomes less than expected, making the comparison inconspicuous. Therefore, a 51-group single assembly control rod withdrawal problem is presented and discussed in the next sub-section.

Problem 2: A Single Assembly Problem With 51-Group Cross-Sections

Problem 2 is present here to further analyze the accuracy and efficiency performance of the AML-PCQM. Problem 2 is a typical 17 × 17 type fuel assembly problem. The C5G7-TD benchmark is used to design the pin cell architecture (Boyarinov et al., 2016). The macro cross-sections are given in a 51-group structure, and the kinetics parameters are also provided in an 8 delayed group neutron representation (Kim, 2016). The MOC ray spacing is fixed to 0.05cm, and 64 azimuthal angles and a Tabuchi-Yamamoto polar quadrature (Yamamoto et al., 2007) with 2 polar angles per half-space are used. The active pin cells are divided into 40 flat source regions consisting of 3 fuel rings, 1 clad ring, and 1 moderator ring with 8 azimuthal divisions, and the axial reflector cells also use the 6 × 6 type grid. The axial configuration includes 24 fuel layers of 5cm thickness and 2 axial reflector layers of 10cm thickness both at the top and bottom of the reflector areas, as shown in Figure 6. The boundary condition of Problem 2 is the same as Problem 1. The transient event in Problem 2 is a prompt withdrawal of the central control rod at the beginning of the 0.12s transient. The assembly is partly rodded for the initial condition, where the central control rod stays at 55cm of insertion.

FIGURE 6

Before further verification for the AML-PCQM, one must explain how the multilevel PCQM scheme works. In the standard PCQM, the flux distribution is separated into an amplitude and shape function, solved by transient transport and EPK, respectively. Then, the EPK replaces the transport and captures the amplitude varies on the fine time step, while the transient transport is moved to a coarse time step for the slower changes of the shape function. Therefore, the overall transient simulation has a much lower computational expense with reasonable accuracy. The first coarse time step in the TML technique becomes the medium time step, on which the transient transport is replaced again by MG CMFD. As a result, the transport time step grows coarser. As a result, the overall computational burden decreases again.

Although the MPACT team has implemented the 1G CMFD into the TML system referred to as TML-4 (Shen et al., 2021), the 1G CMFD was only used to remit the computing expense the MG CMFD level. As a result, the transient transport time step in TML-4 was unchanged, which means the largest computational burden is still not improved. However, it is completely reasonable and possible to extend the transport time step once again through 1G CMFD level between MG CMFD and EPK, which is exactly the primitive purpose of the AML-PCQM.

Cases 2.1 to 2.6 are presented to observe the influence of the transport step in TML cases. The fractional core power history and relative errors are shown in Figure 7. In Table 3, the transport time step grows coarser through cases 2.2–2.6, where the MG CMFD step and the EPK step stay unaltered to capture the flux varies. Case 2.1 with a transport step of 0.2ms is used as the reference case for accuracy. The RMSE of the fractional core power history indicates that the TML scheme in the HNET code does maintain an excellent agreement of about 0.55% to the reference case. The transient transport also requests less execution time when the time step was growing coarser, and the over computational burden is reduced, as it represents in Figure 8. Unfortunately, though, the MG CMFD solution demands a larger ratio of the transport run time, even nearly 35%, although the MG CMFD level maintains the original expense on solutions.

FIGURE 7

TABLE 3

FIGURE 8

In this circ*mstance, the 1G CMFD level can be applied to reduce the MG CMFD run time like the cases 1.5 and 1.6 in Problem 1. Figure 9 provides the fractional core power history and the relative errors for the AML-PCQM case 2.7 to 2.10. Table 4 and Figure 10 show that the 1G CMFD substitutes the original 1ms time step for the MG CMFD. As a result, the RMSE of the fractional core power history slightly increases to about 0.6%, but it is still an excellent agreement to the reference case. This is because the total burden of CMFD solvers keeps decreasing when the MG CMFD time is step getting large, and the tendency to be limited appears after the MG CMFD time step is larger than 5ms. Even so, approximately 70% of the CMFD solver’s run time can be depressed again because of the AML-PCQM scheme.

FIGURE 9

TABLE 4

FIGURE 10

The results here indicate that the 1G CMFD level makes an essential contribution to the intermediate time step in capturing variations in the whole-energy-space amplitude magnitude in front of EPKE. Also, the 1G CMFD level is more computationally efficient in predicting changes in the pin-wise amplitude function on the whole energy space. Therefore, it can be practical to minimize the MG CMFD solver’s computational expense, as shown in Table 5, and presents a considerable potential to handle the circ*mstances with the larger geometry modeling.

TABLE 5

Conclusion

A new multilevel predictor-corrector quasi-static method for pin-resolved neutron kinetics simulation named AML-PCQM is proposed to implement a scheme in the HNET code.

The transient formulation for the multi-group neutron transport equation is given first, followed by two gCMFD TFSP solutions and a summary of the EPKE scheme. Following that, the AML-PCQM algorithm is presented, in which the PCQM iteration technique is used to couple the transport/MG CMFD level, the MG CMFD/1G CMFD level, and the 1G CMFD/EPKE level. In each level, the initial TFSP is solved using a coarse time step in the predictor process, and the flux distribution is factorized to the amplitude and shape functions in the subsequent corrector process, where the predicted solution is rectified using numerous fine time steps. For example, in the transport/MG CMFD level coupling, the spatial shape functions of the angular sub-pin flux are assumed to change slowly over time, and the MG CMFD pin-wise amplitude function is calculated using a multi-step MG CMFD transient equation. The MG CMFD scalar flux calculated in its time step is then corrected by the 1G CMFD scalar flux in the second level. For the last level, the predictor 1G CMFD scalar flux is then corrected by a core-wise amplitude magnitude generated by the EPKE. Finally, two prompt rod withdrawal problems are chosen for the primary verification of the accuracy and efficiency performance of the AML-PCQM solution and to compare different quasi-static modes. The numerical results indicate that the speedup results of AML-PCQM cases reach over 15, and the errors remain less than 0.6% with the reference case of a pure transient transport solution. Further, the AML-PCQM scheme performs remarkable overall efficiency advantages compared to the TML method and shows a considerable potential to handle the larger geometry modeling circ*mstances.

In general, the preliminary numerical results for the prompt rod withdrawal problems show that the AML-PCQM scheme in HNET has successfully utilized several TFSP solvers in a multilevel quasi-static calculation framework, the good agreement of the fractional core power and reactivity to the reference cases has been achieved, and the balance between the accuracy and efficiency can be adjusted through different quasi-static mode or time step set. More verifications are required in the future, particularly for simulations with non-smooth reactivity insertion or explicitly modeled rod movement. Another ongoing research emphasis is the development of the thermal-hydraulic feedback module.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

LK: Develop the transient transport code and the AML-PCQM algorithm; Transient analysis of the benchmark problems. Writing. CH: Propose ideas to AML-PCQM scheme; Develop the AML-PCQM algorithm and verification. Writing - Review and Editing, Supervision. QZ: Verification of the AML-PCQM implementation. YX: Ideas to conduct the AML-PCQM scheme and transient transport. Provide details of the theory and advices to solve problems in the process of the research. Supervision.

Funding

This study is supported by the Chinese National Natural Science Foundation Project 12075067 and the National Key R and D Program of China (2018YFE0180900).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Boyarinov, F., Fomichenko, P., Hou, J., Kostadin, N., Jason, J., Ivanov, V., et al. (2016). Deterministic Time-dependent Neutron Transport Benchmark without Spatial hom*ogenization (C5G7-TD). Paris, France: Nuclear Energy Agency Organisation for Economic Co-operation and Development (NEA-OECD).

Caron, D., Dulla, S., and Ravetto, P. (2016). New Aspects in the Implementation of the Quasi-Static Method for the Solution of Neutron Diffusion Problems in the Framework of a Nodal Method. Ann. Nucl. Energy 87, 34–48.

Hou, J. J., Ivanov, K. N., Boyarinov, V. F., and Fomichenko, P. A. (2017). OECD/NEA Benchmark for Time-Dependent Neutron Transport Calculations Without Spatial hom*ogenization. Nucl. Eng. Des. 317, 117–189.

Joo, H. G., Cho, J. Y., Kim, K. S., Lee, C. C., Zee, S. Q., et al. (2004). Methods and Performance of a Three Dimensional Whole-Core Transport Code DeCART. Chicago, Illinois, US: PHYSOR 2004–The Physics of Fuel Cycles and Advanced Nuclear Systems.

Wu, W. (2014). Research on 2D/1D Coupled 3D Whole-Core Transport Calculations Based on Parallel Computation. Beijing: Tsinghua University. Ph.D. thesis.

Xu, Y., and Downar, T. (2012). Convergence analysis of a CMFD method based on generalized equivalence theory. Knoxville, TN: PHYSOR 2012–Advances in Reactor Physics-Linking Research, Industry, and Education.

Xu, Y., and Hao, C. (2019). A Novel and Efficient Hybrid RSILU Preconditioner for the Parallel GMRES Solution of the Coarse Mesh Finite Difference Equations for Practical Reactor Simulations. Nucl. Sci. Eng. 194 (2), 104–119.

Zhu, A. (2016). Transient Methods for Pin-Resolved Whole Core Transport. Ann Arbor: University of Michigan. Ph.D. thesis.