Optimization of transonic low-Reynolds number airfoil based on genetic algorithm

Abstract

A 2-D airfoil shape optimization in transonic low-Reynolds number regime is conducted. A Navier–Stokes flow solver with a transition model (k-ω SST γ-Re_θ) is used to evaluate the fitness function. Single-point and multi-point formulations of the optimization results are compared. In addition, the effects of Mach number and angles of attack on aerodynamic characteristics of the optimized airfoils are investigated under low Reynolds number (Re = 17,000) and high-subsonic-flow ( $M_{a,\infty }=0.6-0.9$ ) conditions. The results show that the corresponding drag divergence Mach number curves of the conventional airfoil present almost a parallel shifting at the entire Mach number range. By contrast, the unconventional airfoil starts showing a significant drag reduction when Mach number is greater than 0.75. Besides, the maximum lift-to-drag ratio is highly influenced by the Mach number because of the formation, movement, type, and strength of a shock wave. In addition, the distinguishing difference in the conclusion between two airfoils is that the lift fluctuation of the conventional airfoil amplifies with the increase of the Mach number. However, the unconventional airfoil shows an opposite trend.

Introduction

Since 1997, seven rovers have successfully landed on the surface of Mars. However, only 0.3% of the planet's circumferential paths have been explored, which is mainly due to the hilly terrain and limited ground visibility. Rovers can significantly increase their range and maneuverability by providing an aerial view of upcoming paths by planning their route more effectively and avoiding obstacles on the route accurately. Such advantages and functionalities guide scientists to continue to improve and optimize the design of Mars aerial vehicles.

The design of airfoil’s high-Mach low-Reynolds numbers has gained scientific attention due to their application to Mars aerial vehicles. The low density of the Martian atmosphere and the relatively small Mars rotor result in very low chord-based Reynolds numbers (10³ <Re_c< 10⁴) ¹ which often introduce problems of laminar boundary layer separation on the suction side due to the adverse pressure gradient. Additionally, the low temperature and largely CO2-based atmosphere of Mars compound the overall aerodynamic problem by resulting in a lower speed of sound, indicating a limitation on the rotor’s maximum tip Mach number. As a result, a typical flying vehicle design needs to pay attention to: (1) the flow separation without reattachment occurred at the low Reynolds number condition, (2) relatively high blade rotational speed to ensure effective performance of Mars rotor, and (3) low Reynolds number compressible aerodynamics of which very little knowledge is available.

Subsonic low-Reynolds number airfoil design is now fairly well understood, ² the key issue is that the flow field has both strong unsteady and viscous effects, leading to a formation of complex vortical flow, among which a special phenomenon is the formation of laminar separation bubbles. ³ According to Carmichael’s survey, ⁴ the Reynolds number based on the mean chord length is around $5\times {10}^4$ , indicating the boundary-layer state is subcritical. Further decreaseing the chord Reynolds number, that is, $Re\hspace{0.17em}<5\times {10}^4$ , will experience laminar separation with no subsequent reattachment. However, the formation of a laminar separation bubble, as shown in Figure 1, is expected when the chord Reynolds number is higher than 5×10⁴. Laminar separation bubbles usually can get separation, transition, and reattachment occurred within a short distance. ⁵ In contrast, at chord Reynolds number below 30,000, complete laminar flow can occur for small angles of attack. The adverse gradients become more severe and laminar separation occurs as incidence increases, limiting the lift coefficient and significantly increasing the drag. At even lower Reynolds number, that is, ${Re}_c\approx 15000$ , the separation may occur over the entire rear of the airfoil, extending into the wake, which de-cambering effects can be expected. ⁶ Figure 2 shows the maximum section lift-to-drag ratio varies with Reynolds number. A sudden drop on maximum lift-to-drag ratio can be found when chord Reynolds number at approximately ${Re}_c\approx {10}^5$ .

Figure 1.

The mean flow structure of a laminar separation bubble.

Figure 2.

Maximum sectional lift-to-drag ratio versus Reynolds number. ³⁵

The fundamental research on transonic flow in the field of wing/airfoil aerodynamics was carried out in the early second half of the 20th century. Nieuwland and Spee ⁷ reviewed early work on transonic wing/airfoil designs. The supercritical airfoil design is based on the concept of local supersonic flow with isentropic recompression, and the relevant shape was characterized by a large leading-edge radius, reduced curvature over the middle region of the upper surface, and substantial aft camber. ⁸ The concept of aft shift of maximum camber is also utilized for ultra-low Reynolds number airfoil design, resulting higher attainable lift coefficients and lift-to-drag ratios by a less severe reduction of lift past the linear range. The aft shift of camber experiences reduces adverse pressure gradients and hence delays the separation. ⁹ For a given airfoil shape, the thickness ratio (t/c) and lift coefficient are often the most critical parameters influencing the drag divergence Mach number, as shown in Figure 3(a). In general, Reynolds number, relative thickness and location, leading-edge shape, and airfoil camber influence on transition location, and hence the aerodynamic characteristics, that is, Figure 3(b) shows the effect of a blunt trailing edge on the drag coefficient development at various Mach numbers. Recent work by Sudhi ¹⁰ and Badrya ¹¹ showed that the optimum hybrid laminar flow control airfoil has achieved a 25% lower drag than the NLF (natural laminar flow) airfoil.

Figure 3.

Transonic airfoil performance, (a) effect of thickness ratio on drag divergence ³⁶ and (b) effect of blunt trailing edge on drag coefficient. ³⁷

The high-Mach low-Reynolds number combination is unique in aeronautical applications and is unfortunately distinguished by a complete lack of suitable airfoil data, as shown in Figure 4. Early work on this topic was conducted by Drela. ¹² The computational study highlighted the relevant characteristics of laminar separation at a low Reynolds number, where “A”-shocks/boundary layer separation interaction was observed. As the Mach number increases, transition in the separation bubble will be delayed due to a reduction in growth rates of disturbance, resulting in an increased loss. Anyoji et al. ¹³ studied the effects of low Reynolds numbers ( $Re=0.43\times {10}^4–4.1\times {10}^4$ ), Mach number ( $Ma=0.1–0.6$ ), and gas species(air and CO₂) on the aerodynamic characteristics of a thin plat and a NACA 0012-34 airfoil experimentally. Their results indicated that the affection of the Reynolds number was on the formation, shift, and burst of the separation bubble, whereas the compressibility effects were on a delayed laminar-turbulent transition and reattachment of the separated shear layer. More recently, Şugar-Gabor and Koreanschi ¹⁴ presented a family of single-point and multi-point optimized airfoils at a high-lift condition in the Martian atmosphere. Their results indicated that significant differences in the airfoil shapes between transitional and fully turbulent flow were realized, and the optimized airfoil shows improved performance at both on-design and off-design conditions when the transition effects are considered.

Figure 4.

Flight condition on Mars.^12,38

Propellers for Mars UAVs operate under various operating conditions, ranging from vertical take-off to hovering conditions ¹⁵ and cruising. In addition, the angle of attack, Reynolds number, and Mach number for a rotor blade design varies from the root to the tip. Relatively poor performance would occur at the blade root when the Reynolds number is below 10⁵. On the other hand, the presence of transonic flow will, of course, influence the quantitative performance of the airfoils and their off-design behavior, such as a delayed drag divergence Mach number. These objectives result in design conflicts for new airfoils, as improvements in divergence Mach number require thin airfoils with limited camber, while improved maximum lift requires thicker airfoils with greater camber.

Increased payload, vertical take-off and landing capabilities lowered stall speeds, and extended endurance can all be derived from the beneficial effects of improved high-lift airfoil aerodynamics. However, working under the Martian atmosphere condition, the presence of a separated flow, together with a shock wave forming, introduces an interaction between the separation region and the shock wave, leading to pressure fluctuations, structural vibrations, and noise. It is, therefore, an essential task to gain a proper understanding of the phenomenon and seek possible ways to reduce the above-stated adverse impacts.

This paper aims to present an airfoil with good performance over a specified Mach range for the increasingly important transonic low-Reynolds number regime in which Mars rotor aerial vehicles operate. The work contributes to the verification and validation of state-of-the-art RANS-based aerodynamic shape optimization of a transonic low-Reynolds number airfoil. First, a validation case, NACA0012-34 airfoil, was carried out by comparing previous experimental results to current computational results. Second, the Genetic Algorithm^14,16 was adopted for airfoil optimization because of its global nature. The Reynolds number is 17,000, and the Mach numbers are 0.75 and 0.8 for optimization cases. The turbulence intensity (T_i) sets to 0.01% for all simulation cases. The pitching axis is 0.25c downstream of the leading edge on the chord line, as depicted in Figure 7. Finally, the aerodynamic performance of optimized airfoils under off-design conditions is studied with the Mach number varying from 0.6 to 0.9. In addition, unique features of transonic low-Reynolds number flows are revealed, including the evolution of different types of the shock wave, wake elongation, shock foot on the separation region, and strong shock wave on separation suppression.

Aerodynamics in Mars environment

The Martian environment has significantly impacted Mars rotor UAV design.^1,14,17 A comparison of atmosphere properties on Earth at sea level and Mars at ground level are listed in Table 1. However, it is worth mentioning that vehicles flying in the stratosphere would operate under similar conditions to those on Mars, including a combination effect of the low Reynolds number and air compressibility, although with higher gravity force on Earth.

Table 1.

Properties of Mars and Earth atmosphere.^1,17,38

	Mars (ground level)	Earth (sea level)	Ratio
Gas	95% CO2, 2.7% N2 1.6% Ar, 0.03% H2O	78% N2, 21% O2 1% Ar	─
Gravity(m/s2)	3.72	9.78	0.38
Density(kg/m3)	0.0175	1.225	0.0143
Pressure(Pa)	720 1	101300	0.0071
Temperature(K)	210 (–63.1°C)	288 (14.9°C)	─
Gas constant(J/g∙K)	0.189	0.287	0.66
Speed of sound(m/s)	223	340	0.66
Viscosity(kg/m∙s)	1.13×10−5	1.789×10−5	0.63

The low Reynolds number affects the lift and the drag characteristics, while the compressibility effects shear layer stability, that is, delaying the transition to turbulence. ¹⁸ The low temperature and low speed of sound limit the Mars rotor tip speed, indicating airfoils require satisfactory performance over the given set of Mach numbers.

Optimization framework

Airfoil shape parameterization

The selection of design variables and the parameterization of these variables is one of the most important factors in aerodynamic shape optimization. On the premise of ensuring the geometric characteristics, such as smoothness or continuity of an airfoil. Aerodynamic shapes typically are represented using Splines (e.g., B-spline and Bézier curves), ¹⁹ FFD, ²⁰ class-shape transformations (CSTs), ²¹ parametric sections (PARSEC), ²² Bézier-PARSEC,^23,24 and Bézier-GAN (generates smooth shapes that conform to Bézier curves). ²⁵ In the present work, the Bézier curve is used to describe the airfoil geometry (as shown in Figures 5 and 6). Clearly, the parameterization methods, Bezier curve fitting, could successfully approximate the shape of the airfoil profile. The advantage of this choice is the possibility of conjugating the properties of Bézier functions in terms of the regularity of the curve and easy usage with a piecewise structure that also allows to make a local modification to the geometry.

Figure 5.

Sensitivity analysis of parameterized control points for NACA airfoil.

Figure 6.

Sensitivity analysis of parameterized control points for UltraThin airfoil.

A Bézier parameterization is determined by its control points which are physical points in the plane, and is defined in parametric form by the following:

$\mathbf{P}\left(\mathrm{t}\right)={\displaystyle \sum }_{i=0}^n{B}_i^n\left(t\right)\cdot {\mathit{P}}_i$

(1)

where

$B_i^n\left(t\right)=\left(\begin{array}{l}n\\ i\end{array}\right)\cdot t^i{\left(1-t\right)}^{n-i}\cdot {\omega }_i$

is the ith Bernstein basis function in degree n.

$\left(\begin{array}{l}n\\ i\end{array}\right)=\frac{n!}{i!\left(n-i\right)!}$

The points P _i is called control points for the Bézier curve. The polygon formed by connecting the Bézier points with lines, starting with P₀ and finishing with P_n, is called the Bézier polygon (or control polygon). $B_i^n\left(t\right)$ is the Bernstein polynomial and ${\omega }_i$ their associated weights. The sensitivity study of the choice of 10, 12, and 14 control points for NACA are studied, and the results are shown in Figure 5(d and f). Finally, the integrated forces did not show a noticeable discrepancy between different airfoil parameterization expressions, as shown in the lift and drag coefficients in Figure 5 (d and e). The control points sensitivity study of UltraThin airfoil also agrees with this as shown in Figure 6 (c and d).

Flow solver and mesh

The flow solver used is a structured multiblock, finite-volume, cell centered scheme solving the compressible Reynolds-averaged Navier–Stokes equations, and it is expressed in Cartesian coordinates as follows:

$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_j}\left(\rho U_j\right)=0$

(2)

$\frac{\partial }{\partial t}\left(\rho U_i\right)+\frac{\partial }{\partial x_j}\left(\rho U_j{U}_i\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[{\mu }_{eff}\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right)\right]$

(3)

The viscous flux terms are evaluated using second-order central differences in space, and inviscid fluxes are evaluated with a third-order-accurate MUSCL scheme, and turbulent viscosity is modeled by the SST $\gamma -{Re}_{\theta }$ model, in which involves local-correlation–based transition modeling, including two extra transport equations, see equations (4) and (5). ²⁶ The deferred correction approach was implemented for the advantages of easier treating the higher-order approximations, ²⁷ and grid non-orthogonality. ²⁸ The transport equation for the transition-onset momentum-thickness Reynolds number is defined as:

$\frac{\partial \left(\rho {\stackrel{\sim }{Re}}_{\theta t}\right)}{\partial t}+\frac{\partial \left(\rho U_j{\stackrel{\sim }{Re}}_{\theta t}\right)}{\partial x_j}=P_{\theta t}+\frac{\partial }{\partial x_j}\left[{\sigma }_{\theta t}\left(\mu +{\mu }_t\right)\frac{\partial {\stackrel{\sim }{Re}}_{\theta t}}{\partial x_j}\right]$

(4)

The transport equation for intermittency is used to trigger the transition process ( $\gamma >0$ ) and is defined as:

$\frac{\partial \left(\rho \gamma \right)}{\partial t}+\frac{\partial \left(\rho U_j\gamma \right)}{\partial x_j}=P_{\gamma }-E_{\gamma }+\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_f}\right)\frac{\partial \gamma }{\partial x_j}\right]$

(5)

Furthermore, incoming freestream velocity, freestream turbulent intensity ( $T_i$ ) and the turbulence length scale are imposed at the inlet, whereas pressure is prescribed at the outlet boundary. The non-slip boundary condition is applied on the wall. The turbulent kinetic energy is set to zero and the pressure on the wall has zero normal gradients.

Mesh sensitivity study

Global optimization based on a genetic algorithm requires large numbers of solver calls. Therefore, a mesh sensitivity study is essential to balance minimizing numerical drag and minimizing the run time. A NACA0012-34 airfoil at a very low Reynolds number of $1.1\times {10}^4$ and a subcritical Mach number of 0.61 was selected to evaluate the mesh sensitivity study, and a wind tunnel work was performed by Anyoji et al., ¹⁸ using CO2 as the working fluid. These experimental results were obtained as part of a study aimed at understanding airfoil aerodynamics in the Martian atmosphere and are considered a suitably challenging verification case. As shown in Figure 7, a C-type structured mesh was generated around the NACA0012-34 airfoil. The faces of the control volumes form a curvilinear, body-fitted grid with coordinates (ζ, η). The same mesh topology ²⁹ but with different grid levels, as shown in Table 2, is applied for this study. A dense mesh with 50 grid points and an edge expiation ratio of 1.05 was applied inside the boundary layer in order to capture small, turbulent vortices and large flow gradients. The first cell distance was set to 2 X 10⁻⁵c for the present study to meet the requirement of a y⁺ value, $y^{+}<1$ , in which a large value of y⁺ may introduce a noticeable influence on leading-edge vortex formation and development. ³⁰ A computational domain size of 15c upstream and 30c downstream distance is considered in order to capture the potential pressure field around the airfoil, as shown in Figure 7(a).

Figure 7.

Mesh topology for NACA0012-34 airfoil. (a) C-type computational domain and (b) mesh in near wall region.

Table 2.

Grid-sensitivity study for NACA0012-34 at $\mathit{M}\mathit{a}=\mathbf{0.61}$ , $\mathit{R}\mathit{e}=\mathbf{1.1}\times {\mathbf{10}}^{\mathbf{4}}$ , and α = 4°.

Case	Size	Method	C l	C d
G1	300×45	Menter’s transition SST model 26	0.1836	0.0632
G2	300×90		0.1683	0.0599
G3	460×90		0.1293	0.0590
G4	460×180		0.1283	0.0588
G5	920×180		0.1280	0.0573
Wind 18	——	——	0.1179	0.0566

Mesh sensitivity analysis and verification of the numerical methods are carried out for NACA0012-34 high-Mach low-Reynolds airfoil in order to ascertain whether the selected grid density was of sufficient resolution and that the spatial discretization errors were minimal. The first grid, G1, listed in Table 2, was generated by estimating resolution requirements based on previous studies at lower Reynolds numbers. The second grid, G2, was subsequently produced to improve the simulation’s resolution. The maximum grid size of case G5 has 920 nodes streamwise, and 180 nodes are in the normal direction. The gird independence study shows little difference in aerodynamic coefficients, for example, lift coefficient, between G3, G4, and G5. All the RANS results in Figure 8 show a comparable solution in terms of aerodynamic coefficients. Case G5, with doubled grid points in streamwise and normal direction compared to G3, shows a much better agreement with the wind tunnel data, especially in drag coefficient as shown in Figure 8(d). Therefore, for the subsequent analysis, a grid size of 920×180 was considered a sufficient resolution for the present investigations.

Figure 8.

Mesh sensitivity study for NACA0012-34 airfoil at $Ma=0.61$ and $Re=1.1\times {10}^4$ in CO₂.

Airfoil shape optimization for Mars rotor blade

In our previous studies, ²⁹ the laminar separation bubble significantly influences the lift-curve slope and drag characteristics at a low Reynolds number. However, contrary to the Earth’s atmosphere, Mars atmosphere shows a strong interaction between viscous and compressibility effects, indicating a highly complicated flowfield would be expected around the wing or propeller. This paper further investigates the limits to airfoil performance in the unexplored high-Mach low-Reynolds number in Mars environment.

As already stated, a single-point design generally leads to unacceptable off-design performances; therefore, it is necessary to specify multi-point designs. The optimization problem presented here is the improvement of the drag-rise characteristics of high-Mach low-Reynolds number airfoils.

Geometry and design variables

A desired high-Mach low-Reynolds number airfoil should hold good aerodynamic performance and low sensitivity to the leading-edge roughness in unpredictable Mars environments. A third-order Bézier is used to describe the airfoil’s geometry. The Bézier curve parameters are tabulated in Table 3 for NACA and UT airfoils, of which the boundaries are carefully pre-selected to avoid impractical shapes.

Table 3.

Bézier parameterization with upper and lower bounds of their control points.

Control points	NACA8804		UltraThin airfoil
	x boundslower, upper	y boundslower, upper	x boundslower, upper	y boundslower, upper
CoP1	0.88, 0.92	0.070, 0.100	1.2e-3, 1.5e-3	3e-3, 5e-3
CoP2	0.72, 0.88	0.060, 0.080	0.056, 0.074	0.013, 0.017
CoP3	0.60, 0.80	0.010, 0.020	0.13, 0.18	0.050, 0.072
CoP4	0.15, 0.30	0.100, 0.150	0.17, 0.25	0.060, 0.075
CoP5	1.2e-3, 1.5e-3	0.10, 0.13	0.35, 0.56	0.018, 0.028
CoP6	1.0e-3, 3e-3	-0.015, -0.025	0.6, 0.75	0.010, 0.016
CoP7	0.20, 0.28	0.030, 0.050	0.67, 0.78	0.055, 0.070
CoP8	0.65, 0.76	0.050, 0.080	0.78, 0.85	0.050, 0.068
CoP9	0.70, 0.84	0.070, 0.120	0.88, 0.95	0.020, 0.030
CoP10	0.85, 0.89	0.070, 0.085	1.0e-3, 3e-3	-1e−2, -8e−3
CoP11	/	/	0.07, 0.086	4e-3, 7e-3
CoP12	/	/	0.12, 0.15	0.02, 0.04
CoP13	/	/	0.17, 0.20	0.05, 0.065
CoP14	/	/	0.23, 0.35	-9e−3, -1.4e-2
CoP15	/	/	0.45, 0.63	0.008, 0.015
CoP16	/	/	0.65, 0.75	0.015, 0.020
CoP17	/	/	0.76, 0.83	0.055, 0.075
CoP18	/	/	0.82, 0.90	0.018, 0.023

Although up to 10 (18 for UltraThin airfoil) control points have been adopted for airfoil’s parameterization, as shown Figure 5 (d and e) and Figure 6 (c and d) and each control point can be varied in the x and y direction. The designed bounds for the x and y direction are listed in Table 3. These four points primarily control the geometric shape at and near the leading edge of the airfoil, thus significantly influencing the aerodynamic performance.

Genetic algorithm

Genetic algorithms (GA) are based on Darwin’s evolution theory which is one of the gradient-free methods.^16,25,31 In some cases, the optimization problem is ill-conditioned, or the gradient of the objective is in-accessible. The gradientless method can avoid this problem and increase the possibility of finding the global optimal. The GA process includes mutation, recombination, and reproduction of new generation or designs. For a single-point genetic algorithm, the GAs starts an initial population production (i.e., airfoils), and real number encoding was utilized to reduce the storage compared to binary representations. The fitness of each aerodynamic performance is evaluated though a CFD solver. Survival of the fittest translates into discarding the chromosomes with the highest cost. The formulation for single-point and for multi-point will be different, as shown in Table 4. The roulette-wheel sampling method is utilized for pairing. The probability of selecting some individual p _i for reproduction is given by:

$p_i=\frac{fitness\left({Chrom}_i\right)}{{\displaystyle \sum }_{i=1}^{N_p}fitness\left({Chrom}_i\right)}$

(6)

Table 4.

Objectives for single-point and multi-point optimization problem.

Single-point formulation: $\min C_d\left(\mathit{C}\mathit{o}\mathit{P},\alpha ,{Ma}_{\infty }\right)$ , and $\mathrm{s}.\mathrm{t}.:\{\begin{array}{l}{C}_l\left(\mathit{C}\mathit{o}\mathit{P},\alpha ,{Ma}_{\infty }\right)\ge C_l^{*}\\ \begin{array}{l}t\left(\mathit{C}\mathit{o}\mathit{P}\right)\ge {0.98t}^{*},\text{for}\hspace{0.17em}\text{NACA}\hspace{0.17em}\text{airfoil}\\ t\left(\mathit{C}\mathit{o}\mathit{P}\right)\ge t^{*},\text{for}\hspace{0.17em}\text{UT}\hspace{0.17em}\text{airfoil}\end{array}\end{array}$

(8)

Multi-point formulation: ${\mathrm{C}}_{d,multi}={\displaystyle \sum }_{i=1}^k{\omega }_i\left[\frac{C_{d,i}\left(\mathit{C}\mathit{o}\mathit{P},{\alpha }_i,{Ma}_{\infty ,i}\right)}{C_{d,i}^{*}}\right]$ , and $\{\begin{array}{l}\min {\displaystyle \sum }_{i=1}^k{\omega }_i\left[\frac{C_{d,i}\left(\mathit{C}\mathit{o}\mathit{P},{\alpha }_i,{Ma}_{\infty ,i}\right)}{C_{d,i}^{*}}\right]\\ \min \hspace{0.17em}\left({\alpha }_i\right)\end{array}$ $\mathrm{s}.\mathrm{t}.:\{\begin{array}{l}{C}_l\left(\mathit{C}\mathit{o}\mathit{P},\alpha ,{Ma}_{\infty }\right)\ge C_l^{*}\\ \begin{array}{l}t\left(\mathit{C}\mathit{o}\mathit{P}\right)\ge {0.98t}^{*},\text{for}\hspace{0.17em}\text{NACA}\hspace{0.17em}\text{airfoil}\\ t\left(\mathit{C}\mathit{o}\mathit{P}\right)\ge t^{*},\text{for}\hspace{0.17em}\text{UT}\hspace{0.17em}\text{airfoil}\end{array}\end{array}$

(9)

where $\mathit{C}\mathit{o}\mathit{P}$ are control points, as shown in Figure 5(a) and Figure 6(a), for the airfoil with Bezier curve fitted, the asterisks(*) in the above equations refer to the parameters or aerodynamic properties associated with the original airfoil, ${\omega }_i$ is the adaptive weighting factor.

The pure-mutation operator takes one chromosome from the total number of populations with a random value, as shown following:

${Chrom}_{new}=RAND\left({Chrom}_{\mathit{max}}-{Chrom}_{\mathit{min}}\right)+{Chrom}_{\mathit{min}}$

(7)

After the mutations take place, the costs associated with the offspring and mutated chromosomes are calculated. The process described is iterated. The maximum number of generations depends on when an acceptable solution is obtained.

A weakness of GAs is their poor computational efficiency, which prevents their practical use when the cost function evaluation is computationally heavy. Oyama et al. ³² used an EA to redesign the NASA rotor67. The blade shape has 56 design parameters, with a population size of 64 and a generation size of 100 are utilized. However, Timnak, N. et al. ³³ applied the GA method for the optimum design of a transonic airfoil, the PARSEC method with nine design variables was applied for airfoil parameterization. The total population of each generation is set to be 20. In addition, the tournament operator is used with an elitist strategy where the first- and second-best chromosomes in each generation are directly transferred to the next generation. Flows over multi-element airfoils have been analyzed using GA methods by Benini, E. et al. ³⁴ The population size was chosen by multiplying the number of free variables (equal to 22) by a factor of 3 or 4. The value of 70 individuals per generation felt within the cited interval. The total number of generations performed was set to 30 artificially, corresponding to 2100 evaluations of the fitness function. Therefore, in the present case, the total number population, Np, is 100, and the probability values are 0.8 and 0.05 for crossover and mutation, respectively.

Results and discussion

This section shows the optimization results of the NACA8804 and UltraThin airfoils. In particular, starting from the original airfoil coordinates, the aerodynamic shape optimization results will be presented in the following arrangement, including results comparison of single-point and multi-point optimization; then, the airfoil performance is discussed, including effects of Mach numbers and Reynolds numbers.

Single-point Vs multi-point optimization

Two single-point optimizations, applying single-point formulation (8), for two different flight conditions, ${\mathrm{M}}_{a1,\infty }=0.75$ and ${\mathrm{M}}_{a2,\infty }=0.80$ at $R_e=\mathrm{17,000},$ are investigated separately. The constraints parameters for NACA 8804 and UltraThin (UT) airfoil are design lift coefficient ( $C_l^{*}$ ) and maximum thickness ( $t^{*}$ ): $C_l^{*}=0.581$ , $t_{NACA}^{*}=0.04$ , and $t_{UT}^{*}=0.02$ , indicating the lift of the original airfoil at ${\mathrm{M}}_{a1,\infty }=0.75$ , ${\alpha }_{NACA}\approx 0.982°$ , and ${\alpha }_{UT}\approx 1.926°$ for NACA8804 and UT airfoil, respectively. On the other hand, a two-point design is arranged for both airfoils, with the objective function obtained as a weighted factor summation of $C_{d,i}/C_{d,i}^{*}$ at the design points: ${\mathrm{C}}_{d,multi}=w_1{\left(\frac{C_{d,i}}{C_{d,i}^{*}}\right)}_{M_{a1,\infty }=0.75}+w_2{\left(\frac{C_{d,i}}{C_{d,i}^{*}}\right)}_{M_{a2,\infty }=0.8}$ , $w_1$ and $w_2$ are 0.6 and 0.4 for different design conditions. Four different runs have been carried out, with $\mathrm{\alpha }=0.5°,1.0°,1.5°,and\hspace{0.17em}2.0°$ , to obtain the drag values at the designed lift condition.

Figure 9(a) shows the geometries of the optimized airfoils. In all cases, the three NACA optimized airfoils are very smooth, and the percentage of aft camber is reduced compared to the original airfoil at design conditions. The optimized airfoils have a maximum thickness of ${\mathrm{t}}_{\mathit{max}\hspace{0.17em},NACA}=0.0392,0.0396,and\hspace{0.17em}0.0402$ with single-point $M_{a\infty }=0.75,0.8$ and multi-point, respectively. The UT airfoils in Figure 9(b) show that all optimized airfoils have a maximum thickness above $t_{UT}^{*}=0.02$ ( ${\mathrm{t}}_{\mathit{max}\hspace{0.17em},UT}=0.0202,0.0213,and\hspace{0.17em}0.0207$ with single-point $Ma=0.75,0.8$ and multi-point, respectively), indicating the optimizer could handle the design constraints. Figure 9(c and d) presents the drag-rise curves with the Mach number at a fixed C_L for viscous optimizations. The noticeable difference in drag reduction obtained from the UT airfoil between the optimal and the initial airfoils is when the Mach number is beyond 0.75. However, the NACA airfoil shows a drag reduction over the entire Mach number range. The NACA airfoil optimization history for the different formulation is plotted, see Figure 9(c), in terms of objective function. For the single-point optimization, more generations are required converge the solution as compared to lower Mach number case, that is, $M_{a\infty }=0.75$ . In fact, the optimizer suffering more resistance in reducing the shock wave drag. In addition, the multi-point shows the slowest convergence during the optimization as it needs to compromise all the conditions.

Figure 9.

Initial airfoil and single-/multi-point optimized airfoils at $M_{a1,\infty }=0.75$ and $M_{a2,\infty }=0.8$ .

Table 5 and Table 6 report the drag coefficients of the optimized airfoil and the improvement in percentage as compared with the initial airfoil. The aerodynamic characteristics of initial airfoils at corresponding flight conditions are evaluated first, as shown in the first line. Notice that the performance improvement of optimized transonic low-Reynolds airfoils is quite limited. However, one can observe that better improvement was achieved at higher Mach numbers for both airfoils, indicating the gain in total drag is mainly due to the reduction of wave drag.

Table 5.

$C_d$ at different $M_{a,\infty }$ for initial and the optimized airfoils: NACA8804.

	$M_{a1,\infty }=0.75$		$M_{a2,\infty }=0.80$		$M_{a,\infty }=0.75\hspace{0.17em}and\hspace{0.17em}0.80$
Design case	$C_d$	Improved (%)	$C_d$	Improved (%)	$C_d$	Improved (%)
Initial airfoil (NACA 8804)	0.0645	/	0.0736	/	0.0681	/
Single-point ( $M_{a1,\infty }=0.75$ )	0.0613	4.96%	0.0680	7.61%	0.0639	6.05%
Single-point ( $M_{a2,\infty }=0.80$ )	0.0595	7.75%	0.0650	11.68%	0.0617	9.40%
Multi-point	0.0572	11.3%	0.0636	13.6%	0.0598	12.25%

Table 6.

$C_d$ at different $M_{a,\infty }$ for initial and the optimized airfoils: UltraThin airfoil.

	$M_{a1,\infty }=0.75$		$M_{a2,\infty }=0.80$		$M_{a,\infty }=0.75\hspace{0.17em}and\hspace{0.17em}0.80$
Design case	$C_d$	Improved (%)	$C_d$	Improved (%)	$C_d$	Improved (%)
Initial airfoil (UltraThin)	0.0338	___	0.0405	___	0.0365	___
Single-point ( $M_{a1,\infty }=0.75$ )	0.0328	2.96%	0.0361	10.86%	0.0341	6.521%
Single-point ( $M_{a2,\infty }=0.80$ )	0.0329	2.66%	0.0352	13.10%	0.0338	7.34%
Multi-point	0.0331	2.07%	0.0356	12.10%	0.0341	6.58%

Figures 10 and 11 show the time-averaged plots of the multi-point optimized airfoil at the design conditions (i.e., C _l = 0.581). For clarity, Figure 10 a–e shows the Mach contour with a trailing-edge separation, a thickening in the wake separation is observed as the Mach number increases. Periodic lift fluctuations, such as the NACA airfoil in Figure 10(h), indicate the enhanced unsteadiness regarding a strengthened shock wave/separation interaction as the Mach number increases. When the Mach number reaches 0.84, the amplitude of lift fluctuation almost doubles compared with that at Ma = 0.8, which shows an earlier separation in Figure 10(f) marked as an arrow on the pink curve. In addition, the time-averaged pressure coefficient distribution in Figure 10(f) pink curve further highlights an increase in shock wave strength (i.e., a larger pressure rises). The corresponding two fundamental drag components, pressure drag C _dp, and friction drag C _df , have been expressed in Figure 10(g). A sharp increase in pressure drag coefficient is observed at Ma > 0.75. However, the friction drag coefficient is almost constant. This increase in drag is virtually due to the formation of the shock wave on the upper surface.

Figure 10.

Time-averaged plots of the multi-point optimized NACA airfoil at the design conditions of C _l = 0.581 and Re = 17,000.

Figure 11.

Time-averaged plots of the multi-point optimized UltraThin airfoil at the design conditions of C _l = 0.581 and Re = 17,000.

In contrast, Figure 11 a–e shows the airfoil wake variation under the influence of the Mach number, forming a pattern similar to that of a vortex-shedding in the far downstream region. As the Mach number reaches 0.84, a stretched and thickened wake pattern is formed. If we compare the periodic lift fluctuations between two airfoils as shown in Figure 10(h) and Figure 11(h), we immediately see that the UT airfoil shows a decreasing of amplitude in lift coefficient fluctuation as the Mach number increases, while the NACA airfoil shows the opposite trend. A possible reason is that the shock wave thickens the wake separation and hence tend to move the shock foot away from the airfoil surface (i.e., enhanced shock/separation interaction, which is a unique feature of low-Reynolds high-Mach number flow). On the other hand, the UT airfoil shows that as the Mach number increases a shock wave starts to form at the leading edge, and a shock-induced boundary layer separation is observed when Ma = 0.84. The corresponding separation location on the pink curve (i.e., Ma = 0.84) is indicated by a black arrow in Figure 11(f). Finally, the UT airfoil shows the pressure drag does not rise as fast as for the NACA airfoil for the entire Mach number range, as shown in Figure 11(g), indicating an attached flow was retained on upper surface of the UT airfoil.

Aerodynamic performance

To help understanding on the performance of the optimized airfoil under off-design conditions, the influences of the Mach number, from 0.6 to 0.9, on the aerodynamic characteristics of the optimized airfoils for a given Reynolds number are investigated, as shown in Figures 12 and 13. First, NACA and UT airfoils have shown improvements in aerodynamic efficiency at given constraints. A maximum lift-to-drag ratio of 13.9 and 20.5 are obtained for NACA and UT optimized airfoils, respectively. Second, the NACA (4% thickness) and UT airfoils (with 2% thickness) both shown a less variation on the $C_l-\alpha$ curves, that is, NACA airfoil shows from $5°\hspace{0.17em}to\hspace{0.17em}10°$ , and UT airfoil is from $1°\hspace{0.17em}to\hspace{0.17em}10°$ . Third, it can be observed that the angle of attack corresponding to the maximum lift-to-drag ratio decreases with Mach number, and the sharp reduction in ${\left(C_l/C_d\right)}_{\mathit{max}}$ at moderate incidences is believed to be the result of the substantial wave drag of the airfoils. Finally, the lift-curve slope $C_{l\alpha }$ is close to linear between $1°<\mathrm{\alpha }\hspace{0.17em}<\hspace{0.17em}8°$ for UT airfoils. In contrast, the gradient of the lift curve is slightly smaller than the theoretical inviscid lift slope ( $d{C}_l/d\alpha =2\mathrm{\pi }$ ) for NACA airfoils at moderate incidences.

Figure 12.

Effect of Mach number on aerodynamic characteristics of NACA-based airfoils at off-design conditions: optimized airfoil (1) at $M_{a1,\infty }=0.75$ , optimized airfoil (2) at $M_{a2,\infty }=0.8$ , and optimized airfoil (3) at $M_{a1,\infty }=0.75\hspace{0.17em}and\hspace{0.17em}0.8$ .

Figure 13.

Effect of Mach number on aerodynamic characteristics of UT airfoils at off-design conditions: optimized airfoil (1) at $M_{a1,\infty }=0.75$ , optimized airfoil (2) at $M_{a2,\infty }=0.8$ , and optimized airfoil (3) at $M_{a1,\infty }=0.75\hspace{0.17em}and\hspace{0.17em}0.8$ .

Figure 12 through (a) to (h) shows the effect of the Mach number on the lift and lift-to-drag ratio characteristics of NACA and UT airfoils for a fixed Reynolds number. In regards to lift, it can be pointed out that the Mach number has a significant influence at $\mathrm{\alpha }\hspace{0.17em}<\hspace{0.17em}5°$ and $\mathrm{\alpha }\hspace{0.17em}>\hspace{0.17em}10°$ for NACA airfoils, including the initial case. Furthermore, the NACA airfoils show less nonlinear behavior in lift than UT airfoils at moderate angles of attack, roughly $5°\le \mathrm{\alpha }\hspace{0.17em}<\hspace{0.17em}10°$ . For $\mathrm{\alpha }\hspace{0.17em}>\hspace{0.17em}8°$ , the lift slope is clearly below $2\mathrm{\pi }$ , which is caused by the turbulent separation on the airfoil’s upper surface. As further increases incidence, the $C_l-\alpha$ curve starts to level off at a rough incidence of $15°$ and stays near the flat until $\mathrm{\alpha }=30°$ without rising or failing. In contrast to the results for the UT airfoil as described in Figure 13 (a, and c to e), a higher lift coefficient is obtained by NACA-based airfoils at $\mathrm{\alpha }\hspace{0.17em}>\hspace{0.17em}15°$ for all Mach numbers, that is, $C_{l,at\hspace{0.17em}\alpha >15°}\hspace{0.17em}>\hspace{0.17em}1.5$ .

Figures 12 and 13 (b, and f to h) show the effect of the Mach number on each airfoil (2% t/c UltraThin, and 4% NACA airfoils) separately in terms of the lift-to-drag ratio ( $C_l/C_d$ ) as a function of the angle of attack. A noticeable difference between NACA and UT airfoils in the peak value of $C_l/C_d$ was observed. Compared with the NACA airfoil, the UT airfoil shows a higher lift-to-drag ratio at a lower angle of attack, $C_l/C_d=20.5\hspace{0.17em}at\hspace{0.17em}\alpha =4°$ in Figure 13(g), whereas the NACA airfoil showed a peak value of $C_l/C_d=13.9\hspace{0.17em}at\hspace{0.17em}\alpha =5°$ in Figure 12(g). It has been shown that for a circular cambered plate airfoil, 6.35% camber and 1% thickness with sharp leading-edge, the peak value of $C_l/C_d$ was around 8 at $\alpha =5°$ for $M_{a,\infty }=0.7$ and $R_e=3000$ . ¹⁷ As the Mach number increases from 0.6 to 0.9, the peak value of $C_l/C_d$ for UT airfoil in Figure 13 (g) reduces from 20.5 to 13.17 with $\alpha =5°\hspace{0.17em}and\hspace{0.17em}4°$ , respectively. In addition, the peak value of the lift-to-drag ratio for UT airfoil, in Figure 13 (b, f to g), is more sensitive to the change of angles of attack. This is caused by an abrupt separation of the laminar flow near the leading edge as the angle of attack increases to modest values. The NACA airfoil, on the other hand, shows a smooth reduction in $C_l/C_d$ as incidence changes, associated with a trailing-edge stall.

To provide insight into the mechanisms responsible for the trends in the lift and the lift-to-drag ratio, as shown previous, Figure 14 displays shock wave evolution associated with changes in angles of attack, which were obtained for Mach numbers 0.6, 0.75, 0.8, 0.85, and 0.9 at Reynolds number 17,000. At low incidence $\alpha =2°$ in Figure 14, the laminar boundary layers on the upper surface of airfoils are almost entirely attached to the surface, except a trailing-edge separation has occurred on the NACA airfoil, that is, most of the attached flow is shown before 70% of the chord length. In contrast, the UT airfoil shows an unsteady flow pattern in the near wake region. The vortices shed alternatively from the upper shear layer of the separated boundary layer and the lower shear flow. Moreover, the compressibility for NACA airfoil affects mainly the pressure suction peak value, and an increase in the slope of the favorable pressure gradient is obtained, as shown in Figure 15(a), which indicates the flow accelerates continuously into the shock. As the Mach number increased to 0.85, a $\mathrm{\lambda }$ -shock was formed above the separation region on the NACA airfoil upper surface, followed by a secondary shocklet. Since a pressure rising in the flow was introduced by the main shock, as shown in Figure 15 (a, blue line), which propagates upstream in the subsonic part of the boundary layer and hence thickens the trailing-edge separation region in Figure 14 (a, $M_{a,\infty }=0.85\hspace{0.17em}and\hspace{0.17em}\alpha =2°$ ). As the Mach number further increases to $M_{a,\infty }=0.9$ , the shock has strengthened and moved rearward, where the compressibility strongly effects the trailing-edge separation, leading to a stretched wake pattern, shown in Figure 14 (a, $M_{a,\infty }=0.9\hspace{0.17em}and\hspace{0.17em}\alpha =2°$ ).

Figure 14.

Mach contours for optimized airfoil at $M_{a,\infty }=0.6,0.75,0.8,0.85,and\hspace{0.17em}0.9$ , Re = 17,000, $\mathrm{\alpha }=2°,5°,and\hspace{0.17em}10°$ , showing shock pattern varies as angles of attack and $M_a$ changes.

Figure 15.

Time-averaged pressure distribution for optimized airfoil at $M_{a,\infty }=0.6\hspace{0.17em}to\hspace{0.17em}0.9$ , showing (a–c): NACA airfoil, (d–f): UT airfoil.

By contrast, in Figure 15(d), the UT airfoil shows a pressure suction peak followed by an adverse pressure gradient, and the associated suction peak location moves rearwards as the Mack number increases. In addition, trailing-edge separation, x/c = 0.872, is observed on the upper surface of the UT airfoil at $M_{a,\infty }=0.6$ . As the Mach number increases to 0.75 and 0.8, a weak shock is formed on the nose of the airfoil, as shown in, as shown in Figure 14 (b, at $\alpha =2°$ ), and the induced separation occurred at x/c = 0.293 and x/c = 0.261, respectively. However, it should be noticed that a $\mathrm{\lambda }$ -shock formed near the leading edge, and the shock-induced separation occurred at $x/c=0.41$ , which eventually could determine the initial deflection of the separated shear layer.

At moderate angles of attack $\alpha =5°$ in Figure 14, the UT airfoil shows a leading-edge separation at $M_{a,\infty }<0.8$ , whereas a trailing-edge separation is found on the upper surface of the NACA airfoil. An A-shock structure have been captured on the up surface of the airfoil when the $M_{a,\infty }=0.8$ , this is because the laminar boundary layer ahead of the first shock is separated due to the formation of a precompression wave, which impinges on the separated free shear layer and reflects as an expansion wave, indicating a formation of subsequent A-shocks, as shown in Figure 14 (b, $M_{a,\infty }=0.8\hspace{0.17em}and\hspace{0.17em}\alpha =5°$ ).

When the airfoil rises at $10°$ of angles of attack, as shown in Figure 14, the separation point moves toward the leading edge, resulting in a pronounced shear layer. With increasing the Mach number, $M_{a,\infty }>0.8$ , the pressure distribution of the NACA airfoil, in Figure 15(c), indicates a significant effect of compressibility on the leading-edge flow attachment, demonstrating the appearance of a favorable pressure gradient. However, at $M_{a,\infty }=0.85\hspace{0.17em}and\hspace{0.17em}\alpha =10°,$ the UT airfoil shows a separation location at $x/c=0.22,$ where an associated adverse pressure gradient is existed, as shown in Figure 15 (f, marked as blue line). When the Mach number further increases at $M_{a,\infty }=0.9$ , almost fully attached flow, except at the trailing edge, was observed on the up surface of the NACA and UT airfoil, as shown in Figure 14 (a, and b, at $\alpha =10°$ ).

Conclusion and outlook

A Reynolds-averaged Navier–Stokes solution implemented with the κ-ω SST γ-Re_θ turbulence-transition model approach was applied to the study of transitional flows on transonic low-Reynolds number airfoils in the Martina atmosphere. First, a genetic algorithm based on Bezier curve parameterization combined with single-/multi-point objective function formulations was conducted to investigate the drag reduction characteristics. In addition, the current work also revealed a distinct flow characteristic: shock wave patterns vary above the flow separation region on the upper surface of airfoils as the Mach number and angles of attack increase. Finally, compared to the NACA airfoil, the UltraThin airfoil forms a shock wave near the leading edge. The shock movements will bring local but very severe changes in loading, hence bring important significance in the contra-rotating propeller design. The following specific conclusions were drawn from this study:

1) Based on the drag divergence Mach number results of NACA (4% thickness) and UltraThin airfoil (2% thickness), the single-/multi-point optimized airfoil shows some Mach number dependency. The corresponding ${\mathrm{C}}_d$ vs. ${\mathrm{M}}_{a,\infty }$ curves of UT airfoil present almost identical at $M_{a,\infty }<0.7$ . However, the drag reduction effect is significantly improved for higher Mach numbers, that is, $M_{a,\infty }\ge 0.75$ . Therefore, a further assessment for multi-point formulation on drag reduction of transonic low-Reynold number airfoils at higher Mach numbers is necessary.

2) For the performance at off-design, the Mach number shows a negligible effect on the lift coefficient at moderate incidences. On the other hand, the lift-to-drag ratio significantly reduces as Mach number increases, indicating an increase in the drag coefficient. Besides this, the lift-curve slope is nearly linear at moderate angles of attack for NACA and UltraThin airfoils, that is, $\mathrm{\alpha }=5°$ . As angles of attack increases, that is, $\mathrm{\alpha }>\hspace{0.17em}15°$ , the $C_l-\alpha$ curve tends to level off, and the effect of compressibility on the lift coefficient is clearly observed at high angles of attack.

3) The Mach number shows a significant impact on the flow past airfoils, including of flow separation, trailing-edge wake patterns, and shock wave types. A stretched trailing-edge separation pattern is clearly observed from both NACA and UltraThin airfoils, which indicates the suppression of compressible flow on early separated shear layer roll-up motion, hence resulted in a smoother ${\mathrm{C}}_l-\alpha$ curve at beyond stall angle, that is, preventing the vortex-shedding.

4) The noticeable differences between NACA (conventional design) and UltraThin airfoil (unconventional design) are shown as following: a) At $M_{a,\infty }\le 0.8$ , NACA airfoil shows a trailing-edge separation, and a shock wave starts to appear at moderate to high incidences. By contrast, UT airfoil shows a leading-edge separation as incidence increases. As angles of attack rise at moderate, an “A-type” shock appears on the up-airfoil’s surface, which alters the response of the outer flow pressure to displacement surface perturbations, including the influence on the growth, curvature, and even unsteadiness of the separated shear layer. b) At $M_{a,\infty }=0.85$ , a $\mathrm{\lambda }$ -shock appears above the separation region near the trailing edge for the NACA airfoil, and shock with increased strength were found at moderate to high incidences. Whereas UltraThin airfoil shows a leading edge $\mathrm{\lambda }$ -shock, the shock-induced separation occurred at 41% of the chord at $\mathrm{\alpha }=2°$ . c) At $M_{a,\infty }=0.9$ , the increased size of the shock wave was found, which covers almost the entire chord length, including the squeezing and stretching of the wake.