Numerical Investigation on the Transition Flow around NLF Airfoil

Abstract

A natural laminar flow (NLF) airfoil is designed to reduce drag by expanding laminar flow areas. In-depth knowledge of transition performance is essential for its aerodynamic design. The k-ω-γ-Re_θ framework, which consists of the SST k-ω turbulence model and γ-Re_θ transition model, is employed to simulate transitional flows around an NLF wing RAE5243 airfoil. The transition performances of the RAE5243 airfoil under various values of turbulent intensity, temperature, angle of attack, and Mach number are simulated and compared. The results show that the rise of inflow turbulent intensity will promote an earlier transition on both the suction and pressure sides. The influence of wall temperature on transition is limited. The rise of angle of attack will lead to an earlier transition on the pressure side but a later transition on the suction side. With the rise of Mach number, the transition happens earlier under a zero and positive angle of attack but later under a negative angle of attack. In addition, the correlation of transition onset locations with respect to turbulent intensity, surface temperature, angle of attack, and Mach number is established based on numerical results.

1. Introduction

The laminar-turbulent transition is a typical flow phenomenon inside boundary layers, and it is related to flow instability, friction resistance, and heat transfer [1,2]. Because of the significant influence of transition on flow characteristics, it has become a hotspot and difficulty in aerodynamic design [3,4,5]. The natural laminar flow (NLF) and hybrid laminar flow (HFL) aerodynamic designs have attracted many scholars’ attention, since the flow friction, aerodynamic heat, and aerodynamic noise can be greatly reduced under a laminar state. These types of aerodynamic designs rely on knowledge of the transition, especially its onset, and therefore various studies have been carried out to predict the transition process.

The Reynolds averaged Navier–Stokes (RANS) method is the most commonly used computational fluid dynamics (CFD) approach in engineering practice [6,7], since it can seek a balance between simulation accuracy and calculation cost [8,9]. Therefore, most transition models are coupled with the RANS framework, and the transition model based on intermittency mode is the most popular type. The intermittency γ is defined by the occurrence probability of turbulent fluctuation, which is equal to 0 in a laminar boundary layer, and equal to 1 in fully turbulent flow. Dhawan established a prediction model of intermittency along the streamwise direction according to experimental data [10], which was an early attempt to estimate the transition process using intermittency. To improve the accuracy and generality of intermittency mode, the transport equations of intermittency have been constructed by different scholars [11,12]. The intermittency mode is effective in estimating the transition process from laminar to turbulent flow, but it cannot be employed as the criterion to predict transition onset directly. Other parameters are applied to serve as onset criteria, including laminar kinetic energy k_L and momentum thickness Reynolds number Re_θ. The concept of laminar kinetic energy k_L was first proposed by Bradshaw [13]. According to experimental data, it is considered that the energy in the transition region can be split into laminar kinetic energy and turbulent energy, and the laminar kinetic energy plays a vital role in amplifying perturbations and bursting transitions before the transition. The transport equation of laminar kinetic energy has also been constructed in different forms by various scholars [14,15]. The momentum thickness Reynolds number is another commonly used parameter as transition criterion. The relation between Re_θ and turbulent intensity is established empirically according to experiments with a flat plate under a varying pressure gradient [16,17]. However, nonlocal variables are required in the established empirical relation, and the transport equation of momentum thickness Reynolds number has also been constructed to solve its local distribution accurately [18]. So far, based on the combination of RANS framework, intermittency mode, laminar kinetic energy, and momentum thickness Reynolds number, different types of transition models have been developed and implemented into CFD codes, such as k-ω-γ mode [19,20], k-ε-γ-k_L mode [21], and k-ω-γ-Re_θ mode [7,18].

With the development of a transition model with sufficient accuracy and efficiency, the aerodynamic design of an NLF airfoil has been greatly improved using numerical techniques. The transition modes based on linear stability theory [22,23,24] and e^N [25,26] method are widely employed in the numerical simulations of transitional flow around an NLF airfoil. Since these methods always require non-local boundary-layer parameters, and these are difficult to implement in RANS codes, the k-ω-γ mode [27,28] and k-ω-γ-Re_θ [30,31,32] mode have become more widely used in RANS solvers to investigate the flow structure around an NLF airfoil. By means of these numerical techniques, various optimization algorithms are employed to improve the aerodynamic performance of an NLF airfoil, including the adjoint method [27,28], particles swarm optimization [29], and genetic algorithm [30,31]. However, despite the various studies on aerodynamic design, comprehensive investigations on the aerodynamic characteristics of NLF airfoils are relatively rare.

According to the best knowledge of the authors, the k-ω-γ-Re_θ transition model, as a typical representative of the transition model, has proven its effectiveness in a large number of numerical simulations, making it an important and efficient tool for investigating the transition characteristics of airfoils. However, systematic analysis of the influence of different operation conditions on the transition characteristics is relatively limited. It is of great significance to clarify the correlation between the transition position and operation conditions. To this end, in the present paper, the k-ω-γ-Re_θ transition model is employed to investigate the transitional flow around an NLF wing RAE5243 airfoil [32] under various conditions. The influences of turbulent intensity, temperature, angle of attack, and Mach number on transition location and aerodynamic characteristics are discussed and analyzed, which can provide a guide for the aerodynamic design of an NLF airfoil under wide operating conditions.

2. Numerical Method

2.1. Transition Model

The k-ω-γ-Re_θ transition model [7], which is coupled with the SST k-ω turbulence model [6], is applied in the present study. The governing equations of turbulent kinetic energy k, turbulent frequency ω, intermittency γ, and momentum thickness Reynolds number Re_θ are listed as follows:

$$\begin{array}{l}\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{u}_i\right)\\ =\frac{\partial }{\partial x_j}\left(\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right)+{\gamma }_{eff}{\mu }_t{S}^2-\min \left(\max \left({\gamma }_{eff},0.1\right),1.0\right)\rho \beta k\omega \end{array}$$

(1)

$$\begin{array}{l}\frac{\partial }{\partial t}\left(\rho \omega \right)+\frac{\partial }{\partial x_j}\left(\rho \omega u_j\right)\\ =\frac{\partial }{\partial x_j}\left(\left(\mu +\frac{{\mu }_t}{{\sigma }_{\omega }}\right)\frac{\partial \omega }{\partial x_j}\right)+{\alpha }_{\omega }\frac{\omega }{k}{\gamma }_{eff}{\mu }_t{S}^2-\rho \beta f_{\beta }{\omega }^2+2\left(1-F_1\right)\frac{\rho }{\omega {\sigma }_{\omega ,2}}\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}\end{array}$$

(2)

$$\begin{array}{l}\frac{\partial }{\partial t}\left(\rho \gamma \right)+\frac{\partial }{\partial x_j}\left(\rho u_j\gamma \right)\\ =\left(1-C_{e1}\gamma \right)C_{a1}{F}_{length}\rho S{\left(\gamma F_{onset}\right)}^{C_{\gamma 3}}+\left(1-C_{e2}\gamma \right)C_{a2}\rho \Omega \gamma F_{turb}+\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\gamma }}\right)\frac{\partial \gamma }{\partial x_j}\right]\end{array}$$

(3)

$$\begin{array}{l}\frac{\partial }{\partial t}\left(\rho {\mathrm{Re}}_{\theta }\right)+\frac{\partial }{\partial x_j}\left(\rho u_j{\mathrm{Re}}_{\theta }\right)\\ =c_{\theta }\frac{{\rho }^2{u}^2}{500\mu }\left({\mathrm{Re}}_{\theta t}-{\mathrm{Re}}_{\theta }\right)\left(1.0-F_{\theta }\right)+\frac{\partial }{\partial x_j}\left[{\sigma }_{\theta }\left(\mu +{\mu }_t\right)\frac{\partial {\mathrm{Re}}_{\theta }}{\partial x_j}\right]\end{array}$$

(4)

The details of each source term in the above governing equations can be checked in the previous literature [6,7], and the default values of model constants are employed in the present study. In the k-ω-γ-Re_θ transition model, the governing equation of intermittency γ is solved to predict the development from laminar flow to turbulent flow, and the governing equation of momentum thickness Reynolds number Re_θ is employed to evaluate the transition onset. The main advantage of the γ-Re_θ transition model lies in the fact that only local flow variables, rather than the integral of the boundary layer flow, are required during simulation, and therefore it can be easily coupled with other turbulence models.

2.2. Physical Model and Computational Domain

The NLF wing RAE5243 airfoil [32,33] with chord length C = 0.24 m is investigated in the present study. A three-dimensional computational domain with a size of 10C × 2C × C is established, as shown in Figure 1. The inlet of the computational domain is set at the location 3C away from the leading edge of the airfoil, and the outlet is set at the location 7C downstream the leading edge. The height of the computational domain is two times the chord length, and the spanwise length of the computational domain is set the same as the chord length. In the present numerical simulation, a three-dimensional compressible transitional flow around the airfoil under steady state is considered. The density of the fluid is subject to the ideal gas equation of state, and the Sutherland model is employed to determine fluid viscosity.

The commercial software ANSYS ICEM 19.0 is employed to generate the hexahedral structured mesh in the computational domain. To account for the mesh requirement for simulating the transitional flow, mesh refinement is employed around the airfoil surface. The maximum value of y+ on the airfoil surface is within 2, and most values of y+ are at the level of 1. The gird growth ratio near the surface is equal to 1.05. The mesh distribution around the airfoil surface is illustrated in Figure 2.

2.3. Numerical Settings

The commercial CFD software ANSYS Fluent 19.0 is employed to solve the steady transitional flow around the RAE5243 airfoil. The k-ω-γ-Re_θ turbulence and transition coupled model described in Section 2.1 is applied in the present simulation. The pressure-based solver is used to simulate these subsonic and transonic flows, and the SIMPLEC algorithm is used for pressure-velocity coupling. For the boundary conditions, a pressure far-field condition is given at the inlet, with total temperature of 303 K and total pressure of 1 atm, and a pressure outlet condition is given at the outlet. The free slip and adiabatic conditions are imposed on the top and bottom walls of the domain, and the periodic condition is applied on the left and right boundaries. The airfoil surface is treated as a no-slip and isothermal wall in the present simulation. Since steady states are considered in the present study, the simulations start from freestream values in the entire domain.

2.4. Independence Test of Mesh

Mesh quantity can significantly influence the simulation results. In order to validate the independence of the employed mesh, four sets of mesh are employed in the present study, and the simulation results of the lift coefficient c_L and drag coefficient c_D are compared. The mesh quantity and corresponding simulation results are listed in

Table 1. Independence test of mesh.

No.	Mesh Quantity	Streamwise Grids	Lift Coefficient	Drag Coefficient
1	2.39 × 106	282	0.2672	0.0053
2	3.07 × 106	376	0.2781	0.0051
3	3.75 × 106	470	0.2797	0.0051
4	4.43 × 106	564	0.2800	0.0050

.

During the independence test of the mesh, the mesh distribution along the spanwise direction and normal direction remains the same, while the grids along the streamwise direction gradually increase. Therefore, the mesh quantity rises from 2.39 × 10⁶ to 4.43 × 10⁶, with streamwise grids rising from 282 to 564. From Table 1, it can be seen that the lift coefficient and drag coefficient hardly vary when the mesh reaches 3.75 × 10⁶. To find a tradeoff between computation cost and simulation accuracy, Mesh 3, with 3.75 × 10⁶ elements, is employed for the following simulations.

2.5. Verification of Numerical Methods

The experimental results of an NLF(1)-0416 airfoil [34,35] are referred to in the present study to verify the accuracy of the numerical methods. The same size of computational domain and the same topology of mesh as the RAE5243 airfoil are established. The numerical methods mentioned in Section 2.3 are also applied for the NLF(1)-0416 airfoil. The freestream incoming flow condition for the NLF(1)-0416 airfoil includes the Mach number Ma = 0.1, the Reynolds number Re = 2.0 × 10⁶, the static temperature T_∞ = 215.91 K, and turbulent intensity of 0.03%. The experimental and numerical results of the pressure coefficient c_p = (p − p_∞)/(0.5ρu_∞²) under these conditions are illustrated in Figure 3, where p, ρ, and u denote pressure, density, and velocity, respectively, and subscript ∞ represent freestream.

From Figure 3, it can be seen that the experimental and numerical pressure coefficients along the streamwise direction agree well, and the sudden rise of pressure coefficient due to transition onset can be captured. Therefore, the accuracy of the numerical methods can be validated.

In addition, experimental measurement for the RAE5243 airfoil is also carried out by means of the technique of temperature-sensitive paint (TSP) [33], and the only difference lies in the surface of the airfoil being treated as an adiabatic wall rather than a constant-temperature wall. The airfoil is tested in a wind tunnel with channel walls. Therefore, it is reasonable to assume free-slip boundary conditions on the side walls to avoid the interaction between boundary layer and shock wave under no-slip boundary conditions. In addition, the experimental measurement is conducted within a time in the order of minutes, so the wall of the airfoil can be treated as adiabatic. Following these experimental settings, the boundary condition for the surface of the airfoil is revised into an adiabatic wall, and the other numerical methods remain the same as those mentioned in Section 2.3. The RAE5243 airfoil with three angles of attack (α = −3°, 0°, 3°) under a Mach number of 0.4 is tested in the present research. Figure 4 illustrates the distribution of temperature along the streamwise direction on the suction surface. According to Figure 4, the sudden variation of temperature due to transition onset can be predicted accurately, which verifies the accuracy of the employed numerical methods again. It can be seen from Figure 4 that the main difference between the experimental and numerical results appears at the leading and trailing edge of the airfoil. In the present experiment, the TSP technique is employed to measure the distribution of temperature on the airfoil surface. Although the insulation paint is coated on the airfoil surface to achieve adiabatic conditions, a certain heat flux still appears at the leading edge and trailing edge of the airfoil due to the large temperature gradient and thin thickness. Therefore, the distribution of temperature at both ends is more significantly affected by heat conduction effects and no longer meets the conditions of adiabatic walls, which is still the main shortcoming of the TSP technique. Despite this, the agreement on the transition onset location between the experimental and numerical results validates the reliability of the present numerical method.

3. Result and Discussion

Using the verified numerical methods, the aerodynamic performance of an RAE5243 airfoil, with a special emphasis on the transition onset location under various conditions, is investigated in the present research, including turbulent intensity, temperature, angle of attack, and Mach number. The angle of attack is defined by the angle between the chord line of the airfoil and the inflow direction. In the present study, we built models of airfoils under various angles of attack and generated the corresponding meshes. The operating conditions for the present study are listed in

Table 2. Operating conditions.

Turbulent Intensity ε/%	Surface Temperature Tw/K	Angle of Attack α/°	Mach Number Ma
0.1, 0.5, 1	293, 313	−6, −3, 0, 3, 6	0.4, 0.6, 0.7, 0.75

. These conditions are determined according to the range of conditions of interest and the range of parameters in which our experiments can be carried out. For turbulent intensity and angle of attack, the conditions are set almost uniformly within the range. For surface temperature, the conditions are chosen at equal intervals on both sides of the freestream temperature. For Mach number, the conditions are set with decreasing intervals within the range due to more complex flow under high Mach number. Detailed results of the simulated transition onset location are listed in Appendix A.

3.1. Influence of Turbulent Intensity

To investigate the influence of turbulent intensity on transitional flow around an RAE5243 airfoil, the other conditions are fixed as Ma = 0.4 and α = 0°. The transition onset locations under different operating conditions are depicted in Figure 5. According to Figure 5, it can be seen that a rise of inflow turbulent intensity will promote earlier transition on both the suction and pressure sides, and the influence on the suction side is larger than that on the pressure side.

The distribution of intermittency on the midplane of the RAE5243 airfoil under different levels of inflow turbulent intensity is depicted in Figure 6. It can be clearly seen that, with the rise of turbulent intensity, the non-zero region of intermittency gradually expands, and the thickness of the high-intermittency region around the airfoil also increases. This distribution is related the earlier transition onset under higher inflow turbulent intensity.

Figure 7 shows the distribution of intermittency along the streamwise direction under various turbulent intensities. An obvious sharp rise of the intermittency appears on both the suction and pressure sides, which indicates the onset of the transition. After the sharp rise, the intermittency gradually decreases on the suction side, while it remains at a relatively stable level on the pressure side. Finally, the intermittency rises rapidly near the trailing edge of the airfoil.

3.2. Influence of Surface Temperature

The influence of the RAE5243 airfoil’s wall temperature on the transition is discussed in this section under various angles of attack and Mach numbers, while the turbulent intensity is fixed at a value of 0.5%. The simulation results of the transition onset location are depicted in Figure 8, Figure 9, Figure 10 and Figure 11, with respect to different Mach numbers.

The influence of surface temperature on transition onset locations is relatively weak on both the suction and pressure sides. Under low-Mach-number conditions (Ma = 0.4, 0.6), the simulated transition onset location hardly varies under different temperatures. When it comes to higher Mach numbers (Ma = 0.7, 0.75), some discrepancies can be observed at different surface temperatures, especially on the suction side.

Figure 12 depicts the distribution of intermittency under two wall temperatures. The distributions under T_w = 293 K and T_w = 313 K are almost the same, and therefore the transition onset locations under these two temperatures are quite similar.

The distribution of intermittency along the streamwise direction under two surface temperatures is shown in Figure 13. The distribution patterns of the intermittency on both the suction and pressure sides are basically the same as those we discussed in the previous section. Moreover, it can be seen that larger oscillations appear under higher surface temperature.

3.3. Influence of Angle of Attack

The transition performances of the RAE5243 airfoil under different angle of attacks are compared in this section. The angles of attack vary from −6° to 6°, and the inflow turbulent intensity and surface temperature are fixed at 0.5% and 293 K, respectively. The simulated results are depicted in Figure 14.

Under all four values of Mach number, with the rise in the angle of attack from negative values to positive values, the transition on the suction side appears later, while the transition on the pressure side appears earlier. The variation of the transition onset location with respect to the angle of attack decreases as the Mach number rises. In addition, the transition onset location on the suction side is earlier than that on the pressure side under negative angles of attack, and the transition onset location on the pressure side is earlier than that on the suction side under positive angles of attacks. By comparison, under a zero angle of attack, the transition onset locations on the suction and pressure sides are relatively close.

The distribution of intermittency under various angles of attack is plotted, as shown in Figure 15. It can be found that the thickness of the non-zero region at the leeward side is always larger than that at the windward side. As the angle of attack rises from negative values to positive values, the suction side of the RAE5243 airfoil converts from windward side to leeward side, while the pressure surface is the opposite. Since the transition onset happens earlier on the suction side under a negative angle of attack and on the pressure side under a positive angle of attack, it can be concluded that the transition onset always happens earlier on the windward side.

Figure 16 illustrates the distribution of intermittency on the suction and pressure sides under various angles of attack. Under negative angles of attack, the rapid rise of the intermittency appears near the leading edge on the suction side. Similarly, under positive angles of attack, the rapid rise of the intermittency exists near the leading edge but on the pressure side. Such a rapid rise near the leading edge can also lead to a more oscillating distribution afterwards, especially under large angles of attack (α = ±6°).

3.4. Influence of Mach Number

In this section, the influence of Mach number on transition onset location is depicted in Figure 17. The values of inflow turbulent intensity and surface temperature are also fixed at 0.5% and 293 K, respectively.

From Figure 17, it can be seen that the variation law of the transition onset location with respect to Mach number is related to the angles of attack. Under a negative angle of attack, the transition happens later with the rise of the Mach number. By comparison, under a zero or positive angle of attack, the transition happens earlier with the rise of the Mach number.

To illustrate the influence of the Mach number on intermittency, Figure 18 and Figure 19 show the distribution of intermittency with various Mach numbers under negative and positive angles of attack, respectively. Under both negative and positive angles of attack, the non-zero regions of intermittency gradually expand as the Mach number rises. Under the relatively high Mach number (Ma = 0.75), there exists a sharp rise of intermittency on the suction side near the trailing edge. This is due to the appearance of a shock wave under a high Mach number.

The influence of Mach number on the streamline distribution of the intermittency under negative and positive angles of attack is illustrated in Figure 20 and Figure 21, respectively. As the Mach number increases, the intermittency can rise to a higher level, associated with a more complicated distribution of larger oscillations. This complicated distribution is even more significant on the suction side under negative angles of attack, and on the pressure side under positive angles of attack.

3.5. Correlation

On the basis of the above numerical results, the following correlation of transition onset locations with respect to turbulent intensity, surface temperature, angle of attack, and Mach number is established.

$$x_{tol,s}=260.52{\mathrm{Ma}}^{-0.6896}{\theta }^{11.5421}{\tau }^{-0.0438}{\epsilon }^{-0.1875}$$

(5)

$$x_{tol,p}=2.27\times {10}^{-6}{\mathrm{Ma}}^{-0.6514}{\theta }^{-15.2313}{\tau }^{0.1680}{\epsilon }^{-0.1823}$$

(6)

where x_tol,s and x_tol,p denote the transition onset location on the suction side and the pressure side, respectively. θ is non-dimensional angle of attack normalized within the range of [−π/2, π/2], as follows:

$$\theta =\frac{\alpha -\left(-\pi /2\right)}{\pi /2-\left(-\pi /2\right)}$$

(7)

and τ is the non-dimensional surface temperature normalized by freestream total temperature T_∞, as follows:

$$\tau =\frac{T_w}{T_{\infty }}$$

(8)

The above relations are determined by multiple linear regression based on the least squares method in logarithmic form. The R-square values for transition onset locations on the suction and pressure sides are 0.8621 and 0.8097, respectively. By means of the above correlations, the average errors of predicted transition onset locations on suction and pressure sides are 15.44% and 25.72%, respectively. Figure 22 shows the distribution of simulation–prediction pairs. It can be found that most prediction results are located within the error range of ±20%.

4. Conclusions

In the present study, a numerical investigation based on the k-ω-γ-Re_θ framework is carried out on transitional flows around an RAE5243 airfoil. The numerical methods are validated by NLF airfoil results and experimental measurement. The influences of inflow turbulent intensity, wall temperature, angle of attack, and Mach number on transition performance are analyzed.

The numerical results show that the rise of inflow turbulent intensity will promote an earlier transition on both the suction and pressure sides. The influence of wall temperature on transition is limited. The rise of angle of attack will lead to an earlier transition on the pressure side but a later transition on the suction side. With the rise of Mach number, the transition happens earlier under a zero or positive angle of attack but later under a negative angle of attack.

Based on the numerical results, the correlation of transition onset locations with respect to turbulent intensity, surface temperature, angle of attack, and Mach number is established. The average prediction errors on the suction and pressure sides are 15.44% and 25.72%, respectively.