Thermal hydraulics models

The governing equation for the solid temperature is as follows: (1)where ρ, c_p and λ represent the density, specific heat, and thermal conductivity, respectively. t denotes the time, and T is the temperature. Subscript s and g represent solid and gas, respectively. q_s is the heat source of the solid. q_gs is the heat transferred between gas and solid, and it is calculated by Eq. (2) as below: (2)where A_gs and h_gs represent the heat transfer area and convective heat transfer coefficient between the solid and gas, respectively; V_mesh is the mesh volume.

DAYU3D code has incorporated various empirical correlations for the thermal conductivity and specific heat of commonly used materials in HTGRs [17], such as the pebble bed, matrix graphite of fuel pebble, reflector, carbon bricks, concrete, etc. In the pebble bed region, the Zehner-Schlunder formula [18], the Robold formula [19], and the Barthels formula [20] are employed to calculate the pebble-bed effective thermal conductivity.

Radiation heat transfer is also considered in DAYU3D. In order to simplify the calculation, the anisotropic effective thermal conductivity is utilized to simulate the radiation heat transfer, which is a proven practice in the TINTE program [21]. Moreover, DAYU3D incorporates a more accurate heat radiation model based on the view factor [22].

As for the gas flow and gas temperature calculation, the porous media model is employed. The mass, momentum and energy conservation equations are as follows: (3) (4) (5)where φ is the porosity, u represents the superficial velocity, S is the mass source, p denotes the gas pressure, W represents the resistance coefficient, and g is the gravity acceleration. The subscript ‘g’ represents the gas phase. q_gs is calculated by Eq. (2). q_src denotes the energy carried by the mass source and is calculated by the following equation: (6)where the subscript ‘src’ represents the mass source.

The properties of helium are calculated by KTA formulas [23]. The flow resistance and convective heat transfer inside the pebble bed are calculated according to KTA standards as well [24]: (7) (8) (9)where Re, Nu and Pr represent the Reynolds number, the Nusselt number and the Prandtl number, respectively; d is the hydraulic diameter.

In the pebble-bed region, the temperature distribution inside the fuel pebble is calculated. The governing equation for fuel pebble temperature is the one-dimensional (1D) heat conduction equation in spherical coordinate: (10)where r is the radius, subscript ‘fe’ represents the fuel element, i.e. the fuel pebble, and q_fe denotes the nuclear power density inside the fuel pebble.

Since the online refueling is used in pebble-bed HTGR, there may be different batches of fuel pebbles in one porous media mesh, as shown in Fig. 1. Because different batches of fuel elements have different power, they have different surface temperatures. Therefore, radiation heat transfer would occur between different batches of fuel spheres within the same mesh. The amount of radiation heat received by a single sphere of jth batch from other batches can be calculated by [21, 25]: (11)where R is the pebble radius; ε and T denote the surface emissivity and surface temperature of fuel pebble, respectively; subscript n and j represent the n^th and j^th batch of fuel pebble, respectively; is the average surface emissivity. pn is the fraction of nth batch.

Fig. 1Full-screen View

Meshing in pebble-bed region and side reflector

The calculation of fuel particle temperature is conducted through lumped-heat-capacity method, considering the size of fuel particle is relatively small. Once the temperature of each shell layer in fuel sphere has been determined, the temperature of representative particle in each shell layer (as shown in Fig. 1) can be quickly calculated by the following equation [26]: (12)where subscript ‘fp’ represents the fuel particle, T_fe is the fuel pebble temperature in one shell layer, and α_f represents the thermal resistance between the fuel particle and the matrix graphite.

Neutronics models

The neutronics kinetics calculation for HTGR is implemented in DAYU3D code system, employing a 3D cylindrical (r-z-θ) geometry to model the coupled neutronic/thermal-hydraulic (N/TH) transient behavior.

The predictor-corrector quasi-static method (PCQS) is utilized in solving space-time neutronics kinetics equation, which is written as: (13)where is the g^th group (similarly hereinafter) neutron flux, is averaged velocity, is diffusion coefficient, is removal cross section, is scattering cross section from the g′^th group to the g^th group, is fission production cross section, is prompt fission spectrum, is concentration of the i^th group delayed precursor, is decay constant of the i^th precursor group, and is production fraction of the i^th precursor group.

The equation includes six groups of delayed neutron precursor groups. The quasi-static approach decomposes the neutron flux into two components: (14)where N(t) and Ti(t) are amplitude functions, representing the rapidly varying global population of flux and precursor concentration, respectively. and are shape functions describing the spatial distribution of flux and precursor concentration, which changes slowly over time. In order to ensure the uniqueness of the decomposition, the normalization condition is imposed on the shape function by using the adjoint flux: (15)The adjoint flux is obtained by solving the eigenvalue adjoint diffusion equation.

PCQS consists of predictor and corrector phase, which solve shape and amplitude function, respectively. In predictor phase, the space-time equation is temporally discretized into large steps and is solved by implicit scheme, as shown in Eq. (16). (16)Then, by using the normalization condition of Eq. (15), the shape function is obtained as follows: (17)In corrector phase, by substituting Eq. (14) and Eq. (15) into Eq. (13), the point-kinetics equation for amplitude function is obtained as: (18)where the point-kinetics parameters are defined by the following equation: (19)In order to precisely capture the rapid varying of the amplitude function, the point-kinetics equation is discretized and solved over the minor time steps, which also requires the simultaneous updating of point-kinetics parameters. Considering the inclusion of the shape function in the expression of point-kinetics parameters as is shown in Eq. (19), it is necessary to interpolate the shape functions obtained over large steps into the nested minor steps. Fig. 2 summarizes the time-step division and solution process of PCQS.

Fig. 2Full-screen View

Schematic of the time-step division and solution process of PCQS

Compared with the traditional implicit method used in TINTE, PCQS enables the use of large time steps for solving high-dimensional shape equation and minor time steps for low-dimensional amplitude equation, which remarkably improves the computation efficiency while maintaining the accuracy.

Besides, the neutronics calculation module of DAYU3D employs two-level meshing system. The coarse-level mesh is used to characterize the spatial distribution of macroscopic cross sections and cover some of mesh grids for solid temperature computation, which are user-defined as neutronics zone. The fine-level mesh is generated by imposing a refined division on the coarse-level mesh and is used to spatially discretize the space-time kinetics equation by finite-volume method.

The multi-group macroscopic cross section, corresponding to each coarse-level grid, is in the representation of four-orders polynomial expansions parameterized by temperatures of fuel particles and graphite moderators, which is provided as follows: (20)where, T_fuel and T_mod is the temperatures of fuel particle and graphite moderator; and is the basis temperatures of fuel particle and graphite moderator; is the macroscopic cross section at basis state; en and fn is the higher order expansion coefficients of parameters.

The polynomial expansions are generated by the PANGU program [27, 28], which is developed for HTGR physics design. During the coupling calculation of neutronics and thermal-hydraulics, the macroscopic cross section will be dynamically generated by fitting the real-time temperatures into the polynomial expansions.

Code framework

The DAYU3D code employs a modular design to facilitate software development and code maintenance. The code structure and modular composition is shown in Fig. 3. The top-level modules of the code include the input module, calculation module, and output module. The input and output modules primarily handle the XML-formatted input files, the restart files, and the result outputs. The calculation module manages thermal hydraulics and neutronics calculations. The thermal-hydraulic calculation comprises several sub-modules: 3D solid temperature module, multi-batch fuel pebble and particle temperature module, 3D fluid flow module, and 3D fluid temperature module. The calculation module also contains the sub-modules for initialization, matrix solving, and material property calculations.

Fig. 3Full-screen View

Diagram of DAYU3D code structure and modular composition

The computation flowchart of DAYU3D is shown in Fig. 4. The Picard iteration framework is employed to address the N/TH coupling problems. Different time steps are allowed for neutronics and thermal-hydraulics calculations to accommodate their respective time scales and numerical stability requirements. A nested Picard iteration is incorporated within the inner iteration of thermal hydraulics fields to enhance both convergence and computational efficiency.

Fig. 4Full-screen View

Computation flowchart of DAYU3D

Modelling of control rod’s continuous movement

In transient simulations, the insertion of control rods leads to drastic localized neutron flux changes, particularly near the control rod tip, known as the cusping effect. Due to the limitation of mesh division, the rod is often partially inserted into a coarse mesh in the axial direction, thereby making it necessary to homogenize mesh grid. The difficulty lies in the determination of homogenized macroscopic cross sections of the mesh grid. The simple volume-weighted scheme cannot take into consideration the drastic gradient of the local flux distribution. In order to address this issue, TINTE code introduces the ROMO model and uses Eq. (21) to calculate the weighting factor [21]. However, since the parameter S in Eq. (21) is an empirical parameter, users need to continuously adjust its value to obtain a smooth control rod worth curve [21]. (21)where, l is the insertion length of control rod. L is the axial length of mesh. refers to the absorber concentration when the rod is partially inserted. c is the absorber concentration when the rod is fully inserted. S is the interpolation factor.

In order to overcome the limitations of ROMO model, the adaptive axial splitting method is developed into DAYU3D code, which fundamentally circumvents the homogenization process. This method is summarized as follows: in the movement of the control rod, an extra splitting is imposed on the axial layer of fine mesh grids, at where the bottom of control rod is in real-time located. While, once the bottom leaves this layer, the adaptive axial merging will be imposed to restore the previously split layer, while another splitting will be implemented according to the new position. Figure 5 presents this process.

Fig. 5Full-screen View

(Color online) Schematic of adaptive axial splitting method. (a) Splitting at current step. (b) Merging and splitting at the next step

It should be noted that the adaptive axial splitting method is imposed on the fine-level computation mesh, rather than the coarse-level mesh, as is mentioned in Sect. 2.2. This means the additional axial division does not have impact on the mesh mapping of neutronics and thermal-hydraulic fields.

Rigorous radiation heat transfer model

In the initial stage of DAYU3D’s development, the gap radiation heat transfer model primarily based on the TINTE code, in which the radiation heat transfer was treated as anisotropic heat conduction. The idea of this method is to equate the radiation heat transferred in each direction to the corresponding conductive heat transfer, thereby establishing the relationship between thermal conductivities in different directions [17]. Therefore, the principal components of thermal conductivity (λ₁, λ₂, λ₃) can be determined. The schematic of TINTE’s heat radiation model is shown in Fig. 6. However, this method introduces an assumption based on heat radiation between two infinite surfaces when calculating λ₁ [21], as shown in Eq. (22). This assumption may lead to significant computational errors in wider gaps, such as the cavity at the top of the reactor. (22)where, λ₁ refers to one of the principal components of thermal conductivity; T₁ and T₂ are the temperatures of two surfaces, respectively; d indicates the width of the gap; σ refers to the Stefan-Boltzmann constant; εs is the system emissivity of two surfaces.

Fig. 6Full-screen View

Schematic of TINTE’s heat radiation model

To address the problem of TINTE’s heat radiation model, a more accurate method based on 3D view factor and the net radiation method is being developed in DAYU3D. Due to the complex geometry in HTGRs, a numerical integration method is employed to compute view factors between surface elements. The calculated view factors agree well with the reference results from the Computational Fluid Dynamics (CFD) program [22].

Improved break mass flow model

In the TINTE code, the mass flow model for break conditions does not account for the effects of pipe length and inlet pressure, leading to an overestimation of the flow rate. To address this issue, a new set of formulas developed by Dou et al. [29] are incorporated into DAYU3D to consider the impact of different pipe lengths and inlet pressures on critical flow. The critical mass flow rate is calculated by Eq. (23). The numerical comparisons with CFD results indicates that these formulas provide accurate calculations [29]. It will effectively improve the accuracy of break flow simulations in the DAYU3D code. (23)where, A is the cross-section area of break; κ refers to the adiabatic index; p₀ and ρ₀ are the pressure and density of gas at the inlet of pipe; f is a correction factor calculated as follows:

If L < 2 m: (24)If 2 m≤L≤10 m: (25)If L > 10 m: (26)where L refers to the length of the pipe.

Efficient methods for solving flow field in HTGRs

In traditional TINTE code, the fluid pressure equation is typically solved using TINTE matrix ordering and Gauss elimination [17]. However, this fixed solution strategy may be inefficient under complex conditions. To address this, the DAYU3D code introduces an efficient flow solver based on advanced matrix ordering, symbolic factorization, and block matrix techniques.

In the DAYU3D matrix ordering method, the longest path in one-dimensional (1D) flow regions is prioritized and placed at the front of the matrix. This matrix ordering method is termed as “One-Dimensional Flow-path Priority Ordering” (ODFPO). The ODFPO reduces the number of non-zero fill-ins during the Gaussian elimination, and makes the reordered matrix exhibit a block-structured characteristic [30]. Therefore, the ODFPO method can be effectively combined with symbolic factorization and block matrix solving, respectively enhancing the computational efficiency for 2D and 3D simulations. The computational time of DAYU3D is significantly reduced compared to that of TINTE [30].

Global and parallel calculation of multi-batch fuel temperature

In the calculation of multi-batch fuel temperature, the TINTE code adopts a batch-by-batch and grid-by-grid iteration approach, which results in low computational efficiency. In the DAYU3D code, the temperatures of all batches of fuel pebbles within one grid are solved simultaneously, to avoid the iteration between batches and improve the calculation efficiency. Additionally, the DAYU3D code supports parallel computation of fuel pebble temperatures across all pebble-bed meshes, so the efficiency can be further enhanced.

Constructed cases

The verification matrix of thermal-hydraulic calculation is shown in Table 1. First, verification was carried out for individual calculation modules covering 1D, 2D, and 3D conditions, different materials, diverse flow regions, distinct boundary conditions, and both steady state and transient states. After the verification of individual modules, the verification of some coupled modules was conducted. Reference solutions for the verification include analytical solutions, as well as the calculated results by the TINTE code using the same input parameters. The verification results of DAYU3D are summarized in Table 2. The relative deviation of DAYU3D compared to the reference solution is below 1%, and mostly below 0.2%.

Among all the cases, a fully coupled example is presented in more detail to demonstrate DAYU3D’s accuracy. It is a simplified test case simulating the steady-state condition of a pebble-bed HTGR with 15 batches of fuel elements. The simulation results were compared against those from TINTE, as shown in Fig. 7. The solid temperature, helium temperature, and particle temperature exhibit good agreement with the results from TINTE.

Fig. 7Full-screen View

(Color online) Relative deviation between DAYU3D and TINTE results. (a) Solid temperature. (b) Core gas temperature. (c) Central particle temperature of the 1st batch pebbles. (d) Central particle temperature of the 15th batch pebbles

SANA experiments

SANA experiments were conducted in 1990s in Germany with the objective of illustrating the mechanisms of heat transfer in the pebble bed [31]. It consists of a pebble bed, some heating elements, insulating materials, and a pressure vessel. Some thermocouples are installed in the experimental facility to measure the temperature.

In SANA experiments, a comprehensive range of conditions have been tested [32, 33]. It includes not only 2D heating conditions but also 3D heating conditions, making it highly suitable for validating 3D thermal-hydraulic codes for HTGRs. The DAYU3D code has performed calculations for all the conditions involving helium filled and 60 mm graphite pebbles. Some example cases are provided to demonstrate the accuracy of DAYU3D.

The comparison between the calculated results of the DAYU3D code and the experimental data for the 2D conditions of SANA is shown in Fig. 8. For the 2D conditions including full-length heating, upper-section heating, and lower-section heating of the central heating element, the results of DAYU3D accurately reflect the temperature distribution within the pebble bed.

Fig. 8Full-screen View

Comparison of DAYU3D to SANA in 2D conditions. (a) Central full-length heating: 10 kW. (b) Central full-length heating: 30 kW. (c) Central upper-section heating: 20 kW. (d) Central lower-section heating: 20 kW

In 3D conditions, taking the cases of central 0 kW with radial 3×9 kW (0 + 3×9 kW) heating power and central 10 kW with radial 3×5 kW (10 + 3×5 kW) heating power as examples, a comparison between DAYU3D and SANA experimental data is shown in Fig. 9 and Fig. 10. The results illustrate that the DAYU3D code accurately calculates the 3D temperature distribution caused by the radial heating elements, showing good agreement with the experimental data.

Fig. 9Full-screen View

(Color online) Comparison of DAYU3D to SANA at 0 + 3×9 kW heating power. (a) Sectors near the radial heating elements (0°<θ<8.11°). (b) Sector far from the radial heating elements (57°<θ<60°)

Fig. 10Full-screen View

(Color online) Comparison of DAYU3D to SANA at 10 + 3×5 kW heating power. (a) Sectors near the radial heating elements (0°<θ<8.11°). (b) Sector far from the radial heating elements (57°<θ<60°)

3D HTR-PM case

Figure 11 and Fig. 12 shows the 3D calculation model of HTGR, with the meshing and material settings. The helium inlet and outlet of the core are located only in the first sector. It is assumed that the control rod channels are cooled by helium, whereas the absorber ball channels remain uncooled. The coolant mass flow rate is 96 kg/s, with the temperature at the core inlet set to 250 ℃ and the core outlet pressure maintained at 7 MPa. All the control rods are inserted at z = 275 cm, while the absorber balls are not dropped.

Fig. 11Full-screen View

(Color online) Meshing of the reactor inside pressure vessel in r-θ plane

Fig. 12Full-screen View

(Color online) Meshing and material setting in the first sector (0<θ<6°) in r-z plane. (Material setting: 1--Core; 2--Water panel; 6,8,12,13,18,24,34,36,42--Reflector; 7--Material between pebble bed and reflector; 9--Cold helium channel; 10--Helium inlet; 11,15,19,20,21,25,28,29,37--Gap; 14,45,46--Carbon brick; 16,47,48--Reactor pressure vessel; 22--Core top cavity; 23--Cold helium plenum; 26--Hot helium plenum; 27--Core bottom cavity; 31--Pressure vessel bottom cavity; 32--Cavity concrete wall; 35,41,52,53--Metal internals and core barrel; 38--Air boundary; 40--Leakage flow; 43--Bottom metallic supporter flow passage; 55--Helium outlet; 81--Small absorber sphere channel)

The solid temperature distribution within the pebble bed is shown in Fig. 13. As helium flows from top to bottom of the pebble bed, the temperature at the bottom of the core is higher, while the temperature at the top is lower. Additionally, due to the influence of the control rod and absorber ball channels, the temperature at the periphery of the pebble bed is higher near θ = 0°, θ = 60° and θ = 120°, while the temperature in other sectors is lower. The solid temperature near 180° is generally lower than that near 0°, primarily due to the influence of the helium inlet in the core.

Fig. 13Full-screen View

(Color online) Solid temperature distribution within the pebble bed. (a) Distribution in r-z plane of the first sector (0°<θ<6°). (b) Distribution in r-θ plane at z = 192.5 cm. (c) Distribution in r-θ plane at z = 522.5 cm. (d) Distribution in r-θ plane at z = 907.5 cm

The temperature distribution on the sidewall of pressure vessel is shown in Fig. 14. The vessel temperature distribution exhibits two characteristics. First, the temperature is lower in the upper region and higher in the lower region. Second, the temperature is higher near the helium inlet and lower farther from it. The first characteristic is mainly influenced by the temperature distribution in the pebble bed region, where the upper core temperature is lower and the lower core temperature is higher. The second characteristic is mainly because the 250 ℃ helium entering the reactor pressure vessel heats the vessel. Therefore, the helium inlet affects not only the temperature distribution of the pressure vessel but also that of the pebble bed region.

Fig. 14Full-screen View

(Color online) Temperature distribution on the sidewall of pressure vessel

Performance analysis

Typical test cases on HTR-PM, including steady state and depressurized loss of forced cooling accident (DLOFC), are calculated to examine the computational efficiency of the DAYU3D code. In the DLOFC case, the coolant pressure is assumed to drop instantaneously to 1 atm at the onset of the accident. To ensure a fair comparison, all the codes employ identical model setups, and all the cases are executed sequentially without using parallel acceleration. The calculation time of TINTE, MGT-3D and DAYU3D are compared in Table 3. The total calculation time of DAYU3D is reduced by about 90% compared to TINTE in the 2D cases. As for the 3D cases, the time reduction is also remarkable compared with MGT-3D. The efficiency advantage of the DAYU3D code can be attributed to two main factors. First, at the algorithm level, DAYU3D employs more advanced linear system solvers, such as the efficient flow field solving method. Second, at the implementation level, extensive optimizations have been applied to the data structures and memory management, reducing the overhead associated with dynamic memory allocation and improving overall computational performance.