Three-dimensional neutronics/thermal-hydraulics model of reactor core

The reactor core is a key component of a nuclear reactor system, and studying its internal neutronic/thermal-hydraulic behavior is crucial for analyzing the safety and performance of the core. Therefore, the 3D coupled neutronics-thermal-hydraulics behavior was considered for analyzing the reactor core.

The steady-state multigroup neutron diffusion equation was used to describe the neutronic behavior within the reactor core. (1)where ϕg represents the neutron flux corresponding to g-th energy group, ; here, we omitted the independent variable spatial vector r, Dg, , are the neutron diffusion coefficient, macroscopic removal cross-section, and macroscopic neutron production represents the macroscopic scattering cross section from group g′ to group g. χg is the neutron fission spectrum and k_eff is the effective multiplication factor. When considering the neutronic/thermal-hydraulic coupling, the diffusion coefficient and cross-sections can be rewritten as the functions of temperature: D=D(T_u, T_m, D_m) and , where T_u and T_m are the fuel temperature and moderator temperature, respectively; D_m is the moderator density. Because of the enormous computational effort required for full-core detailed neutron diffusion simulations and because the average power and temperature distribution of pin cells or assemblies are usually sufficient for engineering purposes, this study employs the nodal method [56] to solve the reactor core neutron diffusion equation. For the steady-state thermal–hydraulic calculations, the single-channel method [57, 58] was used for each moderator channel. The hydraulic calculation of each channel was achieved based on the one-dimensional mass, momentum, and energy conservation as follows: (2)where G_m, f, g, ρ_m, H_m, p_c, D_e, A_z, P_h are coolant mass flux, channel friction factor, gravitational acceleration, coolant enthalpy, core pressure, flow channel equivalent diameter, flow channel area, and heated perimeter, respectively; θ is flow channel inclination angle; q_l is fuel rod line power density, which is calculated as: (3)where q_v is the volume power density; F_u is the proportion of heat released from the pellets to the total heat release, which is generally taken as 97.4% in engineering; E_f is the energy released from a single fission; A_u is the cross-sectional area of the fuel pellets. Under the single-phase assumption, the pressure drop and fluid enthalpies of the reactor were calculated by considering the inlet temperature, inlet pressure, and inlet mass flux. After the enthalpy and pressure were calculated, the temperature and density distributions along the channel were obtained using the IF97 thermophysical properties table [59].

The thermal calculations mainly involve the radial one-dimensional thermal conductivity calculation of the fuel rod, including the fuel pellet, gap, and cladding. A fuel pellet can be described as: (4)where r represents the position of the radial coordinates, T_u and κ_u are the temperature and heat conductivity of the fuel pellet, respectively, and for the gap and cladding, this equation is generally rewritten as (5)where T_g and κ_g denote the general temperature and heat conductivity of the gap and cladding, respectively.

This is the coupled neutronics/thermal-hydraulics model that is used in this study. As an example, a 3D NEACRP reactor core was used in this study, whose design references the literature [60], as shown in Fig. 1. The neutronic behavior of the reactor core was simulated using the open-source software ADPRES [60]. The 3D thermal–hydraulic process was calculated using our self-developed code, which was coupled with ADPRES for neutronics/thermal-hydraulics coupling. In the calculation, vacuum and reflective boundaries were used for the core boundary and internal symmetry plane in neutron diffusion, respectively. The inlet temperature, inlet pressure, and inlet mass flux were provided for the thermal-hydraulic calculation. In the radial direction, each assembly was divided into two ×2 nodals, and the axial direction was divided into 18 layers. The thermal-hydraulic and physical calculations used a similar mesh configuration.

Fig. 1Full-screen View

(Color online) Geometry of NEACRP reactor [60], including (a) layer 1, (b) layer 2, (c) layer 3, and (d) axial distribution

Three-dimensional CFD calculation of steam generator

Another carefully considered component is the steam generator (SG), which is an important component of a nuclear reactor’s primary circuit and must always be optimized for design. To this end, this study used the three-dimensional CFD method to analyze the detailed thermal-hydraulic behavior inside the steam generator. A schematic of the steam generator is shown in Fig. 2, in which the red and blue parts represent the primary and secondary sides, respectively. To reduce computational cost, an equivalent U-tube model was used instead of a complete steam generator [61], as shown in Fig. 3. In this simulation, the governing equations for the primary side are expressed as follows: (6)where , u, p, k, ε, T, C_p, μ, μ_t, κ are density, velocity, pressure, turbulent kinetic energy, turbulent dissipation rate, temperature, specific heat, viscosity, turbulent viscosity, and thermal conductivity, respectively; Q is heat source; F is the fluid force term; , = 1.92, Cμ = 0.09, σ_k = 1.0 and σε = 1.3 are model constants. Neglecting the stress, strain, and movement of the heat transfer tubes, the governing equation of the heat transfer tube wall can be expressed as (7)where κ_w, c_pw, and T_w denote the thermal conductivity, specific heat capacity, and temperature of the heat transfer tube wall, respectively.

Fig. 2Full-screen View

(Color online) Schematic of steam generator (Red and blue represent the primary and secondary sides, respectively)

Fig. 3Full-screen View

Schematic of equivalently SG model

For the secondary side, refer to Sun et al. [62]. The liquid is solved as a continuous phase (denoted by l), whereas the vapor is solved as the dispersed phase (denoted by v). The ensemble-averaged conservation equations in the two-fluid model for each phase can be written as (8)where i, j={l, v} for the liquid and vapor phases, αi, ρi, u_i, pi, , λi, and Ti are the volume fraction, density, velocity, pressure, effective dynamic viscosity, thermal conductivity, and temperature of i-th phase, respectively. represents the mass transfer rate from phase i to j, F_i,j represents the force acting on phase i relative to phase j, represents the total heat transfer from phase i to j.

On the secondary side, the vapor phase is regarded as a laminar flow and the liquid phase is described using the turbulence model: (9)where d_b represents the bubble detachment diameter, C_μb = 0.6 is a constant. All the other variables follow the same definitions as those in Eq. (6). The heated boiling process on the secondary side was simulated using a thermal phase transition model that considers the mass, momentum, and energy transfer across the interfaces between the two phases. The details of these models can be found in Refs. [62, 63].

In this study, the detailed thermal-hydraulic behavior within SG was calculated using the finite-element method (FEM) with the COMSOL software [64]. Figure 4 illustrates the mesh configuration of SG. After grid-independent verification, 1478800 meshes were used for the calculation. To calculate the stability, the inlet boundaries were applied to the temperature and velocity distributions, whereas the outlet boundary was applied to the pressure distribution.

Fig. 4Full-screen View

Mesh configuration of equivalently SG model

The detailed CFD model used for the steam described above. For further applications, further research can be conducted, such as full-geometry and porous-medium approaches [65].

System model

Other components, including the pipelines, pressurizer, and main pump, were simulated using system models. Because the detailed temperature and pressure distributions inside the pipeline are not as important for the entire core system circuit, a one-dimensional thermal-hydraulics model is used for the pipelines: (10) (11)where the subscript p indicates the physical quantity of the pipeline, θ represents the pipe tilt angle, and d represents the pipe inner diameter (m); f_p represents the Darcy friction factor in the pipes, which can be expressed as follows [66]: (12)where Re_D is the pipe Reynolds number, Re_c is the critical Reynolds number, is the pipe roughness (m).

In the discrete heat-transfer pipe model, an energy equation must be incorporated into the continuity and momentum equations. The energy equations can be expressed as (13)where q_p denotes the volumetric heat loss (W/m³). Using the above equations, the thermal-hydraulic process within different types of pipelines can be simulated, and each pipe is divided into 100 grids for calculation. The pressurizer is taken as a pressure containment boundary of the primary circuit, through which the circuit pressure will be adjusted to the set pressure, such as 15.5 MPa [67]. The effect of the main pump is reflected in the change in pressure (head) in the primary circuit, and its model relationship can be obtained from the head characteristic curve [68]. Therefore, in this system model, these two components are approximated using zero-dimensional (0D) models.

multi-scale multi-physics primary circuit model

The multiscale multi-physics primary circuit system investigated in this study is illustrated in Fig. 5, which includes the 3D reactor, 3D steam generator, 1D pipelines, 0D pressurizer, and main pump. For plug-and-play, all components were calculated using independent codes. In this study, these calculations are referred to as the full-order model (FOM).

Fig. 5Full-screen View

(Color online) Schematic of equivalently SG model

Normally, this multiscale model can be established by coupling the above component models with inlet and outlet values. The calculation is implemented by sequentially solving the models of each component, transferring their output values to the input values of the downstream components, and iterating until the loop converges. Through these coupling processes, the primary circuit system is simulated in a multiscale manner, and the behavior of the entire system is obtained by aggregating the behavior of its individual components. However, the necessity for coupling and transferring data across various software platforms presents significant challenges. These challenges are compounded by the substantial computational time required for integrated simulations. In addition, different methods for solving different components cause greater difficulties in this coupling. As mentioned earlier, in solving the steam generator using the CFD method, to calculate the stability, the inlet boundaries are applied to the temperature and velocity distributions, while the outlet boundary is applied to the pressure distribution. Under these conditions, the classical forward sweep method is not applicable because the upstream pressure is unknown, whereas the downstream pressure is known. To ensure the computational stability of the circuit model, a forward sweep-backward sweep method is proposed, in which the velocity and temperature can be calculated using the forward sweep method (e.g., clockwise), whereas the pressure should be calculated using the backward sweep method (e.g., counterclockwise). Because of this requirement, existing system programs must be significantly modified, and the universality of the programs will be reduced. For this phenomenon, a surrogate model that directly obtains the outlet physical quantity from the inlet conditions without considering the numerical stability of the sweep process may be a good alternative. The inlet pressure of the steam generator can be obtained by setting different outlet pressures and establishing a mapping relationship. To this end, the reduced-order method (ROM) is effective for constructing surrogate models.

In the following section, we present a plug-and-play multiscale reduced-order model for the primary circuit system, in which surrogate models for all components are established and connected. The external file coupling and integration processes in multiscale coupling calculations can be eliminated.

Creation of snapshots

Snapshots of each component were obtained by conducting system or 3D calculations within the corresponding component state parameter ranges. For each component, we used the inlet or outlet parameters as the state parameters according to the requirements of the independent component calculation. For the reactor core, the core power P_n, inlet mass flux G_in, inlet temperature T_in, and inlet pressure p_in are considered as its state parameters, which can be written as μ_c = {T_in, G_in, p_in, P_n}. For pipelines, the inlet velocity U_in, inlet temperature T_in, and inlet pressure p_in are considered as their state parameters, which can be written as μ_p = {T_in, U_in, p_in}. The above state parameters are consistent in both the snapshot calculations and reduced-order circuit calculations.

Specifically, for the steam generator, two types of state parameters are used. Because the primary and secondary circuits are strongly coupled under normal operating conditions, each set of SG primary circuit state parameters corresponds to a specific secondary circuit state. If the primary circuit state parameters are specified, the secondary circuit state can be uniquely determined. Therefore, only the primary loop-state parameters must be specified under normal operating conditions. For snapshot calculation, the inlet temperature T_in, velocity U_in and outlet pressure p_out are considered for computing stability [61], which can be written as . For reduced-order calculation in circuit coupling, to achieve sequential connection of different components between circuits, the inlet parameters are still used as the state parameters of the reduced-order surrogate model, namely μ_s = {T_in, U_in, p_in}. Because of the absence of a CFD calculation process in the reduced-order surrogate model, no stability issue is detected with this state parameter setting.

After defining the state parameters of each component, the ranges of state parameters for each component to calculate snapshots should be determined. This was achieved by conducting one-dimensional system calculations on typical circuits under different core power conditions and summarizing the inlet and outlet parameter ranges of each component.

Calculation of reduced bases

One of the main processes of the reduced-order surrogate model is to extract reduced bases from snapshots. In this study, the proper orthogonal decomposition (POD) method was used because of its maturity and reliability. For self-consistency, the basic properties of the POD were introduced. Taking the reactor core as an example, we consider a collection of N temperature and power snapshots corresponding to the discrete parameter set for power and temperature distributions. The snapshot for k-th parameter μ_k is defined as: (14)where u_P(μ_k) and u_T(μ_k) represent the power and temperature distributions within the reactor core corresponding to k-th state parameter combination μ_k, respectively; M is the total mesh number, the subscript i indicates the mesh index, Ti and Pi are the coolant temperature and power distribution corresponding to ith mesh node, respectively. The correlation matrix S is then formed by computing the inner product of each pair of snapshots as follows: (15)The set of eigenvalues and the corresponding eigenvector of the matrix S_α are computed to obtain the POD bases. The set of POD bases can be determined according to (16)where P_α,i is ith order POD basis vector of the αth variable. The magnitude of i-th eigenvalue describes the relative importance of the i-th order POD basis vector P_α,i. The first N_p orders of the POD bases are used to establish the reduced-order model, which is defined by the following regulation: (17)where ϵ is set to 10^-8 to ensure efficient information. Using this approach, the collection of POD bases is defined as . The snapshots can be expressed using a linear combination of POD bases as follows: (18)where and are ith basis coefficients of the temperature and power distributions corresponding to μ_k.

These were the reduced bases calculated for the reactor cores. Using a similar approach, the reduced bases for other components can be obtained. In contrast, for the steam generator, the velocity and temperature were used as snapshots: (19)where u, v and w denote the velocity components in x, y, and zdirections, respectively. For the pipe, only axial velocity was considered. (20)where U denotes axial velocity.

Establishment of reduced-order surrogate model

The surrogate models for each component, including the computational and reconstruction models, are established in this section.

Computational surrogate models were established for reduced-order circuit calculations in which the physical quantity at the outlet of the component was calculated using the inlet variables. As in the previous section, we consider a collection of N output variable snapshots . The snapshot for k-th parameter μ_k is defined as: (21)where T_out, U_out, and p_out are the average outlet temperature, velocity, and pressure, respectively. By establishing a mapping relationship between the inlet and outlet parameters, a computational surrogate model can be established. (22)In this study, both linear and RBF interpolation methods were used. A computational model was used for system circuit calculations. This complex system calculation can be simplified by connecting various component computational surrogate models.

After circuit calculation using the computational surrogate model, the reconstruction models were used to reconstruct the one-dimensional or three-dimensional distribution of the kernel physical variables, including the temperature and power distribution of the reactor core and the temperature, velocity, and pressure distribution for SG and pipes. This reconstruction model is based on the mapping relationship between the inlet parameters μ and basis coefficients c: (23)Subsequently, the variable distribution can be reconstructed as (24)

Multi-physics multi-scale circuit reduced-order model

Similar to component surrogate models, the circuit surrogate model is divided into computational and reconstruction models. The computational model connects the inputs and outputs of the computational surrogate models for different components, as shown in Fig. 6. For instance, the output parameters of the reactor core, including the outlet temperature, flow rate, and pressure, served as the input parameters for Heat Section 1. Execution of the circuit ROM computational model begins with the reactor core. Initially, the core power was specified as required and the core inlet parameters were assumed for the calculation. Subsequently, core outlet parameters were determined using a computational surrogate model. Subsequently, the computational surrogate model for each component was executed individually and iterated until the entire circuit model converged. After the ROM computation, the detailed distributions of each component were reconstructed using the reconstruction model, as shown in Fig. 7. For each component, by taking the inlet parameters as state parameters to determine the POD coefficients , the 1D or 3D distributions can be reconstructed according to Eq. (24).

Fig. 6Full-screen View

(Color online) ROM computation model for circuit system

Fig. 7Full-screen View

(Color online) ROM reconstruction model for circuit system