ILES of an array of three subsonic counter-flow jets issuing from a wing leading edge exposed to hypersonic aerodynamic heating

The flow of an active thermal protection system exploiting subsonic counter-flow jets for wing leading edges of hypersonic vehicles is numerically studied on the basis of the three dimensional Navier-Stokes equations. The coolant air issuing from around the stagnation point as an array of three jets spreads over both the upper and the lower sides of the cylinder surface and about 40 ~ 60% cooling effectiveness is achieved in the range up to 5 degrees angle of attack despite the occurrence of various three-dimensional fluid-dynamic instabilities. The numerical scheme is second order accurate but simple inclusion of high order polynomial approximation in the reconstruction enables the capturing of finer structure of the flow field.


Introduction
The development of active thermal protection system (TPS) exploiting counter-flow jets for nose caps and wing leading edges of hypersonic vehicles is an old and new problem, which has been studied since the early days of space program [1][2][3][4][5][6][7][8]. Still fresh in our minds is the Antipode (https://edition.cnn.com/travel/article/antipodehypersonic-concept-plane/index.html), a concept design of a hypersonic commercial airplane, which is equipped with supersonic counter-flow jets in the long penetration mode (LPM) for reducing drag and heat and flies from New York to London in 11 min. Most of the existing works are for the case of supersonic jets and such an active TPS, however, still remains a matter of proposal. One of the crucial reasons lies in the simple fact that supersonic jets consume a large amount of coolant gas even in the case of LPM. This paper is written as a newsletter about our latest study concerning the feasibility of active TPS exploiting subsonic counter-flow jets for wing leading edges of hypersonic vehicles. In our previous study [9] a uniform hypersonic flow past a cylinder with a subsonic counter-flow jet issuing from around the stagnation point was numerically studied on the basis of the two-dimensional compressible Navier-Stokes equations. The findings in the previous study are summarized as follows; (i) high cooling effectiveness, about 90%, is achieved in the case of zero angle of attack despite the occurrence of the Kelvin-Helmholtz (KH) instability around the interface between the high temperature air coming from behind a bow shock wave and the coolant air; (ii) when the angle of attack exceeds 2.5 degrees, the coolant gas over-concentrates in the upper side of the cylinder and the lower side is not cooled at all; (iii) the vulnerability of a single subsonic jet against transverse flows can be overcome without an increase in the total mass flow rate of the coolant air by dividing the single jet into an array of three jets; (iv) both sides of the cylinder are covered with the coolant air and comparable high cooling effectiveness is achieved in the range up to 5 degrees angle of attack. Although these results are seemingly promising towards the realization of active TPS exploiting subsonic counter-flow jets, however, they do not necessarily assure such high performance in actual situations, where the occurrence of three-dimensional instabilities such as the secondary KH instability and turbulence transition is expected; the mixing of the high temperature air and the coolant air is enhanced, and consequently, the cooling performance decreases. Some results of our latest three-dimensional numerical simulations are reported here.
Since the flow field involves a strong bow shock wave, the employment of a shockcapturing scheme in the numerical simulation, i.e. the implicit large eddy simulation (ILES), seems reasonable. It is well known that numerical solutions with strong shock waves often exhibit anomalous behaviours such as carbuncle phenomenon and postshock oscillations and the prevention of such pathologies requires a lot of numerical dissipation, which, on the other hand, is very harmful for the clarification of detailed behaviors of flows with instabilities. The numerical scheme employed in the present ILES should yield sufficiently large dissipation around the bow shock wave. At the same time it should be less dissipative in smooth regions, although numerical dissipation is needed to regularise the numerical solution at the mesh scale. High-order accurate shock-capturing schemes with advanced reconstruction techniques such as WENO are typically preferred in ILES. One of the authors (T.O.) and his students have developed a simple high-order accurate shock-capturing scheme, the formal accuracy of which is 5th order in space and 4th order in time [10]. It satisfies the abovementioned seemingly antipodal properties and its performance is largely comparable to that of WENO5-Rusanov-RK4. Although its computational cost is 30~40% less than that of the advanced scheme, the computation is still costly when the resolution is high. There is room for argument about whether a low-order accurate scheme with a fine grid or a high-order accurate scheme with a coarse grid is more appropriate in ILES. The aim of the present study is the reliable prediction of the cooling effectiveness in more realistic situation and it requires a fine grid especially near the cylinder surface, where turbulence transition is expected to take place. And yet demanding is the use of a fine grid sufficient for the achievement of the expected asymptotic behavior of the numerical solution towards convergence in this engineering problem; it is recommended for reliable DNS's of wall turbulent flows that the mesh size in the wall normal direction is smaller than the viscous thickness in near wall region but we are not able to carry out the simulation with such a high resolution because of budget limitation. Under the situation, we adopt the former approach in the present study. The original high-order accurate scheme is simplified and the computational cost is drastically reduced in exchange for the decrease in the formal accuracy to second order. But the use of high order polynomials in the reconstruction is retained to capture detailed structure of turbulent behavior of the flow.

Problem and assumptions
Consider a uniform hypersonic flow of air past a circular cylinder. The radius of the cylinder is 120 mm and its axis is taken as the z axis of the xyz Cartesian coordinate system; the upstream uniform flow is in the direction of (x, y, z) = (cosΩ, sinΩ, 0); the coolant gas (air) is injected through three channels of a height of 4 mm, which are embedded in the cylinder with a spacing of 4 mm; and the three jets are in the negative x direction (Ω corresponds to the angle of attack, see Fig. 1). We carry out the numerical simulation of the flow on the basis of the three-dimensional compressible Navier-Stokes equations for an ideal gas; no real gas effect is taken into account; the heat capacity ratio is 1.4; the Prandtl number is 0.72; and the temperature dependence of the viscosity is assumed to obey Sutherland's law. The no-slip boundary condition and nojump condition are applied to the flow velocity and the temperature, respectively, at the surface of the cylinder and the side-walls of the channels. The isentropic condition is imposed around the inlets of the jets and simple extrapolation is employed as the outflow condition.

Numerical method
In the shock-capturing scheme of Ref. [10], the inviscid part of the numerical flux is computed as a convex combination of three types of fluxes, F A , F D and F C . The F A is dissipative and F D and F C are less dissipative; F D can capture a stationary contact discontinuity without inner points. The arguments of F A and F D are computed by MUSCL with the van-Leer slope limiter and will be hereafter denoted by h L MUSCL and h R MUSCL , where h stands for the primitive variables, i.e. h = ρ, u, v, w, P (density, flow velocity, and pressure). The arguments of F C are computed from fifth order accurate Lagrange's polynomial approximations of primitive functions for conservative variables; the scheme is developed for structured grids. The F A is for shock capturing, F D is for contact-discontinuity capturing and F C is for smooth regions. The reader is referred to Ref. [10] for the definitions of these numerical fluxes. The weights in the combination are controlled accordingly. The numerical flux for the Navier-Stokes equations is simply made by the addition of the physical dissipation, which is approximated by a conventional second-order accurate central finite difference formula.
In order to reduce computational cost, we introduce several simplifications into the original high-order accurate shock-capturing scheme. The reconstruction is based on cell-averaged values in the simplified scheme. The number of Gauss points needed for the integration of the numerical flux over cell-face is reduced (from four in 3D case) to one and the formal accuracy of the time marching is reduced to second order (RK-2). The inviscid numerical flux is computed as a convex combination of F A and F D .
where α is a weight in the range of 0 ≤ α ≤ 1. We prepare two schemes. The first scheme, (S1), the arguments of F A and F D are common and are h L MUSCL and h R MUSCL . The second scheme, (S2), also uses h L MUSCL and h R MUSCL as the arguments of F A but those of F D are modified. Let the arguments of F D in (S2) be denoted by h L S2 and h R S2 and let h k be the value of h at the cell-centre of the index k, which is computed from the cell-average of the conservative variables. A quartic function generated from h k where β is a weight in the range of 0 ≤ β ≤ 1. The weights α and β are computed at each cell-face according to the physical situation and the smoothness of solution. The guiding principle is as follows: (i) F A should be dominant at cell-faces around shock waves and F D should be so except around them; (ii) h H MUSCL should be dominant in h H S2 at cellfaces in regions where h changes abruptly in the scale of cell width, e.g. contact discontinuities, and h H Quartic should be so in regions where h changes slowly in the scale of cell width. In the present study the following formulas for α and β are adopted.
Shimada and Ohwada Advances in Aerodynamics (2020) 2:12 where h H (h = ρ, u, v, w, P, H = L, R) are computed by MUSCL (the subscripts MUSCL is omitted for concise expression), T Ã ¼ T L þT R 2 (T is the temperature, computed from ρ and P by EOS, e.g. T = P/ρ), and C k (k = 1, ⋯, 6) are positive constants [C k = 10 (k = 1, ⋯, 6) in the present study]. Since the weight β is determined from the magnitude of the jump j h L MUSCL −h R MUSCL j, the numerical dissipation of F D could be insufficient for the elimination of waves with high wave numbers such as those with zigzag wave forms when their amplitudes are very small; the jump jh L S2 −h R S2 j could be much smaller than their amplitudes for such waves. This is in contrast to the original scheme.

Results
The physical condition for the present numerical simulation is as follows. The upstream corresponds to a Mach 5 cruise at an altitude of 30 km; the density, pressure and temperature are 0.015kg/m 3 , 1000 Pa and 230 K, respectively; the Reynolds number based on the upstream condition and the radius of the cylinder (120 mm) is 2 × 10 5 . The wall temperature (the cylinder and the channels) is fixed at 400 K. The reservoir condition of the jets is (3.29 × 10 4 Pa, 400 K); the reservoir pressure is slightly larger than the stagnation pressure of the ambient air in the case without TPS; the temperature of the ambient air is about 1500 K. The computational condition is as follows. The region for z is limited to 0 < z < 30 [mm] and the periodic boundary condition is imposed at the boundaries. Multi-block structured curvilinear grids are employed. Three different resolutions, namely Coarse, Midst and Fine, are prepared; the number of total cells is about 2.7 × 10 7 , 1.3 × 10 8 and 4.3 × 10 8 in Coarse, Midst and Fine, respectively. The minimum cell size in the direction normal to the cylinder surface is 4.8 × 10 −2 mm, 2.9 × 10 −2 mm and 1.9 × 10 −2 mm in Coarse, Midst and Fine, respectively, and it was about three times larger than the viscous length scale at θ = 30 even in the case of Fine (and thus the present simulation is not regarded as DNS). The distribution of the time-averaged (mean) shear stress over the cylinder surface was nearly the same as that in the case without jets except around the exit; that of the mean heat flux was not very sensitive to scheme and resolution, although the statistical convergence was not confirmed clearly in the present limited computation. Figure 2 shows a snapshot of the overall temperature field in the plane of z = 15 mm; figures of the numerical results in this paper are for Ω = 0 (Scheme-S2, Midst) unless otherwise stated. The lower and middle jets bend downward and the upper jet does upward; the bend direction of the middle jet depends on the initial data of the simulation. The Mach numbers of the jets were about 0.1 (upper), 0.07 (middle) and 0.15 (lower) and the jet Reynolds numbers based on the channel height were in the range of 15003 300. The total mass flow rate of the coolant air per unit length was about 0.13 kg/ms. These three jets were flapping at the frequency around 2 kHz while the bow shock was almost stationary. The enlarged view around the exits of the three jets is shown in Fig. 3. It is seen that mushroom spikes are formed around jet boundaries. It is considered that the formation of these spikes are due to the occurrence of the Rayleigh-Taylor (RT) instability and the subsequent KH instability; the jets are bent abruptly and the inertial force is acting to the jets in the direction from the heavy (cold) side to the light (hot) side. These mushroom spikes are advected downstream and consequently form longitudinal vortex rolls.
Plotted in Fig. 4 is the time history of the cooling effectiveness η defined by η = 1 − Q J /Q 0 , where Q 0 is the heat transmitted to the cylinder in the range (−60 < θ < 60 [deg], 0 < z < 30 [mm]) per unit time in the case without TPS and Q J is the one in the case with jets. It is about 0.47 at the time of Point A and 0.6 at the time of Point B; irrespective of scheme and resolution, η varied in this range. This is in contrast to the case of the two-dimensional simulations, where η~0.9 is achieved. Figures 5 and 6 Fig. 7. It is seen in the right part of the figure that the high temperature air  Fig. 8. A longitudinal vortex roll is located on the left side of the middle jet and high temperature air is dragged between them; three separation bubbles are formed just downstream of the three jets and the corresponding jet boundaries act as the separated shear layers of these bubbles. The structure of these separation bubbles is not two dimensional; the formation of twodimensional KH vortex rolls was not confirmed although some spanwise vortices were locally observed around the separated shear layer of the upper bubble. A stagnation region is formed inside the middle bubble. Early breakdown to turbulence (bypass transition) was clearly confirmed inside the lower separation bubble. Turbulent boundary layers were developed downstream of the upper and lower bubbles. Figure 9 shows the distribution of the circumferential velocity in the vicinity of the upper cylinder surface.  The formation of near-wall streaks is clearly seen. The radial component of the flow velocity was positive around stripes with smaller circumferential velocity and vice versa around stripes with larger circumferential velocity. The air in the former stripes came from the cold side (the coolant air) and the latter did from the hot side (the high temperature air behind the bow shock) and this resulted in the formation of the stripe pattern of the heat flux distribution (Figs. 5 and 6). The mixing in the turbulent boundary layer also allowed occasional penetration of the high temperature air into the bottom of the boundary layer, which caused an abrupt increase in the heat flux. Figure 10 shows the radial distributions of the averaged circumferential velocity at θ = 30, 60 [deg], which is depicted using the wall scaling. Incidentally, the results of the present simulations are appreciably above the standard plot in the log-law region. The discrepancy is not resolved in the case of θ = 60 even after the transformation of Trettel and Larsson [11], although it yields better results than Van Driest transformation [12]. On the other hand, instantaneous flow field was very sensitive to resolution and scheme. Figures 11, 12 and 13 show the comparison of Q-criterion. It is seen that finer turbulent flow structure is captured with higher resolution and (S2) yields much finer results than (S1). The fineness of the result of (S2, Coarse) or (S2, Midst) looks comparable to that of (S1, Fine), while the ratio of the computational cost of (S2, Coarse) to that of (S1, Fine) is 0.06 and that of (S2, Midst) is 0.4. In other words, the simple inclusion of high order polynomials in the reconstruction enables the capture of a finer instantaneous flow structure at a cheaper cost.
Finally, the results for other cases are briefly mentioned. The cooling effectiveness increased to about 70% when the mass flow rate of the jets increased to 2.6 times. The upper and middle jets bent upwards and the lower jet bent downwards and overconcentration of the coolant gas was not observed in the case of Ω = 5 [deg]. The cooling effectiveness, however, appreciably decreased (η = 0.4~0.5) though. Incidentally, the computations were performed on the Cray XC40 of Academic Center for Computing Media Studies of Kyoto University.

Concluding remarks
The present 3D simulation of the active TPS exploiting an array of three subsonic counter-flow jets exhibited various fluid-dynamic phenomena such as jet flapping, RT and KH instabilities, bypass transition in separation bubbles, turbulent boundary layers and heat flux distribution of stripe pattern caused by near-wall streaks. The cooling effectiveness was about 50~60%, which was in contrast to the case of the twodimensional simulations, predicting high cooling effectiveness of about 90%. The turbulent mixing was one of the major reasons for the decrease in the performance. Besides, Fig. 9 The distribution of near-wall circumferential velocity (Point A)

Shimada and Ohwada Advances in Aerodynamics
(2020) 2:12 longitudinal vortex rolls generated by the RT and KH instabilities occasionally wiped out the coolant air on a large scale and dragged the high temperature air there, which also caused the performance decrease. The mean velocity field was not very sensitive to resolution and scheme, which was in contrast to the instantaneous flow field. By including high order polynomial approximation in the reconstruction, a finer structure of instantaneous flow field was captured. Even after the drastic simplification of the numerical scheme, the computation of the present problem is still huge. Because of the limited budget, we could not reveal the statistical property of the turbulent flow field fully. The correlation between the velocity fluctuation and the mean temperature field was not clarified because the statistical convergence was quite slow except the mean Fig. 10 The radial distribution of the averaged circumferential velocity. MUSCL means (S1) and MUSCL+Quartic does (S2) Fig. 11 Isosurface of Q-criterion. Left: (S1, Coarse), Right: (S2, Coarse). The surface is colored according to the magnitude of the circumferential velocity for reference velocity field. Nevertheless, the present study provides some enlightening information about the feasibility of the active TPS exploiting subsonic counterflow jets. It also demonstrates that ILES approach is very promising as an engineering tool for analyses of turbulent flows especially with strong shock waves, although further investigation of the validity of the less effortful simulation method in standard test cases such as homogeneous isotropic turbulence and channel flows seems to be necessary and is left as our future subject.