Direct numerical simulation of hypersonic boundary layer transition over a lifting-body model HyTRV

Direct numerical simulation (DNS) of transition over a hypersonic lifting body model HyTRV developed by China Aerodynamics Research and Development Center is performed. The free-stream parameters are: the free-stream Mach number is 6, the unit Reynolds number is 10000/mm, the free-stream temperature is 79 K, the angle of attack is 0, and the wall temperature is 300 K. Weak random blowing-and-suction perturbations in the leading range are used to trigger the transition. A high order finite-difference code OpenCFD developed by the authors is used for the simulation, and grid convergence test shows that the transition locations are grid-convergence. DNS results show that transition occurs in central area of the lower surface and the concaved region of the upper surface, and the transition regions are also the streamline convergence regions. The transition mechanisms in different regions are investigated by using the spectrum and POD analysis.


Introduction
Hypersonic boundary layer transitions are of great importance for aerodynamic issues in hypersonic vehicle design [1]. The prediction of transition location is the basics for thermal protection design, which is a popular research problem currently. To obtain the transition locations, American National Aeronautics and Space Administration (NASA) has experimented in wind tunnels and flights, and then the empirical formula for the transition location is obtained [2]. In the transition prediction methods, the e N method is the most used transition prediction method. At present, the theory of transition prediction is still far from perfect.
Generally, the transition of the boundary layer is caused by the modal instability of the boundary layer, and its instability mechanism has been studied for a hundred years. The unstable modes mainly include the first/second modes, the cross-flow instability modes, the streamwise vortex instability modes, the Görtler modes and the attachment line instability modes, etc. [3] The cross-flow instability often occurs on the surface with cross-flow, and it plays an important role in the transition. Recent researchers more and more focus on the cross-flow instability mode. Dong et al. [4] gave a research about the receptivity of the roughness elements to the cross-flow instability. Lu et al. [5] studied the receptivity mechanism of cross-flow instability modes induced in the three-dimensional boundary layer. Lee et al. [6][7][8][9] also have a lot of research about the mechanism of the transition of the boundary layer using various experimental methods.
In the early days, research on hypersonic transitions was limited to simple shapes such as flat plates, followed by blunt cone shapes. Li et al. [10] carried out the earliest direct numerical simulation of the transition of the hypersonic blunt cone boundary layer with a small angle of attack, and gave a non-monotonic transition curve on the cone. As a typical high-speed aircraft head shape, the study of blunt cone/sharp cone hypersonic transition has attracted a lot of attention in recent years, and wind tunnel experiments and DNS research have also been gradually developed. By using the e N method, Su et al. [11] found that the Mack second modes dominate the transition of the Mach 6 blunt cone boundary layer with an isothermal wall boundary. Lots of wind tunnel and flight experiments are performed to measure the transition locations on supersonic/hypersonic blunted or sharp cone boundary layers [2,12].
With the development of the design concept of hypersonic aircraft, lifting bodies are important shapes of recent high-speed aircraft, and the prediction, mechanism, and control methods of the boundary layer transition of hypersonic lifting bodies have become the requirements of aircraft design [13]. There are some researches about lifting bodies. Wang et al. [14] used the Reynolds Averaged Navier-Stokes (RANS) equations to study the aerodynamic and aerothermal character of an X-33-like lifting body. TJ Horvath et al. [15] introduced a 0.02-scale model of a lifting body concept for possible application to the Assured Crew Return Vehicle from Space Station Freedom and the model was tested at Mach 6 and 10 in air. Y Minami et al. [16] gave an introduction about LIFLEX (Lifting-body Flight Experiment), which is a flight experiment using a small and law-cost vehicle. Takagi et al. [17] performed research about the influences of different modes on the transition. China Aerodynamics Research and Development Center completed the MF-1 flight test of hypersonic transition research [18][19][20], which mainly studied the transition problem caused by the second mode instability of the cone boundary layer. Compared with the blunt cone, the shape of the lifting body is more complicated, and its boundary layer has multiple complicated factors such as cross flow, reverse pressure gradient, separation, entropy layer interference, and the transition mechanism of the lifting body is also more complicated. Due to the lifting body's complex shape, high Reynolds number and Mach number, and difficulty of calculations, there is currently little literature that reports on the transition of the shape of hypersonic lifting bodies, especially, within the knowledge of the authors, there is still no DNS study on hypersonic lifting bodies.
In this article, to get more information of the transition mechanism for boundary layers with complex shapes, we performed a direct numerical simulation (DNS) on a lifting body model Hypersonic Transition Research Vehicle (HyTRV) designed by China Aerodynamics Research and Development Center [6], and simulated the wind tunnel test conditions of the lifting body. HyTRV has a lifting-body configuration with analytical profiles, which is designed for the studies of hypersonic transition. By DNS of HyTRV, we expect to obtain more complex mechanisms than normal flat-plates or blunt cones. The disturbance that triggers the transition was the random wall blowing-and-suction perturbations in the lifting body's leading. The skin friction and heat flux distributions of the lifting body are given by DNS, and the transition region of the lifting body surface is also shown. The hypersonic transition database of the DNS lays the foundation for the further analysis of the transition mechanism of the lifting body.

Configuration of the lifting body model HyTRV
The lifting body model HyTRV was designed by China Aerodynamics Research and Development Center [3]. Figure 1 shows the shape of the lifting body model. The whole length of the lifting body model is 1600 mm. Its head is a 1:2:1 ellipsoid, and its major axis and minor axis are 80 mm and 40 mm. The bottom half of the bottom section of the lifting body is a 1:4 ellipse, and the major axis and minor axis are 600 mm and 150 mm respectively.
The upper part of the lifting body is controlled by the Class/Shape Function Transformation (CST) curve and elliptic curve. The function curve of the upper part is: Where z is the vertical coordinate, α = 0.72, z E is the coordinate given by the 1:4 ellipse, z c is the coordinate given by the CST function, and the analytical formula of the CST function is as follows (where y is the horizontal direction and the z-axis is the vertical direction): Where Z x = 270 mm is the height of the upper surface curve, and Y x = 600 mm is the width of the bottom surface of the aircraft.

DNS settings
The incoming flow Mach number is 6, the unit length Reynolds number is 10000/mm, the free-steam flow temperature is 79 K, and the wall temperature of the lifting body is 300 K. The calculation is carried out in two steps: 1) Steady-state flow field calculation.

2) Boundary transition DNS calculation.
In order to save computing resources, we use half of the domain when calculating, and symmetric boundary conditions are used at the symmetry plane. The schematic diagrams of the computational grid and computational domain are shown in Figs. 2 and First, the steady-state (laminar flow) calculation is performed, and the calculation adopts the finite volume code OpenCFD-EC. The Van Leer splitting with MUSCL limiter is used for the convection terms, and the LU-SGS method is applied for the time advance. After the calculation reaches the time convergence, a steady flow field is obtained, which can provide initial values and boundary conditions for unsteady direct numerical simulation.
The schematic diagram of the DNS calculation domain in the symmetry plane is shown in Fig. 3. The calculation domain and grid of DNS are shown in Fig. 4. The boundary conditions of the entrance and upper boundary of the calculation domain are given by the calculation result of the first step (steady simulation). Isothermal non-slip boundary conditions are used on the wall, and non-reflective outlet boundary conditions are used for the outlet boundary.
A high-order finite difference program OpenCFD-SC [10] developed by the authors is applied for DNS calculation. The convection term is discrete by using Steger-Warming splitting and an optimized 6th-order monotonicity-preserving scheme (OMP6) developed by Li et al. [21] and the viscosity term is discrete by using a 6th-order central finite difference scheme, and the 3rd-order TVD type Runge-Kutta scheme is used for time advance. In order to enhance the stability of the calculation, the flow field is filtered every 10 steps [22]. Numerical tests show that this filter can effectively remove high-frequency non-physical oscillations in the flow field to ensure the stability of the simulation. In order to trigger the transition, a random blowing and suction disturbance was added to the wall surface of the aircraft in the range of 50 mm-60 mm from the leading, and the amplitude was 2% of the free-stream flow velocity. This disturbance simulates the perturbations caused by the wall roughness. The author's previous calculation results for blunt cones show that for wall blowing and suction disturbances with a small amplitude, the transition position is not sensitive to the specific form of the disturbance, and the transition position has a certain universality [2].

Verification and analysis of calculation results
In order to verify the grid convergence of DNS calculation, four sets of grids are used in the calculation, and the grid parameters are shown in Table 1. Among them, Case 1 is a coarse grid, Case 2 doubles the circumferential grid, and Case 3 and Case 4 refine the flow grid based on Case 2. In addition, Case 3 also adjusts the flow grid layout to ensure a high-resolution mesh in the transition region. Case 4 refines the grid in the wall normal direction, and the space of the first grid (to the wall) keeps the same as that in other cases. The table shows that for Case 4, the total number of grids has reached 1.4 billion grids, the streamwise grid space is 10.0 wall units (y + ), and the spanwise grid space has 4.1 wall units, reaching the grid resolution standard of DNS. Figure 5 shows the distribution of the skin friction coefficients at the center-line of the bottom surface for the four cases. It can be seen from this that near x = 1000 mm, the friction coefficient rises sharply, indicating the occurrence of transition. The friction  showing the grid convergence. The following analysis is based on the data of Case 4. Figure 6 shows the distribution of skin friction coefficients on the surface of the lifting body. It can be seen from the figure that most of the upper surface keeps laminar, and the transition area only occurs in the concave area near the bottom. On the lower surface, transition occurs in a large area. Figures 7 and 8 are respectively partial enlarged views of the distribution of friction coefficients in the transition regions of the upper and lower surfaces. It can be seen from Fig. 8 that there is a significant transition in the recesses on the upper surface, and streamwise streak structures are formed, which shows that the transition occurs here. In the area close to the side at the bottom, obvious spanwise stripes are formed and in the following study, we find this region is the fully developed turbulence. Then the spanwise stripes appeared twisted and locally deformed, indicating the occurrence of secondary instability. Figure 8 is a partial enlarged view of the skin friction coefficient of the lower surface. It can be seen that the transition appears near the centerline of the lower surface. The transition appears firstly in the form of turbulent spots, and then the turbulent spot area expands and forms streamwise streak-like structures, indicating that the flow has reached turbulence. The figure shows that the disturbance waves converge to the central line from the left and the right side. In the area where the waves intersect, a complex stripe-like perturbation structure appears, and then an obvious turbulence spot structure is formed and developed to turbulence. Figure 9 shows the wall heat flux on the surface of the lifting body. It can be seen that the distributions of the heat flux and the skin friction have very strong similarity. Figure   Fig Figure 10 shows the distribution of instantaneous temperature at different streamwise positions. Along the streamwise direction, the structures become chaotic in four regions and it can be claimed that there are four transition regions. Figure 11 shows the temperature profile at x = 1450 mm. From the figure we can see that there are four transition regions marked as zone 1, zone 2, zone 3 and zone 4. Figure 11 shows that To further study the transition routines of the four zones, the vortical structures and instantaneous pressure fluctuations are investigated. Figure 12 shows the streamwise vortex at x = 1200 mm. We can find the streamwise vortex is formed in zone 1 and zone 4. Figure 13 shows the velocity vector illustration of zone 1 and 4 at x = 1400 mm. From Figs. 12 and 13, we can find the shear stress near the wall. So, we infer that the shear stresses destabilize the vortex in the streamwise direction and then the transition occurs. We also give the instantaneous pressure fluctuations of zone 4 in a streamwise wall-normal plane, shown in Fig. 14, which shows the streamwise fast increasing process of the secondary instability waves. Different from zones 1 and 4, there are no obvious separation bubbles in zones 2 and 3, and from the mean surface streamlines (see Fig. 9), it is deduced that crossflow instability plays important roles in zones 2 and 3, and it needs more study like LST to analyze the disturbance to give more information in the future.

Frequency spectrum and perturbation wave analysis
To further study the mechanism of the transition, the frequency spectrum and perturbation wave analysis is performed. Figure 14 shows several measuring points on the lower surface. Figure 14(a), Fig. 14(b), Fig. 14(c) and Fig. 14(d) give different settings of the measuring points for zone 1, zone 2, zone 3 and zone 4.  Figure 16(a) shows that the perturbation contains two frequency peak regions, the low-frequency peak is about 2KHz and the high-frequency one is about 30KHz5 3KHz. The authors conjecture that the low-frequency relates to the oscillation of the large-scale separation bubble (see Fig. 10b), and the high-frequency relates to the secondary instability vortex (see Fig. 15). Figure 16(b) shows a wide frequency band with a peak frequency at 49KHz, which indicates that the secondary vortex is stronger enough to excite the nonlinear perturbations (see Fig. 10c). There is a wide frequency in Fig. 16(c) and it shows that the transition occurs. In Fig. 17, the pressure spectrum exhibits exponent behavior with − 5/3 law at x = 1450 mm of zone 1, which indicates that the flow at x = 1450 mm of zone 1 is the fully developed turbulence. Figure 18 shows spectra of pressure fluctuations of zone 2. In the early range, the spectrum of pressure is typical double-peak spectrum, which contains a low-frequency  band (20-50KHz) and a high-frequency band (150-250KHz), and then, the spectrum fuses into wide spectrum. Since there are strong crossflows in zone 2 (see Fig. 9), it is conjectured that the low-frequency band is related with the instability of crossflow and the high-frequency band is related with the secondary instability or the Mack 2nd modes, because the typical frequency range of the Mack 2nd modes is similar to the high-frequency band (150-250KHz). In order to clarify the mechanism of instable waves, further research is necessary, especially combined with a stability theory. Figure 19 shows spectra of pressure fluctuation of zone 3. The transition mechanism in zone 3 is somehow like that in zone 2, where crossflow instability plays an important role. Like zone 2, the spectrum is also double-peak in the early ranges (see Fig. 19b, c), which contains a high-frequency band with the central frequency of 200KHz and a low-frequency band with the central frequency of 20KHz. Since the flow in this zone is like that in zone 2, the low-frequency band is conjectured as the effect of crossflow instability, and the high-frequency band is conjectured as the effect of the secondary instability or the Mack 2nd waves. Figure 20 shows spectra of pressure fluctuations of zone 4. The transition mechanism in zone 4 is similar to that in zone 1, but the transition in this zone is much later than that in zone 1. Figure 20b and Fig. 20c also show a double-peak frequency feature. Compared with zone 1, the frequency of the low-frequency band is lower, which is about 1KHz, indicating the oscillations of the separation bubble in zone 4 (see Fig. 10af) are at a lower frequency.

Proper orthogonal decomposition (POD) analysis
To further study the transition mechanism of zone 4, we perform the POD analysis. Collect a set of transient information for a certain flow field variable, {u 1 , u 2 , ⋯, u N }, and it can be described as where u i is instantaneous variables in i-th moment, and v i is the pulsating quantity.  The POD method expresses v i by a linear combination of a set of optimal orthogonal basis functions, and where p j is the POD modal basis function, a j (t i ) is the modal coefficient at time t i corresponding to the modal p j . We define matrix C = V T V, where V = {v 1 , v 2 , ⋯, v N }, and then the eigenvalue can be solved as where λ j is the eigenvalue, A j = [a j (t 1 ), a j (t 2 ), ⋯, a j (t N )] T is the corresponding eigenvector matrix. Rearrange the eigenvalues according to their size, then the first-order POD mode corresponds to the maximum eigenvalue, and the rest of the modes can be deduced by analogy. Further, the POD mode can be solved as POD can use energy to measure the contribution of each mode to the flow, and the energy is defined as We perform POD decomposition on the temperature distribution of zone 4. Figure 21 gives the percentage of the first eight POD modes, which shows that the amplitude of the first mode is much stronger than other modes. The first mode is related with the mean flow, and since we are more interested in the perturbation waves, we pay more attention to the other modes. Figure 22 shows the images of POD modes 2-6, and these are the second inability modes. This figure indicates that mode 2 is approximately symmetry, and mode 3 is approximately anti-symmetry, which shows that both symmetry and anti-symmetry modes play important roles in the transition. Figure 22 also shows that the strongest perturbations occur in the head and the shoulder of the mean separation bubbles, where strong shear occurs.

Conclusions
In this paper, direct numerical simulation (DNS) of hypersonic boundary layer transition of a lifting body model HyTRV is performed, and the perturbations are the small random blowing-and-suctions in the lifting body's head. The grid convergence tests are performed to validate the DNS. The analysis of DNS data came to the following conclusions: 1) The calculation results show that the lower surface of the lifting body has a large range of transition. While the upper surface only has a transition in the concave area of the bottom wall. 2) The DNS shows four transition zones on the surface of the lifting body, and the transition mechanism is different for each zone. Through near-wall streamline analysis, it is found that the transition is easier to occur in the region where streamlines converge.
3) The POD analysis shows both symmetry and anti-symmetry modes in zone 4 (the concave region of the lower surface). 4) The secondary instability of the separation bubble (streamwise vortex) plays an important role in zones 1 and 4, and the crossflow instability plays an important role in zones 2 and 3.