Implementation of high-fidelity neutronics and thermal-hydraulic coupling calculations in HNET

Implementation of high-fidelity neutronics and thermal-hydraulic coupling calculations in HNET 1 first-author https://orcid.org/0000-0001-8487-8272 1 1 corresp * haochen.heu@163.com * haochen.heu@163.com haochen.heu@163.com 2 3 Fundamental Science On Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Harbin 150001, China Fundamental Science On Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Harbin 150001, China Institute of Nuclear and New Energy Technology (INET), Tsinghua University, Beijing 100084, China Institute of Nuclear and New Energy Technology (INET), Tsinghua University, Beijing 100084, China Purdue University, West Lafayette, IN 47907, USA Purdue University, West Lafayette, IN 47907, USA To perform nuclear reactor simulations in a more realistic manner, the coupling scheme between neutronics and thermal-hydraulics was implemented in the HNET program for both steady-state and transient conditions. For simplicity, efficiency, and robustness, the matrix-free Newton/Krylov (MFNK) method was applied to the steady-state coupling calculation. In addition, the optimal perturbation size was adopted to further improve the convergence behavior of the MFNK. For the transient coupling simulation, the operator splitting method with a staggered time mesh was utilized to balance the computational cost and accuracy. Finally, VERA Problem 6 with power and boron perturbation and the NEACRP transient benchmark were simulated for analysis. The numerical results show that the MFNK method can outperform Picard iteration in terms of both efficiency and robustness for a wide range of problems. Furthermore, the reasonable agreement between the simulation results and the reference results for the NEACRP transient benchmark verifies the capability of predicting the behavior of the nuclear reactor. Coupling calculation High-fidelity neutronics Thermal-hydraulics Matrix-free Newton/Krylov method Transient simulation Introduction 1 With the growing computing power and requirement of safety analysis [ 1 1 ], there has been interest in the high-fidelity simulation of the nuclear reactor core [ 2 2 – 4 4 2–4 ]. The operational phenomena in the reactor to be simulated are very complex owing to the close interaction among various physical fields, particularly neutronics and thermal hydraulics [ 5 5 , 6 6 ]. More specifically, the power density from the neutron flux significantly affects the thermal properties, which in turn has a powerful influence on neutronics through cross sections. Hence, to achieve a realistic prediction of reactor behavior under steady-state and transient conditions, the strong bond between neutronics and thermal-hydraulics must be considered. Because the steady-state solution is generally regarded as the initial condition for transients [ 7 7 ], the first priority is to address the steady-state coupling of neutronics and thermal hydraulics. For steady-state coupling calculations, the Picard iteration has been studied and used extensively in high-fidelity simulation codes. However, convergence may be difficult and inefficient when the nonlinearity of the problem increases significantly [ 8 8 ]. Hence, some approaches have emerged to address the limitations of Picard iteration, such as relaxation, partially converging physics [ 8 8 ], and Anderson acceleration [ 9 9 ]. Alternatively, the fully implicit Jacobian-free Newton-Krylov (JFNK) method [ 10 10 ], with fast convergence and good robustness, has come into focus. However, most high-fidelity neutronics codes do not form the residuals required to implement JFNK. Therefore, the application of JFNK requires significant modifications to the algorithms of these well-developed codes. For simple implementation and good convergence performance, the MFNK method was proposed [ 11 11 ]. In contrast to the JFNK method, MFNK constructs residuals through the solutions of individual subsystems within the Picard iteration framework. Thus, MFNK can possess the advantage of implicit coupling while avoiding the reformulation of the coupled subsystems into a larger nonlinear system. After the steady-state coupling calculation, the transient coupling simulation was emphasized. Although the implicit coupling scheme mentioned above has been applied to the nodal method for transient coupling calculation [ 12 12 ], most transient coupled systems still adopt the operator splitting (OS) method to compromise between computational burden and accuracy. Numerical instability may occur because of the explicit exchange of coupled information. Some improvements to the OS method have been developed, and the two main strategies are using staggered time mesh and higher-order time discretization of the nonlinear terms [ 12 12 ]. The high-fidelity Neutron Transport Program (HNET) is a 2D/1D transport code based on multilevel acceleration theory [ 13 13 , 14 14 ]. The accuracy and efficiency of HNET have been verified against various steady-state and transient tests [ 15 15 , 16 16 ], but only for problems without thermal-hydraulic feedback. To consider the interaction between neutronics and thermal hydraulics, a simplified internal thermal-hydraulic module was developed in HNET to provide thermal-hydraulic feedback. The MFNK method was adopted for the steady-state coupling calculation. Moreover, the OS method with a staggered time mesh was applied to the transient simulation. The remainder of this paper is organized as follows. In Sect. 2, the MFNK method and the estimation of the associated optimal perturbation size are introduced. Section 3 describes the implementation of the coupling neutronics with thermal hydraulics. An assessment of the MFNK method for the VERA core physics benchmark problem 6 [ 17 17 ] and the numerical results for the NEACRP transient benchmark [ 18 18 ] are presented in Section 4. Finally, conclusions are drawn in Sect. 5. Matrix Free Newton/Krylov method 1 Derivation of MFNK method 2 The MFNK method is based on the block Gauss-Seidel (BGS)-style Picard iteration. For the nuclear reactor nonlinear system considering neutronics and thermal hydraulics, the BGS-style Picard iteration alternatively solving each physical field can be written as in Eq. (1) and Eq. (2). <math display="block" id="M1"><mrow><msup><mi>x</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msup><mo>=</mo><msub><mi>F</mi><mrow><mtext>TH</mtext></mrow></msub><mo stretchy="false">(</mo><msup><mrow><mover accent="true"><mi>y</mi><mo stretchy="true">¯</mo></mover></mrow><mi>k</mi></msup><mo stretchy="false">)</mo></mrow></math> <math display="block" id="M2"><mrow><msup><mi>y</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msup><mo>=</mo><msub><mi>F</mi><mtext>N</mtext></msub><mo stretchy="false">(</mo><msup><mrow><mover accent="true"><mi>x</mi><mo stretchy="true">¯</mo></mover></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msup><mo stretchy="false">)</mo></mrow></math> where vector x denotes the fuel temperature, coolant temperature, and coolant density in the thermal-hydraulic model, and y denotes the power density derived from the neutron fluxes. F TH and F N represent the thermal-hydraulic solver and neutronics solver, respectively. Superscript k represents the iteration number. Because the coupled problem is implicit, the overline in Eq. (1) and Eq. (2) indicates that the corresponding variables should be transformed further. In this sense, F TH and F N should be regarded as abstract solvers involving the total required operations, rather than only an actual solver. For brevity, the tilde is omitted from the expressions in the subsequent analysis. When the coupled problem converges, the solutions to Eq. (1) and Eq. (2) must be consistent, which implies that the nonlinear system can be expressed in terms of a residual form, as in Eq. (3). <math display="block" id="M3"><mrow><mi>F</mi><mo stretchy="false">(</mo><msup><mi>y</mi><mi>k</mi></msup><mo stretchy="false">)</mo><mo>=</mo><msup><mi>y</mi><mi>k</mi></msup><mo>−</mo><msub><mi>F</mi><mtext>N</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>F</mi><mrow><mtext>TH</mtext></mrow></msub><mo stretchy="false">(</mo><msup><mi>y</mi><mi>k</mi></msup><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow><mo>=</mo><mn>0</mn></mrow></math> Subsequently, Newton method was adopted to resolve the nonlinear equation. By applying Taylor’s expansion to the right-hand side of Eq. (3) around y k , and ignoring the higher-order terms yields Eq. (4). <math display="block" id="M4"><mrow><mi>F</mi><mo stretchy="false">(</mo><msup><mi>y</mi><mi>k</mi></msup><mo>+</mo><mi>δ</mi><mi>y</mi><mo stretchy="false">)</mo><mo>=</mo><mi>F</mi><mo stretchy="false">(</mo><msup><mi>y</mi><mi>k</mi></msup><mo stretchy="false">)</mo><mo>+</mo><msup><mi>F</mi><mo>′</mo></msup><mo stretchy="false">(</mo><msup><mi>y</mi><mi>k</mi></msup><mo stretchy="false">)</mo><mi>δ</mi><mi>y</mi><mo>=</mo><mn>0</mn></mrow></math> Then, Eq. (4) can be rearranged as shown in Eq. (5), which is referred to as Newton iteration. <math display="block" id="M5"><mrow><mi>J</mi><mo stretchy="false">(</mo><msup><mi>y</mi><mi>k</mi></msup><mo stretchy="false">)</mo><mi>δ</mi><mi>y</mi><mo>=</mo><mo>−</mo><mi>F</mi><mo stretchy="false">(</mo><msup><mi>y</mi><mi>k</mi></msup><mo stretchy="false">)</mo></mrow></math> where J ( y k ) represents the Jacobian matrix. The linear system Eq. (5) can be solved using the Krylov subspace method, where the Jacobian matrix is only utilized for the action of the matrix-vector product. In this scheme, the Jacobian matrix-vector product can be approximated by the finite difference approximation to avoid assembling the Jacobian matrix explicitly, as shown in Eq. (6). <math display="block" id="M6"><mtable><mtr><mtd><mi>J</mi><mrow><mo>(</mo><mrow><msup><mi>y</mi><mi>k</mi></msup></mrow><mo>)</mo></mrow><mtext>v</mtext><mo>=</mo><mfrac><mrow><mi>F</mi><mo stretchy="false">(</mo><msup><mi>y</mi><mi>k</mi></msup><mo>+</mo><mi>ε</mi><mtext>v</mtext><mo stretchy="false">)</mo><mo>−</mo><mi>F</mi><mo stretchy="false">(</mo><msup><mi>y</mi><mi>k</mi></msup><mo stretchy="false">)</mo></mrow><mi>ε</mi></mfrac></mtd></mtr><mtr><mtd><mo>=</mo><mtext>v</mtext><mo>−</mo><mfrac><mrow><msub><mi>F</mi><mtext>N</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>F</mi><mrow><mtext>TH</mtext></mrow></msub><mo stretchy="false">(</mo><msup><mi>y</mi><mi>k</mi></msup><mo>+</mo><mi>ε</mi><mtext>v</mtext><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow><mo>−</mo><msub><mi>F</mi><mi>N</mi></msub><mrow><mo>(</mo><mrow><msub><mi>F</mi><mrow><mtext>TH</mtext></mrow></msub><mo stretchy="false">(</mo><msup><mi>y</mi><mi>k</mi></msup><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow></mrow><mi>ε</mi></mfrac></mtd></mtr></mtable></math> where ε is the perturbation size and v is the basis vector of the Krylov iteration. Upon the δ y , the power density is updated: <math display="block" id="M7"><mrow><msup><mi>y</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msup><mo>=</mo><msup><mi>y</mi><mi>k</mi></msup><mo>+</mo><mi>δ</mi><mi>y</mi></mrow></math> The MFNK iteration is repeated until the convergence criteria are met. Subsequently, thermal-hydraulic results are obtained based on the converged value y . From the above analysis, it can be found that only the solution and not the residual of individual subsystems are necessary to implement MFNK. Optimal perturbation size 2 The finite-difference approximation implemented in Eq. (6) results in truncation and rounding errors, as well as convenience. The truncation error is proportional to the perturbation size, whereas the rounding error is the opposite. Therefore, a balance should be found between the truncation error and the rounding error through the selection of the perturbation size. Xu [ 11 11 ] proposed a method to estimate the optimal perturbation size. Assuming that y , v , and ε are given at the machine precision, the error e d introduced by the finite-difference approximation can be arranged into Eq. (8). <math display="block" id="M8"><mtable><mtr><mtd><msub><mi>e</mi><mtext>d</mtext></msub><mo>=</mo><mfrac><mrow><mi>F</mi><mo stretchy="false">(</mo><mi>y</mi><mo>+</mo><mi>ε</mi><mtext>v</mtext><mo stretchy="false">)</mo><mo>−</mo><mi>F</mi><mo stretchy="false">(</mo><mi>y</mi><mo stretchy="false">)</mo><mo>−</mo><mi>ε</mi><msup><mi>F</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mi>y</mi><mo stretchy="false">)</mo><mtext>v</mtext></mrow><mi>ε</mi></mfrac></mtd></mtr><mtr><mtd><mo>+</mo><mfrac><mrow><mover accent="true"><mi>F</mi><mo stretchy="true">¯</mo></mover><mo stretchy="false">(</mo><mover accent="true"><mrow><mi>y</mi><mo>+</mo><mi>ε</mi><mtext>v</mtext></mrow><mo stretchy="true">¯</mo></mover><mo stretchy="false">)</mo><mo>−</mo><mi>F</mi><mo stretchy="false">(</mo><mi>y</mi><mo>+</mo><mi>ε</mi><mtext>v</mtext><mo stretchy="false">)</mo><mo>−</mo><mover accent="true"><mi>F</mi><mo stretchy="true">¯</mo></mover><mo stretchy="false">(</mo><mi>y</mi><mo stretchy="false">)</mo><mo>+</mo><mi>F</mi><mo stretchy="false">(</mo><mi>y</mi><mo stretchy="false">)</mo></mrow><mi>ε</mi></mfrac></mtd></mtr></mtable></math> where <math id="M9"><mrow><mover accent="true"><mrow><mi>y</mi><mo>+</mo><mi>ε</mi><mtext>v</mtext></mrow><mo stretchy="true">¯</mo></mover></mrow></math> is the finite precision representation of <math id="M10"><mrow><mi>y</mi><mo>+</mo><mi>ε</mi><mtext>v</mtext></mrow></math> and <math id="M11"><mrow><mover accent="true"><mi>F</mi><mo stretchy="true">¯</mo></mover></mrow></math> is the function F evaluated for finite precision arithmetic. The terms on the right-hand side of Eq. (8) are referred to as truncation and round-off errors, respectively. After quantifying the error caused by the finite-difference approximation [ 11 11 ], the error upper bound can be obtained. Subsequently, the optimal perturbation size can be estimated by minimizing the upper error bound, as shown in Eq. (9). <math display="block" id="M12"><mrow><msub><mi>ε</mi><mrow><mtext>opt</mtext></mrow></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mrow><mo>‖</mo><mtext>v</mtext><mo>‖</mo></mrow></mrow></mfrac><msqrt><mrow><msub><mi>ε</mi><mrow><mtext>mach</mtext></mrow></msub><mfrac><mrow><mrow><mo>‖</mo><mrow><msup><mi>F</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mi>y</mi><mo stretchy="false">)</mo></mrow><mo>‖</mo></mrow><mo>+</mo><mn>2</mn><mi>η</mi></mrow><mi>γ</mi></mfrac><mrow><mo>‖</mo><mi>y</mi><mo>‖</mo></mrow></mrow></msqrt></mrow></math> where ε mach is the machine precision, γ represents the Lipschitz constant, η is the relative error in computing F N ( F TH ( y )), and <math id="M13"><mrow><mrow><mo>‖</mo><mo>⋅</mo><mo>‖</mo></mrow></mrow></math> denotes the L-2 norm in this study, except where otherwise noted. Note that the numeric values of <math id="M14"><mrow><mrow><mo>‖</mo><mrow><msup><mi>F</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mi>y</mi><mo stretchy="false">)</mo></mrow><mo>‖</mo></mrow></mrow></math> and γ are unavailable during the MFNK scheme process. To overcome this difficulty, the following approximation of Eq. (10) and Eq. (11) were adopted. <math display="block" id="M15"><mrow><mover accent="true"><mi>M</mi><mo stretchy="true">¯</mo></mover><mo>=</mo><munder><mrow><mi>max</mi></mrow><mrow><mn>1</mn><mo>≤</mo><mi>j</mi><mo>≤</mo><mi>m</mi></mrow></munder><mfrac><mrow><mrow><mo>‖</mo><mrow><mi>F</mi><mo stretchy="false">(</mo><mi>y</mi><mo>+</mo><mi>ε</mi><mtext>v</mtext><mo stretchy="false">)</mo><mo>−</mo><mi>F</mi><mo stretchy="false">(</mo><mi>y</mi><mo stretchy="false">)</mo></mrow><mo>‖</mo></mrow></mrow><mrow><mrow><mo>‖</mo><mrow><mi>ε</mi><mtext>v</mtext></mrow><mo>‖</mo></mrow></mrow></mfrac></mrow></math> <math display="block" id="M16"><mrow><mover accent="true"><mi>γ</mi><mo stretchy="true">¯</mo></mover><mo>=</mo><mfrac><mrow><mn>2</mn><mrow><mo>‖</mo><mrow><mi>F</mi><mo stretchy="false">(</mo><msup><mi>y</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msup><mo stretchy="false">)</mo></mrow><mo>‖</mo></mrow></mrow><mrow><msup><mrow><mrow><mo>‖</mo><mrow><msup><mtext>s</mtext><mi>k</mi></msup></mrow><mo>‖</mo></mrow></mrow><mn>2</mn></msup></mrow></mfrac></mrow></math> where, s k is the increment of the k th Newton iteration. By replacing <math id="M17"><mrow><mrow><mo>‖</mo><mrow><msup><mi>F</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mi>y</mi><mo stretchy="false">)</mo></mrow><mo>‖</mo></mrow></mrow></math> and γ in Eq. (9), with that given in Eq. (10) and (11), the optimal perturbation size is computed using Eq. (12). <math display="block" id="M18"><mrow><msub><mrow><mover accent="true"><mi>σ</mi><mo stretchy="true">¯</mo></mover></mrow><mrow><mtext>opt</mtext></mrow></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mrow><mo>‖</mo><mtext>v</mtext><mo>‖</mo></mrow></mrow></mfrac><msqrt><mrow><msub><mi>ε</mi><mrow><mtext>mach</mtext></mrow></msub><mfrac><mrow><mover accent="true"><mi>M</mi><mo stretchy="true">¯</mo></mover><mo>+</mo><mn>2</mn><mi>η</mi></mrow><mrow><mover accent="true"><mi>γ</mi><mo stretchy="true">¯</mo></mover></mrow></mfrac><mrow><mo>‖</mo><mi>y</mi><mo>‖</mo></mrow></mrow></msqrt></mrow></math> From the above derivation, it can be observed that there is almost no additional computing cost in the estimation of the optimal perturbation size. Additionally, the ratio of <math id="M19"><mrow><mover accent="true"><mi>M</mi><mo stretchy="true">¯</mo></mover><mo>/</mo><mover accent="true"><mi>γ</mi><mo stretchy="true">¯</mo></mover></mrow></math> is the most important, not the value of <math id="M20"><mrow><mover accent="true"><mi>M</mi><mo stretchy="true">¯</mo></mover></mrow></math> or <math id="M21"><mrow><mover accent="true"><mi>γ</mi><mo stretchy="true">¯</mo></mover></mrow></math> itself. Implementation of coupling neutronics with thermal-hydraulics 1 High-fidelity Neutron Transport Program (HNET) 2 HNET is a 3D whole-core transport code capable of performing high-fidelity neutron transport calculations for both steady-state and transient problems. The 2D/1D method with multilevel generalized equivalence theory-based CMFD acceleration (gCMFD) was developed in HNET for transport calculations. Specifically, the 2D method of characteristics (MOC) was utilized for the radial solver, and the two-node nodal expansion method is used for the axial dependency. The 2D and 1D solutions were coupled using transverse leakage. To decrease the computational cost of the CMFD, an innovative self-developed RSILU-preconditioned GMRES solver was utilized for the CMFD linear system. Regarding transient simulation, the multi-level predictor-corrector quasi-static scheme, together with the multi-level gCMFD acceleration, was developed to accelerate the time-dependent neutron transport calculation. Additionally, HNET employs the 47-group HELIOS library to provide cross sections, and the subgroup method was used for resonance self-shielding calculations. A detailed description of the transport calculations in HNET can be found in the References [ 14 14 , 15 15 , 19 19 , 20 20 ]. Thermal-hydraulic model 2 To provide the fuel temperature, coolant temperature, and coolant density to the neutronics field, a simplified internal thermal-hydraulic module was developed in the HNET. For pressurized water reactors, the module consists of a 1D radial heat conduction model and a 1D axial heat convection model, which have been proven to be effective for a wide range of applications [ 8 8 ]. The heat-transfer calculation began with a cylindrical conduction equation. Because of the large ratio of axil length to radius and the axisymmetric geometry, axial and circumferential heat transfer can be ignored; thus, only radial heat conduction is considered for the fuel rods, as in Eq. (13). <math display="block" id="M22"><mrow><mi>ρ</mi><msub><mi>c</mi><mtext>p</mtext></msub><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mfrac><mn>1</mn><mi>r</mi></mfrac><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>r</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mi>k</mi><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow><mi>r</mi><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>r</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><msup><mi>q</mi><mo>‴</mo></msup><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow></math> where T , ρ, c p , k , and r denote temperature, density, specific heat, thermal conductivity, and radius, respectively. The volumetric power density <math id="M23"><msup><mi>q</mi><mo>‴</mo></msup></math> contributes to the fission reaction. In addition, it is necessary to provide boundary conditions between heat conduction and convection. The heat flux transferred from the pellet to the fluid is determined by the following boundary condition as in Eq. (14). <math display="block" id="M24"><mrow><msub><mi>h</mi><mtext>w</mtext></msub><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>M</mtext></msub><mo>−</mo><msub><mi>T</mi><mtext>f</mtext></msub><mo stretchy="false">)</mo><mo>=</mo><mo>−</mo><msub><mrow><mrow><mrow><mi>k</mi><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>r</mi></mrow></mfrac></mrow><mo>|</mo></mrow></mrow><mrow><mi>r</mi><mo>=</mo><msub><mi>r</mi><mi>M</mi></msub></mrow></msub></mrow></math> where T M and T f are the temperatures of the cladding outer wall and coolant, respectively, and h w is the convective heat transfer coefficient, which is determined by the Dittus-Boelter correlation. In terms of the flow field, a single-channel model with the following reasonable assumptions is applied: (a) Each assembly is defined as a channel region. In addition, no cross-flow between the adjacent channels was considered. (b) The coolant was always in the liquid phase. (c) The coolant pressure was approximately constant during the steady state and was fast transient. Under these assumptions, the momentum conservation equation was removed. Thus, only the remaining mass and energy conservation equations were resolved in the flow field, and their general formulas in the axial direction are given by Eq. (15) and Eq. (16). <math display="block" id="M25"><mrow><mfrac><mrow><mo>∂</mo><mi>ρ</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>ρ</mi><mi>v</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mn>0</mn></mrow></math> <math display="block" id="M26"><mrow><mfrac><mrow><mo>∂</mo><mi>ρ</mi><mi>h</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>ρ</mi><mi>h</mi><mi>v</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><msub><mi>L</mi><mi>H</mi></msub></mrow><mrow><msub><mi>A</mi><mtext>c</mtext></msub></mrow></mfrac><msub><msup><mi>q</mi><mo>″</mo></msup><mtext>w</mtext></msub><mo>+</mo><msub><msup><mi>q</mi><mo>‴</mo></msup><mi>c</mi></msub></mrow></math> where ρ , v , and h denote the density, velocity, and enthalpy of coolant, respectively. z is the height in the axial direction and t represents time. L H and A c are the heated perimeter and the cross-sectional area of the channel, respectively. <math id="M27"><mrow><msub><msup><mi>q</mi><mo>″</mo></msup><mtext>w</mtext></msub></mrow></math> is the heat flux at the cladding surface, and <math id="M28"><mrow><msub><msup><mi>q</mi><mo>‴</mo></msup><mtext>c</mtext></msub></mrow></math> is the heat source generated directly inside the coolant. To solve the 1D radial heat conduction equations and 1D axial heat convection equations, the finite volume method is used to carry out spatial discretization, the theta method is adopted to discretize the time derivative terms, and the discrete formulas for the numerical calculation can be obtained. Coupling implementation 2 The implementation of the MFNK method in the steady-state coupling calculation is shown in Fig. 1 Fig. 1 . Notably, two coupling configurations can be arranged according to the choice of the initial guess. Because the dimension of power density is small in comparison with thermal-hydraulic solutions, the power density was selected as the initial estimate in this study. Next, to initiate the MFNK scheme with an appropriate initial power density y , the Picard iteration was performed first. Specifically, according to a primary y , thermal-hydraulic solutions are calculated and then used to update the macroscopic cross sections. Then, based on the results of the transport calculation, a new y Picard was obtained. The MFNK method began after the initial guess calculation. Based on the given y , the Picard iteration was performed to form the residual, as described in Eq. (3). Subsequently, the Newton iteration was solved using the GMRES method, where the required matrix-vector product was approximated through a finite difference approximation with the optimal perturbation size. In each Newton step, a new solution y was updated. The iterative scheme was repeated until convergence was achieved. In addition, the maximum Krylov subspace dimension was chosen as two, and the transport solver was fully converged in each calculation. In the transient simulation, the OS method was used to implement the transient coupling scheme. To improve stability, a staggered time mesh was applied to the OS method [ 21 21 ]. In this scheme, the two physical fields are advanced for different time step sizes ( Fig. 2 Fig. 2 ). The thermal-hydraulics problem was solved first with one-half of the neutronics time step. Subsequently, based on the updated thermal-hydraulic conditions, the neutronics solver was invoked to calculate the neutron flux. Then, the thermal-hydraulics solver was advanced again using the new power density with one-half of the neutronics time step. Numerical results and analysis 1 VERA Problem 6 [ 17 17 ] with different power levels and boron concentrations was first employed to study the convergence performance of the MFNK method. A comparison of the numerical results for MFNK, Picard, and Picard with relaxation is presented in this section. Subsequently, the NEACRP transient benchmark [ 18 18 ] was used to preliminarily test and verify the transient coupling calculation. The corresponding results were compared with reference solutions for verification. Assessment of MFNK performance 2 VERA Problem 6 describes a single fuel assembly problem under typical full-power and nominal thermal-hydraulic conditions, which are presented in Table 1 Table 1 . The assembly consisted of 264 fuel rods, 24 guide tubes, and an instrument tube at the center. Meanwhile, there were six intermediate grids and two end grids in the assembly along the axial direction. More detailed geometric and material specifications are provided in [ 17 17 ]. 1 1 Parameter 1 1 Value 1 1 System pressure (psia) 1 1 2250 1 1 Rated power (100%) (MW) 1 1 17.67 1 1 Rated coolant mass flow (Mlbs/hr) 1 1 0.6823 1 1 Inlet coolant temperature (K) 1 1 565 1 1 Boron concentration (ppm) 1 1 1300 For VERA Problem 6, the transport equation was solved using the 2D/1D method via gCMFD acceleration. Based on the 47-group HELIOS library, nuclide density, and temperature, the cross-sections were obtained through resonance treatment. Fig. 3 Fig. 3 illustrates the discretization of MOC mesh in radial direction, the heterogeneous fuel pins are divided into 56 flat source regions (FSR), and the reflector cells use the 2 × 2 Cartesian sub-mesh. In addition, a ray spacing of 0.03 cm with 64 azimuthal angles and a Tabuchi-Yamamoto polar quadrature using 3 polar angles per half-space were utilized for MOC calculation. The active fuel region was divided into 49 axial layers, similar to VERA Problem 3 [ 17 17 ]. Furthermore, the core plates, nozzles, and grids were explicitly modeled for the transport domain. The radial boundaries were set to be reflective, and the vacuum boundary condition was applied axially, owing to the top and bottom water reflectors. For the thermal-hydraulic domain, each fuel pin, including the fuel pellet, gap, and cladding, was explicitly modeled and calculated individually using the radial power density inside the fuel pin. The discretization of the coolant channel in the axial direction is the same as that in the transport calculation. In addition, the inlet coolant conditions were adopted for the bottom structural material below the active fuel, whereas the average exit conditions of the coolant at the top of the active fuel region were used as the top structural material. Generally, a fission source is used to measure the convergence of a coupled system. To evaluate the convergence behavior of MFNK, a strict convergence criterion using the power density residual norm was set for VERA Problem 6, as in Eq. (17). <math display="block" id="M29"><mrow><mtext>RMS</mtext><mo>=</mo><msqrt><mrow><mfrac><mrow><msup><mrow><mrow><mo>‖</mo><mrow><mi>F</mi><mo stretchy="false">(</mo><mi>y</mi><mo stretchy="false">)</mo></mrow><mo>‖</mo></mrow></mrow><mn>2</mn></msup></mrow><mi>n</mi></mfrac></mrow></msqrt><mo>≤</mo><mn>1</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>6</mn></mrow></msup></mrow></math> where RMS is the root mean square of the residual vector and n is the linear system dimension. First, according to Eq. (6), the influence of perturbation size on the MFNK method was investigated through VERA Problem 6. To test the behavior of the optimal perturbation size defined by Eq. (12), three empirical perturbation sizes were considered for the comparison. The first was a simple empirical perturbation with a constant value of 0.001. The second method is based on Pernice’s method [ 22 22 ]. <math display="block" id="M30"><mrow><mi>ε</mi><mo>=</mo><mfrac><mrow><msqrt><mrow><mo stretchy="false">(</mo><mn>1</mn><mo>+</mo><mrow><mo>‖</mo><mi>y</mi><mo>‖</mo></mrow><mo stretchy="false">)</mo><msub><mi>ε</mi><mrow><mtext>mach</mtext></mrow></msub></mrow></msqrt></mrow><mrow><mrow><mo>‖</mo><mtext>v</mtext><mo>‖</mo></mrow></mrow></mfrac></mrow></math> The third empirical perturbation size proposed by Knoll [ 11 11 ] is given by Eq. (19). <math display="block" id="M31"><mrow><mi>ε</mi><mo>=</mo><mfrac><mn>1</mn><mrow><mi>n</mi><mrow><mo>‖</mo><mtext>v</mtext><mo>‖</mo></mrow></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><mi>b</mi><mrow><mo>|</mo><mrow><msub><mi>y</mi><mi>i</mi></msub></mrow><mo>|</mo></mrow><mo>+</mo><mi>b</mi></mrow></mstyle></mrow></math> where b is a constant whose magnitude is within a few orders of magnitude of the square root of the machine round-off. Fig. 4 Fig. 4 illustrates how the perturbation size affects the convergence behavior of the MFNK method for VERA Problem 6. It can be observed that MFNK with a constant perturbation size converges very slowly. The performances of the other two empirical perturbation sizes are relatively consistent, and their residual norms seem to stagnate around 0.1. By contrast, MFNK with an optimal perturbation size can converge further and satisfy the convergence criterion. The primary reason is that the optimal perturbation size can reduce the error introduced by the finite-difference approximation even when the residual norm becomes small, whereas the empirical ones cannot. Hence, the optimal perturbation size is used in the following analysis. Subsequently, a convergence analysis was performed by comparing the results between the MFNK and Picard methods. Notably, in order to make an appropriate comparison between the MFNK and Picard methods, the convergence history was monitored by the iteration index of the MOC, which is a time-consuming procedure, rather than their individual outer iterations. The convergence history of k eff and the maximum fuel temperature for VERA problem 6 are shown in Fig. 5 Fig. 5 and Fig. 6 Fig. 6 , respectively. The results indicate that MFNK achieved a stable state quickly, whereas Picard initially exhibited oscillatory behavior. For instance, the oscillation of k eff in MFNK was a maximum of 1 pcm, whereas the maximum variation of k eff in Picard is approximately 38 pcm. In addition, as illustrated in Fig. 5 Fig. 5 and Fig. 6 Fig. 6 , the solutions resulting from MFNK and Picard match well. To further evaluate the performance of MFNK, VERA Problem 6 with a modification of the power level and boron concentration was simulated. Traditionally, relaxation has been introduced to Picard for better convergence behavior. Hence, Picard and Picard with relaxation applied to the temperature solution were performed for comparison. Because the optimal choice of relaxation factor depends strongly on the problems, a relaxation factor of w = 0.5 was adopted in this paper [ 8 8 ]. The overall convergence performances of MFNK, Picard, and Picard with relaxation are summarized in Table 2 Table 2 and Table 3 Table 3 . In addition to the MOC iterations required for convergence, the total number of resonance calculations is also presented because of the relatively high computational cost. Notably, the resonance calculation was performed in the coupling scheme only when the maximum temperature change between two adjacent resonance calculations was larger than 1 K. 1 1 Couple Scheme 1 1 Iterations 5 1 Power level (%) 1 1   1 1   1 1 60 1 1 80 1 1 100 1 1 120 1 1 140 1 1 MFNK 1 1 MOC 1 1 85 1 1 105 1 1 107 1 1 126 1 1 122 1 1 Picard ( w = 0.5) 1 1 Resonance 1 1 10 1 1 12 1 1 11 1 1 10 1 1 11 1 1   1 1 MOC 1 1 103 1 1 113 1 1 117 1 1 126 1 1 137 1 1   1 1 Resonance 1 1 11 1 1 11 1 1 11 1 1 12 1 1 12 1 1 Picard 1 1 MOC 1 1 72 1 1 122 1 1 197 1 1 divergence 1 1 divergence 1 1   1 1 Resonance 1 1 5 1 1 10 1 1 26 1 1 divergence 1 1 divergence 1 1 Couple Scheme 1 1 Iterations 5 1 Boron concentration (ppm) 1 1   1 1   1 1 0 1 1 500 1 1 900 1 1 1300 1 1 1700 1 1 MFNK 1 1 MOC 1 1 112 1 1 115 1 1 111 1 1 107 1 1 113 1 1 Picard ( w = 0.5) 1 1 Resonance 1 1 11 1 1 12 1 1 12 1 1 11 1 1 13 1 1   1 1 MOC 1 1 119 1 1 116 1 1 117 1 1 117 1 1 122 1 1   1 1 Resonance 1 1 12 1 1 12 1 1 12 1 1 11 1 1 11 1 1 Picard 1 1 MOC 1 1 177 1 1 209 1 1 223 1 1 197 1 1 216 1 1   1 1 Resonance 1 1 20 1 1 27 1 1 30 1 1 26 1 1 29 As shown in Table 2 Table 2 , MFNK and Picard with relaxation can considerably reduce the MOC iterations compared with Picard over all except for the 60% power level case. Although the three methods displayed a dependence on the power magnitude, MFNK and Picard with relaxation were less susceptible to the power level than Picard. Furthermore, it was observed that MFNK and Picard with relaxation could achieve convergence criteria in all cases, whereas Picard performed poorly and even diverged as the power level increased. This is primarily because increasing the power enhances coupling strength. However, Picard is generally unable to handle this stronger coupling. Additionally, the number of MOC iterations for MFNK was the same as or less than that for Picard with relaxation. More specifically, MFNK can decrease the number of MOC iterations from Picard with a maximum relaxation of 15. In terms of boron perturbation, Table 3 Table 3 demonstrates that varying the boron concentration does not affect the convergence behavior. As shown in Table 3 Table 3 , the performance of MFNK is slightly better than that of Picard with relaxation, and poor convergence can be found in the Picard results. Table 2 Table 2 and Table 3 Table 3 show that the number of resonance calculations for MFNK and Picard with relaxation are comparable, and both are less than that for Picard as a whole. In other words, the cross sections can converge faster in MFNK and Picard with relaxation. Fig. 7 Fig. 7 compares the power density residual norm of MFNK against Picard with and without relaxation with respect to the variation in the power level. In contrast to Picard with relaxation, the residual norm of MFNK dropped more rapidly in the 80% and 140% power cases. In addition, Picard converged at a markedly slower rate and divergence behavior occurred in Picard at a 140% power level. Based on the above analysis, it can be concluded that although MFNK does not perform as efficiently as Picard at the 60% power level, it can display a better convergence performance than Picard for a wider range of problems. Moreover, although the performance of MFNK is better than that of Picard with a relaxation factor of 0.5 in given cases, the opposite results may occur if the optimum values of the relaxation factor are selected. Nevertheless, it should be easier to improve the convergence rate and enhance the numerical stability when utilizing MFNK than when applying a relaxation scheme, because the optimal relaxation factor is problem-dependent and generally chosen empirically [ 21 21 ]. Transient simulation of control rod ejection 2 The NEACRP benchmark is a 3D PWR core transient problem constructed to simulate the ejection of a control assembly (CA) from the core under operational conditions. The assembly arrangement is shown in Fig. 8 Fig. 8 . Moreover, 2-group macroscopic cross sections were set for 11 homogenized compositions. However, the NEACRP benchmark was originally specified for diffusion-based codes, resulting in a lack of total and self-scattering cross sections. To obtain an equivalent transport version, the total cross sections can be approximated to transport cross sections, and then the self-scattering cross sections are calculated using Eq. (20). <math display="block" id="M32"><mrow><msub><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mtext>g</mtext><mo>→</mo><mtext>g</mtext></mrow></msub><mo>=</mo><msub><mi>Σ</mi><mrow><mtext>t</mtext><mo>,</mo><mtext>g</mtext></mrow></msub><mo>−</mo><msub><mi>Σ</mi><mrow><mtext>a</mtext><mo>,</mo><mtext>g</mtext></mrow></msub><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mtable columnalign="left"><mtr><mtd><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow></mtd></mtr><mtr><mtd><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>≠</mo><mi>g</mi></mrow></mtd></mtr></mtable><mi>G</mi></munderover><mrow><msub><mi>Σ</mi><mrow><mtext>s</mtext><mo>,</mo><mtext>g</mtext><mo>→</mo><msup><mtext>g</mtext><mo>′</mo></msup></mrow></msub></mrow></mstyle></mrow></math> where Σ t , Σ a , and Σ s are the total, absorption, and scattering cross-sections, respectively. The NEACRP benchmark comprises six cases with different combinations of control arrangements, original positions, and power levels. A detailed description of the geometry, macroscopic cross-sections, and initial configuration of the control assemblies can be found in [ 18 18 ]. In this study, cases A2 and B2 were considered to verify the accuracy of the transient simulation. Both cases were initiated by the rapid ejection of the CA in 0.1 s at hot full power state; the CA to be ejected is marked by a circle in Fig. 8 Fig. 8 . Specifically, the symbols A and B represent the ejected CA of cases A2 and B2, respectively. To resolve the NEACRP benchmark using the transient 2D/1D transport solver in HNET [ 15 15 ], the geometry was uniformly discretized into an FSR mesh size of 0.7202 cm × 0.7202 cm for the MOC calculation. Additionally, the ray spacing for the NEACRP benchmark is the same as that for VERA Problem 6. From bottom to top, the active zone was divided into 26 layers with heights of 7.7, 11, 15 (21 layers), 12.8 (2 layers), and 8.0 cm, which differed from the division of the NEACRP benchmark for high accuracy of results. A reflective boundary was used to the west and north because of symmetry, whereas the rest used vacuum boundary conditions. Concerning the thermal-hydraulic calculation, an average fuel pin, representing all fuel pins of the associated assembly, was modeled in detail to obtain the radial temperature distribution by applying the average assembly power density. In the axial direction, the separation of the coolant channel was consistent with that of the transport calculation. The steady-state thermal-hydraulic conditions are listed in Table 4 Table 4 . The relation between thermal properties, such as thermal conductivity and specific heat capacity, and temperature is derived from [ 18 18 ]. Moreover, to consider the feedback effects from thermal hydraulics and boron concentration, the macroscopic cross-sections for each composition are defined by Eq. (21). <math display="block" id="M33"><mrow><mi>Σ</mi><mo>=</mo><msub><mi>Σ</mi><mn>0</mn></msub><mo>+</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>x</mi><mo>=</mo><msub><mi>ρ</mi><mrow><mtext>m</mtext><mo>,</mo><msqrt><mrow><msub><mtext>T</mtext><mtext>f</mtext></msub></mrow></msqrt><mo>,</mo><mtext>C</mtext><mo>,</mo><msub><mtext>T</mtext><mtext>m</mtext></msub></mrow></msub></mrow></munder><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><mi>Σ</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac></mrow><mo>)</mo></mrow></mrow></mstyle><mo stretchy="false">(</mo><mi>x</mi><mo>−</mo><msub><mi>x</mi><mn>0</mn></msub><mo stretchy="false">)</mo></mrow></math> where ρ , T m , T f , and C are the moderator temperature, moderator density, fuel temperature, and boron density, respectively. The subscript 0 represents the reference values. 1 1 Parameters 1 1 Value 1 1 Core thermal output (MW) 1 1 2775 (HFP) 1 1 Core pressure (Mpa) 1 1 15.5 1 1 Coolant inlet temperature (℃) 1 1 286 1 1 Mass flow (kg/s) 1 1 12893 A comparison of steady-state and transient solutions for cases A2 and B2 is given in Table 5 Table 5 . The reference results listed in Table 5 Table 5 were derived from [ 23 23 ]. As shown in Table 5 Table 5 , there is good agreement between the simulation and reference results. The critical boron concentration at the steady state deviated slightly from the reference value, and the maximum difference was approximately 31 ppm, which is considered acceptable. In addition, the representative parameters in the transient process coincided well with the reference. 1 1 Parameters 3 1 A2 3 1 B2 1 1   1 1 Reference 1 1 Calculation 1 1 Deviation 1 1 Reference 1 1 Calculation 1 1 Deviation 1 1 Critical boron concentration (ppm) 1 1 1156.63 1 1 1187.66 1 1 31.03 1 1 1183.83 1 1 1214.79 1 1 30.96 1 1 Power peak 1 1 1.0800 1 1 1.0821 1 1 0.0021 1 1 1.0630 1 1 1.0651 1 1 0.0021 1 1 Peak time (s) 1 1 0.0100 1 1 0.0105 1 1 0.0005 1 1 0.1200 1 1 0.1200 1 1 0 1 1 Final power at 5 s 1 1 1.0350 1 1 1.0357 1 1 0.0007 1 1 1.0380 1 1 1.0399 1 1 0.0019 1 1 Coolant outlet average temperature at 5 s (℃) 1 1 324.600 1 1 325.049 1 1 0.449 1 1 324.700 1 1 325.164 1 1 0.464 Fig. 9 Fig. 9 depicts the evolution of the core-normalized power and reactivity for case A2. Here, the adopted reference solutions are computed using the PARCS code [ 24 24 ]. As indicated, the HNET results for case A2 are in agreement with that of PARCS. The slight deviation in the normalized power and reactivity curves can be explained by the difference in the neutronics methodology. The NEACRP benchmark is resolved by neutron transport calculations in HNET, whereas the nodal method is used in PARCS. In addition, as illustrated in Fig. 10 (a) Fig. 10 (a) , there is no significant difference in the responses of the coolant exit temperature to the CA ejection between these two codes. Because the result of the Doppler temperature from the PARCS is not provided in [ 24 24 ], Fig. 10 (b) Fig. 10 (b) plots only the solution of the HNET. From Fig. 10 (b) Fig. 10 (b) , an increase in the average Doppler temperature owing to power increments is observed, which behaves as expected. As shown in Fig. 11 Fig. 11 and Fig. 12 Fig. 12 , to further test the level of agreement, a closer analysis of the axial relative power and radial normalized power for case A2 was performed at 0 s, 0.1 s, and 5 s. The relative error e is defined by Eq. (22). <math display="block" id="M34"><mrow><mi>e</mi><mo>=</mo><mfrac><mrow><mi>c</mi><mi>a</mi><mi>l</mi><mi>c</mi><mo>−</mo><mi>r</mi><mi>e</mi><mi>f</mi></mrow><mrow><mi>r</mi><mi>e</mi><mi>f</mi></mrow></mfrac><mo>×</mo><mn>100</mn><mi>%</mi></mrow></math> where calc represents the values calculated by HNET, and ref is the reference value derived from PARCS [ 24 24 ]. The normalized power in the pin and assembly levels are calculated using Eqs. (23) and (24),. <math display="block" id="M35"><mrow><msub><mrow><mover accent="true"><mi>P</mi><mo stretchy="true">¯</mo></mover></mrow><mrow><mtext>pin</mtext></mrow></msub><mo>=</mo><mfrac><mrow><msub><mi>P</mi><mrow><mtext>pin</mtext></mrow></msub><msub><mi>N</mi><mrow><mi>p</mi><mi>i</mi><mi>n</mi></mrow></msub></mrow><mrow><msub><mi>P</mi><mrow><mtext>core</mtext></mrow></msub></mrow></mfrac></mrow></math> <math display="block" id="M36"><mrow><msub><mrow><mover accent="true"><mi>P</mi><mo stretchy="true">¯</mo></mover></mrow><mrow><mtext>assem</mtext></mrow></msub><mo>=</mo><mfrac><mrow><msub><mi>P</mi><mrow><mtext>assem</mtext></mrow></msub><msub><mi>N</mi><mrow><mtext>assem</mtext></mrow></msub></mrow><mrow><msub><mi>P</mi><mrow><mtext>core</mtext></mrow></msub></mrow></mfrac></mrow></math> where P pin and P assem denote the axial-integrated power per fuel pin and per fuel assembly, respectively; N pin and N assem represent the total number of fuel pins and fuel assemblies, respectively; and P core is the core power. With respect to the axial relative power, no significant difference between HNET and PARCS is observed in Fig. 11 Fig. 11 . Fig. 12 Fig. 12 shows that the results of assembly level normalized power at the considered temporal points from HNET agree well with PARCS. Moreover, because the NEACRP benchmark was performed using a transient transport solver, the pin-wise power distribution can be obtained directly without power reconstruction. Fig. 13 Fig. 13 depicts the pin-wise radial normalized power for case A2 at 0.0 s and 0.1 s. The comparison shows that the pin power in the center increases significantly after the withdrawal of the control rod, which coincides with the sharp increase in power observed in Fig. 9 Fig. 9 . The above results prove that the evolution of the power distribution can be effectively captured when the reactivity is being inserted. Finally, because the thermal-hydraulics feedback is provided at the assembly level for the NEACRP benchmark, Fig. 14 Fig. 14 shows the local temperatures in assemblies C-3 and B-8, where the highest and lowest powers appear, respectively. It can be observed that there are obvious discrepancies for fuel centerline temperatures between the two assemblies. As shown in Fig. 14 (b) Fig. 14 (b) , the increase in temperatures between 0 s and 5 s in assembly B-8 is larger than that in assembly C-3, particularly the fuel centerline temperature. Conclusion 1 The primary objective of this study was to couple neutronics with thermal hydraulics in HNET. First, a simplified internal thermal-hydraulic module was developed for thermal properties. Subsequently, a steady-state coupling calculation is achieved by implementing the MFNK method. Next, the OS method with a staggered time mesh was used to accomplish the transient coupling simulation. The investigation of the perturbation size shows that the optimal perturbation size can significantly enhance the convergence of the MFNK method compared with the empirical perturbation size. With respect to VERA problem 6 at different power levels and boron concentrations, the MFNK method is slightly sensitive to the feedback intensity and can considerably reduce the computational cost compared with the Picard iteration in most cases. It can be proven that the MFNK method is more stable and efficient than the Picard iteration for a wider range of problems. Notably, the performance of MFNK is as good as or better than that of Picard with relaxation, and there is no requirement for input data in MFNK compared with the relaxation approach. Furthermore, the numerical results for the NEACRP transient benchmark are in agreement with the reference values, indicating the validation of the transient coupling calculation in the hot full-power state. In conclusion, the results demonstrate that the neutronics and thermal-hydraulics coupling simulation is effectively performed in HNET under both steady-state and transient conditions. References R.R. Yang , Y. Yuan , C. Hao et al. , keff uncertainty quantification and analysis due to nuclear data during the full lifetime burnup calculation for a small-sized prismatic high temperature gas-cooled reactor . Nucl. Sci. Tech. 32 ( 11 ), 127 ( 2021 ). doi: 10.1007/s41365-021-00969-w http://doi.org/10.1007/s41365-021-00969-w keff uncertainty quantification and analysis due to nuclear data during the full lifetime burnup calculation for a small-sized prismatic high temperature gas-cooled reactor J.A. Turner , K. Clarno , M. Sieger et al. , The Virtual Environment for Reactor Applications (VERA): Design and architecture . J. Comput. Phys. 326 , 544 - 568 ( 2016 ). doi: 10.1016/j.jcp.2016.09.003 http://doi.org/10.1016/j.jcp.2016.09.003 The Virtual Environment for Reactor Applications (VERA): Design and architecture D. She , J. Guo , Z.H. Liu et al. , PANGU code for pebble-bed HTGR reactor physics and fuel cycle simulations . Ann. Nucl. Energy 126 , 48 - 58 ( 2019 ). doi: 10.1016/j.anucene.2018.11.005 http://doi.org/10.1016/j.anucene.2018.11.005 PANGU code for pebble-bed HTGR reactor physics and fuel cycle simulations D. She , B. Xia , J. Guo et al. , Prediction calculations for the first criticality of the HTR-PM using the PANGU code . Nucl. Sci. Tech. 32 , 90 ( 2021 ). doi: 10.1007/s41365-021-00936-5 http://doi.org/10.1007/s41365-021-00936-5 Prediction calculations for the first criticality of the HTR-PM using the PANGU code F. Gou , Y. Liu , F.B. Chen et al. , Thermal behavior of the HTR-10 under combined PLOFC and ATWS condition initiated by unscrammed control rod withdrawal . Nucl. Sci. Tech. 29 , 123 ( 2018 ). doi: 10.1007/s41365-018-0472-3 http://doi.org/10.1007/s41365-018-0472-3 Thermal behavior of the HTR-10 under combined PLOFC and ATWS condition initiated by unscrammed control rod withdrawal B. Deng , Y. Cui , J.G. Chen et al. , Core and blanket thermal-hydraulic analysis of a molten salt fast reactor based on coupling of OpenMC and OpenFOAM . Nucl. Sci. Tech. 31 , 85 ( 2020 ). doi: 10.1007/s41365-020-00803-9 http://doi.org/10.1007/s41365-020-00803-9 Core and blanket thermal-hydraulic analysis of a molten salt fast reactor based on coupling of OpenMC and OpenFOAM R. Vadi , K. Sepanloo , An improved semi-implicit direct kinetics method for transient analysis of nuclear reactors . Nucl. Sci. Tech. 30 ( 11 ), 162 ( 2019 ). doi: 10.1007/s41365-019-0690-3 http://doi.org/10.1007/s41365-019-0690-3 An improved semi-implicit direct kinetics method for transient analysis of nuclear reactors Q. Shen , S. Choi , B. Kochunas , A robust , relaxation-free multiphysics iteration scheme for CMFD-accelerated neutron transport k-eigenvalue calculations-II: Numerical results . Nucl. Sci. Eng. 195 ( 11 ), 1202 - 1235 ( 2021 ). doi: 10.1080/00295639.2021.1906586 http://doi.org/10.1080/00295639.2021.1906586 relaxation-free multiphysics iteration scheme for CMFD-accelerated neutron transport k-eigenvalue calculations-II: Numerical results A. Facchini , J. Lee , H.G. Joo , Investigation of Anderson acceleration in neutronics-thermal hydraulics coupled direct whole core calculation . Ann. Nucl. Energy 153 , 108042 ( 2021 ). doi: 10.1016/j.anucene.2020.108042 http://doi.org/10.1016/j.anucene.2020.108042 Investigation of Anderson acceleration in neutronics-thermal hydraulics coupled direct whole core calculation A. Hajizadeh , H. Kazeminejad , S. Talebi , High-order fully implicit SIMPLE-based model for fully implicit simulation of upward two-phase flow . Nucl. Sci. Tech. 29 ( 8 ), 114 ( 2018 ). doi: 10.1007/s41365-018-0456-3 http://doi.org/10.1007/s41365-018-0456-3 High-order fully implicit SIMPLE-based model for fully implicit simulation of upward two-phase flow Y. Xu , A matrix free Newton/Krylov method for coupling complex multi-physics subsystems . Dissertation, Purdue University , 2004 . A matrix free Newton/Krylov method for coupling complex multi-physics subsystems O. Zerkak , T. Kozlowski , I. Gajev , Review of multi-physics temporal coupling methods for analysis of nuclear reactors . Ann. Nucl. Energy 84 , 225 - 233 ( 2015 ). doi: 10.1016/j.anucene.2015.01.019 http://doi.org/10.1016/j.anucene.2015.01.019 Review of multi-physics temporal coupling methods for analysis of nuclear reactors C. Hao , L.X. Liu , L. Kang et al. , An efficient hybrid multi-level CMFD in space and energy for accelerating the high-fidelity neutron transport calculation . Ann. Nucl. Energy 161 , 108446 ( 2021 ). doi: 10.1016/j.anucene.2021.108446 http://doi.org/10.1016/j.anucene.2021.108446 An efficient hybrid multi-level CMFD in space and energy for accelerating the high-fidelity neutron transport calculation C. Hao , Y. Xu , T.J. Downar , Multi-level coarse mesh finite difference acceleration with local two-node nodal expansion method . Ann. Nucl. Energy 116 , 105 - 113 ( 2018 ). doi: 10.1016/j.anucene.2018.02.002 http://doi.org/10.1016/j.anucene.2018.02.002 Multi-level coarse mesh finite difference acceleration with local two-node nodal expansion method L. Kang , C. Hao , Q. Zhao et al. , Two-level time-dependent GET based CMFD acceleration for 3D whole core transient transport simulation using 2D/1D method . Ann. Nucl. Energy 142 , 107405 ( 2020 ). doi: 10.1016/j.anucene.2020.107405 http://doi.org/10.1016/j.anucene.2020.107405 Two-level time-dependent GET based CMFD acceleration for 3D whole core transient transport simulation using 2D/1D method L.X. Liu , C. Hao , Y. Xu , Equivalent low-order angular flux nonlinear finite difference equation of MOC transport calculation . Nucl. Sci. Tech. 31 , 125 ( 2020 ). doi: 10.1007/s41365-020-00834-2 http://doi.org/10.1007/s41365-020-00834-2 Equivalent low-order angular flux nonlinear finite difference equation of MOC transport calculation A.T. Godfrey , VERA core physics benchmark progression problem specifications , CASL-U-2012-0131-004. Consortium for Advanced Simulation of Light Water Reactors ( 2014 ). VERA core physics benchmark progression problem specifications H. Finnemann , NEACRP 3-D LWR core transient benchmark final specifications, NEACRP-L-335 . OECD Nuclear Energy Agency ( 1991 ). NEACRP 3-D LWR core transient benchmark final specifications, NEACRP-L-335 L. Kang , C. Hao , Q. Zhao et al. , The advanced multilevel predictor-corrector quasi-static method for pin-resolved neutron kinetics simulation . Front. Energy Res. 9 , 747148 ( 2021 ). doi: 10.3389/fenrg.2021.747148 http://doi.org/10.3389/fenrg.2021.747148 The advanced multilevel predictor-corrector quasi-static method for pin-resolved neutron kinetics simulation C. Hao , L. Kang , Y. Xu et al. , 3D whole-core neutron transport simulation using 2D/1D method via multi-level generalized equivalence theory based CMFD acceleration . Ann. Nucl. Energy 122 , 79 - 90 ( 2018 ). doi: 10.1016/j.anucene.2018.08.014 http://doi.org/10.1016/j.anucene.2018.08.014 3D whole-core neutron transport simulation using 2D/1D method via multi-level generalized equivalence theory based CMFD acceleration Q. Shen , Robust and efficient methods in transient whole-core . Dissertation, University of Michigan , 2021 . Robust and efficient methods in transient whole-core D.A. Knoll , D.E. Keyes , Jacobian-free Newton-Krylov methods: a survey of approaches and applications . J. Comput. Phys. 193 ( 2 ), 357 - 397 ( 2004 ). doi: 10.1016/j.jcp.2003.08.010 http://doi.org/10.1016/j.jcp.2003.08.010 Jacobian-free Newton-Krylov methods: a survey of approaches and applications H. Finnemann , H. Bauer , A. Galati et al. , Results of LWR core transient benchmark, NEA/NSC/DOC (93)25 . OECD Nuclear Energy Agency ( 1993 ). Results of LWR core transient benchmark, NEA/NSC/DOC (93)25 A2: central rod ejection at HFP . https://engineering.purdue.edu/PARCS/Code/TestSuite/CalculationMode/StandAloneMode/NEACRPPWR https://engineering.purdue.edu/PARCS/Code/TestSuite/CalculationMode/StandAloneMode/NEACRPPWR A2: central rod ejection at HFP All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Yanling Zhu, Xingwu Chen, Chen Hao, Yizhen Wang and Yunlin Xu. The first draft of the manuscript was written by Yanling Zhu, Chen Hao and Yizhen Wang, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-022-01120-z This study was supported by the China Postdoctoral Science Foundation (No. 2021M703045), the National Natural Science Foundation of China (No. 12075067), and the National Key R&D Program of China (No. 2018YFE0180900). 2022-07-11 2022-08-16 2022-09-12 2022-11-18 2022-12-16 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2022 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-022-01120-z http://dx.doi.org/10.1007/s41365-022-01120-z . Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law. 2022 Nuclear Science and Techniques 11 33