AERODYNAMIC CHARACTERISTICS OF LOW REYNOLDS NUMBER AIRFOILS

Original scientific paper Assessment of airfoil aerodynamic characteristics is essential part of any optimal airfoil design procedure. This paper illustrates rapid and efficient method for determination of aerodynamic characteristics of an airfoil, which is based on viscous-inviscid interaction. Inviscid flow is solved by conformal mapping, while viscous effects are determined by solving integral boundary layer equations. Displacement thickness is iteratively added to the airfoil contour by alternating inviscid and viscous solutions. With this approach efficient method is developed for airfoil design by shape perturbations. The procedure is implemented in computer code, and calculation results are compared with results of XFOIL calculations and with experiment. Eppler E387 low Reynolds number airfoil and soft stall S8036 airfoil are used for verification of developed procedure for Reynolds numbers 200000, 350000, and 500000. Calculated drag polars are presented in this paper and good agreement with experiment is achieved as long as small separation is maintained. Calculated positions of laminar separation, reattachment, and turbulent separation closely follow experimental measurement. The calculations are performed in relatively short time, which makes this approach suitable for low Reynolds number airfoil design.


Introduction
This paper presents an efficient and reliable approach for assessment of aerodynamic characteristics of low Reynolds number airfoils.This approach is used for design of optimum airfoils by systematic modification of airfoil shape.Since large number of airfoil aerodynamic characteristics computations are required by this optimization approach it is necessary to have quick and reliable determination of aerodynamic parameters without human intervention.The used method for determination of airfoil aerodynamic characteristics combines incompressible potential flow solution about airfoil and boundary layer correction for assessing airfoil drag and maximum lift coefficient.Comparison of theoretical methods for predicting airfoil aerodynamic characteristics performed by Maughmer [1] and Coder [2] showed that the airfoil drag calculated by computer codes incorporating boundary layer methods generally agrees better with experimental results than did Navier-Stokes solvers, while the maximum lift coefficient is frequently over predicted by all methods.This is also supported by Morgado [3].Furthermore, potential flow coupled with boundary layer equations solvers requires about 100 times less variables compared to N-S solvers without necessity to regenerate computational grid, with runtime enormously faster for similar accuracy as argued by Drela [4].This fact makes inviscid-viscous approach more suitable for airfoil optimizations with shape perturbation where large number of configurations needs to be analysed.
Eppler's [5,6] and Drela's [7] XFOIL codes, however, use inverse airfoil design approach for which a target pressure distribution about airfoil is specified and corresponding airfoil shape is determined by combining inviscid flow solutions and boundary layer viscous correction.Inverse airfoil design approach is successfully used to design new airfoils for variety of applications [8,9,10].Human intervention is necessary by this approach since specified pressure distribution often results in open or intersected airfoil shapes.Airfoil design by shape perturbation, on the other hand, is widely used as airfoil design methodology since calculations always start and end with realistic airfoil shapes.Lissaman and Carmichael [11,12] pointed out that airfoils at low Reynolds number have poor lift to drag ratio mainly due to the presence of laminar separation bubble.Based on this fact, low Reynolds number airfoil design methods should be capable to reliably estimate location of laminar separation, reattachment and transition.Reference [13] provides experimental measurements of upper surface boundary layer features of selected airfoils at Reynolds numbers of 200 000, 350 000, and 500000.These data provide suitable source for validation of aerodynamic computations.Experimental data from reference [13,14] are used in the validation of the presented results.
Satisfactory and fast boundary layer calculations still rely on correlations from experimental measurements.For proper determination of airfoil aerodynamic characteristics laminar separation location, reattachment, and transition locations, have to be carefully determined.
For high Reynolds numbers, flow is essentially turbulent and is more resistant to separation, so separation occurs at the rear part of the airfoil, contrary to laminar flows, as pointed out in references [11,12], where transition from laminar to turbulent flow normally occurs near the minimum pressure point, at the first onset of the adverse pressure gradient.If initial pressure gradient is too high, laminar separation occurs before transition to turbulent flow took place.If the separated flow reattaches a region of separated flow, it is called laminar separation bubble (LSB).This region is characterized by zero skin friction coefficient [8,16].LSB extends over the airfoil surface between the point of laminar separation to the point of reattachment.As discussed in [17] In this paper a conformal mapping to unit circle, by applying Fast Fourier Transformation (FFT) to determine mapping coefficients, is combined with integral boundary layer solution in direct way, where airfoil shape is specified and velocity distribution is found contrary to inverse airfoil design formulation used by Eppler, because it is more suitable for optimization by shape perturbations.Inverse design approach is still possible by our method when suitably defined goal function in the form of velocity or pressure distribution is specified, resulting into airfoil shape which is always realistic.
Transition from laminar to turbulent becomes more important as Reynolds number is decreased because transition location determines LSB length and hence airfoil drag.Thus, when LSB appears a transition criterion based on e  method is utilized instead of Eppler modified criteria which is used to predict natural transition.Laminar separation model is incorporated in the code so drag increase due to bubble is computed without need for user intervention.As a result, the goal of the current work is to validate this computational procedure.This goal is accomplished by comparison with XFOIL results and published experimental measurements in low Reynolds number regime that is important for UAVs and wind turbines.Airfoils selected are typical low Reynolds number airfoils with drag polar having a laminar bucket.The computational procedure is illustrated in Fig. 1.Solution starts with finding inviscid pressure distribution for given aerodynamic shape, angle of attack, and Reynolds number based on conformal mapping technique.This pressure distribution is then used to calculate boundary layer development by integral method.Boundary layer displacement thickness is then added to the initial airfoil coordinates, and solution is repeated until there is no significant change in displacement thickness.Drag is calculated utilizing modified Young-Squares formula applied at the trailing edge.Viscous corrections are applied to lift and pitching moment by modifying the effective angle of attack.

Conformal mapping
Mapping technique is used to generate exact solutions for potential flow problems efficiently.It is also widely used in inverse design problems [7,18].Development here closely follows [22] and [23].There are three basic steps involving transforming an airfoil into a circle (Fig. 2).Firstly, a given airfoil coordinates in physical domain z is transformed to near circle shape via Karman-Trefftz inverse transformation using Eq.(1).
where z 0 is defined midway between leading edge and its center of curvature, and z1 is trailing edge singularity point, as shown in Fig. 2a where  = 2 − /π, and τ is airfoil trailing edge angle.Near circle coordinates  1 can be obtained from Eq. ( 1) as Eq. ( 2a) is valid when near the trailing edge upper surface is above real axis and lower surface is below it.For the points of the lower surface which are above real axis Eq. (2b) is used: The second step is to translate the near circle to the centroid of the coordinate system, using Eq. ( 3), Fig. 2c.where . After taking logarithms of both sides, and equating real and imaginary parts the following two equations are obtained.
The circle radius is given by . The coefficients   , and b n can be determined in iterative manner using Fast Fourier Transform.The procedure starts by choosing 2N equally spaced points on the true circle beginning from trailing edge point image.
If the value of the near circle angle at trailing edge is substituted in Eq. ( 7), an expression for b 0 is obtained as given by Eq. (8).

Determination of airfoil velocity distribution
The conjugate complex velocity V in z-airfoil plane is obtained by the equation where V is conjugate velocity in airfoil plane, d/d 3 is conjugate complex velocity around circle, and the rest of the terms are derivatives of the transformations, as given by Eqs.(10÷12).
Conjugate complex velocity distribution around circle centred at coordinate origin is given by Eq. (13).According to Kutta condition velocity at airfoil trailing edge must be finite.Since the derivative d 1 /d at trailing edge is infinite while other mapping derivatives are finite, to achieve finite velocity at airfoil trailing edge velocity in circle plane which corresponds to airfoil trailing edge (r=R,  = ) must be equal to zero, as given by Eq. ( 14): e e e 0 2 (14) From this equation circulation is determined as: Circulation Γ necessary to satisfy Kutta condition thus depends on angle of attack α and position of trailing edge image in circle plane μ.The Pressure distribution is obtained by Eq. ( 16).
where d/d is given by Eq. ( 9).Figs. 3 to 5 show derivatives of transformations for E387 airfoil which constitute components of Eq. ( 9).Fig. 6 shows corresponding inviscid pressure distribution at an angle of attack of 2°.Calculated pressure distribution obtained by XFOIL code is also shown.The pressure distributions match favourably.

Boundary layer computations
Once the pressure distribution is obtained using the previously described method, the boundary layer equations are then separately solved along upper and lower sides of the airfoil.Boundary layer thicknesses δ 2 and δ 3 are determined by solving integral momentum and energy boundary layer equations Eq. ( 17), and Eq. ( 18): Following Eppler [5] and [6], the closure correlations used to solve the above system of equations are given for laminar and turbulent flows as function of the shape factor H 32 and local Reynolds number based on boundary layer momentum thickness Re δ2 .The initial values for δ 2 and δ 3 are calculated using first step solution starting from stagnation point as discussed in [5].
Boundary layer features are specified in terms of the shape factor H 32 .For instance, Laminar separation is indicated at H 32 value of 1,51509, while turbulent separation is reached when H 32 =1,46.When Eppler modified transition criterion Eq. ( 20) is satisfied [6], switching from laminar to turbulent closure correlations is triggered.
This approximate criterion based on local values of the shape factor H 32 , surface roughness r, and Re δ2 is empirically derived from fitting through wind tunnel data and flight test results [20].When bubble appears Drela method is used as transition criterion.This method assumes transition when an amplification factor for Tollmien-Schlichting waves in the boundary layer has grown to e 9 19.
Airfoil shape is modified by adding calculated upper boundary layer displacement thicknesses δ 1ui to the upper airfoil coordinates ui y and subtracting displacement thickness for the lower side as given by Eq. ( 21) and Eq. ( 22) ui ui ui y y Upper surface boundary layer development charts for low Reynolds number airfoil E387 at an angle of attack of 1,5° is shown in Fig. 7, where logarithm of Re δ2 is plotted against the boundary layer shape factor H 32 .Solution of boundary layer equations starts at the stagnation point with H 32 =1,62, which decreases toward laminar separation limit criterion is satisfied, and thus laminar separation bubble is formed and bubble computations are evoked.The bubble due to reattachment.It is worth to note that H 32 decreases rapidly and extends on the airfoil surface from laminar separation point in the region just before laminar separation as seen in Fig. 7.This behaviour influences solution of Eqs. ( 17) and ( 18) in laminar part of the bubble.Gaster [21] studied laminar separation bubbles experimentally and pointed out the importance of pressure gradient parameter.Horton [22] proposed a semiempirical bubble model with constant velocity profile in laminar part of the bubble, followed by linear decrease in velocity from transition to reattachment.LSB for low Reynolds number airfoils is studied by O'Meara [23].
Dini [24] has developed semi-empirical model intended for airfoils at low Reynolds numbers.Dini model uses improved velocity plateau function in laminar bubble length  1 Eq.( 23), making use of Gaster pressure gradient parameter, which depends on local boundary layer parameters at separation point and average velocity gradient along the bubble.
where the subscript s refers to separation conditions and the derivative is with respect to surface distance s.DU is a quadratic function of Gaster pressure parameter.
Transition is calculated by e n method modified by Drela [19].In the turbulent part of the bubble however, a shape factor distribution for H 32 was developed to derive the solution of the boundary layer integral equations in inverse mode from transition up to reattachment [24].Reattachment location, and thus turbulent bubble length l 2 is then calculated using correlations for spreading angle of turbulent shear layer and separation angle γ as given by Eq. (24).Using this method separation point location, transition, and reattachment points can be calculated.Following Drela, turbulent separation after bubble reattachment is assumed when Eq.( 25) is satisfied.Lift and moment coefficients are obtained from potential flow theory.A correction is then applied to take care of separation effects Eq. ( 26).Drag coefficient is calculated using modified Squire-Young formula applied at the trailing edge Eq. ( 27).This formula is applied at trailing edge for upper and lower surfaces separately [5], and was found to have good agreement with experimental measurement.

R , , l l
S sep is airfoil surface distance for which flow is separated, and δ us is slope of airfoil upper surface at trailing edge.

Results and validation
Current calculation results are compared with published experimental data for two selected airfoils.Eppler low Reynolds number airfoil E387 is used as benchmark for validating low Reynolds number aerodynamic computations.It is extensively tested in NASA Langley Low Turbulence Pressure Tunnel (LTPT, where drag polar, and pressure measurements at low Reynolds numbers are published [13].Recently, E387 airfoil is tested in the University of Illinois at Urbana-Champaign (UIUC) subsonic wind tunnel [14,15], which is intended to validate and refine airfoil low Reynolds number computation methods.The second airfoil is Selig S8036 low Reynolds number airfoil designed for soft stall characteristics.Experimental measurement data for these two airfoils at flow Reynolds numbers are 200.000,350,000 and 500,000 are used in the validation of current computations.These measurements include drag polar and location of upper surface boundary layer flow features.
Figs. 8 to 10 show comparisons of measured [13] and calculated pressure distributions over E387 airfoil at Reynolds number of 300 000 and at angles of attack of 2, 4, and 6 degrees.The location of the separation bubble is clearly observed on the upper surface.Calculated pressure distribution agrees with experimental data and XFOIL results.
The bubble location is calculated with acceptable accuracy for optimization computations.The general observation is that the bubble moves upstream as angle of attack increases, with length being shorter.Fig. 11 shows comparisons of locations of upper surface features of the two airfoils at different angles of attack and Reynolds numbers of 200.000, 350.000, and 500.000.The computed laminar separation, Reattachment, and turbulent separation locations on upper surface are compared to experimental measurements.Laminar separation and reattachment locations from XFOIL are also shown for E387 at Reynolds number 350,000.A laminar separation bubble extends on the upper surface starting approximately at mid chord.As angle of attack increases the bubble moves toward the leading edge, and its length decreases.When the bubble length close to leading edge is very short, it could be interpreted as a transition without bubble.Current calculations follow the general trend of both experimental measurement and XFOIL predictions.As Reynolds number increases the laminar separation bubble tends to shorten in length, which is in agreement with the general fact that laminar separation bubble is more dominant in low Reynolds number range.The results of XFOIL and current calculations seem to underestimate the reattachment point location, this is also noted in [14].For low angles of attack turbulent separation takes place at or very close to the trailing edge.When angle of attack further increases, turbulent separation moves forward causing high increase in drag and loss in lift.In all cases turbulent separation point assessed by current computations and XFOIL code at high angles of attack is more aft than the measured locations.This misspredictions have the consequence of over estimating the angle of maximum lift, and thus the value of maximum lift coefficient.It also limits the method capability to predict lift and drag at high lift coefficients, which is clearly seen from drag polar curves shown in Figs. 12 to 17. Calculation results closely follow typical laminar airfoil aerodynamic characteristics.At low angles of attack the agreement with experimental data is noticeably close.When angle of attack increases further turbulent separation point moves away from trailing edge toward the leading edge limiting method accuracy.

Figure 12
Comparisons between calculated and experimental drag polar for E387 at Re=200,000.

Figure 13
Comparisons between calculated and experimental drag polar for E387 at Re=350,000.

Figure 14
Comparisons between calculated and experimental drag polar for E387 at Re=500,000.

Figure 15
Comparisons between calculated and experimental drag polar for S8036 at Re=200,000.

Figure 16
Comparisons between calculated and experimental drag polar for S8036 at Re=350,000.

Objective functions
Objective function for airfoil design may vary from one application to another.Optimization algorithm manipulates airfoil shape parameters in systematic manner to satisfy the objective.For instance, maximizing range can be formulated in terms of maximizing lift to drag ratio at specified range of angles of attack α and Reynolds numbers as given by Eq. (28).Minimizing drag or optimizing airfoils for specific pressure distribution can be achieved by similar function Eq. (29) and Eq.(30).Objective function can also be combined to fulfil multiple objectives, as illustrated by Eq. (31) where the factor   allows different weights being given to each component at the i th point.

Conclusion
This paper has illustrated a systematic approach to the assessment of aerodynamic characteristics of flow over airfoils at low Reynolds numbers.The computational procedure starts with mapping a given airfoil shape into a true circle in three subsequent steps.Multiplication of derivatives of these transformations with velocity distribution around circle results in the inviscid velocity distribution at a specified angle of attack.Circulation is fixed by applying Kutta condition at trailing edge of the true circle.Boundary layer integral equations solution enables the assessment of lift viscous corrections, total drag, and laminar separation bubble location.The calculation procedure is repeated by adding boundary layer displacement thickness, until change in airfoil shape is negligibly small.This requires only few iterations, making this approach very efficient for airfoil design by systematic airfoil perturbation.
In order to validate this method of calculation a comparison with experimental data for E387 and S8036 airfoils is performed.Curves of boundary layer flow features on upper surface and drag polar show satisfactory agreement with measurement.In low angle of attack range where airfoil optimization is expected, both lift and drag are computed with reasonable accuracy.Separation bubble location can be also assessed in consistence with measurements.Weak laminar separation bubble is not captured by this procedure, however this weak bubble often causes small drag penalties and can be tolerated.When angle of attack is high, and when turbulent separation occurs on the upper surface, maximum lift coefficient is overestimated.Bubble length predicted by current computations is shorter than that obtained from experimental measurements, this may lead to underestimation of bubble effect or to estimate transition without bubble in cases when laminar separation bubble experimentally exists on airfoil surface.Turbulent separation point locations obtained from current computations are between experimental and XFOIL results.
Although computed lift and drag coefficients deviate from measured data at higher angles of attack, the predicted aerodynamic data allows the use of current procedure in design of airfoils for variety of applications without human intervention utilizing shape perturbation approach.

Figure 1
Figure 1 General computational procedure 2 Computational procedure 2.1 Conformal mapping

Figure 2
Figure 2 Conformal Mapping of airfoil to pure circleA continuous representation of logarithm of near circle radius as function of its angle θ is then obtained by cubic spline interpolation.Finally, near circle shape ς 2 is transformed to true circle ς 3 using the general transformation in the form of Eq. (4), Fig.2d.

Figure 3 Figure 4 Figure 5 Figure 6
Figure 3 Potential flow velocity around circular cylinder

Figure 11
Figure 11 Comparisons of locations of upper surface flow features forE387 and S8036 at Re 200 000, 350 000, and 500 000.(Solid lines represent experimental data, dashed lines is XFOIL, and filled symbols are current calculations).

Figure 17
Figure 17Comparisons between calculated and experimental drag polar for S8036 at Re=500,000.

δ 2 :
Boundary layer momentum thickness δ 3 : Boundary layer kinetic energy thickness H 12 : Boundary layer shape factor, δ 1 /δ 2 H 32 : Boundary layer shape factor, δ 3 /δ 2 c f : Friction coefficient C D : Dissipation coefficient r : Roughness factor U : Potential flow velocity V ∞ : Free stream velocity W : Complex potential s : Surface distance measured from stagnation point Re δ2 : Reynolds number based on momentum thickness α : Angle of attack c d , c l : Airfoil drag and lift coefficients respectively.Cp : Pressure coefficient 6 References