The effect of variations in first- and second-order derivatives on airfoil aerodynamic performance

ABSTRACT The geometric factors which influence airfoil aerodynamic performance are attributed to variations in local first- and second-order curvature derivatives. Based on a self-developed computational fluid dynamics (CFD) program called UCFD, the influence of local profile variations on airfoil aerodynamic performance in different pressure areas is investigated. The results show that variations in first- and second-order derivatives of the airfoil profiles can cause fluctuations in airfoil aerodynamic performance. The greater the variation in local first- and second-order derivatives, the greater the fluctuation amplitude of the airfoil aerodynamic coefficients. Moreover, at the area near the leading edge and the shock-wave position, the surface pressure is more sensitive to changes in first- and second-order derivatives. These results provide a reference for airfoil aerodynamic shape design.


Introduction
The aerodynamic performance of an airfoil is determined by its aerodynamic shape. However, the aerodynamic shape of the airfoil can be influenced in various ways, such as the connection of different types of curves during the airfoil design process, the defects or errors caused by imprecision in the manufacturing process, and the deviation of the airfoil profile while dispersing into polylines during the process of aerodynamic performance prediction. The sensitivity of geometric parameters to aerodynamic performance of an airfoil is an important issue, and a number of studies have been conducted in this area. Walraevens and Cumpsty (1995) experimentally proved that the aerodynamic performance of compressor airfoils is quite sensitive to the geometry of the leading edge. In a low-speed wind tunnel, H. X. Liu, Jiang, and Chen (2004) investigated the flow around an airfoil with a circular leading edge and an elliptic leading edge; they found that the geometry of the leading edge had a large influence on the flow separation and transition near the leading edge and also on the later development of the boundary layer. Song, Gu, and Xiao (2014) used Reynolds-averaged Navier-Stokes (RANS) computations to show that a continuous-curvature leading edge reduces the minimum loss and enlarges the working range because the leading edge spike is removed CONTACT Diangui Huang dghuang@usst.edu.cn when the separation bubble is eliminated or reduced. Using a numerical simulation method, Lu, Xu, and Fang (2000) studied the phenomenon that a suction peak can be formed due to excessive expansion of the flow at the surface of a circular leading edge. Elmstrom, Millsaps, Hobson, and Patterson (2005) used a computational fluid dynamics (CFD) method to study the aerodynamic performance of a compressor airfoil while applying uniform or non-uniform coats with different thicknesses at the leading edge. With a change of shape at the leading edge, the pressure gradient near the leading edge shows significant variation. In addition, under certain conditions, laminar flow separation and premature turbulent transition can be caused. Based on the automatic differential principle and RANS with a finite volume method, Xu, Wang, Wu, and Ye (2014) showed that on the suction side of the airfoil surface, the flow field in the supersonic area is more sensitive to geometric error than in the subsonic area. In the area of the leading and trailing edges of an airfoil, the surface coordinates have a large impact on the flow field and vice versa. C. G. Li, Li, and Chen (2012) proved that the attachment on a blade surface can change the geometric factors of an airfoil, and negatively affect the aerodynamic characteristics of a wind turbine in operation. In numerical simulations, using a coarse grid near the geometry might cause deviations in the first-order derivative between two adjacent grid points, which can lead to inaccurate results. Therefore, it is of value to research the influence of these deviations, whether they are caused by the grid resolution or attachments on the surface.
Relatively few studies in the existing literature provide a systematic examination of the relationship between some important geometric parameters that describe an airfoil's shape and its aerodynamic behavior. These geometric factors -including the geometric continuous connection among different types of curves during the airfoil design process, the allowances for flaws or errors incorporated during the manufacturing process, and the deviation due to the conversion of a continuous curve into polylines for the numerical prediction of aerodynamic performance -can all result in variations of local firstand second-order derivatives. In this paper, NACA 0012 and RAE 2822 airfoils are used as a baseline. The influence of variations in local first-and second-order derivatives is studied in detail and observed in several different positions, as well as under different pressure gradients. Accordingly, the main purpose of this paper is to identify the effect that continuity in the curvature smoothness of the profile of an airfoil has on the accuracy of numerical predictions of its aerodynamic performance.

Geometric model design
The geometric factors influencing airfoil aerodynamic performance are attributed to variations in local firstand second-order derivatives. In this paper, the following definition is provided: as shown in Figure 1, Curve 1 and Curve 2 are both parts of the rough airfoil profile, placed on the left and right sides of a reference node O, respectively. At point O, the first-order derivative for Curve 1 is y 1 and the second-order derivative is y 1 . Similarly, the first-order derivative at point O for Curve 2 is y 2 and the second-order derivative is y 2 . As such, the differences between the local first-and second-order derivatives of Curves 1 and 2 are defined as m = y 1 − y 2 and n = y 1 − y 2 , respectively.
In order to study the influence of the local first-order derivative m and second-order derivative n on the aerodynamic performance of the studied airfoils, corresponding profiles of the baseline airfoils (the NACA 0012 and the RAE 2822) were designed in accordance with the following criteria.
2.1.1. Criterion 1: m = 0, n = 0 In order to satisfy the above requirement, the leading edge of the symmetrical NACA 0012 airfoil is replaced with a circular arc and an elliptic arc, as illustrated in Figures 2 and 3, respectively. It is also required that the drawn arc is exactly tangential to the rest of airfoil at the points where they merge. In this way, m = 0 and n = 0 at the join points for the designed airfoil can be guaranteed.

Criterion
2: m = 0, n = 0 (y 1 = y 2 = 0) In order to parameterize the influence of the surface attachment on geometric factors of the airfoil, a small protruding equilateral triangle is placed on the surfaces of both airfoils at different chordwise positions. Considering the machining accuracy during the manufacturing process, the degree of geometric error is 63 μm (Xu et al., 2014) under an IT7 grade of tolerance while the length of a chord is 1 m. Thus, the side length of the small protruding equilateral triangle was set to 50 μm, which is smaller than the geometric error. The values of m and n for two sides of the triangle at its vertex located outside the body of airfoil can be calculated and the influence of their variations on the airfoil performance can then be examined.
Furthermore, a discontinuity in the first-and secondorder derivatives at the points of connection between the sides of the triangle and the curve of the surface of the airfoil should be avoided. Hence, a non-uniform rational B-spline (NURBS) spline curve was adopted to improve the smoothness of the airfoil surface with a protruding triangle. As shown in Figure 3(a), four points on the original airfoil surface and three points on one side of the triangle are selected. A smooth NURBS spline curve is then drawn to pass through the given set of points. The purpose of using a NURBS spline curve is not only to  connect all the ordered data points, but also to ensure that the first-and second-order derivatives at all the selected data points have identical values. In this way, as shown in Figure 3(b), a protruding small equilateral triangle is smoothly joined to the airfoil. Thus, the drawn NURBS spline curve is characterized by the continuity and uniform curvature variation. The continuity of the first-and second-order derivatives at the connecting points is satisfied (Qin & Zhang, 1995). It is the case that m = 0, n = 0 (y 1 = y 2 = 0) occurs only at the third vertex of the added triangle outside of the designed airfoil.
For the symmetrical NACA 0012 airfoil shown in Figure 4, straight lines with several different slopes (m 1 ∼ m 5 ) are designed to close the rear part of the airfoil in order to meet the requirements of m = 0, n = 0 (y 1 = y 2 = 0) at the trailing edge.

Criterion 3
: m = 0, n = 0 (y 1 = y 2 = 0) To meet the above criterion, the middle regions of the upper and lower surfaces of the NACA 0012 airfoil are replaced with two elliptical segments that are symmetrical with respect to each other ( Figure 5). In addition, the elliptical segment keeps completely tangential to the surface of the airfoil at the point where two curves meet and thereby m = 0, n = 0 (y 1 = y 2 = 0) is realized for the modified airfoil.

Numerical schemes and grid strategy
An in-house CFD program called UCFD developed by the authors' research team was utilized to provide a numerical prediction of the aerodynamic performance of the modified airfoils. The numerical results were then compared to those obtained by the variable domain variational finite-element method (FEM).
The UCFD program (Huang, 2011) uses the finitevolume method and adopts the velocity, pressure, and enthalpy as original variables (Huang, Zhuang, & Cai, 2007). A dual-time procedure is implemented into a new preconditioning technology based on the primal variables of pressure, velocity and enthalpy, wherein the physical time layer is used to track the physical variety of the flow and the pseudo time layer is used to iterate every physical time step (Huang, 2006;X.-S. Li & Gu, 2008, 2013X.-S. Li, Gu, & Xu, 2009). The convective and pressure terms are differenced using either the upwind flux-difference splitting technique of Roe (1981) or the flux-vector splitting technique of van Leer (1982).
The governing equation can be expressed as: (1) The absolute pressure p can be expressed as p = p 0 + p g , where p 0 is the atmosphere pressure and p g is the gauge pressure. Thus, ∇(p 0 + p g ) = ∇p g .
In Equation (1), the conservation variables of density, momentum, and total energy can be expressed as: where J is the Jacobian of the transformation between Cartesian (x, y, z) and generalized (ξ , η, ζ ) coordinates: The inviscid fluxes can be expressed by: where U, V, and W are contravariant velocities shown as: The viscous fluxes are: , and The shear stress and flux terms are defined in tensor notations as: where e is expressed as: Based on the preprocessing of Equation (1), the following equation can be obtained: Finally the equation can be expressed as For the investigation of a single airfoil flow field, the Menter's k-ω shear stress transport (SST) model is adopted, as it has been widely used and has proven to be numerically well behaved in most cases.
Based on G. L. Liu's (1995) variable domain theory for airfoil aerodynamic airfoil design, a numerical simulation program has been developed using a technique called the potential function method. The flow field around the airfoil can thus be modeled and analyzed relatively easily.
A structured quadrilateral mesh generated by the Pointwise meshing software was adopted in the simulations in this study, and the meshes around all the joint nodes defined in section 2.1 were refined ( Figure 6). An O-type non-uniform structure grid system was adopted to discretize the entire computational domain. The whole calculation domain is circular and the radius of the calculation domain is 20 times the airfoil chord length. The pressure-far-field boundary condition was used in the flow field and a no-slip wall boundary condition was applied at the surface of the airfoil.
The mesh points used for the numerical study of the first-and second-order derivatives were chosen from among the geometric data points in order to consider the deviation caused by the spatial discretization that breaks lines in geometry into short, straight-line segments for a polyline representation during simulation.

Mesh-independence study and code validation
In order to ensure the accuracy of the results, simulations with different grids and domains were carried out with the UCFD program using the Menter's k-ω SST model. A validation test was conducted on the flow field around an RAE 2822 airfoil at a 2.72°angle of attack and a Mach number of 0.75 (see McDonald & Firmin, 1979).
Four different levels of grid (a coarse mesh of 13,664 cells, a medium mesh of 27,840 cells, a fine mesh of 55,552  Figure 7 shows the results of the pressure coefficient distribution along the airfoil surface simulated for different levels of grid compared with the experimental data of McDonald and Firmin (1979). The results are close to each other in most areas except the shock-wave section (0.55 times the chord length on the suction surface). The results of the fine and very fine grids are almost same and fit the experimental data best. The pressure coefficient distribution along the airfoil surface simulated for different levels of grid domain is shown in Figure 8, demonstrating that the data for both the medium and large domains fit the experimental results best.
Tables 1 and 2 show the results of the lift coefficient (Cl) and convergence time (iteration number when monitoring data remain constant). The results of Cl and pressure coefficient distribution show the same regularity that the simulation results of the fine grid and the medium domain are acceptably closed to the experimental results. Balancing the simulation accuracy against the time it takes to run, it was decided that the fine mesh (55,552 cells) and the medium domain (20 times the chord length) were the best combination for use in this study.

The baseline NACA 0012 airfoil
The NACA 0012 airfoil is a typical symmetrical airfoil that is frequently used for benchmarking test cases. The  Gregory and OReilly (1970), which is shown in Figure 9. The lift coefficient under these working conditions is extremely closed to a magnitude of 10 −6 , which causes difficulty in measuring the variation. Therefore, the pressure coefficient along the airfoil surface was chosen as the main reference for research into the influence of first-and second-order derivatives on the NACA 0012 airfoil. It was found that compared to the results obtained with the RANS analysis method, the surface pressure distribution across the airfoil produced using the variable domain variational FEM is more sensitive to variations in local first-order derivatives. As shown in Figure 9, if the trailing edge of the NACA 0012 airfoil is closed with the extension line of the original airfoil profile and the variable domain variational FEM is adopted, the pressure fluctuations on the airfoil surface near the trailing edge are very obvious, with a fluctuation amplitude of about 0.5. However, the fluctuation amplitude of the pressure coefficient near the trailing edge obtained from the RANS method is only 0.06. Nevertheless, it is obvious that the result from the variable domain variational FEM is less accurate than that of the UCFD program when compared with the experimental data -that is to say, the fluctuation of the pressure coefficient can be magnified by using the variable domain variational FEM and this characteristic is very useful for studying the aerodynamic performance of an airfoil, as tiny pressure coefficient fluctuations are more likely to be discovered. Thus, in this study the variable domain variational FEM was only used to find the tiny pressure coefficient fluctuations, and the variable quantities were measured using the UCFD program.

Shape changes in the NACA 0012 airfoil
In addition, the leading edge of the airfoil profile is formed by a circular arc (Figure 2) with a radius of 0.04 and the conditions that the first-order derivative (m = 0) and second-order derivative (n = 1.19) at the connection points are satisfied. The surface pressure coefficient distributions obtained from the UCFD program are shown in Figure 10, illustrated using the results from the original airfoil. Compared to the pressure distribution over the baseline airfoil's surface, the results show that a distinct 'suction peak' in the pressure coefficient occurs due to variations in second-order derivatives.
Subsequently, the case of an elliptical leading edge (Figure 3) which meets the requirements of m = 0 and n = 0.40 was also examined. The major to minor axis ratio of the ellipse is 3:1. The corresponding surface pressure coefficient distribution is presented in Figure 11 with the results from the original NACA 0012 airfoil.  In comparing Figures 10 and 11, it can be seen that variations in second-order derivatives of the airfoil curve have a noticeable impact on the surface pressure coefficient distribution. When a circular arc leading edge is used, there is a significant change in the local secondorder derivative n at the connection point between the circular arc and the airfoil curve which results in a dramatic fluctuation in the airfoil surface pressure coefficient. In contrast to the circular arc leading edge, the n of the elliptic leading edge turns out to be small and this characteristic leads to the disappearance of the peak of the surface pressure distribution. Therefore, m = 0 and thus larger values of n could be a cause of greater fluctuations along the curve of airfoil surface pressure coefficients.
The influence of a protruding equilateral triangle placed at different chordwise positions (0.03 chord length, 0.4 chord length and 0.95 chord length) on the airfoil surface pressure distribution was also studied. Using an equilateral triangle, the conditions m = 0 and n = 0 can be maintained at a joint between one side of the triangle and the airfoil surface. The surface pressure coefficient curves are provided in Figure 12. When the protruding triangle is placed near the leading edge, the surface pressure coefficient experiences the largest fluctuation with an amplitude of 0.15, as shown in Figure 12(a). When the protruding triangle is placed in the middle of the airfoil or at the trailing edge, the surface pressure is found to undergo small-amplitude fluctuations with amplitudes of 0.04 and 0.02, respectively, as shown in Figures 12(b) and 12(c). These results indicate that the leading edge of an airfoil is more sensitive to variations in the first-order derivative m than in the second-order derivative n. Due to the smaller radius of the curvature, a large pressure gradient exists near the leading edge chordwise. Therefore, even tiny changes in the airfoil profile in the leading edge area will cause a large-amplitude pressure fluctuation. This explains the phenomenon which can be seen in Figure 12 wherein the closer the derivative variation position is to the leading edge, the higher the pressure fluctuation becomes, which means that the leading edge of an airfoil is more sensitive to the variation than the middle and the trailing edge.
The airfoil section created from the NACA 0012 airfoil coordinates software is not closed at the trailing edge (Figure 4). However, it must be a closed curve for conducting the following numerical analysis. In Figure 4, the enlarged view shows that straightline segments with slopes of 0.50, 0.29, 0.20, 0.15, and 0.14, and −0.50, −0.29, −0.20, −0.15, and −0.14 are therefore used to close the trailing edge. The values of the first-order derivative m are m 1 = 0.36, m 2 = 0.15, m 3 = 0.06, m 4 = 0.01, and m 5 = 0.00, and the values of the second-order derivative n are 0.00. Figure 13 presents the surface pressure coefficient distributions obtained accordingly. The results indicate that the higher the value of m, the greater the fluctuation amplitude of the surface pressure coefficient at n = 0.

Shape changes in the RAE 2822 airfoil
In transonic flows, one of the most commonly-used test cases is the flow around an RAE 2822 airfoil. Thus, it was selected by 16 European cooperation projects and Advisory Group for Aerospace Research and Development Figure 13. Results for the trailing edge closed by different line segments. (Zhang & Zhang, 2009) as a typical numerical validation case of a two-dimensional transitional turbulence flow model. Due to the variety of the distribution areas of the pressure coefficient under transonic operating conditions, the flow field around an RAE 2822 airfoil at α = 2.72°and Ma = 0.75 was simulated using the UCFD program, and the numerical results are validated by the experimental results of McDonald and Firmin (1979) in Figure 7.
In order to study the influence of local first-and second-order derivatives on the aerodynamic performance of an RAE 2822 airfoil, a protruding equilateral triangle with a side length of 50 μm was designed and placed at different locations on the surface of the RAE 2822 airfoil. As shown in Figure 14, these locations are: the airfoil leading edge where a large favorable pressure gradient exists; the suction surface near the leading edge where there is a transition zone between the favorable and adverse pressure gradients; the suction surface near the leading edge where there is an adverse pressure gradient; the region close to the trailing edge; the pressure surface where there is a favorable pressure gradient; the pressure surface where there is an adverse pressure gradient; and the region near the shock wave.
When the protruding equilateral triangle is placed at 0.002 of the chord length downstream from the leading edge on the suction surface, n = 0 and m = 0 is achieved at the designed points for the following positions: the transition zone between the favorable and adverse pressure gradients near the leading edge on the suction surface (0.043 of the chord length) and the adverse pressure area (0.060 of the chord length) near the leading edge on  Table 3. By comparing the results from different values of Cp of 0.95, 0.23 and 0.21, it is found that the largest fluctuation amplitude occurs in the favorable pressure area near the leading edge on the suction surface, where the largest difference of Cl and Cd from the original airfoil occurs. This means that the closer the position of the first-order derivative variation is to the leading edge, the higher the fluctuation in aerodynamic performanceand this outcome is quite similar to the results found for the NACA 0012 airfoil.
At the transition zone (0.043 of the chord length) between the favorable and adverse pressure gradients where the lowest surface pressure is located, although the fluctuation amplitude caused by this protruding equilateral triangle appears to be smaller than those caused by the triangle when placed on the upper surface in the favorable pressure area, the surface pressure still fluctuates within a wide range chordwise. This phenomenon indicates that variations in the first-order derivative at this position lead to surface pressure fluctuations over a large percentage of the airfoil chord.
Finally, the protruding equilateral triangle was set to the region near the shock wave (0.55 of the chord length). It was found that the surface pressure fluctuations appear in front of the shock wave ( (Table 3). It is noted in Figure 16 that the position of the shock wave changes significantly due to the existence of small triangle.
The vorticity contours near the shock wave are provided in Figure 17 and clearly demonstrate that the presence of the protruding equilateral triangle at the point of the shock wave substantially disturbs the flow patterns around the airfoil. Consequently, a vortex is generated near the position of the shock wave that produces changes in the pressure distribution on the surface of the airfoil. As a result, the position of the shock wave changes as well.
In the following studies, the protruding equilateral triangle is placed on the pressure surface of the airfoil at the positions corresponding to the favorable pressure area (0.20 of the chord length near the leading edge), the transition zone between favorable and adverse pressure gradients (0.38 of the chord length), and the adverse pressure area (0.70 of the chord length). The resulting surface   pressure distributions are shown in Figure 18, and the fluctuation amplitudes of the predicted pressure coefficients are found to be 0.160, 0.150 and 0.045, respectively, as shown in Table 4 with the results of Cp, Cl and Cd. By comparison, it was found that the surface pressure and Cl are more sensitive to variations in the first-order derivative m in the favorable pressure area (0.20 of the chord length) and the transition zone (0.38 of the chord length). However, the influence of the variation of m on the airfoil surface pressure is less apparent when the designed triangle is placed in the region of the adverse pressure area (0.70 of the chord length).
Finally, the small equilateral triangle was placed at the trailing edge on the airfoil suction and pressure surfaces, and the numerical predictions of the surface pressure coefficients are shown in Figure 19, in which the pressure fluctuation amplitudes obtained are 0.005 and 0.045. This result suggests that the effect of variations in the firstorder derivative m at the trailing edge on the pressure surface is more pronounced than the effect it has on the suction surface.

Conclusions
In this paper, the geometric factors that influence an airfoil's aerodynamic performance were ascribed to variations in local first-and second-order derivatives. NACA 0012 and RAE 2822 airfoils were adopted as the baseline onto which multiple small variations of the geometry were applied. The resulting flow field around the airfoil was then simulated in an attempt to systematically study the influence of variations of local first-and second-order derivatives on the aerodynamic characteristics in different regions of the airfoil surface. The main purpose is thus to identify universal design criteria which can be applied to airfoil design for better aerodynamic performance. The results suggest that variations in both the firstand second-order derivatives along the airfoil profile can cause fluctuations in airfoil aerodynamic performance, including in the pressure coefficient, lift coefficient, and drag coefficient. The greater the variation in the firstand second-order derivatives, the larger the fluctuation amplitude with respect to the airfoil aerodynamic performance.
Under the operating condition of m = 0, larger values of n can lead to greater surface pressure fluctuations. Variations in the local second-order derivative n significantly decrease when an elliptic arc is used at the leading edge, and the peak of the pressure coefficient curve can be decreased accordingly. Under the operating condition of n = 0, larger values of m can cause greater surface pressure fluctuations.
For positions around the leading edge, the shockwave on the RAE 2822 airfoil's surface and the trailing edge, the surface pressure distributions are all sensitive to variations in first-and second-order derivatives. Therefore, particular attention should be paid to these areas during the airfoil design process.

Disclosure statement
No potential conflict of interest was reported by the authors.