Numerical Investigation into the Natural Convection of Cryogenic Supercritical Helium in a Spherical Enclosure

: As an ideal pressurized gas, helium, especially supercritical helium, has been widely used in the pressurization system of various launch vehicles and spacecraft. This work mainly focuses on the natural convection of cryogenic supercritical helium in a spherical enclosure. Firstly, a three-dimensional numerical model is established and veriﬁed with experimental data. Then, the effects of inﬂation pressure and heating power on the ﬂow and heat transfer characteristics are simulated. At the same time, the relationship between the Rayleigh number and Nusselt number is studied in detail. Finally, an improved natural convection heat transfer correlation modiﬁed by introducing the density ratio is obtained. The results show that the increase of the inﬂation pressure in the cavity is helpful to enhance the natural convection heat transfer of the cryogenic supercritical helium, and the temperature distribution in the cavity tends to be more uniform when the inﬂation pressure in the cavity increases. As to the improved natural convection heat transfer correlation, the average error between the simulation results and the calculated values is approximately 8%, which can better describe the natural convection heat transfer of cryogenic supercritical helium in the spherical enclosure.


Introduction
Helium is an ideal pressurized gas and has been widely applicated in the pressurization system of various launch vehicles and spacecraft. However, the method of liquid helium pressurization has the drawbacks that it may cause a liquid-vapor two-phase state, which leads to an unstable pressurization process. To overcome this problem, supercritical helium (p ≥ 0.227 Mpa, T ≥ 5.13 K) has been adopted in the rocket pressurization system with its stability and reliability. During the pressurization procedure, the heat transfer between the hot heater and cold supercritical helium dominates by natural convection; hence, a deep understanding of the natural convection heat transfer mechanism of supercritical helium inside an enclosure is vital for the heater electric power control as well as supercritical helium pressure management in the supercritical helium pressurization system. Natural convection in enclosures with internal bodies has been extensively studied due to its range of engineering applications, such as pipes, electronics thermal management, etc. Bishop et al. [1] experimentally investigated the natural convection of air between two isothermal concentric spheres. They developed two Nusselt number-Grashof number correlations for four different diameter ratios. McCoy et al. [2] tested the natural convection heat transfer rates and temperature fields between isothermal vertical cylinders and isothermal spherical enclosures; the temperature fields were found to be similar to concentric sphere flows. Powe and Warrington [3] studied the natural convection heat transfer phenomena between a body (spheres, cylinders and cubes) and its spherical enclosure; a universal heat transfer correlation was developed with an average deviation of less than 12%. Nasrine and Ghaddar [4] numerically investigated the natural convection from a uniformly heated horizontal cylinder in a large rectangular enclosure using the spectral element method. Kim et al. [5] studied the natural convection between a cold outer square enclosure and a hot inner circular cylinder through numerical simulation. Lee et al. [6] studied the internal natural convection heat transfer between a cubic enclosure and spherical internal body and pointed out a critical Rayleigh number beyond which the Nusselt number decreases as the temperature difference increases. Seo et al. [7,8] numerically investigated the natural convection in a long cold rectangular enclosure with a hot cylinder (circular or elliptical) inside employing the immersed boundary method. They addressed the fact that the heat transfer characteristics were affected by the radius of the cylinders. Cho et al. [9,10] studied the natural convection in a cold square enclosure with a vertical array of one circular cylinder and one elliptical cylinder or two elliptical cylinders. In the one circular cylinder and one elliptical cylinder case, the transition of flow regime from unsteady to steady state depends on the variation of aspect ratios of the elliptical cylinder, while, for the two elliptical cases, the transition of the flow regime from steady state to unsteady state depends on the aspect ratio at Ra = 10 6 .
However, compared with ordinary fluids, the thermophysical properties of supercritical fluids change dramatically with temperature and pressure, which makes their flow and heat transfer characteristics more complex [11]. The research on the heat transfer characteristics of supercritical fluid can be traced back to the 19th century, and a literature review was published by Hall [12]. In recent years, with the development of the energy industry, more and more researchers are devoted to the study of heat transfer characteristics of supercritical fluid. Muller and Estevez [13] studied the mass transfer phenomenon caused by the natural convection of supercritical fluid, and pointed out that, in most cases, the heat and mass transfer phenomenon of supercritical fluid is basically the same. Accarcy and Raspo [14] proposed a three-dimensional finite volume numerical method for predicting the natural convection of supercritical fluid in a cavity, and discussed the influence of the linearization of the van der Waals equation of state on convergence. Based on the Redilch-Kwong equation of state, Teymourtash et al. [15] derived a new formula for calculating the coefficient of thermal expansion and studied the natural convection of supercritical carbon dioxide on a vertical plate under variable heat flux. Shan et al. [16] studied the flow pattern and temperature distribution of supercritical carbon dioxide natural convection in the cold wall side of a nuclear reactor through experiments and simulation.
Up to now, the research on the natural convection of supercritical helium is very limited. Hial et al. [17,18] studied the natural convection heat transfer of supercritical helium in a spherical cavity. They measured the natural convection heat transfer coefficients in the temperature range of 4.2-25 K and the pressure range of 0.3-3.5 MPa. Then, they put forward a criterion formula for calculating the heat transfer coefficient of supercritical helium natural convection in the spherical cavity and found the boiling-like phenomenon in the near-critical region [19]. Hideaki and Hisanao [20] studied the natural convection heat transfer of supercritical helium on the vertical plane with or without grooves. The results show that the influence of grooves on the heat transfer coefficient is in the range of ±20%. Irie et al. [21] studied the natural convection heat transfer of supercritical helium in a near-critical state.
According to the literature survey, most of the research focused on the natural convection of normal fluids. The natural convection heat transfer characteristics of su-percritical helium in cryogenic environment is still not clear. Moreover, the lack of effective heat transfer correlation hinders further research on supercritical helium pressurization technology. In order to further understand the natural convection heat transfer characteristics of supercritical helium in spherical enclosures, a three-dimensional steady-state numerical simulation model will be established in this work based on some experimental data. Then, the flow and heat transfer characteristics of supercritical helium in the spherical enclosure will also be studied, aiming at obtaining a high-precision heat transfer correlation.

Physical Model
The physical model of the natural convection of cryogenic supercritical helium in a spherical enclosure is shown in Figure 1. The spherical cavity is immersed in liquid nitrogen with a radius of 200 mm, and a cylindrical heating rod with a height of 100 mm and a diameter of 8 mm is placed vertically at the center of the spherical cavity. The center of the spherical cavity is the origin of the coordinates. In the initial state, the pressure and the temperature of the helium in the spherical cavity exceed its critical point, and as a result are in a supercritical state. Then, the heating rod located in the center of the spherical cavity begins to heat the supercritical helium with a certain power, and the natural convection occurs due to the density of supercritical helium changes. As the spherical cavity is immersed in liquid nitrogen, the wall temperature of the spherical cavity remains constant. When the temperature of the heating rod remains constant, the natural convection of supercritical helium in the spherical cavity is in a stable state.

Governing Equations
The governing equation of natural convection in the spherical cavity is as follows: Momentum conservation equation: Energy conservation equation: where λ eff is the effective thermal conductivity, its value is the sum of the thermal conductivity λ and turbulent thermal conductivity λ t . "ρ" is the density of supercritical helium and varies with temperature [22]. For the natural convection of cryogenic supercritical helium caused by temperature difference, the flow state can be distinguished by the dimensionless Rayleigh number, which is defined as: where g is the gravitational acceleration, β is thermal expansion coefficient of supercritical helium, ν is the kinematic viscosity, α is the thermal diffusion coefficient, L is the characteristic length, and T h , T c are the temperature of the heating rod and the inner wall of the spherical cavity, respectively. The intensity of the natural convection heat transfer of supercritical helium in the spherical cavity is measured by the dimensionless Nusselt number, which is defined as: In which n is direction normal to the wall. In order to identify the intensity of natural convection under different pressures, the dimensionless local Rayleigh number Ra L is introduced, which is defined as: where T L is the local temperature, β L is the corresponding thermal expansion coefficient at local temperature, the characteristic length L is the difference between the radius of the cavity and the heating rod, and ν L and α L are the physical properties corresponding to the local temperature. Generally, radiation heat transfer is also a main form of heat transfer in cryogenic systems. In this work, the surface-to-surface (S2S) radiation model is used to calculate the radiative heat transfer between the heating rod and the inner wall of the spherical cavity. It is assumed in the S2S radiation model that all surfaces are gray bodies with diffuse reflection. The emissivity and absorptivity of the gray body are independent of the wavelength. In addition, according to the Kirchhoff's law, the emissivity of a radiation surface is equal to the absorptivity. The energy flux leaving a given surface is composed of directly emitted and reflected energy. The reflected energy flux depends on the incident energy flux from the surroundings, which can be expressed in terms of the energy flux leaving all other surfaces. Hence, the energy reflected from surface k is: where q out,k is the energy flux leaving the surface, ε k is the emissivity, σ is the Boltzmann's constant, and q in,k is the energy flux incident on the source from the surroundings. However, after considering the radiant heat transfer between the heating rod and the inner wall of the spherical cavity, the simulation analysis found that even if the heating rod power was increased to 100 W, the proportion of radiation heat transfer in the total heat transfer is still less than 1%. Therefore, in order to save computing resources and time, the radiation heat transfer between the heating rod and the wall of the spherical cavity is not considered in the following simulation study.

Numerical Approaches and Boundary Conditions
The structured grid of the computational domain was generated by ICEM. ICEM is the abbreviation of "The integrated Computer Engineering and Manufacturing", which is a piece of software used to generate mesh [23]. The near wall meshes as well as the near heating rod meshes were refined, as shown in Figure 2. After the verification of grid independence, 817,320 grids were selected for the following simulation calculations. ANSYS FLUENT is used in the simulation process. For calculation purposes, it is assumed that the heat generated by the heating rod is evenly distributed on its surface when the heating rod is heated at a certain power. Therefore, the surface of the heating rod is a constant heat flux density boundary in this model, and the heat flux density is related to the heating power. The inner surface of the spherical cavity is an isothermal boundary condition, and the wall temperature is equal to the liquid nitrogen temperature, 77 K.

Model Validation
In order to verify the accuracy of the established model, a supercritical helium natural convection test platform is built and will be introduced in the following section.

Experiment Setup
The supercritical helium natural convection test platform is mainly composed of a cryogenic Dewar, a spherical cavity, a liquid nitrogen filling system, a helium filling system, a heat transfer test unit and a data acquisition system, as shown in Figure 3. During the experiment, the spherical cavity was placed inside the Dewar, which contained liquid nitrogen. The outer side of the Dewar was well insulated to ensure that the liquid nitrogen was at its saturation temperature, 77 K, which can provide a constant wall temperature boundary condition for the spherical cavity. Then, the spherical cavity was filled with helium gas at a certain pressure (above its critical pressure). After the helium in the cavity was cooled to 77 K, both the temperature and pressure of the helium were above their critical points; hence, the helium in the cavity was in a supercritical state. When the temperature and pressure were stable, the heating rod was turned on at a certain heating power and a temperature difference between the heating rod and the spherical cavity wall was generated. As a result, the natural convection of supercritical helium in the spherical enclosure was achieved. The volume of the spherical cavity inside the Dewar is approximately 30 L, which is used to store supercritical helium. The diameter of the spherical cavity is 200 mm, the maximum working pressure is 20 MPa, the leakage rate of the experimental spherical cavity is less than 1 × 10 −8 Pa·m 3 /s. In addition, a heater with adjustable power of 0~150 W is set in the center of the spherical cavity, and the length of the heater is 100 mm, the diameter of the heater is 8 mm. There are 9 temperature measuring points arranged in the test platform; eight are located in the spherical cavity, of which five are used to measure the supercritical helium temperature (T1-T5), two are used to measure the wall temperature in the cavity (T6, T7), and one is used to measure the temperature of heater (T8). The other one is located in the liquid nitrogen Dewar to measure the liquid nitrogen temperature (T9).

Model Validation
In order to verify the accuracy of the model, the steady-state temperature at T1 and T8 under different heating powers was selected for comparison, as shown in Figure 4. It can be seen from the figure that the difference between the simulation temperature and the experimental temperature at point T1 is within 5%, and the simulated temperature of T8 at a lower temperature is also in good agreement with the experimental value. When the pressure in the cavity is low and the temperature of the heating rod is high, there is a certain error between the simulation and the experimental temperature for T8. Through the analysis, the possible reasons are as follows: (1) In the experiment, the heat rod is wrapped with a heat shrink tube, and the heat generated by the heat rod is transferred to the heat shrink tube in the form of heat conduction. While in the simulation process, it is considered that the heating rod is directly in contact with the supercritical helium and the heat flux is distributed evenly. In fact, the thermal conductivity of the heat shrinkable tube is much higher than that of supercritical helium. (2) In the experiment, there is a certain contact thermal resistance between the platinum resistance and the heating rod, and the temperature gradient near the heating rod is very large. Based on the above two reasons, the measured temperature will be lower than the actual temperature. In general, the simulation temperatures are in good agreement with the experimental values, so the model is considered to be effective in calculating the natural convection of supercritical helium in the spherical enclosure.

Results and Discussions
Compared with normal fluids, one of the important characteristics of supercritical fluid is that its thermophysical properties change dramatically with the temperature and pressure. In order to calculate the flow and heat transfer characteristics of cryogenic supercritical helium in natural convection more accurately, the thermophysical properties of cryogenic supercritical helium are analyzed firstly. The data are calculated with the supercritical helium real gas model from NIST (National Institute of standards and Tech-nology), as shown in Figure 5. It can be seen from the figure that the density and specific heat capacity at the constant pressure of cryogenic supercritical helium decrease rapidly with the increase in temperature, while the kinematic viscosity increases slowly at a lower temperature. When the temperature continues to increase, the change rate of density and specific heat capacity at constant pressure gradually slows down, while the kinematic viscosity increases rapidly, while the thermal conductivity changes almost linearly with temperature. In addition, through the comprehensive analysis of density, specific heat capacity at constant pressure and thermal conductivity, it can be obtained that the thermal diffusion coefficient has the same change trend as the kinematic viscosity.

Effects of Inflation Pressure and Heating Power on the Natural Convection
Firstly, the flow state of the supercritical helium in a spherical cavity is analyzed. Taking the inflation pressure of 2 MPa and the heating power of 10 W as an example, the streamline is shown in Figure 6. Since the heating rod is located in the center of the cavity, there is a "flame like" region with the highest velocity just above the heating rod, as shown in Figure 6a. When the heated supercritical helium reaches the top of the cavity, it then flows downward along the wall and forms a symmetrical vortex. The streamline direction of the natural convection in the thermal boundary layer near the heating rod is vertical upward, while it presents a spiral upward shape in the intermediate region. This can also be seen from Figure 6b,c-that is, the supercritical helium flows from the cavity wall with a low temperature to the heating rod with a high temperature in the horizontal direction, while, in the vertical direction, it flows from the bottom of the cavity to the top under the influence of gravity.  Figure 7 shows the temperature distribution of the supercritical helium under different inflation pressures when the heating power is 10 W (Section z = 0). Since only the temperature of the supercritical helium near the heating rod is relatively high, the temperature in the rest of the regions is between 77-78 K. Therefore, the temperature range is set at 77-78 K to show the temperature distribution more clearly, and the black parts are the areas where temperature is higher than 78 K. It can be seen from the figure that there is an obvious temperature gradient under all working conditions, and the "flame like" region with the highest temperature is surrounded by lower temperature supercritical helium. With the increase in the inflation pressure, the temperature gradient in the upper part of the cavity decreases gradually, while the supercritical helium region at 77 K in the bottom of the cavity increases. This shows that the natural convection heat transfer in the cavity is enhanced when the inflation pressure increases.  Figure 8 shows the temperature distribution of the supercritical helium with the heating power at 100 W (Section z = 0). The temperature range in Figure 8 is set at 77-85 K. It can be seen from the figure that the temperature distribution is similar to that when the heating power is 10 W, but the overall temperature in the spherical cavity is higher. When the heating power increases to 100 W, the supercritical helium in the temperature range of 77-78 K is almost in the bottom of the cavity, and the temperature gradient in the spherical cavity is more obvious. Figure 9 shows the velocity contour and streamline of the supercritical helium at section z = 0 under different conditions. It can be seen from the figure that the velocity in the narrow "flame like" region is the highest under all conditions. With the increase in the inflation pressure, the velocity in the narrow regions decreases gradually, and the velocity distribution becomes more uniform. One of the main reasons for this is that when the inflation pressure increases, the supercritical helium density increases. The other reason is that with the increase in the inflation pressure, the heat transfer between the wall and the heating rod is enhanced and the surface temperature of the heating rod is reduced, so the supercritical helium density near the heating rod is also increased. Therefore, the natural convection velocity decreases gradually. In addition, it can also be seen that the velocity boundary layer near the heating rod becomes thinner with the increase in the inflation pressure. From the streamline point of view, the flow state of supercritical helium under different inflation pressures is very similar. In the upper part of the cavity, there are two symmetrical vortices, while in the lower part, supercritical helium flows from the middle to the wall of the cavity. However, with the increase in pressure, there is a tendency for vortex formation on both sides of the wall in the middle of the spherical cavity.   Figure 10 shows the distribution of local Rayleigh numbers Ra L under different inflation pressures when the heating power is 10 W. The contour maps of Ra L are drawn in the range of 10 10 to 10 11 , and it is colored red when greater than 10 11 . It can be seen from the figure that as the inflation pressure increases, the overall local Rayleigh number in the spherical cavity increases significantly, indicating that there is a significant improvement of the natural convection with the increase in the inflation pressure. The reason for this may be that although the temperature difference between the cold and hot ends decreases as the inflation pressure increases, the decrease in viscosity and thermal diffusion coefficient is more obvious. Therefore, the viscosity and thermal diffusion coefficient become the leading factor that affect Ra L . Compared with the temperature distribution in Figure 7, it can be found that Ra L also decreases gradually from the region near the heating rod to the cavity wall with obvious delamination.  Figure 11 shows the variation of the local Nusselt number on the surface of the heating rod in the vertical direction. Overall, the local Nusselt number increases with the increase in the inflation pressure. In other words, the heat transfer performance of the whole region is improved as the inflation pressure increases. In addition, under different pressures, the local Nusselt number shows the same change patterns from the bottom to the top of the cavity, that is, it first increases sharply with the rise in height, then decreases and tends to be stable, and then decreases sharply near the top of the heating rod. This is because the upper and lower surfaces of the heating rod are perpendicular to the direction of gravity, and the heat generated by the upper and lower surfaces can only be diffused in the form of heat conduction, which makes the temperature at both ends of the heating rod increase and the local Nusselt number decrease. Besides, because of the natural convection, supercritical helium with a lower density and higher temperature accumulates above the heating rod, resulting in the decrease of the local Nusselt number. Therefore, the local Nusselt number in the middle part of the heating rod first decreases and then tends to be stable.

Natural Convection Heat Transfer Correlation of Supercritical Helium in a Spherical Enclosure
As mentioned above, the Rayleigh number is used to measure the change of flow pattern in natural convection, and the Nusselt number represents the intensity of convective heat transfer. Both are important dimensionless parameters to investigate natural convection heat transfer. Figure 12 shows the variation of Ra and Nu with the heating power of the heating rod. As can be seen, Ra and Nu decrease with the increase in the heating power when the pressure of supercritical helium is less than 4 MPa. When the pressure is more than 6 MPa, the variation of Ra will change; it first increases and then decreases with the increase in the heating power. This is mainly caused by the change in physical properties of the supercritical helium, as shown in Figure 5. In addition, it can be found from the figure that when the pressure is relatively high, Ra shows the opposite trend under the conditions of high and low heating power. This is mainly because when the pressure is high, the influence of viscosity and thermal diffusion coefficient will be more significant only at high temperatures. The analysis of Nu is similar to that of Ra. The main factors affecting Nu are heat flux, thermal conductivity, and temperature difference between the hot and cold ends. When the heating power increases, the temperature of the heating rod and its surroundings increases, and the thermal conductivity also gradually increases. Therefore, when the pressure is relatively low, Nu will first keep almost unchanged and then decrease with the increase in heating power, while as the pressure increases to a certain value, Nu will increase first and then keep almost unchanged. The variations of Ra and Nu under different heating powers have been analyzed, and the relationship between Ra and Nu will now be further explored. In general, it is considered that the specific heat capacity, thermal conductivity, thermal expansion coefficient and other thermophysical parameters are constant when discussing the relationship between Ra and Nu. As a result, they satisfy the classical heat transfer correlation of Nu = C·Ra a . Figure 13 shows the relationship between Ra and Nu with the thermophysical properties of helium calculated at 77 K and 0.1 MPa. As can be seen, the Ra and Nu of supercritical helium natural convection no longer satisfy the exponential relationship due to the special thermophysical properties of helium. In fact, the Ra calculated by constant thermophysical properties is thousands of times different from that calculated by variable thermophysical properties. Therefore, the thermophysical properties of supercritical helium must be considered when studying the correlation of the natural convection heat transfer of supercritical helium in a spherical enclosure.  Figure 14a shows the relationship between Ra and Nu in Figure 12. It can be seen from the figure that Nu increases with the increase in Ra when the Ra is relatively low (lower than 10 11 ), and they basically meet the power function relationship. However, when Ra increases to larger than 10 11 , there is no strict one-to-one correspondence between Ra and Nu. In order to better explain the special relationship between Ra and Nu when Ra exceeds a certain critical value, the influence of physical properties should be fully considered. Therefore, an improved heat transfer correlation modified by density ratio (∆ρ/ ρ c ) is proposed based on the power function correlation, where ∆ρ = ρ w − ρ s , ρ w and ρ s represents the helium density corresponding to heating rod surface temperature and the cavity wall temperature, respectively; ρ c is the helium density corresponding to the initial temperature. The final heat transfer correlation is as follows: The comparison between the simulation results and the calculated values by the improved heat transfer correlation are shown in Figure 14b. It can be seen from the figure that Nu calculated by heat transfer correlation is in good agreement with the simulation results. In the calculation range, the error can be controlled within 10%, and the average error is approximately 8%. It shows that the heat transfer correlation has high accuracy in calculating the natural convection heat transfer of cryogenic supercritical helium in a spherical enclosure.

Conclusions
In this work, the natural convection of supercritical helium in a spherical enclosure is numerically investigated. Firstly, a three-dimensional steady-state numerical model is established and verified by experimental data. Then, the flow and heat transfer characteristics of supercritical helium in the spherical enclosure are studied, and the effects of different parameters on the heat transfer characteristics are also analyzed. Finally, the natural convection heat transfer correlation of supercritical helium in spherical enclosures is studied. The main conclusions are as follows: (1) The increase in the inflation pressure in the cavity is helpful to enhance the natural convection heat transfer of the cryogenic supercritical helium. With the increase in the inflation pressure in the cavity, the velocity and temperature gradient in the boundary layer near the heating rod decrease, and the temperature distribution in the cavity tends to be more uniform. (2) The influence of the inflation pressure in the cavity on the flow pattern is relatively small, but the velocity in the "flame like" region above the heating rod decreases with the increase in the inflation pressure. (3) The classical heat transfer correlation is modified by introducing the density ratio as Nu = 0.40(∆ρ/ρ c ) 0.12 Ra 0.29 (1 × 10 10 < Ra < 5 × 10 11 ). The average error between the simulation results and the calculated values is approximately 8%, which can better describe the natural convection heat transfer of cryogenic supercritical helium in the spherical enclosure.