Computational methods for investigation of surface curvature effects on airfoil boundary layer behavior

Abstract

This article presents computational algorithms for the design, analysis, and optimization of airfoil aerodynamic performance. The prescribed surface curvature distribution blade design (CIRCLE) method is applied to a symmetrical airfoil NACA0012 and a non-symmetrical airfoil E387 to remove their surface curvature and slope-of-curvature discontinuities. Computational fluid dynamics analysis is used to investigate the effects of curvature distribution on aerodynamic performance of the original and modified airfoils. An inviscid–viscid interaction scheme is introduced to predict the positions of laminar separation bubbles. The results are compared with experimental data obtained from tests on the original airfoil geometry. The computed aerodynamic advantages of the modified airfoils are analyzed in different operating conditions. The leading edge singularity of NACA0012 is removed and it is shown that the surface curvature discontinuity affects aerodynamic performance near the stalling angle of attack. The discontinuous slope-of-curvature distribution of E387 results in a larger laminar separation bubble at lower angles of attack and lower Reynolds numbers. It also affects the inherent performance of the airfoil at higher Reynolds numbers. It is shown that at relatively high angles of attack, a continuous slope-of-curvature distribution reduces the skin friction by suppressing both laminar and turbulent separation, and by delaying laminar-turbulent transition. It is concluded that the surface curvature distribution has significant effects on the boundary layer behavior and consequently an improved curvature distribution will lead to higher aerodynamic efficiency.

Introduction

As the basic element of wings or blades, airfoils are geometric profiles that are designed to generate lift force due to a relative motion in a fluid. Their aerodynamic performance directly determines the overall performance and efficiency of a system. Therefore, airfoil design plays a key role in an aeronautical device design. Early design methods were based on specifying a thickness distribution around a camber line. ¹ These methods did not provide enough flexibility to control both the suction and pressure surfaces to obtain a desirable aerodynamic performance. Modern designers use different definitions for airfoil and blade design methods: direct design, ² semi-inverse design, ³ full-inverse design, ⁴ or full- optimization methods. ⁵ The analysis of aerodynamic performance can be studied by several methods, including Conformal Mapping, ⁶ Thin Airfoil Theory, ⁷ Surface Panel Method, ⁸ or computational fluid dynamics (CFD) solutions.^9,10

The working airflow moves over an airfoil along the surface curves, rather than Cartesian coordinates. When the boundary layer equations are written in curvilinear coordinates, the equations show that local pressure on an airfoil surface has a strong dependence on local radius of curvature. ¹¹ Local curvature is determined by the radius of local curvature, which is directly related to the second-order derivative of local surface equation $y″$ . The importance of $y″$ for aerodynamic performance can also be found in the thin airfoil theory. Van Dyke ⁷ presented the second-order extension of thin airfoil theory showing that the surface pressure and velocity distributions are strongly dependent on the second-order derivative of the airfoil surface function y. Thus, smooth velocity and pressure distributions along the airfoil surface require smooth second-order derivative of the airfoil surface function, and hence smooth surface curvature distributions (i.e. continuous slope-of-curvature distributions) along the surface. A continuous slope-of-curvature requires a continuous third-order derivative of the surface function.

One must distinguish between “smooth surface” and “smooth curvature distribution” as these terms are often confused in the literature. When we say an airfoil has a smooth surface, mathematically it means that the first-order derivative of surface function $y\text{'}$ is continuous, not the second or third order. However, smooth curvature distribution means that it has a continuous slope-of-curvature distribution, which requires that the third-order derivative of the surface function is continuous. In this article, we deal with the latter concept and we address “smooth curvature distribution” as “continuous slope-of-curvature distribution” in order to be consistent. Surface curvature discontinuity often exists at particular areas on airfoils such as the leading edge (LE), e.g. NACA four-digit airfoils, or the trailing edge (TE), e.g. non-symmetrical sharp or blunt TE airfoils. Many blades and airfoils present slope-of-curvature discontinuities, at the point where the LE circle or other shape joins the main part of the airfoil. Whether discontinuity of slope-of-curvature is present depends on the airfoil design methods and airfoil coordinates resources. The discontinuities in the distribution of surface curvature and slope-of-curvature are invisible to the observer in the airfoil shapes and the airfoil surfaces can appear to be smooth. However, the airfoils with these discontinuities often produce unusual loads, higher losses, and thicker wakes. ¹²

To eliminate unnecessary local changes of surface curvature, surface curvature distribution becomes one of the key factors in the design of high-efficiency airfoils and blades.^13–15 Siddappaji et al. ¹⁶ developed a parametric 3D blade design tool for turbomachinery, and they used the definition of splines to modify the blade shapes and obtain the desired blade model. In their research, the curvature as well as the slope-of-curvature distributions of the blade surface are both continuous due to the application of the splines. This design tool was later improved by Nemnem et al. ¹⁷ by the addition of a 2D curvature-defined mean-line blade airfoil geometry generator. They created a fifth-order mean-line by twice integrating the cubic B-spline, which ensures curvature and slope-of-curvature distributions are continuous, and improved the performance of newly designed blades. A design method called the prescribed surface curvature distribution blade design (CIRCLE) method ¹⁸ was proposed to optimize airfoils and blades by ensuring continuous curvature and slope-of-curvature distribution along their surfaces, and showed that the aerodynamic and heat transfer performance strongly depends on surface curvature distribution. Massardo et al.^19,20 used streamline curvature distribution calculations to determine the 3D variation of inlet and outlet flow angles for axial-flow compressor design and improved the compressor efficiency. Based on surface curvature distribution, Song et al.^13,21 showed that improving the curvature continuity at the leading-edge blend point of a compressor blade improved performance by helping to eliminate the separation bubble. They also proved that the main surface curvature continuity improves blade performance. Korakianitis et al. ²² applied the CIRCLE method to two wind turbine airfoils and concluded that more continuous surface curvature and slope-of-curvature distributions improved the aerodynamic performance compared with the original wind turbine airfoils.

Although the effects of surface curvature discontinuities are widely researched in turbomachinery blade and airfoil design, more research attention is required to examine the impact of surface curvature distribution on the behavior of the airfoil boundary layer in external aerodynamics applications. At low Reynolds numbers, the flow on airfoils separates due to a sufficiently large magnitude of the adverse pressure gradient and changes in flow geometry including local surface curvature variations. The flow experiences transition to turbulence and can result in a stall (without turbulent re-attachment) or a laminar separation bubble (LSB; with turbulent re-attachment). ²³ Subjected to inherent effects of flow separation and stall, the airfoils operating at low Reynolds numbers do not always offer good aerodynamic performance. ²⁴ Therefore, more research is needed on the aerodynamic mechanisms causing the improved performance yielded by airfoils with continuous slope-of-curvature distributions.

In this study, the CIRCLE method ²⁵ is used to redesign airfoils by removing their surface curvature and slope of surface curvature discontinuities from LE to TE. A symmetrical airfoil NACA0012 and a non-symmetrical airfoil E387 are judiciously selected for this investigation due to the existence of discontinuities of curvature and slope-of-curvature, their widespread use and the availability of detailed experimental measurements of airfoil performance.^27,28 NACA 0012 is used as an example of an airfoil with a surface curvature discontinuity and E387 is used as an example of an airfoil with slope-of-curvature discontinuities. Reynolds Averaged Navier–Stokes (RANS) computations are applied to both airfoils to numerically investigate the effects of surface curvature on boundary layer behavior, which significantly affects airfoil aerodynamic performance. An inviscid–viscid interaction flow computing scheme of thin airfoil theory and boundary layer models named TATO is used to calculate the external flow velocity and LSB characteristics of the original E387 and the airfoil optimized with the CIRCLE method. The flowchart of this study is presented in Figure 1.

Figure 1.

The main flow chart for the design and computation algorithms in this study. $\mathrm{Curv}$ and $\mathrm{Curv}\text{'}$ represents curvature and slope-of-curvature, respectively.

Methodology

Design approach: The CIRCLE method

In this section, the CIRCLE method for the section side of an airfoil is summarized, and the design procedure for the pressure side of the airfoil is similar. ²⁵ Curvature $\mathrm{Curv}$ and slope-of-curvature $\mathrm{Curv}\text{'}$ are defined as follows:

$$\mathrm{Curv}=\frac{1}{r}=\frac{y″}{(1+y\text{'}2)3/2}$$

(1)

$$\mathrm{Curv}\text{'}=\frac{d(\mathrm{Curv})}{\mathrm{dx}}=\frac{y″\text{'}(1+y\text{'}2)-3y\text{'}y\text{'}\text{'}2}{(1+y\text{'}2)5/2}$$

(2)

where r denotes radius of curvature, $(x,y)$ denotes Cartesian coordinates and $y\text{'}$ , $y″$ , $y″\text{'}$ are the first-, second-, and third-order derivative of y. The sign of the curvature is usually defined as the direction of the unit tangent vector moving along the curve. We define curvature as positive if the vector rotates clockwise, otherwise it is negative. As described in the equation, the sign of the curvature is always the same as the sign of the second derivative of airfoil surface function $y″$ .

Circular geometry is used for LE and TE design to illustrate the design process, yet in principle, any proper geometry (e.g. ellipse) with a continuous slope-of-curvature distribution can be used to design. The CIRCLE method starts from the TE circle and designs the airfoil shape in three line segments sequentially: y3 near the TE (treated as a small circle in this section), y2 in the middle part of the surface, y1 near the LE (also treated as a small circle in this section), and the LE circle, as shown in Figure 2(a). By specification, the method ensures that airfoil surface curvature and slope-of-curvature is continuous from the LE stagnation point to the TE stagnation point. This section describes the approaches for the design of the TE, main body, and LE.

Figure 2.

The CIRCLE method for an airfoil design: (a) Main blade geometry (adapted from Korakianitis et al. ²⁶ ); (b) trailing edge design; (c) leading edge design; (d) suction side surface prescribed curvature distribution.

TE design

The most difficult and complicated part of TE design is joining the TE circle to the airfoil surfaces, as there is a transition from the constant curvature of the circle region to the locally varying curvature of the remaining airfoil surface. The TE circle is located by the TE origin and its radius. On the suction side, the airfoil surface connects to the TE circle at point P _s2 , as shown in Figure 2(b). The position of point P _s2 is specified by the input parameter TE origin and its local airfoil surface angle β _s2 . The line segment y3 extends from point P _s2 to P _sm as in Figure 2(a) and (b), and is specified by an analytic polynomial y = f(x) of the exponential form:

$$y3=f(x)=c_0+c_{1}x+c_{2}x2+c_{3}x3\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}+c_4{k}_1[x-x(P_{s2})]+c_5{k}_2[x-x(P_{s2})]$$

(3)

where k ₁ and k ₂ are exponential functions, c ₀ to c ₅ are thickness coefficients. k ₁ and k ₂ present increasing importance when approaching point P _s2 . Thus, the equation can be considered as a cubic equation near point P _sm ; and as exponential modifications of the basic cubic equation when it approaches the TE circle at point P _s2 .

The six thickness coefficients c ₀ to c ₅ are evaluated from six conditional equations, which are the conditions of the point continuity, the first-, second-, and third-order derivative continuity (four conditions) of the airfoil surface line at P _s2 , and prescribing the point and its first-order derivative continuity at P _sm (two additional conditions). As the TE circle is usually small, the variation in the streamwise curvature from the TE circle to the line segment y3 is usually large. The CIRCLE method enforces curvature and slope-of-curvature continuity through these locally large streamwise changes of curvature. ¹⁸

The circular (or elliptical) form of the TE is necessary. It guarantees the slope-of-curvature continuity near TE region, which provides smooth surface curvature distribution and therefore smooth pressure and velocity distribution, and also avoids a TE mathematical singularity (sharp TE) causing a bursting increase of the skin friction at the TE region, which will be described in the following sections. Moreover, it is easier to manufacture as compared with airfoils with a sharp TE as in reality, it is impossible to produce a zero thickness TE.

Main body design

Bezier splines are used to design the line segment y2 between points P _sm and P _sk . y2 is constructed by mapping the curvature distribution from the original airfoil geometry. The curvature segment from P _s2 to P _sm is evaluated from the analytic polynomial y3 (equation 3) starting from the TE and ending at point C6s in Figure 2(d). The slope of the curvature at point C6s is computed from equation (3) and becomes an input parameter to further calculations. Points C1s to C5s are user specified parameters. Point C1s is specified at an x location corresponding to the point P _sk . The tangency condition at point C6s ensures slope-of-curvature continuity from P _sm to P _sk because the slope of the Bezier curve is tangent to the line of both ends of the curve.

Starting from point P _sm , the design progresses explicitly point-by-point toward the LE to points P _sk and we obtain each y_i + 1 by using central differences equation for curvature at airfoil point i as a function of (x, y) coordinates of points i − 1, i, and i + 1. The Bezier spline is iteratively manipulated until the slope and the y location of the airfoil surface are achieved at points P _sk . ²⁵

LE design

The CIRCLE method used a hybrid method to design the LE region. The suction side surface connects the LE circle at point P _s1 , which is specified by an input parameter LE local surface angle β _s1 . Here, a parabolic construction line is defined and a thickness distribution is added perpendicularly to the construction line, as in Figure 2(c). The construction line starts from the LE of the airfoil, or the center of the LE circle. A thickness distribution is added orthogonally to this parabolic construction line, so that the thickness distribution has continuous point, first, second, and third derivative (continuous y, y′, y″, y″′ and therefore continuous Curv′) at both points P _s1 and P _sk , as discussed below. Point P _s1 is the point connecting the LE thickness distribution to the LE shape. Point P _sk is the point connecting resultant line y1 (from the thickness distribution) to the main part of the airfoil y2, corresponding to curvature point C1s, as in Figure 2(d). The suction side construction line can be of the form: y(x) = Ax² + Bx + C, and the thickness distribution y _t can be expressed as:

$$y_t=c_0+c_{1}x+c_{2}x2+c_{3}x3+(c_4{k}_{11}+c_6{k}_{13})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\times [x-x(P_{s1})]+(c_5{k}_{12}+c_7{k}_{14})[x-x(P_{\mathrm{sk}})]$$

(4)

where functions k ₁₁ , k ₁₂ , k ₁₃ , and k ₁₄ are exponential polynomials, similar to k ₁ and k ₂ as explained in the subsection of TE design. These exponential functions have increasing importance when approaching the end points P _s1 and P _sk of line segment y1. Equation (4) is a cubic polynomial away from these two end points. The eight parameters of the thickness function c ₀ to c ₇ need eight equations and they are derived from the conditions to match: y, y′, y″, and y″′ at point P _s1 and at point P _sk , respectively. This approach ensures continuity of curvature and slope-of-curvature from the TE circle to the main part of the airfoil surface through the leading-edge thickness distribution and finally into the LE circle.

Computational methods

Two computational methods are used in this article. An inviscid–viscid interaction scheme called Thin Airfoil Theory cOde (TATO) is introduced as an engineering algorithm to calculate LSB positions. RANS computations are applied to simulate the flow field of the airfoils. The results of these two methods are compared in the results section.

TATO

An inviscid–viscid interaction scheme is adopted for the LSB and flow transition calculation. The inviscid solution is achieved by thin airfoil theory calculations, providing flow field information including the external velocity of the airfoil. Boundary layer calculations are carried out to determine boundary layer displacement effects, supplying separation, transition, and re-attachment states. The general procedural elements are presented in Figure 3.

Figure 3.

The program flow chart for Thin Airfoil Theory Code (TATO).

Classical thin airfoil theory

The profile of the upper and lower surfaces of an airfoil can be described by the camber and thickness functions of the airfoil. A solution based on potential incompressible theory of thin airfoil of moderate thickness and camber at a moderate angle of attack was proposed by Van Dyke. ⁷ The inviscid external velocity q is determined by two horizontal velocity disturbances on the chord line assuming first- and second-order perturbations.

Boundary layer methods

Thwaites' method ²⁹ is adopted to determine the location of the laminar separation point. Thwaites proposed a laminar separation parameter m:

$$m=-\frac{\theta 2}{\nu }\frac{\mathrm{dq}}{\mathrm{dx}}$$

(5)

where θ denotes boundary layer momentum thickness and ν denotes air kinematic viscosity. He suggested that the laminar flow is on the separation point when m = 0.082. It was later corrected to 0.09 by Curle and Skan. ³⁰ Michel's empirical formula are used to calculate the location of transition according to the momentum thickness of local boundary layer. ³¹ The formula suggests that transition happens when the Reynolds number based on the momentum thickness Re _θ meets the following equation:

$$R{e}_{{\theta }_{\mathrm{tr}}}=1.174(1+\frac{22400}{R{e}_{x_{\mathrm{tr}}}}){\mathrm{Re}}_{x_{\mathrm{tr}}}^{0.46}$$

(6)

where $R{e}_{\theta }=q\theta /\nu$ , $R{e}_x=\mathrm{ql}/\nu$ and l is the curve length between local position on airfoil surface and stagnation point as the reference length. If the separation point is before the transition point, a LSB forms. The turbulence re-attachment can be calculated with Horton's criterion ³² and Roberts' correction ³³ :

$${\Lambda }_R=(\frac{\theta }{q}\frac{\mathrm{dq}}{\mathrm{dx}}{)}_R=-0.0059$$

(7)

where Λ _R is a re-attachment parameter and q is the inviscid external velocity.

RANS equations

RANS computation is used to investigate the aerodynamic characteristics of both the optimized airfoils and original airfoils. ANSYS ICEM CFD 16.0 is used as the mesh generator for all the cases. In the case of NACA0012, an airfoil study with angle of attack (AoA) ranging from 0 ° to 14 ° was conducted for validation and comparison with the experiments done by Lee and Gerontakos. ²⁸ The simulated chord based Reynolds number is $1.35\times 105$ with a turbulence intensity of 0.08%. A structured O grid with 200 nodes perpendicular to the airfoil surfaces and 500 nodes wrapped over each surface was applied to the two airfoils.

In the case of E387 as shown in Figure 4(a), a structured C-grid with 200 nodes perpendicular to the airfoil surfaces and 400 to 500 nodes wrapped over each surface was applied to both airfoils. Because of the circular TE of A7, an O-domain with a diameter of three chord lengths was used in the C-grid to wrap the airfoil, as presented in Figure 4(b). For different cases of AOA and Reynolds numbers, the LSB region has been refined with a minimum of 150 uniformly distributed points to capture the bubble characteristics in detail. To compare with McGhee et al.'s wind tunnel experiment, ²⁷ three of the Reynolds numbers (1, 2, and $3\times 105$ ) carried out in the experiment were selected for the RANS computation. The corresponding incoming turbulence intensities were 0.16%, 0.16%, and 0.2%, respectively.

Figure 4.

The schematic figure of the mesh domains of (a) E387 and (b) A7.

In both cases, the mesh was refined with twice the resolution in the streamwise and wall-normal directions from a coarse mesh, and it was determined to be an independent mesh when lift coefficient predictions showed negligible difference from a finer mesh but significant difference from a coarser mesh. The computational domain extends to 20 chord lengths in all directions for each airfoil, with a total number of approximately 200,000 grid cells and y+ < 1.0. The transition SST turbulence model which is based on the coupling of the k − ω SST transport equations with two other transport equations (in terms of the intermittency and momentum-thickness Reynolds number, respectively) has been selected as the eddy viscosity model in ANSYS Fluent 16.0.

Redesign of the LE of NACA0012

As one of the most common airfoils, NACA0012 is widely studied and applied.^34–36 The thickness distribution equation of NACA0012 is:

$$y/c=\pm 0.6[0.2926\sqrt{x/c}-0.126x/c\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}-0.3516(x/c)2+0.2843(x/c)3\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}-0.1015(x/c)4]$$

(8)

The second-order derivative of this equation is infinite at the LE point and the surface curvature at this point is singular by the definition in equation (1). This is a typical curvature discontinuity, suitable for application of the CIRCLE method. Except the LE point, no curvature or slope-of-curvature discontinuity exists in the NACA0012 profile, because the second- and third-order derivatives of equation (8) are continuous. Consequently, only small improvement in aerodynamic performance can be achieved by fixing only one discontinuous point on the airfoil and keeping a nearly identical profile simultaneously. This case is used to illustrate the effects of curvature discontinuity on aerodynamic performance.

The design resulting from application of the CIRCLE method is shown in Figure 5, and this optimized airfoil is denoted as “QM13F.” The newly designed airfoil profile deviates minimally from the original NACA0012, and some key parameters are shown in Table 1.

Figure 5.

The comparison between (a) profiles and (b) surface curvature distributions of original airfoil NACA0012 and optimized airfoil QM13F. The small figure in (b) is a magnification in the LE area.

Table 1.

Some key design parameters in the CIRCLE method.

TE radius	βs2 (°)	Psm	βsm (°)	βs1 (°)	LE radius	LE position
0.00145	79.5	(0.5, 0.0529)	−3.75	38	0.02	(0, 0)

QM13F is also a symmetrical airfoil and has almost the same thickness distribution as the original airfoil. Since the CIRCLE method has been applied, the newly designed airfoil has a continuous surface and slope-of-surface distribution from TE to LE without the singularity at the LE point. The curvature distributions of the two airfoils are shown in Figure 5(b).

The lift coefficients (Cl) and drag coefficients (Cd) of the airfoils and experimental data of NACA0012 at Reynolds number 1.35 × 105 are presented in Figure 6. The RANS results of NACA0012 match well with the experimental data by Lee and Gerontakos. ²⁸ The computation accurately predicts the stalling AoA and the corresponding Cl and Cd, although it slightly overestimates Cl values from 7 ° to 10 ° in Figure 6(a) and overestimates Cd values when the AoA is low, as shown in Figure 6(b). Both of these figures indicate that the two airfoils have the same lift and drag performance at this Reynolds number until the AoA reaches 10 °. At the AoA > 10 °, however, QM13F has a larger Cl and a reduced Cd compared with NACA0012. This numerical comparison between NACA0012 and QM13F indicates that the LE singularity of NACA0012 significantly improves the airfoils aerodynamic performance near the stalling AoA.

Figure 6.

Comparison between (a) lift coefficients and (b) drag coefficients from RANS results and experimental work by Lee and Gerontakos. ²⁸

As shown in Figure 6, the main differences between the lift and drag performance of the two airfoils are shown at the AoA of approximately 13 °. Therefore, the pressure coefficient (Cp) and skin friction coefficient (Cf ) distributions of these airfoils at 13 ° are examined to further investigate the boundary layer behavior for different aerodynamic performance, as shown in Figure 7.

Figure 7.

The RANS results comparison of (a) pressure coefficients and (b) upper surface skin friction coefficients at 13 ° angle of attack between NACA0012 and QM13F.

The difference in Cp distribution can be observed on the suction side near the LE region. In the magnified view in Figure 7(a), the RANS results of QM13F shows increased negative pressure over the LE suction side compared with the original NACA0012 due to the continuous surface curvature distribution, which results in better lift performance. The differences between the Cf distributions on the suction side are presented in Figure 7(b). This indicates that QM13F has a delayed separation point near the LE region because the Cf value of QM13F before the transition point is generally larger than that of NACA0012. QM13F provides reduced skin friction force without the LE curvature discontinuity as shown in Figure 7(b).

Both TATO and RANS methods are applied to NACA 0012 and QM13F at AoA = 13 °, Re = 1.35 × 10⁵ and the predictions are provided in Table 2. In the experimental work by Lee and Gerontakos, ²⁸ the LSB on NACA 0012 at the same operating conditions appeared near the LE and the length is approximately 0.07 chord length (the bubble is too small to give an accurate length). Both numerical methods reproduced the bubble position and acceptably predicted the bubble length. TATO and RANS predicted that the length of the LSB on QM13F is 0.07 and 0.016 chord-length shorter than that on NACA 0012, resulting in better aerodynamic performance as shown in Figure 6.

Table 2.

The bubble length comparison between TATO and RANS for NACA 0012 and QM13F at Re = 1.35 × 10⁵, AoA = 13 °.

	TATO			RANS
	LS position	TR position	BL	LS position	TR position	BL
NACA0012	0.018	0.079	0.059	0.016	0.099	0.083
QM13F	0.024	0.076	0.052	0.025	0.092	0.067

LS: laminar separation; TR: turbulent re-attachment; BL: bubble length.

It can be concluded that the removal of the LE surface curvature discontinuity results in higher lift and lower drag near stalling AoA due to the higher pressure differences between pressure and suction sides and delayed laminar separation.

Redesign of the airfoil E387

Like many other airfoils, the original geometry of airfoil E387 exhibits slope-of-curvature discontinuities on its surface, observable at the region where the LE circle joins the airfoil surfaces among other places. The airfoil E387 is optimized by the CIRCLE method. The resulting airfoil A7 is presented and compared with the original airfoil. Figure 8(a) shows a comparison between the two airfoil geometries. A7 has a circular TE and E387 has a sharp TE, but the almost the same LE circle is present in both airfoils. A7 is slightly thicker than E387 from LE to TE, and A7 has mildly greater camber from LE to approximately 0.65 chord length and smaller camber for the rest of the chord length than E387. As indicated in Figure 8(b), E387 has slope-of-curvature discontinuities mainly on the suction surface, and an obvious “kink” is observable on each side of the airfoil E387. These slope-of-curvature discontinuities have been removed in the redesigned airfoil A7.

Figure 8.

The comparison between (a) airfoil profiles and (b) curvature distributions of E387 and A7.

Lift and drag coefficients

It is a very straight forward way to compare airfoils performance by comparing their lift and drag coefficients. In this subsection, the simulated lift and drag coefficients of E387 and A7 are compared with show the aerodynamic performance improvement caused by removing slope-of-curvature discontinuities at Reynolds number 10⁵, 2 × 10⁵, and 3 × 10⁵. The experiment carried out in the NASA Langley LTPT ²⁷ is usually regarded as the benchmark to validate results of E387. This experiment has consequently been selected as the reference to validate the results from this article.

The Cl-Cd polar curves at three Reynolds numbers are shown in Figure 9(a) to (c). At all Reynolds numbers, the results of E387 RANS computation match well with the experimental data, and meanwhile the A7 RANS results shows improved Cl and Cd at all the AoA investigated as shown in Figure 9. These results also indicate that the higher the AoA, the better performance airfoil A7 can offer compared with E387.

Figure 9.

Comparison between lift-drag coefficients polar curves from RANS results and experimental work of E387 by McGhee et al. ²⁷ at (a) Re = 1 × 10⁵, (b) Re = 2 × 10⁵, and (c) Re = 3 × 10⁵. The angles of attack of the experimental data range from 0 ° to 10 °.

It is worth to mention that the RANS predictions of Cl for E387 at Re = 10⁵ are smaller than the experimental data when the AoA exceeds 8 °, around which E387 experiences earlier transition to turbulence, while the improved A7 airfoil maintains the increasing trend of the Cl, as shown in Figure 9(a). The earlier laminar-turbulence transition of E387 causes a loss of lift as well as increased drag. The airfoil A7 has reduced Cd practically throughout the range of AoA.

Investigation of LSB

LSB dominates airfoil aerodynamic performance at low Reynolds numbers. ²⁷ The details of the separation bubbles on both E387 and A7 are presented in this subsection to investigate the reasons for the improvements in aerodynamic performance. Figure 10(a) and (b) summarizes the separation points and re-attachment points of LSB at Reynolds number 1 × 10⁵ and 2 × 10⁵, respectively.

Figure 10.

Separation and re-attachment points at Re = 1 × 10⁵ (a) and Re = 2 × 10⁵ (b).

At Re = 10⁵, TATO results provide good predictions of laminar separation positions while the RANS predictions starts to move rapidly toward the LE when the AoA is above 6 °. This LE toward moving phenomenon causes larger separation bubble, resulting in lower Cl and higher Cd results than experimental data as previously shown in Figure 9(a). Meanwhile, RANS computations offer better predictions of turbulent re-attachment positions compared with TATO at Re = 10⁵. At Re = 2 × 10⁵, both RANS and TATO methods provide acceptable results of laminar separation positions of E387 compared with the experimental data, while TATO slightly overestimates the re-attachment points when the AoA is low.

To analyze this table, it is essential to understand the growth and bursting process of LSB on E387 with increasing AoA at low Reynolds number as shown in Figure 11. In Horton's theory, ³² the process can be summarized as three stages. At a constant low Reynolds number, the LSB presents on the airfoil suction side at low angles of attack (Stage 1), and it moves toward the LE region and gets smaller with increasing AoA. When the AoA increases to a moderate value (e.g. E387 at 8 ° and Re = 2 × 10^{5 27}), the laminar-turbulent transition occurs earlier than laminar separation (x _sep  > x _tr in Figure 1), and the airfoil experiences natural transition rather than LSB. When the AoA is higher than these transition angles, the LSB appears again at the LE region as a “short bubble” (Stage 2). The LSB bursts at the stalling AoA resulting in LE stall (Stage 3).

Figure 11.

The growing and bursting process of LSB on an airfoil with increasing AoA at low Reynolds numbers. The predicted bubble lengths at Reynolds numbers 10⁵ and 2 × 10⁵ are summarized in Table 3. The prediction accuracy and computation time are also included in the table. At both Reynolds numbers in the table, the bubble length summarized from the experimental report ²⁷ are acceptably reproduced in all cases except the 7 ° AoA. The TATO predictions provide more accurate results at Re = 1 × 10⁵ than at Re = 2 × 10⁵.

At the AoA of 7 ° in Table 3, the bubble lengths of E387 predicted by both RANS and TATO are more than 20% longer than those observed in the experiment. There are two possible explanations for this phenomenon. First, the sensitivity of the LSB at 7 ° has to be considered. The 7 ° AoA is close to the natural transition AoA as explained in the previous paragraph and the separation bubble appears to be sensitive and unstable. The sensitivity problem of the LSB at AoA around 7 ° has been observed by other researchers. Lin and Pauley ³⁷ found that the freestream turbulence in the experiment may delay the separation and hence make a shorter bubble. They also found that the RANS computation overestimates the bubble length at the AoA of around 7 °. Moreover, the inaccuracy of TATO calculations are mainly from the turbulent re-attachment predictions, which are calculated from equation (7). One possible reason for the delayed re-attachment predictions is that the classical re-attachment criterion does not take turbulence intensity of the experiment into consideration, which plays a key role in the formation of the LSB, especially the location of the re-attachment point. Also, equation (7) is an empirical formula which is based on other experimental results, and it does not take the natural transition AOA of E387 or A7 into account. Consequently, longer bubbles at 7 ° are estimated by both computations.

Table 3.

Bubble length prediction accuracy and computation time comparison between RANS and TATO.

Airfoil		E387	A7	E387	A7	E387	A7	E387	A7	E387	A7	E387	A7
Re × 105		1	1	1	1	1	1	2	2	2	2	2	2
AoA (°)		0	0	4	4	7	7	0	0	4	4	7	7
Bubble length (c)	Experiment	0.43	N/A	0.38	N/A	0.25	N/A	0.26	N/A	0.22	N/A	0.15	N/A
	TATO	0.444	0.419	0.401	0.388	0.316	0.332	0.316	0.305	0.254	0.242	0.188	0.190
	RANS	0.395	0.380	0.359	0.321	0.302	0.313	0.281	0.271	0.230	0.225	0.205	0.210
Error (%)	TATO	3.25	N/A	5.52	N/A	26.4	N/A	21.5	N/A	15.5	N/A	25.3	N/A
	RANS	−8.14	N/A	−5.53	N/A	20.8	N/A	8.08	N/A	4.55	N/A	36.67	N/A
Computation time (s)	TATO (1 core)	3.6	3.5	4.1	4.9	3.8	4.0	4.2	4.8	3.9	3.9	4.2	3.9
	RANS (8 cores)	2305	1735	1669	1352	1169	1563	1712	1533	2428	1666	1093	1111

The experimental data are from Mcghee et al. ²⁷

The predicted bubble lengths of A7 are generally smaller than those of E387 except at the AoA of 7 ° at both Reynolds numbers. This can be explained by the growing and bursting process of an LSB in Figure 11. The bubble sizes of E387 have been efficiently reduced at 0 ° and 4 ° by removing the slope-of-curvature discontinuities. The redesigned airfoil A7 hence has better aerodynamic performance at low AoAs as shown in Figure 9. Compared with the LSBs on A7, the smaller bubbles of E387 at 7 ° indicate that the contraction speed of these bubbles become faster at 7 ° than at lower AoAs. They start moving rapidly toward the LE and getting into Stage 2 with increasing AoAs. However, as the AOA increases, the LSBs on the airfoil A7 contracts at the usual rate and stays in Stage 1. This phenomenon is more obvious at higher AoAs, which will be presented in the next subsection.

A comparison of the two approaches in terms of accuracy and computation time is presented in Table 3. Both methods offer acceptable results compared with the experimental data and RANS gives slightly better accuracy on bubble length prediction. However, the TATO method offers a much quicker response (within 5 s) and is hundreds of times quicker than RANS even if eight cores MPI are applied to the RANS computations. This makes TATO an affordable computing tool for rapid estimation of the LSB on an airfoil.

Table 3 also indicates that when the Reynolds number reaches 2 × 105, the curvature effects become less important for reducing the bubble size because the separation bubble is no longer the key factor affecting the airfoil performance. More specific details of the lift and drag force components are shown in Table 4.

Table 4.

The details of Cl and Cd components for E387 and A7 at Re = 3 × 10⁵, AoA = 4 °.

	Cl			Cd
	Pressure	Viscous	Total	Pressure	Viscous	Total
E387	8.037E-01	4.419E-04	8.042E-01	5.628E-03	5.699E-03	1.133E-02
A7	8.091E-01	5.162E-04	8.096E-01	4.994E-03	5.767E-03	1.076E-02

While Cl is mostly dependent on the pressure forces acting on the airfoil as expected, there is also some effect due to the pressure and skin friction on the drag coefficient Cd. The continuous slope-of-curvature distribution not only increases the pressure component, which is the main source of the total Cl, but also increases the viscous component by 16.81%. The removal of the surface slope-of-curvature discontinuity from the upper surface leads to an obvious rise in pressure difference between the two surfaces and hence leads to an increased Cl. The continuous slope-of-curvature distribution also results in an 11.27% reduction of the pressure Cd and a 5.03% reduction of total Cd.

Cases of higher angles of attack

In “Lift and drag coefficients” section, RANS results indicated that the higher the AoA, the better performance the airfoil A7 can offer compared with E387. This subsection explains these observations by investigating the skin friction coefficient distributions and streamlines between E387 and A7. These are shown in Figure 12.

Figure 12.

Comparison of skin friction coefficient distributions between E387 and A7 at (a) Re = 3 × 10⁵, AoA = 8 °, (b) Re = 1 × 10⁵, AoA = 10 °. The streamline figures are the streamline of E387 (top) and A7 (bottom), respectively, for the corresponding Re and AoA.

In Figure 12(a), E387 exhibits a short bubble at LE, making the airfoil experience earlier transition to turbulence at Re = 3 × 10⁵ and AoA = 8 °. With a turbulent re-attachment at the position x/c = 0.051, as shown in the streamline figure, the boundary layer becomes more energetic and hence significantly increases the skin friction force compared with airfoil A7. The continuous slope-of-curvature distribution results in a lower Cf distribution for the redesigned airfoil A7 with a LSB, which does not dominate the aerodynamic performance at Re = 3 × 10⁵.

As shown in Figure 11, the skin friction force generally increases from Stage 1 to 2 due to the more energetic turbulent boundary layer, and a bursting increase from Stage 2 to 3 due to the LE stall. This explains the difference in the LSB size in the case shown in Figure 12 and the case described previously in Table 3 at 7 °. The RANS results indicate that discontinuous slope-of-curvature distributions make the LSB on E387 move into Stage 2 while the bubble on the airfoil A7 with an improved curvature distribution is still in Stage 1. This leads to an earlier laminar and turbulent separation, as well as an earlier transition process, and therefore higher Cf. The separation bubble on the airfoil A7 is caused only by the unavoidable adverse pressure gradient with no contribution from curvature discontinuity. Consequently, A7 has increased magnitudes of Cl and Cd performance (as the LSB is still in Stage 1) comparing with an earlier laminar and turbulent separation on the airfoil E387.

Both TATO and RANS methods are applied to Figure 12(b) to calculate bubble lengths and compared in Table 5. In the TATO method, equations (5) and (6) are used to calculate laminar separation position x _sep and transition position x _tr based on the inviscid external velocity from thin airfoil theory. As shown in Figure 1, equation (7) is used to calculate the re-attachment point when x _sep  < x _tr . In the RANS method, the separation streamline and the airfoil surface are used to bound the bubble, and the coordinates from the separation point to the re-attachment point are obtained. Both methods predicted a short bubble near the LE of E387 and a longer bubble on the main body of A7. The predictions are consistent with Horton's theory and the authors' inferences in last subsection. The continuous slope-of-curvature distribution delays the LSB's transition to move into Stage 2 and helps to control the LSB contraction speed when the AoA increases.

Table 5.

The bubble length comparison between TATO and RANS at Re = 1 × 10⁵, AoA = 10 °.

		TATO			RANS
		LS position	TR position	BL	LS position	TR position	BL
Re = 1 × 105 AoA = 10 °	E387	0.020	0.133	0.113	0.008	0.101	0.093
	A7	0.051	0.245	0.194	0.039	0.211	0.172

LS: laminar separation; TR: turbulent re-attachment; BL: bubble length.

Figure 12(b) also shows at Re = 1 × 10⁵ and AoA = 10 °, the continuous slope-of-curvature distribution delayed the position of turbulence TE separation by more than 10% for the redesigned airfoil A7. As soon as the flow separates near the TE, the separated region merges with the wake and results in TE stall, causing lift loss and drag increase. Continuous slope-of-curvature distributions make the boundary layer flow attach to the airfoil surface for a longer surface distance and suppress the flow separation in the TE, and consequently improve airfoil performance.

Concluding remarks

This article presented a completed solution for profile design and aerodynamic performance computation methods for the airfoils with discontinuities in their curvature and slope-of-curvature distributions. A symmetrical airfoil NACA0012 and a non-symmetrical airfoil E387 were judiciously selected as examples of airfoils with surface curvature discontinuity and slope-of-curvature discontinuity, respectively. They were redesigned using the CIRCLE method to specifically investigate the effects of surface curvature on aerodynamic performance at different angles of attack and Reynolds numbers.

NACA0012 exhibits a surface curvature discontinuity at the LE stagnation point. The newly designed airfoil removed the LE singularity. RANS computational results indicated that the airfoil surface curvature discontinuity affects pressure and skin friction distribution, and hence affects lift and drag coefficients near the stalling AOA. The airfoil performance near the stalling angle is improved (an average 5% increase in lift and average 10% reduce in drag) by removing the curvature discontinuity.

E387 exhibits slope-of-curvature discontinuities in both surfaces. The newly designed airfoil deviates little in its profile geometry to the original airfoil but with no slope-of-curvature discontinuities. In this case, a rapid-response inviscid–viscid interaction scheme called TATO was introduced to compute the aerodynamic performance. The results are acceptable compared with the results from the RANS computations and experimental work. The TATO scheme provided much faster response than the RANS computation. The flow performance was analyzed for three different Reynolds numbers and it was shown that continuous slope-of-curvature distributions increased the lift coefficient and the reduced drag coefficient for all the Reynolds numbers studied. The separation bubble sizes were computed by RANS results, and compared well against available experimental data in the literature and an inviscid–viscid scheme of the second-order thin airfoil theory and boundary layer methods. For the cases at Reynolds numbers below 2 × 10⁵ and relatively low angles of attack, the airfoil surface curvature affects the size of the LSB which dominates airfoil performance in these cases. The bubble sizes were successfully reduced by removing slope-of-curvature discontinuities, and the aerodynamic performance is consequently improved. For the cases at Re > 2 × 10⁵, the factor that is important for the airfoil performance is no longer the separation bubble size but the inherent profile of the airfoil. At higher angles of attack for all cases, continuous slope-of-curvature distribution delayed LE separation, laminar-turbulence transition, and TE separation. Through these mechanisms, skin friction is reduced and improved aerodynamic performance is achieved.

It has been shown that airfoil surface curvature and slope-of-curvature distributions have significant effects on the boundary layer behavior. Continuous curvature and slope-of-curvature distributions provide better control of flow separation, inherently better profile drag and thereby improve Cl/Cd, leading to improving airfoil aerodynamic performance.