Numerical Solution of Transition to Turbulence over Compressible Ramp at Hypersonic Velocity

: This work deals with the numerical solution of hypersonic ﬂow of viscous ﬂuid over a compressible ramp. The solved case involves very important and complicated phenomena such as the interaction of the shock wave with the boundary layer or the transition from a laminar to a turbulent state. This type of problem is very important as it is often found on re-entry vehicles, engine intakes, system and sub-system junctions, etc. Turbulent ﬂow is modeled by the system of averaged Navier–Stokes equations, which is completed by the explicit algebraic model of Reynolds stresses (EARSM model) and further enhanced by the algebraic model of bypass transition. The numerical solution is obtained by the ﬁnite volume method based on the rotated-hybrid Riemann solver and explicit multistage Runge–Kutta method. The numerical solution is then compared with the results of a direct numerical simulation.


Introduction
The design of hypersonic vehicles is becoming increasingly common nowadays, and it is therefore necessary to develop and use appropriate engineering tools.
Hypersonic flow has a very complex nature. Strong shock waves may be present in the flow field as well as in their interactions with each other or with the boundary layer. Contact waves, large areas of a subsonic flow field, supersonic jets, and other complicated elements may also be present in the flow field. Measuring these flow fields in experimental aerodynamic laboratories can be technically difficult and expensive. Therefore, computational fluid dynamics (CFD) methods are increasingly being used.
Solving hypersonic flow using CFD methods is still a major challenge. It is necessary to use a suitable numerical method that is sufficiently robust, stable, and capable of approximating the flow at low and high Mach numbers. The numerical scheme must also be able to accurately capture all the flow field elements described above. Recent developments in numerical methods, e.g., [1][2][3][4][5], address some of these problems.
Another issue is the appropriate modeling of turbulence, which is still an open problem. Many commonly used models of turbulence neglect terms that may only be significant at hypersonic velocities [6]. A significant part of the problem is modeling the transition from a laminar to a turbulent flow regime. The answer to all these problems may be the use of the direct numerical simulation (DNS), which can accurately predict hypersonic turbulent flows, e.g., [7][8][9][10]. However, the DNS method requires a very precise numerical method and a very fine computational mesh (the number of cells required for the DNS method is equal to Re 9 4 [6]). Therefore, its use for industrial applications in current computers is not possible.
The main objective of this work is to test the possibility of using a relatively common CFD method based on averaged Navier-Stokes equations closed by the modified EARSM model of turbulence and second-order accurate numerical method for hypersonic flow prediction.
The problem selected for this test is flow over a compressible ramp. Although the geometry of the problem is simple, the solution is complex and involves the interaction of shock waves, shock wave interaction with the boundary layer, and separation, reattachment, and even transition from a laminar to a turbulent state. This problem is the subject of current research, e.g., [7,[11][12][13]. These publications use the DNS method, which is suitable for research purposes but very expensive for use in design calculations of re-entry vehicles or engine intakes, which are typical applications of solved problems.
The importance of this test is to determine whether common CFD methods can be used for such cases with sufficient accuracy for engineering designs.
The paper is organized as follows. First, the averaged Navier-Stokes equations together with the modified EARSM model of turbulence are introduced in Section 2. The following section deals with the numerical solution using the finite volume method. After that, the problem of flow over a compressible ramp is formulated and solved numerically. Finally, the obtained results are discussed in the last section.

Governing Equations
Turbulent flow of viscous compressible fluid is modeled by the system of averaged Navier-Stokes equations in vector form (summation over repeated indices will be used throughout the paper), ∂W ∂t where W = (ρ, ρu i , ρE) T is the vector of mean conservative variables (Reynolds averaging is used for density and static pressure, and Favre averaging is used for total specific energy and for components of the velocity vector [14]), ρ is density, u i are components of the velocity vector, E is total specific energy, p is static pressure, and δ ij is the Kronecker delta. The components of the mean stress tensor τ ij are given by relation where µ is the mean viscosity computed using the Rayleigh relation, with µ re f , ρ re f , and p re f being constants. The components of the mean heat flux vector q j are defined as where κ is the specific heat ratio (in this work, the flow of perfect gas (air) is assumed, for which κ = 1.4 holds), and the Prandtl number Pr is assumed to be constant (in this work, Pr = 0.72). The components of the turbulent heat flux vector q t j are analogically ([6]) defined as where the turbulent Prandtl number Pr t is also assumed to be constant (in this work, Pr t = 0.9). The term d t corresponds to the molecular diffusion and turbulent transport of the turbulent kinetic energy k. This term is usually neglected but can be significant in hypersonic flows [6]. The components d t j are given by the relation The model constant σ * , the turbulent viscosity µ t , and components of the Reynolds stress tensor τ t ij will be given in the next subsection. The system of averaged Navier-Stokes Equations (1) is also equipped with the equation of state of a perfect gas in the form Model of Turbulence The system of averaged Navier-Stokes Equations (1) is not closed and therefore needs to be completed by some suitable model of turbulence.
The model of turbulence used in this work is based on the EARSM model, which is further enhanced by the algebraic model of bypass transition [14].
The model consists of the transport equation for the turbulent kinetic energy k, and transport equation for the specific dissipation rate ω, where is the original production term of the turbulent kinetic energy [15], C d represents crossdiffusion given by the formula and β * = 0.09, σ * = 1.01, α = 0.553, β = 0.075, σ = 0.5, and σ d = 0.52 are model constants.
The algebraic model of bypass transition is represented by two functions, namely f SS and γ, where the first is sheer sheltering factor where Ω = 2Ω ij Ω ij is the vorticity tensor magnitude, and constant C SS = 2.75. The latter function represents an intermittency factor defined as where C T = 1.8125 is constant, and the Reynolds number Re kΩ , which serves as the transition onset parameter, is given by the relation Finally, the last two terms in Equations (8) and (9) (τ t ij and µ t ) are the components of the Reynolds stress tensor and turbulent viscosity, respectively. Both terms are parts of the EARSM model's constitutive relations. The components of the Reynolds stress tensor are given as where turbulent viscosity is given by the relation and extra anisotropy is given as The turbulent viscosity µ t and extra anisotropy a (ex) ij are nonlinear terms, which depend on the turbulent time scale dimensionless strain-rate and vorticity tensors and on the related invariants II Ω = Ω * kl Ω * lk and IV = S * kl Ω * lm Ω * mk . The beta coefficients β 1 -β 9 are defined as where the denominator Q is with N representing the turbulent kinetic energy production and dissipation rate ratio defined as where and C 1 is given by relation where and

Numerical Method
The system of averaged Navier-Stokes Equations (1), together with the transport equations of the EARSM model, can be written as

∂W ∂t
where W = (ρ, ρu i , ρE, ρk, ρω) T is the vector of the unknown conservative variables, F j are corresponding inviscid fluxes, R j are corresponding viscous fluxes, and vector Q contains source terms. The numerical solution is obtained by in-house software based on the cell-centered finite volume method [16] in semi-discrete form where W i is an averaged solution over the cell D i , viscous flux R = R · n, inviscid flux F = F · n, ∆S is the area of the cell face between two adjacent cells, and n is the unit vector normal to a face. Since hypersonic flow is considered in this work, a rotated-hybrid Riemann solver [1] is used for approximation of inviscid flux. By using the rotated-hybrid Riemann solver concept, a carbuncle phenomenon caused by strong shock waves can be avoided. Moreover, the rotated-hybrid Riemann solver is based on a combination of the low Mach number variants of the HLLC and the HLL schemes, respectively, which is important for accurate prediction of subsonic parts of the flow field.
Second-order accuracy is accomplished by piece-wise linear WENO reconstruction [17] of the primitive variables. The viscous flux is discretized by the central scheme with aid of diamond-shape dual cells [18]. Time integration of the semi-discrete finite volume method (30) is done by the explicit two-stage TVD Runge-Kutta method [19] with point implicit treatment of the source terms.

Numerical Solution of Hypersonic Flow over a Compressible Ramp
The hypersonic flow over a compressible ramp is characterized by the free-stream Mach number M ∞ = 7.7 and the Reynolds number Re = 8.6 · 10 5 . This setup corresponds to the same problem with total enthalpy h 0 = 1.7 MJ·kg −1 , as in the paper [7]. The relatively low total enthalpy allows a calorically perfect gas to be used as the fluid.
The simulation was performed in the computational domain shown in Figure 1. The lower boundary of the computational domain consists of a flat plate of length L = 100 mm, followed by a ramp of length 220 mm, with an angle of 15 degrees. The left boundary (x = 0) of the computational domain and the upper boundary, which is equidistant by 30 mm from the lower boundary, are considered as inlet, while the right boundary is considered as outlet. The computational domain was covered by a hyperbolically generated grid with 400 × 100 cells, with minimum grid spacing in the y-direction ∆y = 3.8 · 10 −6 m, which corresponds to y + ≈ 0.9.
The free-stream turbulent kinetic energy k ∞ was calculated using relation where the free-stream turbulence intensity Tu is equal to half a percent, which is approximately the lower limit of the bypass transition model used in this work. Therefore, we are in the region between the natural transition and the bypass transition. The free-stream value of the specific dissipation rate ω ∞ is based on the dimensional analysis where L is the characteristic size (in this case the length of the flat plate in front of the ramp), and C is constant of order O(1). In our case, the value C = 7 is chosen because it allows us to correctly capture the transition onset. Finally, the value of the specific dissipation rate ω wall is calculated using relation which is based on the solution of the simplified Equation (9) in the near-wall region. Here, the constant C ω = 1.5, and ∆y 1/2 is the distance between the center of the first cell and the wall. Figure 2 shows the Mach number distribution. It shows a relatively weak shock wave generated by the leading edge of the flat plate and a separation bubble around the corner, with boundary layer separation located at x/L = 0.67 and reattachment at x/L = 1.17. Both the separation and reattachment of the boundary layer generate stronger shock waves, which then interact at the so-called triple point located at x/L = 1.37 (Figure 3). The structure of the flow field is in good qualitative agreement with results obtained by the DNS method in [7].   Table 1 shows a comparison with the results obtained by the DNS method in [7]. There is a relatively large difference in the predicted location of the boundary layer separation compared to the DNS method. Since it is in the laminar part of the boundary layer, this difference is almost certainly due to the much more accurate numerical method and to the much finer grid used by the DNS method. On the other hand, the reattachment location is predicted reasonably well by the presented method. The difference in the predicted location of the triple point is due to above-mentioned shift in the separation location downstream. The shift of the triple point is caused by the change in the angle of the separation shock wave. Finally, the transition onset location is predicted, in good agreement with the DNS method.  Figure 4 shows the turbulent kinetic energy distribution. It shows a rather sharp transition onset located at x/L = 1.67.  Figure 5 shows distributions of the pressure coefficient. A small pressure increase can be observed due to boundary layer separation and to a much higher increase downstream of the reattachment location. The first increase is predicted later compared to both the DNS method and the experiment. This is again due to the much higher accuracy of the DNS method. The second increase is predicted more accurately. A qualitatively different shape can be observed around the point x/L = 1.6. This is due to the typical sharp bypass transition from the laminar to turbulent flow regime, which is notable even in the pressure distribution. The location of the pressure drop in the vicinity of x/L = 2 is predicted in good agreement with the DNS method, although the magnitude of the drop is smaller in the case of the presented method. Finally, the magnitude of the pressure coefficient in the turbulent part of the boundary layer is predicted in very good agreement with the DNS method.
A grid convergence study was also performed. The same problem was solved for three progressively refined meshes. The coarse mesh consists of 200 × 100 cells, the medium mesh of 400 × 100 cells, and fine mesh of 600 × 100 cells. Figure 5 shows all three results, which are almost indistinguishable. Therefore, the results can be considered to be grid-convergent.

Discussion and Conclusions
A numerical simulations of the hypersonic flow of viscous fluid over a compressible ramp with transition from laminar to turbulent flow regime has been performed.
The numerical solution was obtained by the common and much less demanding numerical method used on relatively coarse meshes (40-60 thousand cells) in comparison with the direct numerical simulation (627 million cells). This allows us to obtain solutions in a much shorter time, which is required for efficient analysis and industrial designs. This is clearly the biggest advantage of the presented method. On the other hand, the DNS method provides more accurate results that agree very well with the experiment, but at the cost of very long computation times (it can take months or even years, or it may be nearly impossible on current computers for cases with high Reynolds numbers). This makes the DNS method more suitable for research purposes or as a substitute for experiments.
Another advantage is the relative simplicity of the presented method, which can be used to solve more complex problems with complicated geometry, whereas the DNS method is generally only used for problems with simple geometry. Additionally, since many existing codes are based on the averaged Navier-Stokes equations and two-equation models of turbulence, it is easy to implement the presented method in existing solvers.
The lack of ability to model the natural transition is the biggest weakness of the presented method. The transition from a laminar to a turbulent flow regime is modeled by a very simple algebraic model of bypass transition, whereas the real transition process in the case of flow over a compressible ramp is probably caused by means of the natural one. This leads to a significant difference in the shape of the pressure coefficient distribution in the transition region.
Another weakness is the use of a numerical scheme with only second-order accuracy. There are visible differences (especially the location of the laminar boundary layer separation) in the solution due to the significantly lower accuracy of the presented method compared to the DNS. This could be resolved by using a higher-order reconstruction in space.
However, results obtained by the presented method are qualitatively correct and are reasonable accurate, at least for common industrial design applications. The computational time required to obtain a solution using the presented method is on the order of hours on commonly available computers. Therefore, it can be stated that the presented method has the potential to be used for design calculations of hypersonic vehicles.