A CFD-Compatible Amplification Factor Transport Equation for Oblique Tollmien-Schlichting Waves in Supersonic Boundary Layers

Boundary layer transition is a hot research topic in fluid mechanics and aerospace engineering. In low-speed flows, twodimensional Tollmien-Schlichting (T-S) waves always dominate the flow instability, which has been modeled by Coder and Maughmer from 2013. However, in supersonic flows, three-dimensional oblique Tollmien-Schlichting waves become dominant in flow instability. Inspired by Coder and Maughmer’s NTS amplification factor transport equation for twodimensional Tollmien-Schlichting waves in low-speed flows and Kroo and Sturdza’s linear stability theory (LST) analysis results for oblique Tollmien-Schlichting waves in supersonic flows, a new amplification factor transport equation for oblique Tollmien-Schlichting waves has been developed based on LST. The compressible Falkner-Skan similarity equations are introduced to build the relationships between nonlocal variables and local variables so that all the variables used in the present model can be calculated using local variables. Applications of this new transport equation to the flows over supersonic flat plate, 3% thick biconvex airfoil, and one modified supersonic laminar airfoil show promising results compared with the standard LST analysis results.


Introduction
Since laminar flow has less drag than turbulent flow, laminar flow design technology has been a research hotspot in energy conservation of the green aviation [1]. In the process of laminar flow design, the accuracy of transition prediction plays a crucial role on the design effect. Therefore, it is very important and meaningful for aircraft designers to pay close attention to smart and efficient transition prediction methods. In recent years, there are two main routes to predict transition for airplanes and other complex aerodynamic configurations. One is the local transition models established by experimental data and stability analysis results, such as Menter et al.'s γ − Re θt correlation-based transition model [2][3][4], Walters et al.'s k-k L -ω model based on laminar kinetic energy mechanism [5,6], Fu and Wang's k-ω-γ transition model for high-speed flows [7], and Xu et al.'s physical mode-based transition models [8][9][10][11]. These models play an important role to predict transition for three-dimensional complex aerodynamic flows. The advantages of these transition models are local, convenient, efficient, and compatible with CFD parallel computations. The disadvantage is relying too much on experimental data and empirical parameters.
Compared with the local transition models mentioned above, in the 1950s, a semiempirical method named e N based on linear stability theory, proposed by Smith and Gamberoni [12] and Van Ingen [13], is widely used to predict transition in industry aerodynamic applications [14]. This LSTbased method has been chosen for transition prediction in subsonic and transonic boundary layers by Boeing Inc., Airbus Inc., German Aerospace Center (DLR), National Aeronautics and Space Administration (NASA), France Aerospace Center (ONERA), etc. Subsequently, the e N method was simplified by Drela and Giles [15] based on LST analysis results of similarity solutions for laminar flow. This simplified method built the functions between the most unstable amplification factor and the streamwise boundary layer shape factor H 12 , which has been implemented into the famous airfoil design soft "X-foil" [16]. With the development of CFD technique, Krimmelbein and Krumbein [17], Bégou et al. [18], Pascal et al. [19], and Shi et al. [20] coupled the LST-based e N method with Reynolds-averaged Navier-Stokes (RANS) code. This kind of coupling is reliable, but it is still complex because it also needs to solve boundary layer equations and linear stability theory equations, to integrate the eigenvalues, and to search nonlocal flow variables at the edge of boundary layers.
Based on Drela's idea, Coder and Maughmer [21] established an amplification factor transport equation to solve the amplification factor based on the approximate envelope method, which was extended using new local pressure gradient parameters [22,23] recently. This transition model can predict the two-dimensional Tollmien-Schlichting (T-S) instabilities and laminar separation bubble-(LSB-) induced transition in low-speed flows. It is worth mentioning that this transport equation for N TS factor combines the advantages of e N method and local transition models. All the variables in Coder and Maughmer's transport equation can be calculated using local flow variables so that it can be compatible with modern CFD codes conveniently, especially for unstructured codes.
In 2016, Xu et al. [24] constructed a transport equation for amplification factor of crossflow waves, which is restricted to winglike geometries. In 2019, Xu et al. [25] established a local amplification factor transport equation for crossflow instability in low-speed boundary layers, which performs well in several classical transition prediction cases. Hence, it is time to develop this modeling idea to high-speed flows. As known, in subsonic and low-transonic flows, two-dimensional T-S waves dominate the T-S instabilities. However, in supersonic boundary layers, oblique T-S waves play a dominant role [25]. In this paper, we are trying to build a brand new amplification factor transport equation for the oblique T-S waves in supersonic flows. Because the instability mechanism of two-dimensional T-S waves and oblique T-S waves is different, the present transport equation only has the similar form but different content compared with Coder's transport equation. Note that all the nonlocal variables are fitted using the solution database of compressible Falkner-Skan similarity equations. Since a suitable critical value of amplification factor can be found in the flows below Mach number 3.0 using Mack's relations [14,26,27] with freestream turbulence intensity, the present work is very valuable and meaningful for natural laminar flow (NLF) optimizations of supersonic airfoils and wings.

Modeling of the Transport Equation
With Illingworth transformation, the two-dimensional boundary equations can be written as [28] ρμ subject to the boundary conditions In the equations above, f ′ = u/U e and g = T/T e indicate the velocity profile and temperature profile, respectively. β H = ð2ξ/U e ÞðdU e /dξÞ is the Falkner-Skan pressure gradient parameter, Pr is the Prandtl number, γ H is the ratio of specific heats, M is the Mach number, ρ is the density, and μ is the dynamic viscosity. Note that the subscript "e" means the variables at the edge of boundary layer.

Transport Equation Description for Oblique T-S Waves.
Firstly, the transport equation takes the form where S is the strain rate magnitude, F crit indicates the onset function, F growth stands for the development function of growth rate, and dN TS /d Re θ represents the slop of the amplification factor N TS and momentum thickness Reynolds number Re θ . The source term is mainly established based on the extensive linear stability analysis results by Kroo and Sturdza [29] for oblique T-S waves in supersonic flows. The details can be found in Ref. [29]. Secondly, the onset function, F crit , is given by Here, H k , proposed by Drela and Giles [15] and defined as H k = Ð ð1 − u/U e Þdy/ Ð ð1 − u/U e Þu/U e dy, means the kinematic shape parameter. Thirdly, the function F growth is similar to Drela and Giles's [15] and Coder and Maughmer's [21] formulations: where The functions lðH k Þ and mðH k Þ have the same expressions with Drela and Giles's formulations. The DðH k Þ correlation is developed to modify the behavior of the source term at various kinematic shape parameters through lots of calibrations.
Fourthly, the slope function is modeled as where K 1 , K a , K b , K c , and K d are all the compressibility correction coefficients and T w is the wall temperature. It is worth mentioning that K m means the history effect. Kroo and Sturdza use K m = 11: H k ds to calculate K m , which is very difficult to compute using local variables. For the integral, transport equation with additional source term can be applied. However, it seems very difficult to compute the average of upstream parameter H k using local variables. Therefore, in this paper, the history effect term is set as zero temporarily and this term will be developed in the next step.
Not only two-dimensional T-S waves but also threedimensional oblique T-S waves are stabilized by favorable pressure gradient and increase near the adverse pressure gradient region. It should be mentioned that the H k ranges from 2.5 to 3.0 for these formulations above, which means this transport equation can only be used to describe the development of pure oblique T-S waves in the supersonic flows with moderate favorable pressure gradient and adverse pressure gradient. For the flows with strong favorable pressure gradient, like the stagnation point flows around blunt leading edge, subsonic regions always appear so that H k exceeds the current modeling scope. The present model cannot predict the two-dimensional T-S waves which are different from the oblique T-S waves.

International Journal of Aerospace Engineering
Fifthly, after filtering the existing parameters to describe the pressure gradient, the local parameter H L = Sy/U e , proposed by Coder and Maughmer [21], was adopted to calculate the kinematic shape parameter H k . In the definition, y is the distance to the nearest wall. Figure 1 Sixthly, Xu et al.'s engineering estimation method [9] for Mach number at the edge of boundary layer is chosen, which has the expressions as follows: where a stands for the sound speed, P is the pressure, and the subscript ∞ denotes variables in freestream. Furthermore, ρ e could be given by ρ e = ½ðρ γ H ∞ /P ∞ ÞP 1/γ H and the Mach number at the edge of boundary layer is determined as M e = U e /a e . In some nonlocal models, M e can be obtained through searching operations. However, through the aerodynamic equations (11) and (12), it can get a relatively accurate estimation of M e , which is much more accurate than using the freestream Mach number M ∞ directly. It should point out that the prediction of M e is a CFD issue, which can be computed by the aerodynamic formulations and CFD variables. This part does not belong to the transition modeling content. Finally, the last unknown parameter is the local momentum thickness Reynolds number Re θ = ρU e θ/μ. From the similarity solution database, the relationship between Re θ and the vorticity Reynolds number Re V = ρSy 2 /μ is described as Since there is a small height difference between the maximum value of Re V and H L in favorable pressure gradient boundary layers, a minor height correction is introduced: Re θ = Re θ 1:086 exp 0:3455λ θ ′ + 0:01279 exp 18: where λ θ = ðρθ 2 /μÞðdU e /dsÞ is the Thwaites pressure gradient factor and λ θ ′ can be computed using the following equations [25,30]:

International Journal of Aerospace Engineering
Note that the incompressibility correction for λ θ comes from Ref. [31] and λ θ ′ corresponds to the modified Thwaites pressure gradient factor.
Finally, the effective intermittency factor takes the form where γ is the intermittency factor and the effective intermittency factor γ eff is used to trigger transition in the turbulence model. The coupling way between the transition model and the Menter's Shear Stress Transport (SST) turbulence model from Coder and Maughmer's paper [32] in 2013 is selected in this paper, which is given by Note that γ eff is the switch function for the generation of turbulent kinetic energy by the production term P k,original and the destruction term D k,original in the SST turbulence model. For comparison, the saddle point method, proposed by Cebeci and Stewartson [33], is adopted for the standard LST analysis.

Results and Discussion
In the present work, the open-source structured Reynoldsaveraged Navier-Stokes solver named CFL3D was used as the basic flow solver. Details of this solver could be found in NASA's website (data available online at https://cfl3d.larc .nasa.gov/) and documents [34]. In this work, all of the transition prediction results were obtained by the present transport equation coupling with Menter's k-ω SST turbulence model.  International Journal of Aerospace Engineering employed to conduct the mesh sensitive test, and the results are illustrated in Figure 3(a). The Mach number is 2.2 and Reynolds number is 5:6 × 10 7 /m. The near-wall mesh is fine enough so that y + ð1Þ of the cell next to the wall is smaller than 1.0. It can be seen that as the amount of the grid increases, the results converge gradually, demonstrating the robustness of the new transport equation. Figure 3(b) displays the contour of N TS factor calculated by the present transport equation. Note that the predicted maximum value of N TS factor is extracted from the N TS contour at each streamwise location. With further validations, various freestream conditions were considered. Figure 4 displays the contour of predicted N TS factor and the comparison with the standard LST analysis results under the conditions of M ∞ = 1:5 and Re = 5:4 × 10 7 /m. Meanwhile, the Mach number is 2.2 in Figure 5 and 3.0 in Figure 6. Moreover, the unit Reynolds number in Figure 5 is 3:0 × 10 7 /m and 5:6 × 10 7 /m in Figure 6. Although the Mach numbers and Reynolds num-bers change, the predicted results of N TS factor seem robust and accurate. When a critical value of N TS factor was set as 10.0, the results in the case with M ∞ = 3:5 and Re = 5:89 × 10 7 /m are plotted in Figure 7. Before the threshold value of N TS factor, the development of N TS factor is described accurate and transition occurs near the threshold. Subsequently, in the turbulent region, N TS factor will be quickly dissipated. As a result, the whole process of oblique T-S wave-induced transition is simulated reasonably.

Validation Test Case 2: Biconvex Laminar
Airfoil. Secondly, the 3% thick biconvex airfoil, i.e., the "Airfoil-1" as plotted in Figure Figure 9, compared with the contour of local Mach number. The predicted value of M e agrees well with the contour of local Mach number. The calculated N TS contour and the maximum value at each streamwise location are sketched in Figure 10.
When the angle of attack is 2 degrees, Figure 11 demonstrates the predicted Mach number at the edge of boundary layer which is in accord with the value at the outer edge of the boundary layer in the local Mach number contour. The predicted N TS contour is displayed in Figure 12. Furthermore, the extracted data on the upper surface and lower surface are plotted in Figures 13(a) and 13(b), respectively. It can be seen that most of the N TS factor are simulated accurately.
It is worth noting that the defects without history effect term gradually appear downstream. This is the main reason of the deviations between the predicted value and the standard LST results in the downstream region. Certainly, even though the nonlocal methods are used, there are still deviations compared with the standard LST analysis data [18]. Consequently, the small deviations predicted here can be accepted.

Validation Test Case 3: Modified Laminar
Airfoil. The third case for validation is the "Airfoil-2" as shown in Figure 8, which has a larger radius of curvature near the leading edge. The freestream Mach number is 1.8 and the unit Reynolds numbers contain Re = 3:5 × 10 7 /m, 5:4 × 10 7 /m, and 8:0 × 10 7 /m. The estimated M e is still in good agreement with the reference data  Figure 14. In addition, Figures 15-17 illustrate the computed N TS factor contour and the comparison with standard LST analysis results. Overall, the present transport equation works well. However, due to the lack of history effect term, there are obvious deviations between the predicted value of N TS and the standard LST results. Therefore, for the computations in moderate or strong favorable pressure gradient flows, it seems quite difficult to compute the average of upstream integrated variables using local variables. This problem will be solved in the future research.

Conclusions
In summary, a new amplification factor transport equation for the oblique T-S waves in supersonic flows has been developed based on linear stability analysis results and validated in several typical two-dimensional supersonic flows. The good agreement between the present transport equation and the standard LST analysis results in two-dimensional supersonic airfoils which show that the present transport equation is very promising and encouraging. Obviously, this model has the potential to be extended and applied to predict the oblique T-S waves on three-dimensional wings with sharp leading edge and cones with sharp nose. Moreover, crossflow instability, playing an important role in 3D boundary layers, should be taken into account in the future research. It should be pointed out that the present model can only simulate the pure oblique T-S waves, which means it cannot predict the airfoils with blunt leading edge. Because there is a subsonic region near the blunt leading edge, the instability mechanism near the leading edge may start from two-dimensional T-S waves to oblique T-S waves. Hence, this complex process cannot be captured using the present model.    On the whole, although the application scope is limited to supersonic configurations with sharp leading edge, all the formulations established in this model provide a good idea and a reasonable fundamental framework for the modeling of oblique T-S waves in supersonic boundary layers. The extensive validations at various Mach numbers, Reynolds numbers, and angles of attack are all in accord with the standard LST analysis results, which has proven that the present transport equation is constructed reasonably and correctly. As for the history effect term, it should be taken into account by a new transport equation in the next step.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.