# Numerical Study on Regenerative Cooling Characteristics of Kerosene Scramjets

Authors: Xuan Jin; Chibing Shen; Xianyu Wu
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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–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–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.


## 2. Physical Model and Numerical Methodology

### 2.1. Computational Domain and Boundary Conditions

Geometry of the rectangular cooling channel with the coordinate system is shown in Figure1. Kerosene flows along the positive direction of z at supercritical pressures, and buoyancy is neglected for eliminating the interruption factors.Figure 1
Computational configuration of the rectangular cooling channel.The entire channel was divided into three sections. The middle section withL=2m 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).


### 2.2. 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 n-decane (3.45 MPa~11.38 MPa) was employed [22, 23]:
(1)C10H22→0.151H2+0.143CH4+0.256C2H4+0.126C2H6+0.230C3H6+0.180C3H8+0.196C4H8+0.102C4H10+0.171C5H10+0.124C5H12+0.195C6H12+0.089C6H14+0.169C7H14+0.072C7H16+0.152C8H16+0.012C8H18+0.053C9H18+0.003C9H20A 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 (C3H8) 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 corresponding-state 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–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.Figure 2
Verification of thermophysical property evaluation against experimental data [30–33] at 5 MPa.
(a)
Density and specific heat(b)
Viscosity and thermal conductivity


### 2.3. Numerical Method

Three-dimensional steady-state numerical simulations in this paper concern the supercritical-pressure heat transfer of kerosene flows in the channel and the thermal conduction inside the solid channel region. A second-order, 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:(2)∇⋅ρu=0,∇⋅ρuu=∇⋅τ−∇p,∇⋅ρuet=∇⋅λ∇T−∇⋅pu,∇⋅ρYiu=−∇⋅ρYiud,i+ωiMwi.The thermal conduction equation for the solid channel region is determined in the following form:(3)∇⋅λ∇T=0.As for the turbulent flow of kerosene, the RNG k-ε two-equation 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 near-wall region (Re≤200) [34].
(4)∇⋅ρuk=∇⋅αkμeff∇k+Gk−ρε,(5)∇⋅ρuε=∇⋅αεμeff∇ε+C1εkGk−C2ρε2k.The numerical method and its applicability for supercritical-pressure heat transfer of kerosene flows have been demonstrated by other researchers [14, 35].


### 2.4. Grid Convergence and Method Validation

The Solidworks software was employed to establish the three-dimensional 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 Figure3). 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 (Tw_s) and y+ at the third wall-adjacent node in the fluid domain. Compared with Grid 3, the maximum deviation of Tw_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.Figure 3
Cross-sectional grid distribution of the solution domain.Figure 4
Axial distributions of the interior-wall temperature at the side surface (a) andy+ at the third wall-adjacent node in the fluid domain (b).
(a)(b)Prior to detailed numerical investigations, the above numerical method was validated against published data [36], which has a similar computational configuration with the one shown in 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/m2, and the length L of the middle section is 300 mm.Figure5 gives the axial distribution of bulk temperature (Tb) and interior-wall temperature at the bottom surface (Tw_b) along the axial direction, obtained from the present computation and published data. The maximum discrepancies of Tb and Tw_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.Figure 5
Numerical method validation based on published data [36].


## 2.1. Computational Domain and Boundary Conditions

Geometry of the rectangular cooling channel with the coordinate system is shown in Figure1. Kerosene flows along the positive direction of z at supercritical pressures, and buoyancy is neglected for eliminating the interruption factors.Figure 1
Computational configuration of the rectangular cooling channel.The entire channel was divided into three sections. The middle section withL=2m 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).


## 2.2. 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 n-decane (3.45 MPa~11.38 MPa) was employed [22, 23]:
(1)C10H22→0.151H2+0.143CH4+0.256C2H4+0.126C2H6+0.230C3H6+0.180C3H8+0.196C4H8+0.102C4H10+0.171C5H10+0.124C5H12+0.195C6H12+0.089C6H14+0.169C7H14+0.072C7H16+0.152C8H16+0.012C8H18+0.053C9H18+0.003C9H20A 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 (C3H8) 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 corresponding-state 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–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.Figure 2
Verification of thermophysical property evaluation against experimental data [30–33] at 5 MPa.
(a)
Density and specific heat(b)
Viscosity and thermal conductivity


## 2.3. Numerical Method

Three-dimensional steady-state numerical simulations in this paper concern the supercritical-pressure heat transfer of kerosene flows in the channel and the thermal conduction inside the solid channel region. A second-order, 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:(2)∇⋅ρu=0,∇⋅ρuu=∇⋅τ−∇p,∇⋅ρuet=∇⋅λ∇T−∇⋅pu,∇⋅ρYiu=−∇⋅ρYiud,i+ωiMwi.The thermal conduction equation for the solid channel region is determined in the following form:(3)∇⋅λ∇T=0.As for the turbulent flow of kerosene, the RNG k-ε two-equation 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 near-wall region (Re≤200) [34].
(4)∇⋅ρuk=∇⋅αkμeff∇k+Gk−ρε,(5)∇⋅ρuε=∇⋅αεμeff∇ε+C1εkGk−C2ρε2k.The numerical method and its applicability for supercritical-pressure heat transfer of kerosene flows have been demonstrated by other researchers [14, 35].


## 2.4. Grid Convergence and Method Validation

The Solidworks software was employed to establish the three-dimensional 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 Figure3). 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 (Tw_s) and y+ at the third wall-adjacent node in the fluid domain. Compared with Grid 3, the maximum deviation of Tw_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.Figure 3
Cross-sectional grid distribution of the solution domain.Figure 4
Axial distributions of the interior-wall temperature at the side surface (a) andy+ at the third wall-adjacent node in the fluid domain (b).
(a)(b)Prior to detailed numerical investigations, the above numerical method was validated against published data [36], which has a similar computational configuration with the one shown in 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/m2, and the length L of the middle section is 300 mm.Figure5 gives the axial distribution of bulk temperature (Tb) and interior-wall temperature at the bottom surface (Tw_b) along the axial direction, obtained from the present computation and published data. The maximum discrepancies of Tb and Tw_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.Figure 5
Numerical method validation based on published data [36].


## 3. 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 pressure is 5 MPa; the mass flow rate is 5 g/s; and the heat flux is 1.25 MW/m2.Figure6 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.5m, 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 temperature 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.5m.Figure 6
Temperature distribution in both solid and fluid regions.Figure 7
Thermophysical properties of kerosene at the symmetry plane of the cooling channel.
(a)(b)(c)(d)Figure 8
Distributions of fluid velocity (a) and mass fraction of unreacted kerosene (b).
(a)(b)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 Figure9 that when z<1.5m, qw_b undergoes great variation along the flow direction while qw_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 qw_b.Figure 9
Distributions of the interior wall heat fluxes at the bottom surface (qw_b) and side surface (qw_s).


### 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:
(6)h=qTw−Tb.The cross-sectional averaged bulk temperature (Tb) can be calculated by
(7)Tb=∫STdSS.Figure10 shows the distributions of the heat transfer coefficient and wall temperature with constant heat flux values of 1, 1.25, and 1.5 MW/m2. 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 (Tw-Tb) 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.Figure 10
Distributions of the heat transfer coefficient (a) and wall temperature (b) under different heat fluxes.
(a)(b)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/m2, 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 Figure11.Figure 11
Distributions of the heat transfer coefficient (a) and Reynolds number (b) under different outlet pressures.
(a)(b)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.Figure12 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.Figure 12
Distributions of the heat transfer coefficient under different mass flow rates.


### 3.3. Cooling Performance Estimates

The effects of various operating mode parameters on the cooling performance were 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.Table 1
Design space for design of experiments.
ParametersNomenclatureValue rangeHHeat flux (MW/m2)1~1.5PoutOutlet pressure (MPa)4~6qmMass flow rate (g/s)4~6In 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].
(8)y=c0+∑i=13cixi+∑i=13ciixi2+∑i=12∑j=i+13cijxixj,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.Figure 13
Pareto graphs for the wall temperature (a), bulk temperature (b), and mass fraction (c) of unreacted kerosene at outlet.
(a)(b)(c)It can be concluded from Figures13(a) and 13(b) that main factors have same effects on Tw and Tb at outlet. H and qm devote about 35% positive and negative contribution, respectively, which means that both of the two performance indexes always increase as H increases or qm decreases. Second orders of H and qm also have notable effects on Tw and Tb at outlet, but the interaction effect of H and qm on them is adverse (-6.4% for Tw and 4.8% for Tb, resp.). However, the effects of other factors about Pout are nearly negligible, which is consistent with the results mentioned above. As shown in Figure 13(c), the positive effects of qm and its interaction with H on ηu at outlet are the most remarkable (27.6% and 18.3%, resp.), while H and the second order of qm have the most negative effect (-23.9% and -19.9%, resp.). The remaining factors contribute less than 5% effect on ηu at outlet.


## 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 pressure is 5 MPa; the mass flow rate is 5 g/s; and the heat flux is 1.25 MW/m2.Figure6 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.5m, 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 temperature 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.5m.Figure 6
Temperature distribution in both solid and fluid regions.Figure 7
Thermophysical properties of kerosene at the symmetry plane of the cooling channel.
(a)(b)(c)(d)Figure 8
Distributions of fluid velocity (a) and mass fraction of unreacted kerosene (b).
(a)(b)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 Figure9 that when z<1.5m, qw_b undergoes great variation along the flow direction while qw_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 qw_b.Figure 9
Distributions of the interior wall heat fluxes at the bottom surface (qw_b) and side surface (qw_s).


## 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:
(6)h=qTw−Tb.The cross-sectional averaged bulk temperature (Tb) can be calculated by
(7)Tb=∫STdSS.Figure10 shows the distributions of the heat transfer coefficient and wall temperature with constant heat flux values of 1, 1.25, and 1.5 MW/m2. 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 (Tw-Tb) 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.Figure 10
Distributions of the heat transfer coefficient (a) and wall temperature (b) under different heat fluxes.
(a)(b)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/m2, 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 Figure11.Figure 11
Distributions of the heat transfer coefficient (a) and Reynolds number (b) under different outlet pressures.
(a)(b)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.Figure12 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.Figure 12
Distributions of the heat transfer coefficient under different mass flow rates.


## 3.3. Cooling Performance Estimates

The effects of various operating mode parameters on the cooling performance were 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.Table 1
Design space for design of experiments.
ParametersNomenclatureValue rangeHHeat flux (MW/m2)1~1.5PoutOutlet pressure (MPa)4~6qmMass flow rate (g/s)4~6In 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].
(8)y=c0+∑i=13cixi+∑i=13ciixi2+∑i=12∑j=i+13cijxixj,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.Figure 13
Pareto graphs for the wall temperature (a), bulk temperature (b), and mass fraction (c) of unreacted kerosene at outlet.
(a)(b)(c)It can be concluded from Figures13(a) and 13(b) that main factors have same effects on Tw and Tb at outlet. H and qm devote about 35% positive and negative contribution, respectively, which means that both of the two performance indexes always increase as H increases or qm decreases. Second orders of H and qm also have notable effects on Tw and Tb at outlet, but the interaction effect of H and qm on them is adverse (-6.4% for Tw and 4.8% for Tb, resp.). However, the effects of other factors about Pout are nearly negligible, which is consistent with the results mentioned above. As shown in Figure 13(c), the positive effects of qm and its interaction with H on ηu at outlet are the most remarkable (27.6% and 18.3%, resp.), while H and the second order of qm have the most negative effect (-23.9% and -19.9%, resp.). The remaining factors contribute less than 5% effect on ηu at outlet.


## 4. 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



---

*Source: 8813929-2020-10-28.xml*