A decomposition formula for the wall heat flux of a compressible boundary layer

Understanding the generation mechanism of the heat flux is essential for the design of hypersonic vehicles. We proposed a novel formula to decompose the heat flux coefficient into the contributions of different terms by integrating the conservative equation of the total energy. The reliability of the formula is well demonstrated by the direct numerical simulation results of a hypersonic transitional boundary layer. Through this formula, the exact process of the energy transport in the boundary layer can be explained and the dominant contributors to the heat flux can be explored, which are beneficial for the prediction of the heat and design of the thermal protection devices.


Introduction
The heat transfer prediction is of great importance for hypersonic vehicles. When the transition or turbulence happens, the accurate prediction becomes even more challenging. Therefore, exploring the generation mechanism of the heat flux has drawn wide attention. Many efforts have been put into understanding the generation mechanism of the heat flux, which provides important guidance for the design of the thermal protection system and the thermal management [1]. Some efforts took advantage of the similarities between the generation of friction and heat flux, and constructed the Reynolds analogy [2] Ra = 2St/C f to connect friction with heat flux, where St is the Stanton number, C f the skin friction coefficient. Hopkins and Inouye [3] predicted the surface heat flux of the hypersonic boundary layer with Ra = 1. However, the accuracy of this correlation decreases significantly on a cold wall at high Mach numbers. More intrinsic mechanisms should be considered. Huang et al. [4] derived the formula q w ¼ −u b τ w by the assumption that the heat transfer into the walls equals the total pressure work done across the channel. Chen et al. [5] established the exact relations for the skin friction with other dynamic and kinetic quantities, and they found these relations revealed that the skin friction is intrinsically coupled with the surface temperature through the heat flux. Abe and Antonia [6] proposed a simple relation between the scalar dissipation rate and the wall heat transfer coefficient for the channel flow. Kim et al. [7] proposed a direct approach for the time-dependent heat flux by assuming the temperature approximated as a third-order polynomial of position. Ebadi et al. [8] obtained the wall heat flux by a triple integration of the Reynolds averaged energy equation for the two-dimensional incompressible turbulent boundary layer. Fukagata et al. [9] studied the relationship of the wall heat flux and the Nusselt number in the incompressible channel. And the similar relation was obtained by Liu [10] when studying the compressible pipe flow. Ghosh [11] proposed that the wall heat flux can be expressed by the sum of the integrations of different viscous terms in the incompressible turbulent channel and pipe flow. Zhang and Xia [12] further proposed a formula to assess the contributions of the viscous stresses to the heat transfer in a turbulent channel flow by the method of Fukagata et al. [13], but the turbulent stresses were missing in their expression due to the simplification assumptions used in the channel flow.
Previous studies have not built the relation of the generation mechanism of the heat flux to the dynamic energy transport of the compressible boundary layer. Therefore, a more reliable and accurate method to gain insights into the generation mechanism should be derived. In the paper, we propose a new decomposition formula for the wall heat flux by integrating the conservation equation of total energy, which can explain the energy transport process in the boundary layer and reveal the main factors affecting the wall heat flux. As far as the authors' knowledge, this is the first decomposition formula proposed for the wall heat flux of the compressible boundary layer. The new formula is derived by the method of Renard et al. [14], which is more physical and feasible than those [11,12] based on the method of Fukagata et al. [13] It serves an insight into the complex transport processes of the wall heat flux, and helps us find the key factors.
The work of this paper is organized as follows. In section 2, we describe numerical methods and case setup in brief. In section 3, the detailed derivation of the decomposition formula is presented. In section 4, the validation of the DNS of a hypersonic transitional boundary layer is performed, and the proposed decomposition formula is applied to analyze the heat flux. Finally, some conclusions are drawn in section 5.

Numerical methods and case setup
To analyze the heat flux decomposition of the hypersonic transitional boundary layer, a direct numerical simulation is performed. The compressible Navier-Stokes equations in the curvilinear coordinate are adopted as governing equations.
Where u 1 , u 2 and u 3 are streamwise, normal and spanwise velocities, respectively. And p and ρ are pressure and density. The expressions of ρe, σ ij , and q j are defined as The usual indicial notation is used. Prandtl number Pr is set as 0.72 and the specific heat ratio γ is 1.4. The dynamic viscosity μ is obtained by using Sutherland's law. Ma ∞ is the freestream Mach number and Re ∞ is the Reynolds number. The working fluid is air, and the gas model is assumed to be the perfect gas model.
An in-house code with high-order schemes is employed to perform the DNS. This code has been applied in many DNS simulations of compressible turbulent cases, including compressible homogeneous turbulence [15], turbulent boundary layer [15] and shock wave/boundary layer interaction [16][17][18]. The accuracy and robustness have been well validated. The hybrid optimized WENO scheme is adopted to discretize the inviscid fluxes with a novel discontinuity detector [19]. In the regions with discontinuities, a seventh-order WENO scheme [20] is activated and in the smooth region, a fourth-order bandwidth-optimized upwind-biased scheme [19] is used to resolve the small structures in the turbulent boundary layer. The viscous terms are discretized by the fourth-order central scheme. The thirdorder TVD Runge-Kutta method is adopted as the temporal algorithm. Ma ∞ is set as 6.0, and Re ∞ is 12,000 which is based on the millimetre and freestream parameters. The freestream temperature is 65 K and the wall temperature is 305 K. The millimetre is adopted as the reference length.
The sketch of the computational domain is presented in Fig. 1. The lengths in the streamwise, normal and spanwise direction are L x = 416, L y = 35 and L z = 14, respectively. The computational domain is discretized with N x × N y × N z = 1151 × 320 × 149 points. The grid points are equally spaced in the spanwise direction, and clustered near the wall in the normal direction. The grid spacing normalized by the wall unit in the three directions is dx + = 7.9, dy þ w ¼ 0:36 and dz + = 3.4, respectively. A hypersonic laminar profile is fixed at the inlet. The outlet boundary condition is a supersonic outflow with a sponge layer to further suppress the disturbances originated from the outlet. The upper boundary is a non-reflective condition. And an isothermal nonslip condition is applied at the wall. The periodical conditions are set on both sides in the spanwise direction. To trigger a bypass transition, a blow and suction forcing method is employed, which sets a normal Fig. 1 The sketch of the computational domain velocity disturbance v n = Af(x) cos (2πλ z /L z ) at the wall, where f(x) is used to control the streamwise extent of force [21], λ z is the spanwise wavelength and is set as 4.0 in the present simulation. The forcing region ranges from x = 50 to x = 55.

Derivation of the decomposition formula
Under the assumption of homogeneity in the spanwise direction and nonslip condition at the wall, the time-averaged conservative equations for the specific total energy can be expressed as, when ϕ is an arbitrary variable, ϕ denotes Reynolds average,φ Favre average and ϕ ′′ is the fluctuations concerning the Favre average.Ẽ denotes the specific total energy, q L,j the heat conduction, q T,j the turbulent transport of heat, D j the molecular diffusion, T j the turbulent transport of turbulent kinetic energy (TKE), MS j work by the molecular stresses and RS j work by Reynold stresses. In Table 1, the specific expressions of the terms at the right hand of Eq. (5) are presented. More details about Eq. (5) can be found in the book of Wilcox [22]. The work of Renard [14] and Li [23] about skin friction decomposition provides good hints for building the integration relation of heat flux. In order to isolate the effects of the energy transport terms, the absolute reference frame is employed. The reference absolute frame is attained by assuming that the wall moves at the speed U ∞ . The expressions of the time t a , coordinates x a , y a , velocities u a , v a , pressure p a and the density ρ a in the reference absolute frame satisfy where the subscript 'a' represents 'absolute' variables under the reference frame. By substituting Eq. (6) into Eq. (5), Eq. (5) takes the form Table 1 Expressions for the terms at the right hand of Eq. (5) Work by Reynolds stresses By multiplyingũ a at both sides of Eq. (7), combining the continuity equation and moving the heat fluxes in the normal direction to the left side and all the other terms to the right, the following equation could be obtained, Integrating Eq. (8) from the wall to the farfield, the heat flux on the wall (y = 0) can be expressed as the following form.
½ũ a ∂ ∂y a ðD y;a þ T y;a þ MS y;a þ RS y;a Þdy Finally, Eq. (9) can be transformed into the initial frame and the heat flux coefficient is obtained by introducing the definition of the coefficient In Eq. (10), the heat flux coefficient C h,decomp has been decomposed into seven parts, i.e., C h,1 to C h,7 . The first part C h,1 represents the contribution of heat conduction, the second part C h,2 the turbulent transport of heat, the third part C h,3 the normal component of the molecular diffusion, the fourth part C h,4 the turbulent transport of turbulent kinetic energy (TKE), the fifth part C h,5 the work by the molecular stresses and the sixth part C h,6 the work by Reynolds stresses. And the last part C h,7 contains the streamwise heterogeneity, pressure work and the variation of the specific total energy with time.

Decomposition results of the wall heat flux
In this section, the reliability and accuracy of the new decomposition formula will be demonstrated by the DNS results of a hypersonic transitional boundary layer. The contributions of different energy transport processes such as the heat conduction, the turbulent transport of heat, the work done by molecular stresses will be calculated. Moreover, the key normal locations where the structures impact the wall heat will be assessed by the integrand functions of these contributions.

DNS results of the hypersonic transitional boundary layer
The instantaneous vortical structures displayed by Q criterion [24] are presented in Fig. 2. The transition process can be observed. Three streamwise locations (P1, P2 and P3) are chosen to validate the decomposition results, which are marked by shade planes. In Fig. 3, the mean streamwise velocity after van Driest transformation at P3 is presented and compared with DNS results of hypersonic results by Priebe [25]. Good agreement is obtained. The discrepancy in the outer part of the boundary layer is caused by the different friction velocity at the wall. In Fig. 4, the Reynolds stresses normalized by the friction velocity are presented. Several supersonic and hypersonic results [26][27][28] are also displayed. The information of the setup of reference simulations is presented in Table 2, where δ ref is the nominal boundary layer thickness in the fully developed region. The present data generally agree with those in the reference results. These comparisons demonstrate that the performance of our present direct simulation is acceptable, and the obtained data can be used for decomposition analysis.

Contributions of different energy transport processes
The contributions of different energy transport processes to the wall heat flux at different streamwise locations, e.g. P1-P3 are calculated according to Eq. (10). The ratios of the contributions of the terms in Eq. (10) to the reconstructed heat flux coefficient C h,decom and the relative errors are presented in Table 3. The relative error is defined as Error = (C h,decom − C h,0 )/C h,0 × 100%, where C h,0 is the time-averaged heat flux coefficient. The relative errors are very small in both transitional and turbulent regions. In Table 3, it is found that the work done by the molecular stresses (C h,5 ) and Reynolds stresses (C h,6 ) plays dominant roles in heat production. The sum of the contributions of these two terms takes over 2.7 times of C h,0 in the transition region and 1.9 times in the turbulent region. The contributions of the heat conduction (C h,1 ) and the turbulent transport of heat (C h,2 ) are negative, which indicates both processes act as the transporters of the heat and carry the extra heat to the outer regions of the boundary layer.
The magnitude of C h,1 and C h,5 gets greater in the turbulent region due to the large gradient of the profiles. On the other hand, C h,2 and C h,6 in the transitional region are much greater than those in the turbulent region, because of the quick nonlinear development of the structures. The trends of the contributions can be further assessed when more sample locations are considered. The trends of the contributions in the transition process are presented in Fig. 5. In addition, the instantaneous density gradient and the time-averaged heat flux coefficient are also presented to show the overall evolution procedure. The symbols in the line of C h,i indicate different sample locations (from left to right labelled as X0 to X10). The first sample point (X0) is located in the laminar region; therefore, the work of the molecular stress and the heat conduction is the main factor affecting the wall heat. After entering the transition region, the effects of the Reynold stresses become dominant. Moreover, a rapid variation of C h,6 is observed in the transition region due to the nonlinear development of the coherent structures. While in the turbulent region, an equilibrium state is reached and a smooth variation is obtained. The trends of contributions of the heat conduction C h,1 and the molecular stresses C h,5 are very smooth in both transition and turbulent regions, which indicate they are not sensitive to the change of the turbulent fluctuations. In addition, the effects of the molecular

Key normal locations affecting the wall heat
The key normal locations where the structures affect the wall heat can be attained by analyzing the integrand functions of the terms in Eq. (10). Two main transporters, i.e., the heat conduction (C h,2 ) and the turbulent transport (C h,3 ) and two main contributors, i.e., the work by molecular stresses (C h,5 ) and Reynolds stresses (C h,6 ) are analyzed in this section. The integrand functions of the heat conduction (C h,2 ) and the turbulent transport (C h,3 ) are presented in Fig. 6. The lines in the figure denote the locations in the transitional region and the symbols are the locations in the turbulent region. The contribution of C h,2 is positive when the normal locations are near the wall due to the positive temperature gradient. When y + reaches 10, a valley is observed, which is formed by the decrease of the temperature after the peak of the mean temperature. Compared with the integrand function of C h,2 , two valleys are observed in the integrand functions of C h,3 . The first valley is located at y + ≈ 65, the magnitude of the valley reduces and finally disappears as the locations move downstream. The second valley is located at around y + = 10, which corresponds to the buffer layer of the turbulent boundary layer. The valley will get larger as the locations move downstream. The key normal locations of the work by molecular stresses (C h,5 ) and Reynolds stresses (C h,6 ) are analyzed in Fig. 7. The positive contribution of C h,5 is mainly located in the regions near the wall. As the normal location increases, the integrand function will decrease and a valley is formed around y + = 10. Meanwhile, a peak is observed around y + = 10 in the integrand function of contribution of C h,6 . The peak is located in the buffer layer of the boundary layer. In addition, in the transition process, the magnitude of the peak will get larger as the streamwise distance increases.

Conclusions
In this paper, we proposed a new decomposition formula for the wall heat flux. And the performance of the formula has been well demonstrated by DNS results of a hypersonic transitional boundary layer. Fig. 6 The distribution of the integrand functions of the contributions of (a) the heat conduction C h,2 and (b) the turbulent transport C h,3 (1) Through this formula, the wall heat flux can be decomposed into contributions of seven terms, i.e., the heat conduction, turbulent transport of heat, molecular diffusion, turbulent transport of TKE, molecular stresses, the Reynolds stresses and the streamwise heterogeneity. (2) The contributions of each term can be calculated quantitatively. For the present case, it is found that the heat flux produced by the work done by Reynolds stresses and molecular stresses is much higher than the time-averaged heat flux on the wall. The heat conduction and turbulent convection will carry the extra heat into the outer part of the boundary layer. (3) The structures in the buffer layer (y + = 10) play a dominant role in the production and transfer of the heat flux for the present simulation. Moreover, the contributions of the heat conduction and the molecular stresses are mainly affected by the gradient of the temperature. (4) As the assumption used in the derivation is only the spanwise homogeneity and nonslip wall, this formula can be applied in the analysis of the heat transfer of a hypersonic transitional/turbulent boundary layer at high Mach numbers, which can be employed to identify the main factors affecting the wall heating and provide good guidance for the design of the thermal protection system.