Numerical Study on Regenerative Cooling Characteristics of Kerosene Scramjets

The use of kerosene-based regenerative cooling for scramjet has been found widespread attention due to its inherent nature of high energy utilization efficiency and good thermal protection performance. In order to provide a reference for the later design and experiments, three-dimensional turbulence simulations and sensitivity analysis were performed to determine the effects of three operating mode parameters, heat flux, mass flow rate, and outlet pressure, on the regenerative cooling characteristics of kerosene scramjets. A single rectangular-shaped channel for regenerative cooling was assumed. The RNG k-ε turbulence model and kerosene cracking mechanism with single-step global reaction were applied for the supercritical-pressure heat transfer of kerosene flows in the channel. Conclusions can be drawn that as the kerosene temperature rises along the channel, the decrease of fluid density and viscosity contributes to increasing the fluid velocity and heat transfer. When the kerosene temperature is close to the pseudocritical temperature, the pyrolysis reaction results into the rapid increase of fluid velocity. However, the heat transfer deterioration occurs as the specific heat and thermal conductivity experience their turning points. The higher heat flux leads to lower heat transfer coefficient, and the latter stops rising when the wall temperature reaches the pseudocritical temperature. The same rising trend of the heat transfer coefficient is observed under different outlet pressures, but the heat transfer deterioration occurs earlier at smaller outlet pressure for the reason that the corresponding pseudocritical temperature decreases. The heat transfer coefficient increases significantly along with the rise of the mass flow rate, which is mainly attributable to the increase of Reynolds number. Quantitative results indicate that as the main influence factors, the heat flux and mass flow rate are respectively negatively and positively relative to the intensification of heat transfer, but outlet pressure always has little effects on cooling performance.


Introduction
As the vehicle speed moves into supersonic and hypersonic regimes, the supersonic combustion ramjet (scramjet) system must provide higher combustion pressure and larger efficient thrust [1][2][3][4][5]. Sorts of thermal protection techniques [6], such as the air-breathing cooling and passive thermal protection techniques, are unable to realize effective thermal protection for high-temperature propulsion components. Confronted with this challenge, the kerosene-based regenerative cooling technology has been proposed to sustain more stable operation for the scramjet system [7][8][9]. In order to achieve effective thermal protection, kerosene (coolant) absorbs heat flux from the combustor while flowing through cooling channels at supercritical pressures and undergoes a state change from liquid to supercritical state [10]. Meanwhile, the heat-absorbing capacity of kerosene significantly increases by means of pyrolysis [11], which is a nonlinear transformation process with strong coupling effect between heat and mass transfer.
Many efforts, both experimentally and numerically, have been directed to obtain the fundamental understanding of supercritical-pressure heat transfer characteristics during the above active cooling process. Ulas et al. [12] performed some numerical analyses of regenerative cooling in liquid rocket engines and determined the effect of different aspect ratios and number of cooling channels on the coolant temperature and pressure drop in cooling channels. Pu et al. [13] developed and validated the CFD code for conjugate heat transfer simulation with OpenFOAM software and conducted the coupled thermal analysis of the regenerative cooling structure. Based on the experimental and numerical study, Li et al. [14] and Hu et al. [15] presented new correlations of heat transfer and compared with conventional formulas, such as Dittus-Boelter and Stiegemeier. Zhang et al. [16] experimentally investigated the turbulent flow and convective heat transfer of China RP-3 kerosene in a horizontal straight tube and examined the validity of Chilton-Colburn analogy. The turbulence simulations of Choubey et al. [17] show that compared with other multistrut and wall injection schemes, the collaboration of multistrut combined with 2 wall injection can achieve higher mixing and combustion efficiency for scramjet. Furthermore, they conducted the sensitivity analysis for the effects of various structural parameters on the mixing process [18,19] and combustion enhancement [20,21] of the scramjet. As the actual kerosene pyrolysis process is too complex for numerical simulation, Ward et al. [22,23] and Zhu et al. [24] proposed the global or lumped chemical kinetics to predict kerosene pyrolysis based on the experimental data. Jiang et al. [25] modified the Kumar-Kunzru model to predict pyrolysis behaviors (including secondary reactions and aromatic formation) of China RP-3 aviation kerosene at 5 MPa. The newly developed model contains 18 species and 24 elementary reactions and is applicable at high conversion up to approximately 90%. Ruan et al. [26] proposed an approach for further simplification of the proportional product distribution (PPD) model [22]. The simplified pyrolysis mechanism, with great maneuverability, was validated and applied for numerical studies of supercritical-pressure flows and heat transfer of n-decane in a microcooling tube.
In this paper, three-dimensional turbulence simulations were performed to investigate the supercritical-pressure heat transfer characteristics of kerosene flows in a single cooling channel of kerosene scramjets. Then, the sensitivity analysis was conducted to reveal the effects of heat flux, mass flow rate, and outlet pressure on the heat transfer coefficient and cooling performance.

Computational Domain and Boundary Conditions.
Geometry of the rectangular cooling channel with the coordinate system is shown in Figure 1. Kerosene flows along the positive direction of z at supercritical pressures, and buoyancy is neglected for eliminating the interruption factors.
The entire channel was divided into three sections. The middle section with L = 2 m in length is heated with a constant heat flux imposed on the bottom surface. The unheated sections before and after the heated section are treated as adiabatic just for ensuring a fully developed channel flow and avoiding the outlet effect, respectively. The temperature and mass flow rate of kerosene are specified at the inlet, while the back pressure at the outlet. Other variables at boundary regions are extrapolated from the interior. Steel was adopted as the channel material with constant thermal conductivity of 16.27 W/(m·K).

Kerosene Pyrolysis Model.
In order to reduce the difficulty of kerosene pyrolysis modelling, n-decane was selected as the surrogate of kerosene as they have similar thermophysical properties and cracked product distribution [27,28]. Then, the single-step global reaction for mild cracking of ndecane (3.45 MPa~11.38 MPa) was employed [22,23]:

International Journal of Aerospace Engineering
A key issue in solving the supercritical-pressure heat transfer of kerosene flows is to accurately figure out the kerosene thermophysical properties. The general framework for the same purpose has been established and validated in the previous researches [26,29]. As the reference for thermophysical property evaluation, the density of propane (C 3 H 8 ) was acquired through the NIST program and the modified Benedict-Webb-Rubin equation of state. Then, the density, viscosity, coefficient, and thermal conductivity of kerosene cracked products were obtained by means of the modified correspondingstate method, while specific heat through solving the Soave-Redlich-Kwong state equation and basic thermodynamic equations. As seen in Figure 2, the verification result of thermophysical property evaluation against the published experimental data [30][31][32][33] at 5 MPa shows good agreement except for the specific heat. The deviation of the specific heat between the calculated data and experimental one exceeds 10% as the temperature approaches the pseudocritical temperature (708 K) at 5 MPa. This makes the heat transfer deterioration more serious to a certain extent, but the impacts on the overall performance of kerosene-based regenerative cooling are limited.

Numerical Method.
Three-dimensional steady-state numerical simulations in this paper concern the supercriticalpressure heat transfer of kerosene flows in the channel and the thermal conduction inside the solid channel region. A secondorder, double-precision solver was adopted for spatial discretization with the finite volume method in the commercial code ANSYS, and then variables were calculated through solving the discretized governing equations with the SIMPLEC algorithm. The user-defined function (UDF) of the aforementioned kerosene pyrolysis model was coded and implemented into ANSYS.
The following conservation equations of mass, momentum, energy, and species are numerically solved in the fluid phase:

International Journal of Aerospace Engineering
The thermal conduction equation for the solid channel region is determined in the following form: As for the turbulent flow of kerosene, the RNG k-ε twoequation model (Eqs. (4) and (5)) was adopted in the fully developed turbulent region (Re > 200), while the enhanced wall treatment (Wolfstein one-equation model) in the nearwall region (Re ≤ 200) [34].
The numerical method and its applicability for supercritical-pressure heat transfer of kerosene flows have been demonstrated by other researchers [14,35].

Grid Convergence and Method
Validation. The Solidworks software was employed to establish the threedimensional solution domain. Structured meshes with different grid distributions were, respectively, applied to the fluid and solid domains; then, the two parts of structured meshes were jointed with the coupled wall as mesh interfaces. The concern for the grid convergence analysis is to accurately capture the sharp temperature gradient in the near-wall region. Three sets of grids with different cross-sectional grid number in the fluid domain (361 for Grid 1, 625 for Grid 2, and 1089 for Grid 3) were constructed (presented in Figure 3). Herein, the axial grid number and cross-sectional grid number of the solid domain were, respectively, set to 4800 and 420, which can achieve the required precision by contrast with the prior numerical studies [35,36]. The numerical set-ups consist with the basic case in Subsection 3.1. As shown in Figure 4, the simulation results were compared through the quantitative analysis of the interior-wall temperature at the side surface (T w s ) and y + at the third wall-adjacent node in the fluid domain. Compared with Grid 3, the maximum deviation of T w s for Grid 1 is 5%, while that for Grid 2 is within 1%. Besides, the axial distribution of y + indicates that there are at least three layers of grids in the viscous sublayer region (y + < 5) for Grid 2 and Grid 3, satisfying the meshing requirement of the enhanced wall treatment. Thus, the grid distribution of Grid 2, i.e., 4800 (axial)×(625 +224) (cross-sectional), was used as the grid-independent solution for subsequent simulations in the paper.
Prior to detailed numerical investigations, the above numerical method was validated against published data [36],     Figure 1. The inlet pressure, temperature, velocity, and Reynolds number of kerosene are 4 MPa, 300 K, 4 m/s, and 5520, respectively. The value of heat flux is 3 MW/m 2 , and the length L of the middle section is 300 mm. Figure 5 gives the axial distribution of bulk temperature (T b ) and interior-wall temperature at the bottom surface (T w b ) along the axial direction, obtained from the present computation and published data. The maximum discrepancies of T b and T w b are both within 5%, which is mainly aroused by different evaluation methods of thermophysical properties. The simulation results agree well with the published data [36] in general, proving the accuracy of the present numerical method.

Results and Discussion
3.1. Basic Case. The basic case for illustrating the fundamental characteristics of kerosene-based regenerative cooling was set as follows: the inlet temperature is 300 K; the outlet  International Journal of Aerospace Engineering pressure is 5 MPa; the mass flow rate is 5 g/s; and the heat flux is 1.25 MW/m 2 . Figure 6 illustrates the detailed temperature distribution in both of solid and fluid regions. As presented in Figures 7(c) and 7(d), two turning points, i.e., the peak of specific heat and the bottom of thermal conductivity, are observed in the vicinity of z = 1:5 m, which represents that the kerosene temperature is close to the pseudocritical temperature (708 K) at 5 MPa. Figures 7(a) and 7(b) show that the density and viscosity decrease as the kerosene tempera-ture rises, contributing to the rise of fluid velocity along the channel. Additionally, as clearly seen in Figure 8, the pyrolysis reaction of kerosene results into the rapid increase of fluid velocity when z > 1:5 m.
As the thermal resistance in the solid region is constant, the increase of fluid velocity indicates a more fully developed flow and a thinner viscous sublayer, which does benefit to the intensification of heat transfer at the bottom surface inside the cooling channel. It is observed from Figure 9 that when z < 1:5 m, q w b undergoes great variation along the flow 7 International Journal of Aerospace Engineering direction while q w s only falls slightly. Although the fluid velocity increases continuously when the kerosene temperature at the bottom is close to the pseudocritical temperature (708 K) at 5 MPa, the occurrence of the turning points of specific heat and thermal conductivity increases the thermal resistance in the near-wall region, which holds back the previous trend of q w b : 3.2. Convective Heat Transfer. By comparison with the basic case, the influences of heat flux, outlet pressure, and mass flow rate on the heat transfer coefficient in the rectangular cooling channel were numerically studied. From the principle of Newton's cooling, the heat transfer coefficient (h) is defined as follows: The cross-sectional averaged bulk temperature (T b ) can be calculated by Figure 10 shows the distributions of the heat transfer coefficient and wall temperature with constant heat flux values of 1, 1.25, and 1.5 MW/m 2 . As the bulk temperature on the right hand side of Eq. (6) rises, the change of the heat transfer coefficient depends on the comprehensive effects of wall temperature and heat flux. When the increase of (T w -T b ) exceeds that of the corresponding heat flux, the heat transfer will be weakened. In general, the higher the heat flux is, the lower the heat transfer coefficient is. The heat transfer coefficient stops rising when the wall temperature is close to the pseudocritical temperature (708 K) at 5 MPa, which indicates that the heat transfer deterioration occurs.
The difference of wall temperatures among the three sets of heat fluxes is large at the initial heated section. Then, the difference descends mildly along the flow direction. When the heat flux is greater than 1 MW/m 2 , the wall temperature falls firstly and then increases. Furthermore, the larger the heat flux is, the larger the variation amplitude of the wall temperature is.
For the normal convective heat transfer with fluid pressure and temperature below their critical points, the Reynolds number and heat transfer coefficient almost do not change under the fully developed flow condition. However, the difference is the case at supercritical pressures. Distributions of the heat transfer coefficient and Reynolds number, with different outlet pressures of 4, 5, and 6 MPa, are shown in Figure 11.
The rising trend of the heat transfer coefficient is in accordance with that of Reynolds number before the wall temperature reaches the corresponding pseudocritical temperature at different supercritical pressures. And the maximum relative deviations of the heat transfer coefficient and Reynolds number, among the three sets of outlet pressure, are approximately equivalent (within 6%). The Reynolds number drops slightly as the pressure increases, mainly due to the increase of kerosene density and viscosity. Consequently, the residence time increases with the rising pressure, making heat transfer coefficient decreasing.
When the wall temperature is close to the corresponding pseudocritical temperature at different supercritical pressures, the heat transfer deterioration occurs even if Reynolds number increases sharply. The smaller the outlet pressure is, the earlier the heat transfer deterioration occurs. The main reason for this phenomenon is that the pseudocritical temperature of kerosene drops corresponding to the decrease of supercritical pressure. In general, the outlet pressure produces relatively slight effect on the heat transfer coefficient and Reynolds number. Figure 12 presents the distributions of the heat transfer coefficient with different mass flow rates of 4, 5, and 6 g/s. The difference of the heat transfer coefficient among the three sets of the mass flow rate reaches the maximum when the heat transfer deterioration occurs, which to some extent resembles that among the three sets of heat flux shown in Figure 10(a). The greater the mass flow rate is, the larger the heat transfer coefficient is. This phenomenon is mainly attributable to the increase of Reynolds number, which is consistent with the related findings from above.

Cooling Performance
Estimates. The effects of various operating mode parameters on the cooling performance were  International Journal of Aerospace Engineering studied through sensitivity analysis. Three operation mode parameters, heat flux, outlet pressure, and mass flow rate, were selected as design variables, while the bulk temperature, wall temperature, and mass fraction of unreacted kerosene at outlet as performance indexes. As a statistical method, the Optimal Latin Hypercube method was employed to better explore the design space for design of experiments [37,38] (listed in Table 1); then, 256 sampling points were obtained, which spread evenly and afford to capture higher-order effects.
In order to conduct the sensitivity analysis efficiently, a surrogate model is required to express the relation between performance indexes and design variables. In the present study, the polynomial-based response surface (PRS) model of second-order (shown in Eq. (8)) was chosen to map the given design space and estimates the cooling performance [39].
where y represents different objective performance indexes, x represents the given design variables, and c represents the polynomial coefficients. Through normalization processing and least squares fitting for the polynomial coefficients of Eq. (8), Pareto graphs could fairly display the relative effects of each factor on a performance index using ordered bars. The Pareto graphs in Figure 13 show the contribution of all the factors to three performance indexes.
It can be concluded from Figures 13(a) and 13(b) that main factors have same effects on T w and T b at outlet. H and q m devote about 35% positive and negative contribution, respectively, which means that both of the two performance indexes always increase as H increases or q m decreases. Second orders of H and q m also have notable effects on T w and T b at outlet, but the interaction effect of H and q m on them is adverse (-6.4% for T w and 4.8% for T b , resp.). However, the effects of other factors about P out are nearly negligible, which is consistent with the results mentioned above. As shown in Figure 13(c), the positive effects of q m and its   International Journal of Aerospace Engineering interaction with H on η u at outlet are the most remarkable (27.6% and 18.3%, resp.), while H and the second order of q m have the most negative effect (-23.9% and -19.9%, resp.). The remaining factors contribute less than 5% effect on η u at outlet.

Conclusion
The regenerative cooling characteristics of kerosene flows in a rectangular-shaped channel are investigated numerically in this study. The accuracy of the numerical method was proven by comparison with published data. Based on the analysis of this study, the following conclusions can be made: (1) For the basic case, the kerosene density and viscosity decrease with the increase of the kerosene temperature rises, which leads to the rise of fluid velocity and the intensification of heat transfer. The specific heat and thermal conductivity, respectively, reach the maximum and minimum when the kerosene temperature is close to the pseudocritical temperature (708 K) at 5 MPa, and the heat transfer deterioration occurs although the pyrolysis reaction promotes the rapid increase of fluid velocity. The interior-wall heat flux at the bottom surface undergoes great variation along the channel, while those at other surfaces are basically unchanged (2) The heat transfer coefficient stops rising when the wall temperature reaches the pseudocritical temperature. The heat transfer coefficient decreases as the heat flux increases, and the wall temperature undergoes a greater change. The rising trends of the heat transfer coefficient under different outlet pressures are almost the same before the kerosene is heated to the pseudocritical temperature. However, the heat transfer deterioration occurs earlier at smaller outlet pressure for the reason that the corresponding pseudocritical temperature decreases. The heat transfer coefficient increases significantly along with the rise of the mass flow rate of kerosene, which is attributable mainly to the increase of Reynolds number (3) The quantitative results of sensitivity analysis indicate that the heat flux and mass flow rate have great but adverse influences on the heat transfer coefficient and cooling performance. However, outlet pressure has little effect on regenerative cooling characteristics, and the effect appears to be relatively obvious only when heat transfer deterioration occurs Nomenclature C 1 ,C 2 : Constants in turbulent model c p : Specific heat, kJ/(kg·K) e t : Total internal energy, J G k : Turbulent generation term h: Heat transfer coefficient, kW/(K·m 2 ) k: Turbulent kinetic energy, J/kg L: Length of heated section, m M w : Molecular weight, g/mol p: Pressure, Pa q: Heat flux, MW/m 2 S: Cross-section area, m 2 T: Temperature, K u: Velocity vector, m/s V: Fluid velocity, m/s Y: Species mass fraction y: Distance to the nearest wall, m y + : Dimensionless wall distance, y + = ðy ffiffiffiffiffiffiffiffiffi ffi τ w /ρ p Þ/ν.

b:
Bulk quantity i: Species w: Channel walls w b: Interior-wall at the bottom surface w s: Interior-wall at the side surface.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.