A modified local thermal non-equilibrium model of transient phase-change transpiration cooling for hypersonic thermal protection

Aiming to efficiently simulate the transient process of transpiration cooling with phase change and reveal the convection mechanism between fluid and porous media particles in a continuum scale, a new two-phase mixture model is developed by incorporating the local thermal non-equilibrium effect. Considering the low-pressure and high overload working conditions of hypersonic flying, the heat and mass transfer induced by capillary and inertial body forces are analyzed for sub-cooled, saturated and super-heated states of water coolant under varying saturation pressures. After the validation of the model, transient simulations for different external factors, including spatially-varied heat flux, coolant mass flux, time-dependent external pressure and aircraft acceleration are conducted. The results show that the vapor blockage patterns at the outlet are highly dependent on the injection mass flux value and the exter-nal pressure, and the reduced saturation temperature at low external pressure leads to early boiling off and vapor blockage. The motion of flying has a large influence on the cooling effect, as the inertial force could change the flow pattern of the fluid inside significantly. The comparison of the results from 2-D and 3-D simulations suggests that 3-D simulation shall be conducted for practical application of transpiration cooling, as the thermal protection efficiency may be overestimated by the 2-D results due to the assumption of an infinite width length of the porous plate.


Introduction
During the flight of a hypersonic vehicle, it is subjected to aerodynamic heating from high-speed external flows and high-temperature gases in the engine combustion chamber.A highly efficient thermal protection system is essential to ensure its safe and reliable flying [1,2].Thermal protection can be broadly classified into two categories: passive and active thermal protection.Passive thermal protection relies on materials' heat sink and ablative properties [3] or the use of high-temperature-resistant materials for insulation.Active thermal protection, on the other hand, involves cooling the high-temperature components by absorbing heat through fluid convection or creating gas/liquid films to isolate the high-temperature airflow from the structural components.This category includes methods such as regenerative cooling [4], film cooling [5], and transpiration cooling.The advantages of active thermal protection lie in its robust cooling capability and minimal influence on structural design, but it comes with the complexity of implementing these systems.
Transpiration cooling has been proposed for a few decades, where a coolant passes through a porous wall that undergoes convective heat transfer before injecting out into the high-temperature mainstream flow.This process separates the high-temperature flow from the wall and thickens the boundary layer of the hot mainstream, thereby reducing the amount of heat transferred from the high-temperature mainstream to the solid wall.Theoretically, the maximum cooling capacity can reach up to 1.4 × 10 9 W/m 2 [6].Commonly used materials for the solid wall include layered plates, sintered particle porous materials [7] and sintered wire mesh porous structures [8].When the coolant is a single-phase fluid without phase change within the porous wall, such as using various gases as the coolant, it is referred to as single-phase transpiration cooling, which is the focus of current transpiration cooling technologies [9][10][11][12].When the coolant absorbs heat that undergoes phase change within the porous wall, such as using water as the coolant, it is termed as phase-change transpiration cooling, which performs great potential due to the large phase-change latent heat of the coolant for active thermal protection [13].For numerical simulations of phase-change transpiration cooling, accurately establishing a multiphase heat and mass transfer model is of paramount importance.However, due to the complex multiphase-coupled heat transfer process in the porous media and external flow, there are still many challenges in numerical simulations, such as low and changing pressure in the hypersonic boundary layer, heat exchange in the porous structure between the solid skeletons and the gas/liquid phases, complex motion of aircraft induced body force and 3-D effects of the porous region.
Existing numerical simulation methods can be generally categorized into three approaches: Separate Phase Model (SPM), Semi-Mixed Model (SMM) and Two-Phase Mixture Model (TPMM) [14].The SPM simultaneously computes the governing equations for both liquid and gas phases, allowing for the consideration of two-fluid interactions [15,16].Researchers like Wang [17] and He [18] established SPMs to study the effects of different heat fluxes and coolant mass fluxes on physical quantities such as temperature, pressure and phase saturation in porous media flat plates.Wang et al. [19] employed a SPM to investigate the response of internal gas-liquid phase change to time-varying boundary heat flux and pressure.This kind of model reserves most of the flow information, however, it has higher computation complexity and computing load.To simplify the method, the SMM solves the mass and momentum equations for both fluids, but treats the energy equation as a two-phase mixture.Dong and Wang [20,21] established a SMM and considered modifications to the diffusion coefficient in the fluid energy equation to study the heat transfer deterioration effects caused by local high heat flux.Xin et al. [22] used a SMM to explore the influence of thermal conductivity and porosity of porous media materials on transpiration cooling effects.He et al. [23] employed a transient SMM to investigate the impact of different injection modes and periodic operations.To further reduce the computing load, focusing on only the critical physical quantities, the TPMM solves the governing equations of a single fluid, which represents the constituents of a binary mixture.Shi and Wang [24] established a TPMM in steady and phase change states under constant saturation pressure, considering the influence of various factors, such as particle diameter, porosity and thermal conductivity.Su et al. [25] introduced a modified fluid temperature as the primary variable to solve the TPMM and studied the interaction between supersonic external flow and phase-change transpiration cooling blunt body under steady state and a constant saturation temperature.Liu et al. [26] modified the fluid enthalpy under constant saturation state, and studied the effects of coolant injection mode without gravity.Hu et al. [27] modified the TPMM energy equation considering the change in the saturation temperature in twophase regions, and then investigated the transient effects of phase-change transpiration cooling under low pressure, but assumed the fluid and the solid maintained the same temperature.Cheng et al. [28] utilized the model in [27] to simulate the transpiration cooling effects of porous media with linearly changing porosity and compared it with experimental results.Chen et al. [29] considered variations in cooling properties with temperature and pressure in the TPMM without body forces like gravity, and ignored the influence of saturation temperature on the diffusion coefficient of fluid.
To investigate heat transfer phenomena within porous media at a macroscopic level using the volume-averaging method, the heat transfer model between the solid skeleton and the fluid can be categorized into two main types, known as the Local Thermal Equilibrium (LTE) [30,31] model and the Local Thermal Non-Equilibrium (LTNE) [32,33] model.When using the LTE model, it is assumed that the temperature of the solid skeleton is equal to the temperature of the fluid at the same location.Under this assumption, the heat transfer characteristics within the porous media can be described using a single energy equation.The LTNE model adopts a more realistic representation of heat transfer within porous media, as it accounts for the variations in temperature between the solid skeleton and the fluid.Therefore, it is often used to study convective heat transfer problems.Alomar et al. [34] conducted simulations comparing the gas-liquid phasechange processes within porous media using both LTE and LTNE models.Their findings indicated that the LTNE model provides a more accurate representation of heat transfer effects.The two assumptions can both be applied to the above SPM, SMM and TPMM to establish the heat and mass transfer relationships between the solid and fluid phases.
To take advantage of low computing capacity and low risk of numerical divergence, this work develops a new two-phase mixture model for phase change transpiration cooling by considering the LTNE effect.Based on the convection coefficients between porous media and fluid suggested by [32,33], the governing equations in Ref. [14] are modified into a LTNE form.The primary variable for the fluid energy equation is chosen to be the kinematic enthalpy, as previously proved beneficial for calculation by [26,27].To consider the low-pressure working condition of a hypersonic vehicle, the theoretical equations for the fluid energy in Ref. [14] are established under the sub-cooled state, saturated state and super-heated state, respectively, considering the capillary induced heat and mass transfer under varied saturation pressures and saturation temperatures.In addition, the contributions of the inertial body force to the momentum equation and the mass transfer modeling of the fluid energy equation are considered, as the effect of aircraft overload may be large.Moreover, the iterative update of the kinematic enthalpy is realized by using the relative mobility as a criterion of the fluid saturation state, as it is a good representation of liquid saturation.Then the model is verified by comparing the simulation results with previous experimental and numerical simulation data.This modified model developed based on the previous research can consider more critical factors in the simulation of transpiration cooling.Firstly, the model contains transient terms in every equation, which can provide time-dependent simulation results that are meaningful in applications compared with the steady model.Additionally, the governing equations of the model simultaneously consider local thermal non-equilibrium effects and variations in the boiling point with the saturation pressure.It can simulate the temperature change of the saturated coolant under low-pressure working conditions and temperature non-equilibrium phenomena between the solid medium and the internal fluid.Thirdly, the momentum source and energy source related to the complex variablespeed motion of the aircraft are considered to incorporate their influence to the coolant distribution.Previous models mostly ignored the effect of the body force or only considered the static ground condition.
After the validation of the new model, the transient effects of transpiration cooling under different external influence factors, including time-dependent and spatially-varied external heat flux, different injection mass flux, time-dependent external pressure and aircraft acceleration, are analyzed by comparing the critical properties for thermal protection.Among the selected external factors to be discussed, there is relatively little available information on the influence of the inertial force caused by the acceleration of a porous plate.Finally, the modified model is applied to a 3-D porous plate and the results are compared with the 2-D case.For simplicity, previous studies mostly use 2-D plates; however, the dimension of the model will influence the simulation results of thermal protection under hypersonic operating conditions.

Basic governing equations of the modified model
To develop a new formulation for two-phase flow through capillary porous media that is both physically meaningful and practically useful, Wang and Beckermann [14] converted the separate flow model (SFM) to a multiphase mixture form.
After the mixture density ρ is defined, as in Table 1, the conservation of mass in porous media can be written as: where ε is the porosity of porous media, ρ is the mixture density, and u is the Darcy velocity.
Considering the resistance to flow and gravity of the flow inside porous media, the momentum equation can be written as: where the subscripts l and k represent the variables for liquid phase and kinematic term.p, s l , µ(s l ) , ρ k , K and g represent the mixture pressure, liquid saturation, dynamic (1) viscosity of the mixture, kinematic density, permeability of porous media, and gravity, respectively.
For the conservation of energy, the heat and mass transfer were modeled under the assumption of local thermal equilibrium: where the subscripts f, s, eff and v represent the variables for fluid, solid, average effective value due to porous porosity and vapor phase, respectively.For the enthalpy related to fluid, h f , h k , h v and h l are the specific enthalpy of fluid, kinematic mixture enthalpy, specific enthalpy of vapor and specific enthalpy of liquid, respectively.h s , k eff , T f , j and Q are the specific enthalpy of solid, effective heat transfer coefficient of fluid and solid, temperature of fluid, total mass flux and volume heat source, respectively.The new variables above are defined in [14] and the constitutive relationships are shown in Table 1.

Details of the model modification
In order to conduct the numerical simulation of phase change transpiration cooling under the low-pressure and high acceleration working condition, the above TPMM is modified based on the following assumptions: (1) The saturation temperature of the fluid changes with the mixture pressure when two phases coexist.(2) The temperature of the fluid and the solid at the same location are different, and the convection heat transfer process between them is considered.(3) The liquid phase of the fluid is assumed to be incompressible, while the vapor phase is considered as an ideal gas.(4) The influence of the inertial force on the momentum and energy equations due to the motion of aircraft is considered.(5) The flow inside the porous media is assumed to be at a low velocity laminar condition. (3)

Table 1 Constitutive relationships used in the established model
Variable Formula Relative permeability of liquid phase and vapor phase, Relative mobility of liquid phase and vapor phase, Effective heat transfer coefficient of fluid, Under the above assumptions, considering the transient term and possible acceleration of porous media, the new momentum equation based on (2) can be written as: where, a represents the acceleration of aircraft that causes the inertial force of fluid.
To model the energy equations under the local thermal non-equilibrium assumption, the fluid-to-solid heat transfer coefficients [33] and the specific surface area of porous media based on geometrical considerations [32] are adopted to convert (3) to separate solid energy equation and fluid equation.The solid energy equation can be written as: where, ρ s , c ps , k seff and Qfs are the solid density, solid specific heat capacity, effective solid heat transfer coefficient and convective volume heat source.
The subscript fs means the heat transfer direction is defined as from fluid to solid, where, the subscript sf means the heat transfer direction is defined as from solid to fluid.After eliminating the contribution of solid to the transient term and diffusion term in (3), the fluid energy equation can be written as the following form, which is applicable in every state of the fluid (sub-cooled, saturated and super-heated).
where, k feff is the effective heat transfer coefficient of the fluid.
To implement the numerical simulation, Eq. ( 7) is modified into three forms, which model the fluid energy under the sub-cooled state, the saturated state and the superheated state, respectively.The primary variable for the fluid energy equation is selected as the kinematic enthalpy h k : where, l and 1 − l = v are the relative mobility of the liquid phase and the vapor phase, which represent the relative kinematic viscosity of liquid and vapor to the multiphase mixture, respectively.The detailed constitutive relationship between relative mobility and liquid saturation is in Table 1.
From the second term on the left-hand side of (7), it can be seen that using this variable can maintain the coefficient of the convection term as the mass flow rate ρ u and therefore stabilize the numerical simulation, which had been implemented well in [26,27].
Under the sub-cooled state, the liquid saturation field s l and the relative mobility field l are both the constant value one.Therefore, the mixture density, specific fluid enthalpy and kinematic enthalpy are consequently simplified to variables for the pure liquid phase: where, c pl is the specific heat capacity of water liquid, and T sat,ref is the reference tem- perature of the saturation state.
In addition, when the liquid saturation field is a uniform value one, the total mass flux j on the second term of the right-hand side of (7), including the capillarity-induced dif- fusive flux and the inertial force induced migrating flux, does not have any effect on the mass diffusivity of fluid: where, D(s l ) is the capillary diffusion coefficient, k rv is the relative permeability of vapor, ν v is the kinematic viscosity of vapor.The detailed constitutive relationships are in Table 1.
Therefore, the fluid energy equation under the sub-cooled state can be written as: The convective volume heat source Qsf is modeled by the empirical convective heat transfer coefficient: where, the specific surface area of porous media α sf = 6(1 − ε)/d p .And as suggested in [36,37]: where, d p , k l , Re l and Pr l represent the particle diameter of porous media, the heat trans- fer coefficient, the Reynolds number and the Prandtl number of liquid, respectively.
Under the saturated state, the liquid saturation 0 < s l < 1 and the relative mobil- ity 0 < l < 1 , and the fluid temperature is constrained by the saturation pressure.The relationship between the specific enthalpy of the mixture and the kinematic enthalpy of the mixture can be written as: where, the subscripts sat and fk represent the values under the saturated state and the coefficient transforming the fluid enthalpy to the kinematic enthalpy.h fg (p sat ) is the latent heat under the saturation pressure p sat .
For the implementation of the numerical simulation, the mass diffusion term and the heat diffusion term induced by mass transfer in (7) can be transformed into the (10) diffusion term of the primary variable h k .Considering the contribution of varying satu- ration enthalpy to the liquid saturation field and thermal diffusion, the coefficients used for transformation are calculated as shown in ( 17) and ( 18): Therefore, the fluid energy equation modeled for varying two-phase region temperature can be changed as followed: When the liquid phase and the vapor phase coexist, the volume heat source Qsf is established from the liquid saturation of the mixture [24] and the empirical correlation of surface boiling latent heat [38], and is presented as follows: Under the super-heated state, the liquid saturation field s l and the relative mobility field l are both the constant value zero.Similar to the sub-cooled state, the fluid energy equation under the sub-cooled state can be written as: where, the volume heat source Qsf = h sv α sf (T s − T f ).
Based on the above conservation equations of fluid energy and the constitutive relationships in Table 1, the primary variable h k can be used to represent the critical physi- cal properties under the sub-cooled state, the saturated state and the super-heated state, including the temperature of fluid, the relative mobility of liquid and the liquid saturation.
The relationship between the fluid temperature and the kinematic enthalpy is: The relationship between the relative mobility of liquid and the kinematic enthalpy is: (17) The relationship between the relative mobility of liquid and the liquid saturation is:

Boundary conditions
To solve Eqs. ( 1), (4-24), as explained in Refs.[39,40] with details, the boundary conditions of the physical models are established as follows.
At the inlet boundary Ŵ inlet , the mass flow injection and the heat convection are con- strained as:

Mass flow inlet:
Solid energy:

Fluid energy:
where, ρ c is the density of the inlet coolant, u c is the injection velocity, the convection heat transfer coefficient h c = 0.664Pr c 1 3 Re c 1 2 [41], and h kc is the kinematic enthalpy of the inlet coolant.
At the outlet boundary Ŵ external , the pressure and the heat flux are modeled as follows: Pressure outlet: Solid energy: Fluid energy: where, p external is the external pressure, and qexternal is the heat flux applied by the exter- nal environment.At the wall boundary Ŵ wall , the velocity and the heat flux are modeled as follows: No slip condition: Adiabatic wall for solid and fluid: where, n is the normal vector of the wall boundary.

Solution procedure
The modified LTNE-TPMM is solved under the basic frame of the commercial software ANSYS FLUENT, with the support of User-Defined Memory (UDM) to store the newly defined parameters, User-Defined Scalar (UDS) to implement the newly defined scalar transport equations and User-Defined Functions (UDFs) to establish boundary conditions and update criteria [29].The solution procedure of the numerical method is shown in Fig. 1.After entering a new time step, at the beginning of each iteration, the boundary condition, the relative mobility of liquid, the liquid saturation and the fluid temperature are updated according to the kinematic enthalpy of the last iteration, and the mass and momentum equations are then solved, followed by the fluid and solid energy equations.
After solving all the equations, update the physical properties by the constitutive relationships.At the end of each iteration, check whether the convergence criteria are satisfied.End the overall iteration procedure while reaching the target total simulation time.
For the transient simulation, Pressure-Implicit with Splitting of Operators (PISO) is applied for the pressure-velocity coupling scheme and the second-order upwind scheme is used for the spatial discretization of pressure, momentum, kinematic enthalpy and solid temperature.The time step is adjusted between 0.01 s and 0.1 s, and the number of iterations of each time step is adjusted between 500 and 1000, according to the convergence situation.The convergence criteria are reached when the residual for the energy equations is smaller than 10 -5 .

Model validation
To validate the transient simulation effectiveness of the modified model, simulation is conducted and the results are compared with the experimental results of Wu et al. [42].In this experiment, an alloy of titanium and aluminum was sintered into porous media and heated up under quartz lamps with a heat flux of 500 kW/m 2 .The temperature of the inlet coolant was 293.15K and the mass flow rate was 0.1306 kg/(m 2 •s), (30) corresponding to the 0.5 MPa pressure difference in the experiment.The experimental results indicated that the outer wall temperature, heated directly by the quartz lamp, fluctuates around 373.15 K, demonstrating excellent heat protection performance but with a significant fluctuation range.This was attributed to the rapid increase in the temperature of the solid wall due to the quartz lamp heating, coupled with the heat exchange with liquid water and water vapor causing temperature reduction, resulting in the pronounced fluctuation in the readings of the outer wall thermocouple.For the physical properties in the numerical simulation, the properties of the water coolant used are shown in Table 3.The porosity of the porous media is 0.23 and the particle diameter is 5 μm.From Fig. 2, it can be seen that the numerical method can obtain relatively accurate temperature rising rates and steady temperature levels.
To further validate the transient simulation effectiveness of the modified model, simulation is conducted and the results are compared with the numerical results of Chen et al. [29].Their numerical model had been compared with the steady state experiment  3, the transient value of the liquid saturation at the outlet is in good agreement with that in [29].

Physical models and properties
The physical models used in the following paper are shown in Fig. 4a and b.The water coolant with selected mass flux and temperature is injected from the bottom of the plate and the upper side is exposed to the external environment.The heat flux on the upper surface varies with x coordinates and time.Ŵ inlet , Ŵ external and Ŵ wall stand for the inlet boundary, the outlet boundary and the wall boundary.The 2-D porous plate and the 3-D porous plate have the same dimensions in length and height, while the 3-D model has a finite width of 50 mm compared with the 2-D model.The Start Region and the End Region labeled in Fig. 4 are used for the following result analysis.
Aiming to analyze the transient effect of transpiration cooling under various conditions, including the influence of injection mass flux, external pressure, acceleration and model type, seven numerical cases with different boundary conditions are carried out, as shown in Table 2.For all the cases, the inlet coolant temperature is 300 K, and to represent the flight process during acceleration, the external heat flux distribution decreases with the x coordinate in the form of an inverse proportional function and increases with time.The heat flux reaches a steady value of 1 MW/m 2 after 30 s.In addition, to represent the low-pressure environment of the leading edge of the hypersonic vehicle, the external pressure for Model 1 -Model 3 and Model 5 -Model 7 is selected as 50 kPa, which is largely lower than the standard atmosphere.Three 2D porous plate models (Models 1, 2 and 3) with different mass flux injection are simulated and compared together.Model 4 which is applied with linear decreasing external pressure is discussed with Model 1. Models 1, 5 and 6, which assume that the aircraft has no x-axis direction acceleration, positive x-axis direction acceleration and negative x-axis direction acceleration, respectively, are compared together.The model type of the final model is 3-dimensional, which is simulated to compare with the 2-D plate that has the same boundary conditions.
The physical properties of the inlet water coolant are shown in Table 3.Some critical variables like dynamic viscosity, saturation temperature and latent heat are fitted using proper functions.The selected sintered metal porous medium has the properties shown in Table 4.It has a uniform porosity of 0.3 an average particle diameter of 100 μm.The density of the solid skeleton is 9000 kg/m 3 , the specific heat capacity is 400 J/(kg•K) and the heat transfer coefficient is 30 W/(m•K).

Grid independence test
In order to verify the grid independence of the numerical simulation results, four meshes with the node number 100 × 150, 200 × 250, 400 × 500, and 600 × 750 are generated by the software ICEM CFD.The grids are uniformly refined, with consistent gradients everywhere, gradually becoming denser as the number of nodes along the edges increases.More nodes are arranged along the flow direction of the edges, as shown in Fig. 5.The mass flux of the inlet coolant is 0.5 kg/(m 2 •s) and the coolant temperature is 300 K.The physical properties of the coolant are also provided in Table 3.The external pressure is 1 atm and the applied heat flux is 0.5 MW/m 2 .Their spatial distributions are both uniform.The porous plate with a size of 100 × 20 mm 2 has a heat transfer coefficient of 20 W/(m•K), a porosity of 0.35 and a particle diameter of 500 μm.
The steady state of the four meshes are calculated, and the solid temperature and the liquid saturation at the selected point are plotted.As shown in Fig. 5, the simulation results of the third and fourth meshes are very close.Therefore, the mesh with the node number 400 × 500 will be used for further calculation.For the 3-D physical model, the node number along the flow direction is set as the same with the 2-D model.

The effect of time-dependent and spatially-varied heat flux
Before the result discussion of multiple factors that influence the thermal protection effect of transpiration cooling, the transient effect of Model 1 in Table 2, which is simulated under a time-dependent and spatially-varied heat flux and an injection coolant mass flux of 0.05 kg/(m 2 •s), is analyzed below.As shown in Fig. 6a, the transient change of relative mobility, which represents the liquid saturation, indicates that the phase change takes place first at the start region because of higher heat flux.The liquid at the start region starts to boil at t = 11 s and becomes pure vapor at t = 15 s.At the end region, the liquid saturation begins to increase due to the vapor blockage at the start region that compels more coolant flow to the end region.At t = 30 s, the relative mobility at the end region decreases dramatically because of the step increase of the heat flux.As shown in Fig. 6b, the transient fluid and solid temperature is in consistent with the relative mobility in Fig. 6a.When t = 11 s and t = 17 s, the fluid temperature at the start region and the end region rises to 352 K, which is the saturation temperature of water under a saturation pressure of 50 kPa.The local thermal non-equilibrium phenomenon is obvious between the solid matrix and the fluid inside when the fluid is under a two-phase state.After the fluid becomes pure vapor, the temperature increases rapidly.
To further discuss the distribution of the critical physical properties, the contours of liquid saturation, streamline and fluid temperature are shown in Fig. 7.As shown in Fig. 7a, the start of the phase change will change the flow direction of the coolant, as the kinematic viscosity of water vapor is larger than that of water liquid.Consequently, as shown in Fig. 7c and e, the vapor blockage at the start region will become more severe over time and the fluid temperature will increase as the volume of the liquid phase decreases.In addition, as shown in Fig. 7f, the fluid temperature in the two-phase region is different and it is related to the local saturation pressure.As shown in Fig. 8, with the increase of the mass flux, the solid temperature at the outlet boundary decreases.However, the decline range at different positions is not uniform due to the non-uniform distribution of the liquid coolant.From the results of the velocity distribution, it can be seen that the velocity of water vapor will be much larger than that of water liquid because of much lower density.And the sharp change of velocity can indicate the interface of the liquid phase and the vapor phase.
To analyze the detail of the liquid saturation, flow pattern and solid temperature under different mass fluxes, Fig. 9 shows the contours of them when t = 40 s.As shown in Fig. 9a and b, when the mass flux is 0.55 kg/(m 2 •s), the thermal protection can reach an ideal state, as the coolant maintains the liquid-gas mixture at the outlet without boiling to the pure vapor phase.From Fig. 9c, it can be seen that a second vapor blockage region is formed following the Start Region with the largest heat flux.This phenomenon can be captured under a certain inlet mass flux, when the positions adjoined to the first vapor blockage that has a relatively large heat flux and a relatively large mass flow are more difficult to boil than the rear region that has a relatively low heat flux and a relatively low mass flow.As shown in Fig. 9e and f, when the mass flux is 0.05 kg/(m 2 •s), the thermal protection effect is not desirable.the pressure will decrease the fluid and solid temperature when the fluid is at the twophase region.However, the liquid tends to boil at lower temperature due to the decrease of the overall saturation temperature, which eventually makes the temperature of the solid and the fluid larger after 40 s.
In Fig. 11, the saturation temperature, liquid saturation and streamline contours when t = 35 s under constant and time-dependent external pressure when the injection mass flux is 0.05 kg/(m 2 •s) are shown.As shown in Fig. 11a and c, the saturation temperature is lower when the external pressure is lower, and the distribution is in consistent with the local mixture pressure.As shown in Fig. 11b and d

The acceleration effect
The overload (acceleration) of aircraft can induce an inertial force for the coolant inside the porous media, which will act as a body force to change the velocity pattern.Model 5 and Model 6 in Table 2 have acceleration in opposite directions, and the inertial force induced by them is in the positive and negative direction of the x-axis, respectively.
As shown in Fig. 12a and b, when the inertial force is in the positive direction of the x-axis, the phase change rate of the coolant is increased, while when the inertial force is in the negative direction of the x-axis, the coolant at the start region remains a liquid-gas mixture during the simulation heating process.
In Fig. 13, the contours of pressure, streamline and solid temperature when t = 35 s under different aircraft overload are shown to analyze the influence of inertial forces to the coolant.It can be seen that the change of the flow direction is mainly due to the vapor blockage in Fig. 13a, while the change of the flow direction is mainly due to the inertial force in Fig. 13b and c.The thermal protection effect is therefore dramatically influenced by the motion of aircraft.

The three-dimensional effect on the transient phase-change predictions
In practical engineering applications, the porous plate used is three-dimensional, which means the width of the plate is a finite length.In order to analyze the lateral effect of the transpiration cooling, the following simulation is conducted using a 3-D physical model that has the same length and height with the above 2-D model.The inlet mass flux for both cases is fixed to 0.05 kg/(m 2 •s).
As shown in Fig. 14a, the solid temperature of the 3-D model increases more rapidly during the whole process.The maximum solid temperature is approximately 100 K higher than that of the 2-D simulation results.In addition, the solid temperature of From Fig. 15, it can be seen that the vapor blockage phenomenon is more obvious at the start region under the 3-D model, making more coolant flow through the end region.Therefore, compared with the 2-D simulation results, not all of the coolant boil to pure vapor at t = 40 s.From the above comparison, it can be concluded that the model dimensions have a strong influence on the flow development, especially for the system with finite depth, and therefore change the thermal protection effect.

Conclusions
In this study, a modified model is developed by incorporating the local thermal nonequilibrium effect for transpiration cooling: the governing equations of fluid energy of TPMM are established in a full LTNE form, and the low-pressure and high overload working conditions of hypersonic vehicles are considered by establishing momentum and energy equations under varied saturation pressure, saturation temperature and acceleration conditions.Transient simulations for the effect of external heat flux, injection mass flux, external pressure, aircraft overload and model dimensions are conducted, and the main conclusions are as follows.
(1) The transient simulations of water coolant under time-dependent and spatially-varied heat fluxes suggest that the start of phase change can change the flow distribution of the coolant.Therefore, the vapor blockage at the start region becomes more severe over time, leading to a subsequent rapid increase of the solid temperature.Distributed mass injection should therefore be considered.(2) Increasing the mass flux can obviously reduce the solid temperature, while the decline of temperature at different positions may not be uniform due to the vapor blockage.Moreover, the temperature changing results indicate that when increasing the mass flux to 0.50 kg/(m 2 •s) for the performed case study, the effect of thermal protection is relatively ideal, as the latent heat can be fully utilized.With the reduction of saturation temperature, the vapor blockage becomes more severe at low pressure conditions, leading to further difficulties in thermal protection.(3) The motion of the hypersonic vehicle has a large effect on the cooling process, as the inertial force can significantly change the flow pattern of the fluid inside.A horizontal acceleration e.g., 10 g, can totally change the flow phase distribution of the coolant.Moreover, a 3-D physical model for realistic applications of transpiration cooling is suggested, otherwise the thermal protection efficiency may be overestimated by a 2-D simplification.

Fig. 1
Fig. 1 Flowchart of the numerical solution procedure for the modified LTNE-TPMM of transient phase-change transpiration cooling for active thermal protection

Fig. 2
Fig.2Temperature variations of the top wall using the established model compared with the results of Wu et al.[42]

Fig. 3 Fig. 4
Fig.3The variation of liquid saturation at the outlet using the established model compared with the results of Chen et al.[29]

4. 2
The transient effect of different coolant mass fluxes After discussing the overall transient behavior of the transpiration cooling with phase change, the thermal protection effect of transpiration cooling with different mass fluxes, including 0.05 kg/(m 2 •s), 0.30 kg/(m 2 •s) and 0.50 kg/(m 2 •s) (Models 1, 2 and 3), respectively, are analyzed by comparing the solid temperature distribution at the outlet when t = 40 s.

Fig. 5 Fig. 6
Fig. 5 Comparison of solid temperature and liquid saturation at the selected point obtained by four different meshes

Fig. 7 Fig. 8 Fig. 9 Fig. 10
Fig. 7 Contours of liquid saturation, streamline and fluid temperature under a time-dependent and spatially-varied heat flux when the injection mass flux is 0.05 kg/(m 2 •s) , the case under linear changing external pressure at 35 s has a region at the outlet boundary where the coolant does not change into pure vapor.Model 4 that has the same heat flux distribution with Model 1 results in different final liquid saturation contours because of the time-dependent lowpressure effect.

Fig. 11
Fig. 11 Contours of saturation temperature, liquid saturation and streamline at t = 35 s under constant and time-dependent external pressure when the injection mass flux is 0.05 kg/(m 2 •s)

Fig. 12 Fig. 13
Fig. 12 Transient variations of temperature under different aircraft accelerations when the injection mass flux is 0.05 kg/(m 2 •s)

Fig. 14 Fig. 15
Fig. 14 Comparison of solid temperature of the 2-D and 3-D models under an injection mass flux of 0.05 kg/ (m 2 •s)

Table 3
Physical properties of the coolant used in the numerical simulations

Table 2
Model details for numerical simulation setup in this work

Table 4
Physical properties of the porous plate used in the numerical simulations