OPTIMIZATION STUDIES OF TRANSPIRATION COOLING USING POROUS MEDIUM WITH GRADUALLY-CHANGED STRUCTURE

Non-uniform heating and vapor blockage deteriorate the effectiveness of transpiration cooling, and an optimization method by using porous media with a gradually-changed structure is proposed. The numerical tool based on a two-phase mixture model and local thermal non-equilibrium assumption considering variable properties of coolant is applied. Porous media with linearly-changed porosity or particle diameter is analyzed. The transient results presented that the structure of gradually-changed porosity (or particle diameter) with appropriate parameters can delay the heat transfer deterioration. And it is confirmed that the present theoretical model is an effective tool for optimization design of transpiration cooling system.


INTRODUCTION
Transpiration cooling is considered as an effective active cooling technique to protect the material from extremely high heat flux in different systems (Song et al., 2006;Boehrk, 2014). The excellent cooling performance comes from multiple heat transfer processes. The coolant flows through the tortuous channels in the porous media and convectively exchanges heat with a porous skeleton. Meanwhile, the coolant film is generated which serves as a protective layer at the hot surface and its downstream. When a liquid coolant is applied, the cooling efficiency could be further improved due to the large latent heat during phase change.
Experimentally, the effectiveness of transpiration cooling is validated based on different test facilities and different porous media. Van Foreest et al. (2009) tested the transpiration cooling performance using liquid water based on the wind-tunnel facility. The nose cone porous models made of Procelit 170 were sharped into different radii. The extremely high temperature at stagnation point could be reduced by over 1500 K due to the large cooling capacity of water, which demonstrated the signification of phase change in the cooling process. Huang et al. (2017) experimentally investigated a water-cooled sintered bronze porous plate heated by accelerated hot air. The influence of coolant injection ratio and particle diameter of porous media was studied. Moreover, the vapor-blockage effect is characterized in the hightemperature region. Wu and Zhu (2017) applied a quartz lamp and oxyacetylene flame to generate medium and higher heat flux respectively. The transpiration cooling using water as coolant could protect the porous media below its allowable temperature in both two cases. Meanwhile, the numerical tools are also developed to reveal the detailed heat transfer phenomena in transpiration cooling. Both the enthalpy-based two-phase mixture model (TPMM) (Ramesh and Torrance, 1993) and the separated phase model (SPM) (Wei et al., 2012) could be used. Shi and Wang (2011) developed a coupling method based on the two-phase mixture * Corresponding author. Email: yalinghe@mail.xjtu.edu.cn model and local thermal non-equilibrium model. All the properties of water coolant are assumed as constant. The influence of porosity, thermal conductivity of porous media and coolant injection rate on temperature distribution and area of the two-phase region were investigated. He and Wang (2014) introduced a source term related to phase change in the momentum equation in a separate phase model, and consider the compressibility of vapor by considering it as an ideal gas. The sensitive analysis was conducted to investigate the influence of coolant mass flow rate and heat flux. Chen et al. (2020) investigated the influence of the variable properties of water in transpiration cooling. And it was found that some new phenomena were observed when thermal properties varied with temperature and pressure, which indicated the necessary of this consideration.
In above literatures, only the uniform porous structures were applied. Li and Peterson (2010) studied the different performance of uniform and modulated porous structures. Therefore, it is inspired that one efficient way to improve the transpiration cooling is to modify the local properties of the porous material. On the one hand, researches consider changing the properties along the coolant main flow direction. Shi and Wang (2008) analyzed a two-layered porous medium by local thermal non-equilibrium (LTNE) model. The genetic algorithm was coupled to find out the lowest temperature of the hot surface by a combination of composition, porosity, and thickness of the two layers. Under both constant pressure and mass flow rate conditions, the properties of the porous layer directly heated by the heat flux have more significant effects on the hot surface temperature. The material with high thermal conductivity was recommended to use if the synthetic cost was not considered. Similarly, the two-layer design was accepted by Liu et al. (2013). The layer thickness was fixed in the analysis but the thermal conductivity and porosity in each layer were studied. It was found that higher thermal conductivity and larger porosity near the heating surface are beneficial for reducing the temperature increase, which leads to lower thermal stress.
On the other hand, it is noticed that the heat flux in the porous material surface is non-uniform due to aerodynamic heating. A locally

Frontiers in Heat and Mass Transfer
Available at www.ThermalFluidsCentral.org heated spot may lead to the local generation of water vapor which causes high pressure and flow resistance. Therefore, it is necessary to alleviate the cooling deterioration relating to this vapor block phenomenon by applying a new design of the transpiration cooling process. Zhao et al. (2014) experimentally studied the transpiration cooling of a nose cone with unequal wall thickness. The experiment was conducted based on the arc heated supersonic free jet facility and liquid water was used as a coolant. More coolant could be directed to the stagnant zone by decreasing wall thickness in this region, which resulted in high cooling effectiveness. The gaseous nitrogen was selected as a coolant and local thermal equilibrium model was used in the numerical study of Shen and Wang (2017). A step coolant allocation method was proposed to improve the cooling performance in the leading edge of the nose cone structure. A porous material with non-uniform porosity was further used to distribute the coolant mass flow ratio in different regions. Large porosity in the stagnation zone showed small flow resistance that allowed more coolant to flow through. Wu et al. (2018) applied the pressing and sintering methods coupled with different particle size and pressing pressures to fabricate the non-uniform nose cone with varying permeability. The air was used as coolant to decrease the surface temperature of porous matrix heated in the wind tunnel. It was found that the cooling efficiency could be improved in both the leading-edge region and the entire nose cone under different coolant injection ratios.
From the above literature review, it is noticed that the numerical model for the transpiration cooling with non-uniform porous structure and phase-change process is not complete. Since the non-uniform aerodynamic heat flux at the surface of aircraft will always change during flight, it is meaningful to further study the dynamic response characteristics of the transpiration cooling system. The multi-layer porous model in which abrupt change between layers is used which may leads to the large thermal stress in the interfaces. Moreover, the phasechange of liquid coolant and its real properties are not totally considered in the previous researches. Therefore, a more comprehensive transient numerical model using real coolant thermophysical properties is therefore established. With the application of 3D printing technology, the design and manufacture of gradient-changed porous structure can be realized. And the transpiration cooling performance in a porous media with gradually-changed porous structure is theoretically investigated in this paper.

NUMERICAL MODEL
The TPMM and LTNE assumption are applied to solve the fluid flow and heat transfer in the porous media with gradually-changed porous structure. The TPMM is featured because it avoids the problem of interface tracking. And LTNE is used to capture the temperature of both coolant and porous skeleton.

Governing Equations and Empirical Parameters
The definition of mixture density, velocity and enthalpy are: where ρl and ρv represent the density of liquid water and water vapor, h is specific enthalpy, and hv,sat is the specific enthalpy of water vapor at saturation pressure and temperature. The continuity, momentum, and energy governing equations are: where ρs, cp,s, Ts, and ks,eff represent density, specific heat, solid temperature, and effective thermal conductivity of porous media. ε and K are the porosity and permeability in the porous media. K is evaluated based on ε and dp (particle diameter, characteristic size of porous matrix) (Shi and Wang, 2011): The mixture viscosity is: where krl=s 3 and krv=(1-s) 3 are the relative permeabilities.
In the coolant energy conservation equation, γ and Γh are the advection correction coefficient and diffusion coefficient, which are defined by: where D is the diffusion coefficient, hfg is the latent heat, and kf,eff = εkf is the effective thermal conductivity in the porous matrix. The diffusion coefficient is evaluated based on the capillary pressure J as presented by Eqs (12) and (13).
The energy exchange between the coolant and porous skeleton composes convective heat transfer and boiling. According to the pore size and the Reynolds number, the convective heat transfer coefficient is obtained based on Xu and Jiang (2008), shown in Eq. (14). The energy source due to phase change is presented by Eq. (15).
Moreover, the relationships of coolant temperature and saturation with mixture enthalpy are shown by Eqs. (16) and (17).

Computational Domain and Boundary Conditions
The two-dimensional computational domain of transpiration cooling in porous media with a gradually-changed porous structure is presented in Fig. 1. The coolant goes from the bottom of the water tank and flow through the porous media and then exits from the heating surface.
The mass flow rate at the inlet is constant and the initial temperature is 300 K. In the heating surface, the non-uniform heat flux is applied which increases linearly from 450 to 550 kW/m 2 in the +y direction. The pressure at the outlet is constant which equals 10 kPa. The two-side wall of the porous media and water tank is assumed to be non-slip and adiabatic. The structural parameters of porous media can be changed in each simulation case.

Numerical Treatments
The thermophysical properties of liquid water and water vapor are dependent on both temperature and pressure. The temperature and pressure ranges are 283.15-1003.15 K and 0.01-0.2 MPa. The properties such as density, specific heat, thermal conductivity, and viscosity are obtained by NIST REFPROP (Lemmon et al., 2010) and summarized as correlations as functions of temperature and pressure. The specific enthalpy of saturated water or vapor as well as the latent heat are referred to Glass et al. (2001). The governing equations are solved based on the secondary development platform of the commercial computational fluid dynamic software, ANSYS Fluent. The User-Defined Functions (UDF) are mainly applied to conduct the implementation. The flow model of porous matrix provided by the software is adopted to solve the mass conservation equation and the momentum equation, and User-Defined Scalar (UDS) transport equations are applied to solve the fluid and solid energy equations. More detail ed information is presented in our previous work (Chen et al., 2020).

Model Validation
To verify the reliability of the present model, the calculation results were compared with the simulation results in Shi and Wang (2011) and the experimental results in Dong and Wang (2018).
In Shi and Wang (2011), the thickness of the porous plate is 0.1 m; the porous porosity is 0.35; the average diameter of porous particles is 500 μm; The hot-end heat flux is 10 6 J/(m 2 •s); the inlet mass flow is 0.5 kg/ (m 2 •s). Since the particle diameter of the numerical model studied in Shi and Wang's (2011) paper is on the order of several hundred microns, the empirical correlation of convection heat transfer applied in the model in this section is consistent with that in the literature, which is shown in Eq. (18). Fig. 2 (a) shows the distribution of fluid temperature along the flow direction. It can be seen that the simulation results of the present model are in good agreement with that in Shi and Wang (2011). The main reason for the deviation is that the influence of gravity is considered in the literature while the gravity is ignored in the governing equation of fluid energy in this study. Meanwhile, the authors did not give a specific value of solid thermal conductivity in their case, which will also cause deviations in the results.  (2011) and (b) comparison with experiment data in Dong and Wang (2018).
In Dong and Wang's (2018) experiment study using a hot wind tunnel, the porous plate has a geometric size of 56×32×8 mm, a porosity of 0.347, a thermal conductivity of 18.3 W/(mꞏK) and a permeability of 1.31×10 -13 m 2 . The mass flowrate at the cold boundary is 0.047 kg/(m 2 ꞏs), and the temperature of inlet coolant is 300 K. The hot boundary condition is evaluated using following equations (Dong and Wang, 2018): where T∞ =1073.15 K is the temperature of hot air, h∞(x) is the heat transfer coefficient calculated by the empirical expression suggested by Incropera and DeWitt (1996): where k∞ is the thermal conductivity of hot air, L=0.056 m is the length of the porous plate in x direction, Pr∞=0.695, Re∞=8608, x0 is the length of impermeable part in experiment given in Dong and Wang (2018).
The comparison results of the outlet temperature are shown in Fig.  2 (b). The predicted temperature distribution is quite similar with that in the experiment, with a maximal deviation less than 5.5%. This is acceptable and the reason of the error might be that an empirical formula was adopted to evaluate the outlet heat flux, which was different from the experimental condition.
Therefore, it can be deemed that the present model and corresponding numerical methods are valid in a certain extent.

Optimization Study
Previous studies show that the non-uniform heat flux will lead to inconsistent gasification rates throughout the flow field in a transient transpiration cooling process. Flow resistance is higher in areas with a higher degree of gasification, resulting in a decrease of coolant flow. Then the temperature will rise more rapidly and more vapor will be generated in this area. Thus, a deterioration of heat transfer will finally occur in this area.
Based on the above analysis, the key to solve the problem of deterioration of heat transfer is to redistribute the coolant flow in the porous matrix by taking certain measures, including adjusting porosity distribution and modifying particle diameter distribution. The effect of porosity and particle diameter is mainly reflected in two aspects: permeability (Eq. (8)) and heat convection between the fluid phase and the solid phase (Eq. (14)). The influence of porosity and particle diameter on hsi and K is depicted in Fig. 3. On one hand, the permeability grows with the increment of porosity (or particle diameter), which means the flow resistance is smaller in the higher-porosity (or larger-diameter) area. To optimize the heat transfer of transpiration cooling, higherporosity (or larger-diameter) should be set in the area with higher heat flux. On the other hand, the volumetric heat transfer coefficient decreases with the increment of porosity (or particle diameter), indicating that the energy source term between solid and fluid is smaller and the rise of fluid temperature is delayed. Taken these into consideration, higher-porosity (or larger-diameter) should also be set in the zone with relatively higher heat flux.
Thus, optimization studies are conducted to numerically analyze the effect of the mentioned measures on the deterioration of heat transfer. The parameters of structure and condition in each case are shown in Table 1. Case S (S: standard) is an unoptimizable case under non-uniform heat flux with a value linearly changing from 450 kW/m 2 to 550 kW/m 2 . The non-uniform heat flux of other cases is all the same as that in Case S. Porosity distribution in Case P1 -P3 (P: porosity) is linearly adjusted, while particle diameter distribution is modified in D1 -D3 (D: diameter).
More than 10 cases with different distribution of porosity or particle diameter were simulated but only some typical results are shown in this paper. All the optimization cases are transient. At the time of 0 s, the initial temperature is uniform with a value of 300 K and the field of pressure and velocity is calculated without solving energy governing equations in a steady state in advance.

Influence of porosity
As described in Table 1, the distribution of porosity in Case P1 linearly changes from 0.29 to 0.31 along the y-direction, while the one in Case P2 changes from 0.28 to 0.32 and that in Case P3 changes from 0.27 to 0.33. The average porosity in these three cases maintains at the value of 0.30, which is consistent with that in Case S. The distribution of fluid temperature at the outlet in cases with different porosities at different times is shown in Fig. 4. At the time of 2 s, there is no phase change in Case S, where the fluid temperature is linearly distributed at the outlet, which is consistent with the non-uniform heat flux loaded at the boundary. And it can be observed that the fluid temperature on the left (y < 0) is higher than that on the right (y > 0) in Case P1, P2, and P3 in Fig. 4 (a). Moreover, the phase change happens on the left of the outlet in Case P2 and P3. This phenomenon can be explained by the initial non-uniform flow distribution caused by the gradually-changed porosity. Then at t = 5 s in Fig. 4 (b), the vapor is generated in all the cases. The maximum temperature in Case P3 reaches about 485 K at the left end of the outlet, which is much higher than those in other cases. Later at the t =10 s, the temperature of 4 cases further rises. And the temperature distribution of "higher on the left" changes to "higher on the right" in Case P1. Finally, at t = 15 s, the maximum fluid temperature is located at the right end of the outlet in Case S and P1, while the one in Case P2 and P3 is at the left end of outlet. The maximum temperature in all three optimized cases with changing porosity is lower than that in Case S, as depicted in Fig. 4 (d).
To better display the flow field and temperature field in the dynamic process of transpiration cooling, the contours of fluid temperature, solid temperature, saturation and pressure distribution with added stream lines at t = 9 s are depicted in Fig. 5, Fig. 6, Fig. 7, and Fig. 8, respectively.
It is indicated that the maximum fluid temperature in the computational domain in some cases does not appear at the outlet from Fig. 5. This phenomenon is caused by the variable thermal properties. The saturation temperature is higher in the region with higher pressure, resulting in a reversed temperature distribution along the flow direction. The maximum fluid temperature appears at the right half of domain (y > 0) in Case S and P1, while that in Case P2 and P3 is found in the left half (y < 0). Combining with the saturation results in Fig. 7, it is clearly found that whether the vapor distributions or fluid temperature distribution in Case P1, P2 and P3 is more uniform than that in Case S.
The stream lines showing in Fig. 8 help us better understand the flow distribution. The streamline distribution is similar in Fig. 8 (a) and Fig. 8 (b), where massive coolant floods into the upper left corner. This phenomenon with uneven coolant distribution is not conducive to cooling, resulting in a deterioration of heat transfer happens earlier. In contrast, the flow distributions in Case P2 and P3 are better than those in Case S and P1. More coolant is redistributed into the region with high heat flux.
To quantitatively describe the degree of non-uniformity of the fluid temperature at the outlet, the temperature standard deviation σ is introduced: where N is the total number of equally spaced samples at the outlet, Tave is the average temperature of all samples. A higher σ represents a greater unevenness of the temperature distribution.
To illustrate the dynamic process in these four cases, Fig. 9 shows the transient results at the outlet in cases of porosity optimization, including average fluid and maximum fluid temperatures, average saturation, and temperature standard deviation. It is found that both the average temperature and maximum temperature in Case P2 and P3 are higher than that in Case S between t =5 s and t = 9 s. It is due to the overregulation of coolant flow caused by the large-range porosity distribution in Case P2 and P3. As shown in the Fig. 5 (c) and (d), the maximum fluid temperature appears in the left area (y < 0) in Case P2 and P3, where the heat flux is relatively lower but the coolant flowrate is even much smaller during this period. After t = 10 s, the temperature in all the three modified cases become lower than that in Case S. As given in Fig. 9 (c), the standard deviation of outlet temperature in each case shows a trend of rising fluctuation. The decrease at a certain section of the curve (e.g. 5 s -7 s of Case S) is related to the dynamic process of vapor generation. Firstly, the vapor is formed where the coolant temperature rises to the saturation temperature. Then the vapor region and the two-phase region expand. The distribution of liquid and vapor at the outlet is non-uniform and a gap of the liquid outflow with great mass flowrate (e.g. the upper right corner in Fig. 7 (c)) is possibly formed. As the interface of the two-phase region reaches the wall, the heat and flow will be redistributed and the temperature gradient will therefore change, causing a phenomenon that the previous gap of liquid outflow disappears and then a new gap at some other location will be formed later. The deviation of temperature at the outlet will reduce with the disappearing of the liquid outflow gap.
The curves of average saturation in different cases in Fig. 9 (d) show a similar pattern. At first, the saturation is kept at 1.0 before the fluid temperature reaches the saturation temperature. Then the vapor is produced and the saturation decreases rapidly. After reaching the minimum value, the saturation might rise slightly, indicating that the gap of liquid outflow is forming during this period.
Comparing these porosity-modified cases to Case S, the temperature is lower with better uniformity in the former. These results suggest that adjusting the porosity distribution may help improve the heat transfer deterioration of non-uniform heat flux to some extent. However, it should be noted that the matching process of coolant redistribution and the non-uniform heat flux is dynamic. And the unsuitable distribution of porosity might even worsen the heat transfer of cooling.

Influence of particle diameter
As described in Table 1, the distribution of particle diameter in Case D1 linearly changes from 19 μm to 21 μm along the y-direction, while the one in Case D2 changes from 18 μm to 22 μm and that in Case D3 changes from 16 μm to 24 μm. The average particle diameter in these three cases maintains at the value of 20 μm, which is consistent with that in Case S. Fig. 10 shows the distribution of fluid temperature at the outlet in cases with different particle diameter at different times. At the time of 2 s, there is no phase change in all 4 cases, and the fluid temperature is linearly distributed at the outlet. It can be observed that the fluid temperature on the left (y < 0) is higher than that on the right (y > 0) in Case D2 and D3. Then at t = 5 s, the vapor is generated in all the cases. Later at t =10 s, the fluid temperature remains near the saturation point in Case D2 while the maximum fluid temperature in the other 3 cases exceeds 420 K. Finally, at t = 15 s, the temperature further rises. The maximum fluid temperature is located at the left end of outlet in Case D3, while the one in the other 3 cases is at the right end of the outlet. The maximum temperature in all the optimized cases with changing distribution of particle diameter is lower than that in the Case S.
The contours of fluid temperature, solid temperature, saturation and pressure distribution with stream lines at t = 5 s in cases with different particle diameter are depicted in Fig. 11, Fig. 12, Fig. 13, and Fig. 14, respectively. The maximum fluid temperature appears at the right half of domain (y > 0) in Case S, D1, and D2, while that in Case D3 is found in the left half (y < 0). Combined with the saturation distribution in Fig. 13 and stream lines in Fig. 14, it is found that the distribution of fluid temperature and coolant flow in Case D2 is most reasonable in these 4 cases at the current time.     15 shows the transient results at the outlet for diameter optimization, including average fluid temperature, maximum fluid temperature, average saturation, and temperature standard deviation. As seen from Fig. 15 (a), the average fluid temperature starts rising rapidly at t = 7 s after maintaining a stable value near the saturation point in Case S, while the beginning time in Case D1, D2 and D3 is at t = 9 s, t = 10 s and t = 8 s, respectively. Either the curves of average temperature or the curves of maximum temperature indicate that temperature increase in Case D1 -D3 is slower than that in Case S. Meanwhile, the results of standard deviation in Fig. 15 (c) suggest that the temperature distribution in optimized cases is also more uniform than the one in Case S.
Simulation results suggest that the heat transfer in a non-uniform heat flux condition can be improved by adjusting the porous structure. However, the state of the most uniform distribution of heat or flow in transient transpiration cooling is unstable. In other words, this phasechange transpiration cooling is a positive feedback regulation system: once the equilibrium state of matched coolant flow and heat is broken, the temperature in the area where vapor is firstly generated will keep rising and the coolant flow here will keep decreasing, causing a heat transfer deterioration inevitably. In all optimization cases with graduallychanged porous structures, the deterioration of heat transfer occurs in the final stages of simulations. But a proper structure which is optimized for specific working conditions do will delay the deterioration. These findings are helpful for guiding the later design work of the experimental study and will be beneficial to the practical application of the active thermal protection system.

CONCLUSION
In this paper, optimization studies of transpiration cooling using porous media with a gradually-changed porous structure are conducted by applying the numerical phase-change model combining TPMM and LTNE assumption with variable coolant properties. It is confirmed that this theoretical model is an effective tool for the optimization design of the porous structure in transpiration cooling with phase change. The main conclusions are as follows: (i) The principle of improving the heat transfer deterioration caused by non-uniform heat flux is that changing the pore structure, including distribution of porosity and particle diameter in porous media, can influence the permeability and heat transfer convection coefficient of fluid and solid, resulting the redistribution of coolant flow to match the non-uniform heat flux.
(ii) Comparing with the uniform structure, the structure of gradually-changed porosity (or particle diameter) with appropriate parameters can delay the heat transfer deterioration but cannot prevent this phenomenon because of the unstable characteristics of the matched distribution of coolant flow and heat in transient transpiration cooling.
(iii) An over-adjusted porous structure might lead to heat transfer deterioration occurs in the zone with relatively lower heat flux, indicating that the design of the porous structure should match the real flight condition.