## 2 Numerical Methods
### Numerical Procedure  
In the buoyancy-driven natural convection flow, the intensity of flow can be categorized by the Rayleigh number [36; 37]. When \(Ra<10^{8}\), the flow is laminar. When \(Ra>10^{10}\), the natural convection flow will transform into a turbulent flow. The definition of \(Ra\) number is shown in Equation (9),  
\[Ra=\frac{g\beta\Delta T_{l}L^{3}}{v\alpha} \tag{9}\]  
where \(g\) is the local gravitational acceleration and \(\Delta T_{l}\) is the characteristic temperature difference, which is set as 100 K in this study. \(\beta\) is the thermal expansion coefficient, \(v\) is the kinematic viscosity, and \(\alpha\) is the thermal diffusivity of the fluid. According to Equation (9), the calculated maximum Rayleigh number in this study is \(1.82\times 10^{7}\), which is less than \(10^{8}\). Therefore, the fluid field follows a laminar flow pattern.  
In order to simulate conjugate heat transfer under a subatmospheric environment, the reference pressure of the fluid domain is set in the form of Equation (10). According to the ideal gas law, under a fixed temperature the fluid density is proportional to the pressure. Therefore, the reference density of the air is set to match with the corresponding vacuum level, as shown in Equation (11).  
\[p_{r}=\eta p_{a} \tag{10}\]\[\rho_{r}=\eta\rho_{a} \tag{11}\]  
In Equations (10) and (11), \(\eta\) is the ratio of environmental pressure to the atmospheric pressure, \(p_{r}\) is the reference pressure, and \(p_{a}\) and \(\rho_{a}\) are the atmospheric pressure and air density at the atmospheric pressure, respectively.  
The commercial code STAR-CCM+ is utilized to carry out the simulations. Both separated-flow and coupled-flow solvers have been used in this study. The results from both solvers exhibit little difference. However, the separated-flow solver requires less computational resources than the coupled-flow solver. Therefore, the separated-flow solver is used. The radiation model adapts DO model. Air is assumed to be the ideal gas. The effect of gravity is also taken into account in the numerical model. In the separated-flow solver, the second-order upwind scheme is adopted for the convection term of the equations. The SIMPLE algorithm is used in the overall solution. The second-order precision is used for time discretization. The time step used in this study is calculated by Equation (12) based on the theoretical study by [PERSON] [38]:  
\[\tau=\frac{L}{U}=\frac{L^{2}}{\alpha\sqrt{RaPr}}=\frac{L}{\sqrt{\beta gL\Delta T _{l}}} \tag{12}\]  
where \(\tau\) is the time constant, \(Pr\) is Prandtl Number and \(U\) is the velocity scale. The minimum \(\tau\) value for this study according to Equation (12) is 0.183 s. Based on the results from Equation (12), the iterative time steps of \(\tau\), \(\tau/4\) and \(\tau/16\) have been selected for the convergence test [39]. It is found that when the iterative time step is around 0.18 s, during the heating process, the fluctuation of convection heat transfer coefficient is not negligible. As the time step decreases to a level between 0.04 s and 0.01 s, the numerical stability on the convection heat transfer coefficient can be observed. Hence, the time step of 0.04 s is applied in this study. The stopping criteria for each time step is set as when the residues of continuity equation and velocity equation fall below \(10^{-3}\) and residues of energy equation fall below \(10^{-6}\). The vacuum levels investigated in this study include 0.1 atm, 0.3 atm, 0.4 atm, 0.5 atm and 1.0 atm.