Derivation of MFNK method

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). $x^{k+1}=F_{\text{TH}}({\overline{y}}^k)$ (1) $y^{k+1}=F_{\text{N}}({\overline{x}}^{k+1})$ (2) 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). $F(y^k)=y^k-F_{\text{N}}\left(F_{\text{TH}}(y^k)\right)=0$ (3)

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). $F(y^k+\delta y)=F(y^k)+F^{\prime }(y^k)\delta y=0$ (4)

Then, Eq. (4) can be rearranged as shown in Eq. (5), which is referred to as Newton iteration. $J(y^k)\delta y=-F(y^k)$ (5)

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). $\begin{array}{c}J\left(y^k\right)\text{v}=\frac{F(y^k+\epsilon \text{v})-F(y^k)}{\epsilon }\\ =\text{v}-\frac{F_{\text{N}}\left(F_{\text{TH}}(y^k+\epsilon \text{v})\right)-F_N\left(F_{\text{TH}}(y^k)\right)}{\epsilon }\end{array}$ (6) where ε is the perturbation size and v is the basis vector of the Krylov iteration.

Upon the δy, the power density is updated: $y^{k+1}=y^k+\delta y$ (7)

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

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] 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). $\begin{array}{c}{e}_{\text{d}}=\frac{F(y+\epsilon \text{v})-F(y)-\epsilon F^{\prime }(y)\text{v}}{\epsilon }\\ +\frac{\overline{F}(\overline{y+\epsilon \text{v}})-F(y+\epsilon \text{v})-\overline{F}(y)+F(y)}{\epsilon }\end{array}$ (8) where $\overline{y+\epsilon \text{v}}$ is the finite precision representation of $y+\epsilon \text{v}$ and $\overline{F}$ 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], 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). ${\epsilon }_{\text{opt}}=\frac{1}{\Vert \text{v}\Vert }\sqrt{{\epsilon }_{\text{mach}}\frac{\Vert F^{\prime }(y)\Vert +2\eta }{\gamma }\Vert y\Vert }$ (9) where ε_mach is the machine precision, γ represents the Lipschitz constant, η is the relative error in computing F_N(F_TH(y)), and $\Vert \cdot \Vert$ denotes the L-2 norm in this study, except where otherwise noted.

Note that the numeric values of $\Vert F^{\prime }(y)\Vert$ and γ are unavailable during the MFNK scheme process. To overcome this difficulty, the following approximation of Eq. (10) and Eq. (11) were adopted. $\overline{M}=\underset{1\le j\le m}{\max }\frac{\Vert F(y+\epsilon \text{v})-F(y)\Vert }{\Vert \epsilon \text{v}\Vert }$ (10) $\overline{\gamma }=\frac{2\Vert F(y^{k+1})\Vert }{{\Vert {\text{s}}^k\Vert }^2}$ (11) where, s^k is the increment of the kth Newton iteration.

By replacing $\Vert F^{\prime }(y)\Vert$ and γ in Eq. (9), with that given in Eq. (10) and (11), the optimal perturbation size is computed using Eq. (12). ${\overline{\sigma }}_{\text{opt}}=\frac{1}{\Vert \text{v}\Vert }\sqrt{{\epsilon }_{\text{mach}}\frac{\overline{M}+2\eta }{\overline{\gamma }}\Vert y\Vert }$ (12) 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 $\overline{M}/\overline{\gamma }$ is the most important, not the value of $\overline{M}$ or $\overline{\gamma }$ itself.

High-fidelity Neutron Transport Program (HNET)

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, 15, 19, 20].

Thermal-hydraulic model

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].

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). $\rho c_{\text{p}}\frac{\partial T}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(k\left(T\right)r\frac{\partial T}{\partial r}\right)+q^{‴}\left(r,t\right)$ (13) where T, ρ, c_p, k, and r denote temperature, density, specific heat, thermal conductivity, and radius, respectively. The volumetric power density $q^{‴}$ 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). $h_{\text{w}}(T_{\text{M}}-T_{\text{f}})=-{k\frac{\partial T}{\partial r}|}_{r=r_M}$ (14) 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). $\frac{\partial \rho }{\partial t}+\frac{\partial \rho v}{\partial z}=0$ (15) $\frac{\partial \rho h}{\partial t}+\frac{\partial \rho hv}{\partial z}=\frac{L_H}{A_{\text{c}}}{q^{″}}_{\text{w}}+{q^{‴}}_c$ (16) 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. ${q^{″}}_{\text{w}}$ is the heat flux at the cladding surface, and ${q^{‴}}_{\text{c}}$ 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

The implementation of the MFNK method in the steady-state coupling calculation is shown in 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.

Fig. 1Full-screen View

Solution scheme of MFNK for steady-state coupling calculation

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]. In this scheme, the two physical fields are advanced for different time step sizes (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.

Fig. 2Full-screen View

(Color online) OS method with staggered time step

Assessment of MFNK performance

VERA Problem 6 describes a single fuel assembly problem under typical full-power and nominal thermal-hydraulic conditions, which are presented in 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].

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 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]. 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.

Fig. 3Full-screen View

(Color online) The division of MOC mesh in radial direction: (a) fuel pin; (b) reflector cell

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). $\text{RMS}=\sqrt{\frac{{\Vert F(y)\Vert }^2}{n}}\le 1\times {10}^{-6}$ (17) 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]. $\epsilon =\frac{\sqrt{(1+\Vert y\Vert ){\epsilon }_{\text{mach}}}}{\Vert \text{v}\Vert }$ (18)

The third empirical perturbation size proposed by Knoll [11] is given by Eq. (19). $\epsilon =\frac{1}{n\Vert \text{v}\Vert }{\displaystyle \sum _{i=1}^{n}b\left|y_i\right|+b}$ (19) 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 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.

Fig. 4Full-screen View

The effect of perturbation size in the MFNK method

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 and 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 and Fig. 6, the solutions resulting from MFNK and Picard match well.

Fig. 5Full-screen View

(Color online) Convergence history for k_eff

Fig. 6Full-screen View

(Color online) Convergence history for maximum fuel temperature

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]. The overall convergence performances of MFNK, Picard, and Picard with relaxation are summarized in Table 2 and 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.

As shown in 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 demonstrates that varying the boron concentration does not affect the convergence behavior. As shown in 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 and 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 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.

Fig. 7Full-screen View

(Color online) Convergence behavior of the residual norm for VERA Problem 6 with power perturbation: (a) 80% power level; (b) 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].

Transient simulation of control rod ejection

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. 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). ${\Sigma }_{\text{s},\text{g}\to \text{g}}={\Sigma }_{\text{t},\text{g}}-{\Sigma }_{\text{a},\text{g}}-{\displaystyle \sum _{\begin{array}{l}{g}^{\prime }=1\\ g^{\prime }\ne g\end{array}}^G{\Sigma }_{\text{s},\text{g}\to {\text{g}}^{\prime }}}$ (20) where Σ_t, Σ_a, and Σ_s are the total, absorption, and scattering cross-sections, respectively.

Fig. 8Full-screen View

(Color online) Radial layout of the quadrant core in NEACRP benchmark

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]. 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. 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], 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. The relation between thermal properties, such as thermal conductivity and specific heat capacity, and temperature is derived from [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). $\Sigma ={\Sigma }_0+{\displaystyle \sum _{x={\rho }_{\text{m},\sqrt{{\text{T}}_{\text{f}}},\text{C},{\text{T}}_{\text{m}}}}\left(\frac{\partial \Sigma }{\partial x}\right)}(x-x_0)$ (21) 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.

A comparison of steady-state and transient solutions for cases A2 and B2 is given in Table 5. The reference results listed in Table 5 were derived from [23]. As shown in 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.

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]. 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), 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], Fig. 10 (b) plots only the solution of the HNET. From Fig. 10 (b), an increase in the average Doppler temperature owing to power increments is observed, which behaves as expected.

Fig. 9Full-screen View

(Color online) Transient responses in case A2: (a) normalized power; (b) reactivity

Fig. 10Full-screen View

(Color online) Transient responses in case A2: (a) coolant exit temperature; (b) core averaged Doppler temperature

As shown in Fig. 11 and 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). $e=\frac{calc-ref}{ref}\times 100\%$ (22) where calc represents the values calculated by HNET, and ref is the reference value derived from PARCS [24].

Fig. 11Full-screen View

(Color online) Axial relative power for case A2: (a) time at 0.0 s; (b) time at 0.1 s; (c) time at 5.0 s

Fig. 12Full-screen View

(Color online) Relative error of assembly normalized power for case A2 in quarter core geometry (%): (a) time at 0.0 s; (b) time at 0.1 s; (c) time at 5.0 s

The normalized power in the pin and assembly levels are calculated using Eqs. (23) and (24),. ${\overline{P}}_{\text{pin}}=\frac{P_{\text{pin}}{N}_{pin}}{P_{\text{core}}}$ (23) ${\overline{P}}_{\text{assem}}=\frac{P_{\text{assem}}{N}_{\text{assem}}}{P_{\text{core}}}$ (24) 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. 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 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. The above results prove that the evolution of the power distribution can be effectively captured when the reactivity is being inserted.

Fig. 13Full-screen View

(Color online) Pin normalized power distribution for case A2 in quarter core geometry: (a) time at 0.0 s; (b) time at 0.1 s

Finally, because the thermal-hydraulics feedback is provided at the assembly level for the NEACRP benchmark, 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), 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.

Fig. 14Full-screen View

(Color online) Local temperatures and change in temperatures in the assemblies with the highest and lowest powers: (a) local temperatures; (b) change in temperatures between 0 s and 5 s