**1. Introduction**

The Chinese Spallation Neutron Source (CSNS) generates intensive neutrons due to the impact of 1.6 GeV high-energy pronton on heavy metal targets. Thermal and cold neutrons after deceleration can be used to study the atomic structure and motion of certain objects [1,2]. The current beam power of a CSNS target station is 100 kW; the phase II target is plans to upgrade this to 500 kW. As the core component of the target station, the moderator slows down the leaking neutrons generated in the spall target for neutron scattering experiments [3]. The Decouple Poison Hydrogen Moderator (DPHM) is taken as the research object owing to its complex internal structure. At the same time, the existence of a poisoned plate inside the container leads to uneven distribution of flow and heat transfer.

During the operation of the target station, special working conditions should be considered, such as power off mode, distortion mode and refrigerator failure, etc., which will lead to a sharp increase in the thermal load of the hydrogen circulation system and a rise in temperature and pressure. In order to protect the cryogenic equipment, liquid hydrogen must be discharged in an emergency. In the above process, the ambient pressure

**Citation:** Tong, J.; Zhu, L.; Lu, Y.; Liang, T.; Lu, Y.; Wang, S.; Yu, C.; Dong, S.; Tan, H. Heat Transfer Analysis in Supercritical Hydrogen of Decoupled Poisoned Hydrogen Moderator with Non-Uniform Heat Source of Chinese Spallation Neutron Source. *Energies* **2021**, *14*, 4547. https://doi.org/10.3390/en14154547

Academic Editor: Vladislav A. Sadykov

Received: 15 June 2021 Accepted: 24 July 2021 Published: 27 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

decreases from 1.5 MPa to 0.1 MPa, and the liquid hydrogen in the DPHM changes from a supercritical state to supercooled state [4,5]. Due to the transcritical process, the change of pressure and temperature has a strong influence on the thermophysical properties of liquid hydrogen, for instance density, viscosity and thermal conductivity, which all make the flow heat transfer in the container more complicated [6,7]. In addition, thermal deposition increases significantly with the enhancement of proton beam power at the target station, thus leading to the further rise of the overall temperature of DPHM. To ensure the stable operation of the system, the average temperature of liquid hydrogen and the temperature difference between the inlet and outlet should be maintained in certain extent.

The heat transfer characteristics of supercritical fluid have been investigated extensively, especially the investigation of in tubes have become a focal point of research content [8–11]. Considering the influence of different factors on heat transfer, Wang et al. [12] found that the buoyancy effect under the non-uniform heat source can effectively alleviate the degree of Heat Transfer Deterioration (HTD). Meanwhile, a modified turbulent model, which can accurately predict the influence of semicircular heating condition on the flow and heat transfer, was also proposed. Zhu et al. [13] conducted a numerical study about the effects of gravity on the heat transfer performance of supercritical CO2 flowing within a vertical tube. They found that the effect of gravity on heat transfer is pronounced and closely related to the variations of thermophysical properties, particularly with low mass flux condition. The experiments on turbulent heat transfer via supercritical CO2 in a vertical tube were carried out by Kim et al. [14]. It is indicated that the distribution of wall temperature, which had a noticeable peak value when the wall heat flux was moderated and the mass flux was low, varied with buoyancy effect and flow acceleration. Nevertheless, it seems that issues such as the heat transfer characteristics of supercritical liquid hydrogen have not received sufficient attention. Youn et al. [15] obtained the variable rules of thermophysical properties of supercritical hydrogen by changing inlet conditions, for example the temperature, pressure, mass flow rate, etc. The evaluations of numerous correlations for heat transfer to supercritical hydrogen flowing turbulently in circular tubes were presented by Locke and Landrum [16]. Furthermore, In comparision with other correlations of supercritical hydrogen, by modifying the relevant parameters of correlations about non-hydrogen supercritical fluid with variable properties, this method can be applied to liquid hydrogen and obtain more accurate prediction results. James [17] described the entire validation process for a model of the heat transfer coefficient and the expected content required to complete the validation. By utilizing the above model, the convection heat transfer coefficient between the supercritical, cryogenic hydrogen and the moderator vessel walls was verified. Xie and Zhang [18–20] have carried out a number of studies on enhanced heat transfer of supercritical liquid hydrogen. Specifically, the rib structure was not only added to the wall, but also included some surface grooves and bulges, which were conducive to the enhancement of cooling performance and reducing the influence of the non-uniform distribution on heat transfer. The numerical simulation method of convective heat transfer is also important, which determines the reliability and accuracy of the simulation results. Mosavati et al. [21,22] proposed a novel calculation method for convective heat transfer in a closed cavity, specifically for the distribution factors using backward Monte Carlo method. Moreover, the effects of Rayleigh number, temperature ratio, radiation conductivity and other parameters on heat flux were studied, and the results were compared and analyzed. They conclude that lowering the temperature ratio will make the heat source surface temperature distribution smoother and lead to greater radiation heat transfer flux. At the same time, by increasing the Rayleigh number, the convective heat transfer flux can be significantly improved, and the heat flux on the heat source surface becomes uneven.

Therefore, it is necessary to explore the heat flow characteristics of liquid hydrogen in DPHM under the conditions of changing pressure, heat source and mass flow rate. The study reported here was undertaken to validate the accuracy of the supercritical model by taking the liquid hydrogen in the moderator as the research object. Meanwhile, the corresponding heat sources under different beam power of DPHM were investigated based on the neutron physical-thermal coupling method. Finally, the Computational Fluid Dynamics (CFD) method was used to simulate the flow and heat transfer characteristics of liquid hydrogen by changing the boundary conditions of DPHM.

#### **2. Mathematical and Physical Model**

#### *2.1. Physical Model*

Due to the internal complexity of the DPHM model, it was simplified for the convenience of numerical simulation analysis, as shown in Figure 1. The main structure of DPHM is made of the inlet-pipe of hydrogen, the outlet-pipe of hydrogen, container and poisoned plate, in which the structure of the pipeline is coaxial multi-layer casing. The distance between the exit of the inlet-pipe of hydrogen and the bottom surface of the container is 5 mm. After the liquid hydrogen flows along inlet-pipe of hydrogen for a distance, the poisoned plate separates the flow of it into two uneven sides, which are injected vertically to the bottom surface of the inner cavity under the action of pressure difference to form a circular jet for cooling [23]. The pressure generated by the impact forces the liquid hydrogen to flow radially along the wall, and then, due to the limitation of the internal structure, it is concentrated upward at the neck of the container and finally discharged by the hydrogen return tube. Due to the thermal deposition caused by neutron collisions, liquid hydrogen is mainly heated by pipelines, poisoned plate and container during flow.

**Figure 1.** The Computational Model of DPHM: (**a**) Model Diagram; (**b**) Detailed Size (mm).

The MCNPX is a universal Monte Carlo transport program developed by Los Alamos National Laboratory in the United States, which has been proved to be able to accurately simulate the scattering reaction of high-energy particles inside the moderator [24]. In the current work, the external coupling method was used to modify the common parameters of MCPNX 2.5 and CFX 11.0 software [25,26], the heat source distribution of the moderator calculated by the former was taken as the thermal boundary condition of the input of the latter.

### *2.2. Meshing*

The mesh of DPHM generated by commercial software ICEM 14.0 is shown in Figure 2, where the unstructured mesh with strong adaptability was selected to discretize the fluid domain. The overall quality of the grid is detailed; the maximum and average skewness are 0.652 and 0.223 respectively, and the average aspect ratio is 3.124, which proves the reliability of grid division. Considering the influence of boundary layer on the main flow area, the local mesh of areas with large curvature and complex flow, such as inlet, outlet and wall boundary area, was encrypted so as to accurately capture the flow characteristics. In order to accurately predict the turbulent flow field, the height of the first layer of the boundary layer grid is set to be small enough to meet the criterion of *y* <sup>+</sup> value close to 1. Five inflation layers are divided in the sensitive region, so the estimation of the thermal gradients can be improved. The selection height of the first layer grid in the calculation model is calculated by the following Equation [27],

$$Re = \frac{\rho u L}{\mu} \tag{1}$$

$$C\_f = 0.058 Re^{-0.2} \tag{2}$$

$$
\pi\_w = \frac{1}{2} \mathbb{C}\_f \rho u^2 \tag{3}
$$

$$
\mathcal{U}\_{\pi} = \sqrt{\frac{\overline{\mathbf{r}\_{\text{w}}}}{\rho}} \tag{4}
$$

$$y = \frac{y^+ \mu}{l L\_\tau \rho} \tag{5}$$

where *Re*, *ρ*, *u*, *L*, *μ*, *Cf*, *τ<sup>w</sup>* and *U<sup>τ</sup>* are the dimensionless number, density, bulk fluid velocity, characteristic length, dynamic viscosity of fluid, boundary layer friction coefficient, boundary layer shear stress and shear velocity. Finally, after calculation, the value of the first layer grid height *y* <sup>+</sup> is 2 mm.

**Figure 2.** Generation mesh in the computational domain.

#### *2.3. Governing Equations*

The liquid hydrogen jet cooling in the moderator belongs to high-speed flow, so the influence of gravity on heat transfer flow can be ignored. In the transcritical process, the pressure has a great influence on the density of liquid hydrogen [28]. Therefore, this article assumes that the liquid hydrogen in the moderator is a compressible fluid, and only the steady flow is studied. The specific governing equations are as follows:

The continuity equation [29]:

$$\nabla \cdot (\rho \mathbf{u}) = 0 \tag{6}$$

The momentum equation:

$$\begin{cases} \nabla \cdot (\rho \iota \mathbf{u}) = -\frac{\partial p}{\partial x} + \frac{\partial \mathbf{r}\_{xx}}{\partial x} + \frac{\partial \mathbf{r}\_{yx}}{\partial y} + \frac{\partial \mathbf{r}\_{zz}}{\partial z} \\ \nabla \cdot (\rho v \mathbf{u}) = -\frac{\partial p}{\partial x} + \frac{\partial \mathbf{r}\_{xy}}{\partial x} + \frac{\partial \mathbf{r}\_{yy}}{\partial y} + \frac{\partial \mathbf{r}\_{zy}}{\partial z} \\ \nabla \cdot (\rho w \mathbf{u}) = -\frac{\partial p}{\partial x} + \frac{\partial \mathbf{r}\_{yx}}{\partial x} + \frac{\partial \mathbf{r}\_{yz}}{\partial y} + \frac{\partial \mathbf{r}\_{zz}}{\partial z} \end{cases} \tag{7}$$

The energy equation [30]:

$$\nabla \cdot \left[ (\rho E + P)\mathbf{u} \right] = \nabla \cdot \left( \lambda\_{eff} \nabla T - \left( \mathbf{r}\_{eff} \cdot \mathbf{u} \right) \right) + S\_h \tag{8}$$

where **u**, *τ*, *p*, *P*, *E*, *T*, *λeff*, *τeff* and *Sh* denote the velocity vector, stress tensor, bulk fluid pressure, static pressure, total energy, bulk fluid temperature, effective thermal conductivity, effective stress tensor and volumetric heat source, respectively.

In this paper, the steady-state heat specifies and the material is isotropic, so the Fourier law was employed to describe the heat conduction within the solid wall [31].

$$\nabla \cdot (\lambda \nabla T) = 0 \tag{9}$$

where *λ* is the thermal conductivity of the fluid. The shear stress transport (SST) k-ω turbulence model, which takes into account the transport of turbulent shear stress and has a high accuracy in predicting the complex flow, is used to solve the three-dimensional flow of liquid hydrogen.

#### *2.4. Boundary Conditions*

The boundary conditions for each structure of the computational model, including fluid domain and solid domain, should be set before the simulation begins. Moreover, the non-uniform heat source of DPHM is imported into CFX software by compiling User Defined Function (UDF).



**Table 1.** Comparison of heat deposition calculation results between MCNPX and CFX (W).

#### *2.5. Verification of Grid Independence*

In this paper, the grid independence of DPHM is verified under the conditions of beam power of 100 kW, pressure of 15 bar and inlet flow of 30 g/s. Five groups of grids were selected for comparison, and the three parameters of maximum temperature, maximum pressure of liquid hydrogen and maximum container temperature were detected. The results are shown in Figure 3. It can be seen that the temperature variation tends to be stable when the number of grid elements exceeds 1.68 million, indicating that the independent solution can be obtained for the grid number at this condition. The Grid Convergence Index (GCI) was used to quantify the grid independence [32]. The GCI34 for fine, and medium grids was 1.76%. The GCI45 for medium and coarse grids was 4.61%. The value of GCI45/(*r*pGCI34) was 1.03, which was approximately 1 and indicates that the solutions were well within the asymptotic range of convergence. Therefore, considering the computational efficiency and accuracy, the fourth group of grids was finally employed for calculation.

**Figure 3.** Grid independence validation of DPHM.

#### *2.6. Verification of Physical Parameters*

Pseudo-critical temperature refers to the temperature corresponding to the maximum specific heat capacity of fluid at constant pressure in a supercritical state (Critical pressure and temperature of liquid hydrogen are 12.9 bar and 33.15 K, respectively). In the vicinity of this temperature range, the thermophysical properties of the fluid will change dramatically, thus affecting the flow and heat transfer. As shown in Figure 4a, the density *ρ* and dynamic viscosity *μ* decrease with the gradual increase of temperature under 15 bar but converge to a constant value as the temperature is far beyond the critical value. Meanwhile, the change in specific heat *cp* is most significant, which is further reflected in Figure 4b. There is a prominent phenomenon called the "thermal spike phenomenon", which is not conducive to the stability of heat transfer. This primary peak in *λ* and *cp* of hydrogen disappears with the augment of pressure due to the corresponding Pseudo-critical temperature increasing roughly from 32 to 34 K [33].

**Figure 4.** Variations of thermophysical properties of liquid hydrogen under different pressure: (**a**) Pressure = 15 bar; (**b**) Specific Heat (different pressure).

Accurate prediction of the thermal properties of hydrogen is crucial to the reliability of the simulation results. Based on equations of state (EoSs) theory [34,35], the physical parameters of hydrogen, in which the α parameter and pressure–volume–temperature (PVT) relationships have an important influence, can be accurately predicted. In order to more accurately simulate the flow heat transfer of non-ideal fluids under supercritical conditions, under the premise of considering the calculation accuracy and efficiency, this paper adopts the Peng–Robinson (PR-EOS) equation of state [31], which is universal in engineering, to calculate the physical parameters of working fluids:

$$p = \frac{RT}{v\_{PR} - b} - \frac{a(T)}{v\_{PR}(v\_{PR} + b) + b(v\_{PR} - b)}\tag{10}$$

where

$$a(T) = 0.45724 \frac{R^2 T\_c^2}{p\_c} \left[ 1 + \kappa \left( 1 - \sqrt{(T/T\_c)} \right) \right]^2 \tag{11}$$

$$b = 0.0778 \frac{RT\_c}{p\_c} \tag{12}$$

$$
\kappa = 0.37464 + 1.54226\omega - 0.26992\omega^2 \tag{13}
$$

where *R*, *vPR*, *Tc*, *pc*, and *ω* are the molecular gas constant, the molar volume, the critical temperature, the critical pressure and the Pitzer acentric factor.

Whereas the prediction results of liquid density obtained by the PR-EOS formula have poor accuracy, especially when the pressure is close to the critical region and in a small temperature range. Khashayar [36] evaluated 11 correlations for predicting hydrogen density and found that the modified Redlich–Kwong EoS by Mathias and Copeman (RKMC) [37] was widely applicable to predict liquid hydrogen properties in various ranges. The deviation between the experimental results and the RKMC EoS equation was small in the temperature range of 14~32 K, which proved the rationality and accuracy of the RKMC EoS equation. Since the liquid hydrogen in this paper belongs to the supercritical state, and the temperature change is located in the above range, the RKMC EOS equation is used for the subsequent calculation in this paper. The specific equation are as follows:

$$p = \frac{RT}{v - b} - \frac{a\_c a \left(T\_r\right)}{v(v + b)}\tag{14}$$

$$b = 0.08664 \frac{RT\_c}{p\_c} \tag{15}$$

$$a\_c = 0.42748 \frac{R^2 T\_c^2}{p\_c} \tag{16}$$

where

$$\mathfrak{a} = \begin{cases} \left[ 1 + c\_1 \left( 1 - \sqrt{T\_r} \right) + c\_2 \left( 1 - \sqrt{T\_r} \right)^2 + c\_3 \left( 1 - \sqrt{T\_r} \right)^3 \right], T\_r \le 1\\ \left[ 1 + c\_1 \left( 1 - \sqrt{T\_r} \right) \right]^2, T\_r > 1 \end{cases} \tag{17}$$

where *ac* is the critical attractive parameter, *b* is the molar co-volume, *α* is the temperature dependence function (alpha function), *c1*,*c2*,*c3* are the coefficients of the Mathias and Copeman alpha function, and *Tr* is the reference temperature.

Figure 5 shows the density and viscosity distributions at temperatures of 15~40 K at pressures of 5, 10 and 15 bar. In order to verify the accuracy of the above model, a comparison was made with the data of GASPAK [38], the results showing that the two sets of data were well in agreement and that the error was within the allowable range.

**Figure 5.** Comparison of simulation results with GSPAK data of thermophysical properties: (**a**) Density; (**b**) Viscosity.
