Nonlinear geometric decomposition of airfoils into the thickness and camber contributions

In the thin airfoil theory, the camber line and the thickness distribution of general airfoils are mainly extracted by a linear combination of the upper and lower surfaces, giving rise to geometric distortions at the leading edge. Furthermore, despite the recent eﬀort to obtain analytic expressions for the zero-lift angle of attack and quarter-chord moment coeﬃcient, more analytic generalization is needed for the camber line component in the trigonometric series coeﬃcients. This paper presents a straightforward algorithm to extract the camber line and thickness distribution of general airfoil shapes based on a ﬁnite diﬀerences method and the Bézier curve ﬁtting. Integrals in the thin airfoil theory involving a Bernstein basis are performed, leading to series coeﬃcients involving Gegenbauer polynomials. The algorithm is validated against analytical expressions of the NACA airfoils without introducing or adapting geometric parameters, and the results demonstrate good accuracy. In addition, the present work indicated a signiﬁcantly diﬀerent geometric behavior for the SD7003 airfoil camber slope at the leading edge obtained by the proposed algorithm and the classical linear approximation. Moreover, the method can be coupled conveniently in recent reduced-order models established on the thin airfoil theory to obtain expressions in closed forms for general airfoils.


I. Introduction
A lthough the advancement of computational fluid dynamics (CFD) in the last decades has permitted numerical modeling of complex aerodynamic designs, there is still room for analytical theories and simplified models to shed light on aerodynamic physics.Applying the classical first-order small-disturbance approximation for a thin wing leads to a linear equation for the boundary condition, which can be solved separately as mean-camber line and thickness problems [1].As revealed by Prandtl, the main aerodynamic features of wings with infinite span -the airfoilare significant to aerodynamic investigations of finite-span wings, establishing the relevance of two-dimensional flow analyses around airfoils [2].In this sense, Laplace's equation solutions for inviscid, incompressible, and irrotational flow problems are extensively employed, especially in the study of lifting surfaces [3].
The results obtained in CFD calculations are valid as long as the governing equations and boundary conditions are correctly formulated and implemented [4].In potential flow problems, the linear distribution of singular solutions as sources and vortices on the problem boundaries is an elementary approach.Earlier steady low-speed aerodynamic models have one of its classical thin airfoil theories developed by Glauert [5].By assuming a small-disturbance problem, the model distributes a continuous vortex sheet along the chord and then calculates its vorticity using the no-penetration boundary condition at the camber line.As the solution is not unique, an additional physical condition -known as the Kutta condition -is required [6].Expanding the airfoil vorticity in a trigonometric series gives the relation between the camber line and the main aerodynamic characteristics in Glauert's model.
Despite the proper establishment of steady aerodynamics, the unsteady phenomena represent a branch that still requires substantial investigation compared to the steady one.Examples are the flapping wings of insects and birds and the revolving blades of helicopters and wind turbines [7].Since then, numerous authors in the past century have developed analytical and computational approaches to model and predict unsteady effects, especially for two-dimensional potential flows.Wagner [8] obtained lift and circulatory responses due to a step change in the angle of attack of a flat plate by using the complex analysis of conformal mappings, while Theodorsen [9] solved for small-amplitude harmonic oscillations.Rosenhead [10] introduced the method of approximating the continuous wake's vortex sheet by an array of discrete vortices.Afterward, Katz and Weihs [11] improved Glauert's model by implementing the time-stepping method and discrete vortex shedding for unsteady thin airfoils, and Ramesh et al. [12] have enhanced this model by assuming an intermittent leading-edge vortex (LEV) shedding governed by the maximum instantaneous suction peak.
Recently, an approach approximates the leading edge as a parabola and uses pressure distribution from potential-flow solutions to identify inflow quantities such as airspeed and angle of attack [13,14].
Mostly, these reduced-order methods are established on the thin airfoil theory and obtain the camber line through the airfoil's surface set of points.In this regard, two perspectives exist for the airfoil section construction through the thickness distribution and camber line.While the British convention plots the thickness distribution over the camber line perpendicular to the chord, the American convention plots the distribution perpendicular to the camber slope [15].
The difference between the methods is assumed to be small as traditional airfoils have a maximum camber, usually in the range of 0% (symmetrical sections) to 5%.Nevertheless, much higher camber slopes are used in turbine blades.
Using cubic B-spline curves, Li et al. [16] applied mathematical morphology to reconstruct the blade section from a point-sampled turbine surface.Thus, introduced the method to extract the mean-camber curve, in which the process becomes a nonlinear least-squares problem.
Several applications of curve parametrization, such as the cubic B-spline, were developed for shape design and aerodynamics optimization.Olha et al. [17] constructed a naturally parameterized curve with cubic curvature using an optimization algorithm to solve nonlinear integral equations in aerodynamic problems.Boehm [18] discussed the errors in Bézier and cubic spline presentations of 4-digit NACA airfoils, and Jaiswal [19] studied Bezier curve approximation of airfoils using fourth-order, sixth-order, and eighth-order Bezier curves, where it was solved the linear least squares parametric functional problem for the control points and the nonlinear least squares problem for the nearest points.Ribeiro, Awruch, and Gomes [20] employed genetic algorithms and artificial neural networks coupled with Bézier curves in airfoil optimization for wind turbines.In contrast, Lepine et al. [21] investigate the performance of a nonuniform rational B-spline (NURBS) geometrical representation for the aerodynamic design of wings.
Regarding direct methods of resolving an airfoil as lifting and thickness problems, Barger [22] presented an adaptation of the Theodorsen theory of wing sections of arbitrary shape [23].While Theodorsen uses thin airfoil approximations, Barger exploits the full potential theory to represent the airfoil as a combination of lifting line and thickness distribution.Their model applies a complex analysis of conformal mappings to transform the airfoil potential field into the potential field around a circle.Recently, Cibin and Ranjith [24] have presented analytical expressions of aerodynamic characteristics for arbitrary wing shapes.The airfoil's camber line is obtained as a linear combination of the upper and lower surfaces and is represented in Bernstein polynomials to generate closed-form expressions for the zero-lift angle and the moment coefficient in terms of the beta function.
Despite the precision of the thin airfoil theory in simulating conventional shapes, the thickness effect plays a role in the aerodynamic features.Motta, Guardone, and Quaranta proposed modifying the Theodorsen model to accurately predict the effects of the airfoil thickness on unsteady loads [25].Typically, the classical zero-thickness consideration gives good accuracy for airfoil thickness ratios up to 12% [26].The recent research effort to comprehend the aerodynamics of flapping wings is related to the development of micro-aerial vehicles and power extraction devices, in which thicker shapes generate higher thrust and propulsive efficiency for plunging airfoils [27].Nevertheless, experimental data for plunging maneuvers are historically scarce, which leads to adapting well-known pitching results and methods for predicting flapping airfoil responses [28], where a fundamental concept is the airfoil's effective angle of attack [29].
The proper analysis of the nonlinear geometric decomposition of airfoils into thickness and camber contributions for the thin airfoil theory is deficient in the literature.Numerical computations are focused mainly on obtaining the zero-lift angle of attack and the quarter-chord moment coefficient by a linear approximation and do not present any analytical expressions for the generalized trigonometric series coefficients.The present work aims to develop a straightforward method to extract the camber line and the thickness distribution of general airfoils without resorting to small-angle assumptions.Integrals in the classical thin airfoil theory involving the Bernstein basis are performed, giving rise to series coefficients of Gegenbauer polynomials.Although static aerodynamic validations are conducted, the proposed method is easily coupled to the Fourier series coefficients of many recent unsteady aerodynamic models [30][31][32][33][34][35][36][37][38][39][40].

II. Methodology
This Section is structured as follows.In Section II.A, the first-order small-disturbance approximation is applied to derive the lifting and thickness problems.Section II.B is dedicated to defining the Bernstein polynomials and Bézier curve fitting.In Section II.C, the algorithm based on the finite differences method and Bézier curves is introduced, and its implementation into the thin airfoil theory is explored.

A. Thin airfoil theory
The airfoil's upper and lower surfaces complex functions,   and   respectively, are prescribed for a given thickness distribution   () applied perpendicular to the camber line   () =  +   (),  ∈ [0, ], as illustrated in Fig. 1: where  is the airfoil chord and  () = tan −1  ′  ().Furthermore, Eq. ( 1) can be written as Applying the classical first-order small-disturbance approximation leads to a linear equation for the boundary condition.Considering the two-dimensional airfoil problem, the boundary condition can be transferred to the -axis: where (, ) is the velocity potential,  is the freestream velocity, and  is the geometric angle of attack.One can solve both lifting and thickness problems through the linear boundary condition (4) by distributing singularity elements in the -axis as detailed next.

Lifting problem
Assuming a vorticity distribution () over the -axis (see Fig. 2), the resulting velocity potential is given by and the boundary condition (4) considering the lifting problem becomes The Glauert variable transformation [5] is appropriate, where  and  are related as and the vorticity distribution over the airfoil is taken as a trigonometric series: where the coefficients  0 ,  1 , . . .,   are given by By making use of the Bernoulli equation [3], the pressure coefficient is To obtain the airfoil's lift coefficient, one must integrate the pressure difference over the chord: From Eqs. ( 9), ( 10) and ( 12), the angle of zero-lift  0 is given by The moment coefficient is where From Eqs. ( 10) and ( 14), the quarter-chord moment coefficient is

Thickness problem
The thickness problem is modeled by a source distribution () as shown in Fig. 3.The velocity potential is given by Hence, the induced velocity at  = 0 is where   is the -component of the induced velocity   .The boundary condition (4) leads to the source strength formula From Bernoulli's equation and Eq. ( 17) the pressure coefficient is

B. Bézier curve fitting
Exploring the Bernstein polynomials in describing complex geometries was first made by Casteljau and Bézier in the 1960s [41].Although the slow convergence in approximating continuous functions, Bernstein polynomials found an application in automotive designing.Currently, the so-called Bézier curves often arise in wing design and aerodynamic optimization, and their main properties are detailed next.

Bernstein polynomials
The Bézier curves are described by a sum of control points weighted by a Bernstein basis of degree  and parameter  ∈ [0, 1], defined by Applying the binomial theorem on (1 − ) − gives Therefore, the Bézier curve z() is obtained by the Bernstein basis as: where   ∈ C,  = 0, 1, . . ., , are the control points, and the coordinates x and ỹ are the real and imaginary parts of z(), respectively.

Curve fitting
To find the curve that best fits a given set of points  = [ 1 ,  2 , . . .,   ]  , the following optimization problem should be solved: To solve the linear problem for , an initial set of nodes  is required.In this way, the initial guess  (0)  = (−1)/(−1),  = 1, 2, . . ., , can be taken.Hence, the vector  is calculated as where  ★  =      −1    is the Moore-Penrose pseudo-inverse of   .Therefore, one can apply Newton's method to solve (26).As the Hessian  () of the function  is a diagonal matrix, the method proceeds through the iterations where "⊘" is the Hadamard division (element-wise division).The process restarts with the new control points  ( +1) obtained from Eq. ( 27) until  ()/ ≤   , where   is a suitable error criterion.

C. Finite difference method based on Bézier curves (FDM-B)
Bézier curve shapes the airfoil surface from the given discrete points to proceed with the finite difference method, which starts with an initial guess   =   +   ,  = 1, 2, . . ., , for the camber line as depicted in Fig.

The FDM-B algorithm
Bézier curve fits the surface through the approach described in Section II.B.2 and the set of parameters  = [ 1 ,  2 , . . .,   ]  is obtained.The parameter   , 0 ≤  1 <   <   ≤ 1, related to the leading edge is obtained by solving the root-find problem using a Newton-Raphson iteration: To constrain the problem for the upper and lower surfaces, the change of variable is taken: where  = (  +   )/(  −   ).In this way, we have (  ,   ) = ( 1 ,   ) for the upper surface and (  ,   ) = (  ,   ) for the lower surface.Considering the camber line with  inner points, we have Δ = 1/( + 1) and   = Δ,  = 1, 2, . . ., .Hence, an initial guess  (0)  for the set of points is taken.The normal line at   is given by where a second-order central finite difference is used to approximate the derivative, and its intersections z,  with the airfoil surfaces are obtained solving the root-finding problem: and the parameters    and    are obtained from the transform (30).
The proposed method assumes   (0) =   () = 0 and it is taken into account in Eq. (32).Moreover, the camber slopes at the airfoil edges are calculated from the FDM-B algorithm as From the geometry of the problem shown in Fig. 4 the following root-find problem can be solved: which enforces the camber line to lay at the midpoint of each vector z  − z  .In this way, the Jacobian  ( ) of ( 38) is the diagonal matrix given by where x′ ,  = ℜ( z′ ( ,  )) and ỹ′ ,  = ℑ( z′ ( ,  )), and the camber-line inner points are adjusted at each iteration: until the stopping criteria ||  ( )||/ ≤ ε is achieved.Thus, we have the discrete points   (  ) ≈   = (  ), and the thickness distribution as One can take the generalized thickness distribution function as where   ,  = 0, 1, . . .,   , are shape coefficients -which can be calculated by the Moore-Penrose pseudo-inverse method.The radius of curvature at the leading edge is given by

Implementation in thin airfoil theory
The FDM-B algorithm allows obtaining the discrete points of the camber line and the camber slope (the latter from Eq. ( 32)) of general airfoil sections.In this way, the Bézier curve fitting provides an analytic expression for the camber-line slope.Thus, from Eq. ( 22) we have: where the parameter  is defined from the leading edge to the trailing edge i.e.,   to   .From Eqs. ( 21), ( 30) and (44) and setting (  ,   ) = (  ,   ), we have: From Eqs. ( 10) and ( 45), the following general integral must be solved (see Appendix B) in order to couple the FDM-B method into the thin airfoil theory: where    () is given by Eq. (C27).Therefore, the Fourier coefficients (Eqs.( 9) and ( 10)) become: where are the camber slope coefficients.Thus, different from an infinite Fourier series -which needs truncation of termsthe proposed method uses the exact  terms for the aerodynamic contribution of the camber line.From Eq. ( 13) the angle of zero-lift is and from Eq. ( 14) the quarter-chord moment coefficient is Furthermore, the pressure coefficient (11) becomes Regarding the aerodynamic analysis of the thickness problem, the pressure coefficient ( 19) is obtained from Eqs. ( 17), ( 18) and ( 42) as (see Appendix A): where  and  are related through Eq. ( 7).
Thus, the upper and lower surface pressure coefficients in the first-order approximation problem are the linear superimposition of the lifting and thickness contributions ( 52) and ( 53).The complete methodology is summarized in the flowchart presented in Fig. 5, and its subprocesses are illustrated in Fig. 6.

Yes
No

III. Results & Discussion
During the 1930s, the National Advisory Committee for Aeronautics (NACA) investigated a group of airfoils in the NACA variable-density wind tunnel [42].The related airfoils were developed by combining a thickness distribution and the camber line where the camber-line slope is given by for NACA 4-digit airfoils with the nomenclature "NACA  ", where  = /100 is the maximum camber,  = /10 is the position of the maximum camber, and  = /100 is the maximum thickness.In this way, the referred airfoils become practical in the present analysis as they have closed-form expressions.Thereafter, numerical computations are performed to validate the FDM-B algorithm against analytical results for NACA 4-digit airfoils.
Theoretically, the first-order linear approximation gives good accuracy for airfoil thickness ratios lesser than 12% [26].As the approximation holds as thinner the airfoil is, the present analysis is focused on airfoil thicknesses close to 12%.Then, Bézier curves of degree  = 16 fitted the surface of the NACA 4409 and 2412 airfoils, containing 60 points, with a stopping criterion   = 2 × 10 −9 .The FDM-B algorithm is employed with  = 30 and ε = 1 × 10 −10 for obtaining the main geometric parameters.As shown in Fig. 7, the proposed method demonstrated good accuracy in predicting the camber line, whereas the linear approximation deviates considerably from the exact equation (55), especially in the first quarter of the chord.A reasonable error is found even for small maximum camber airfoils, as in the NACA 2412 airfoil (Fig. 7b), where   (0)/ = 1.52 × 10 −3 in the linear approximation, leading to distortions of the effective angle of attack [22].For the NACA 4409 airfoil, the linear approximation gives   (0)/ = 1.66 × 10 −3 .
One should notice that the proposed algorithm imposes the camber line at the leading and trailing edges to be zero i.e.,   (0) =   () = 0. Overall, the FDM-B algorithm kept the relative errors under 0.5% for the camber line.
After obtaining the camber-line slopes at the inner points, second-order accuracy forward and backward difference approximations give the first derivatives  ′  (0) and  ′  (), respectively, as detailed in Eqs.(36) and (37).Fig. 8 compares the camber slope obtained by the FDM-B algorithm and the linear approximation.As in the camber line analysis, the linear approximation contrasts with the exact formula (56) in the quarter-chord region.The comparison of the camber slope at the leading edge obtained from both methods is presented in Table 1, where the linear approximation gives Although satisfactory accuracy when predicting the thickness distribution along the chord (see Fig. 9), the linear approximation fails at the leading edge, giving   (0)/ = 1.66 × 10 3 and   (0)/ = 1.52 × 10 3 for NACA 4409 and 2412 airfoils, respectively.As the FDM-B algorithm fits the inner points, special attention must be paid to the airfoil edges.In the present analysis, Eq. ( 42) imposes the leading edge to be zero, whereas no assumption is made to the trailing edge.Assuming the generalized thickness distribution function (42) with   = 4 (as in the NACA thickness formula (54)), one can get the shape coefficients   ,  = 0, 1, . . ., 4, by the Moore-Penrose pseudo-inverse method.In this way, the FDM-B algorithm predicts the trailing-edge thickness with 2.9% and 1.5% errors for NACA 4409 and 2412 airfoils, respectively.One should add the equation   () =    to the process of obtaining the shape coefficients to impose a specific thickness    at the trailing edge.Furthermore, the radius of curvature at the leading edge (from Eq. ( 43)) is presented in Table 1.Finally, Fig. 10 shows the Bézier fitting and the resulting composition of the geometric outcomes from the FDM-B algorithm by making use of Eq. ( 1), demonstrating good accuracy when compared against the airfoil data.Even though many more geometric parameter validations would be carried out, the radius of curvature and camber-line slope at the leading edge find applications in recent aerodynamic analysis, as the parabola-airfoil matching [13,14,39].In addition to the geometric verification against analytical results, aerodynamic validations were performed considering the FDM-B algorithm implemented into the thin airfoil theory, as described in Section II.C.2.While Joseph and Mohan [24] proposed analytical expressions for two main aerodynamic characteristics in terms of the beta function, the present work expresses the Fourier coefficients ( 9) and (10) as Gegenbauer polynomials [43].In this way, not only closed-form expressions for the zero-lift angle of attack and the quarter-chord moment coefficient are obtained as for any aerodynamic quantities based on the thin airfoil theory.The analysis method compares the FDM-B algorithm and the linear approximation results against analytic functions obtained from the integration of Eqs. ( 13) and ( 15) by directly using the NACA camber-line slope (56) and the variable transformation (7).Thus, the zero-lift angle of attack and the approximation predicts the camber-line slope at the leading edge 68% greater than the FDM-B algorithm.Furthermore, the linear approximation generates a monotonically non-increasing camber slope -similar to the NACA airfoils.In contrast, the FDM-B algorithm predicts a non-monotonic function, with a leading edge slope  ′  (0)/ = 0.092 and a peak of 0.115 at / = 0.027.This behavior in the camber slope leads to a concave shape in the camber line leading edge, whereas usual airfoils have a convex contour.Bézier curves of degree  = 16 fitted the camber line and camber-line slope points given by the FDM-B algorithm and the linear approximation, with stopping criterion   = 2 × 10 −9 .For the thickness fitting, the Moore-Penrose pseudo-inverse method is applied to obtain the shape coefficients, considering   = 10 and   () = 0. Thus, a composition of these geometric fitted functions is applied in Eq. ( 1) to reconstruct the SD7003 airfoil surface, as depicted in Fig. 13.Despite the good accuracy of the proposed method, one can notice a slight discrepancy in the leading edge composition.Indeed, the FDM-B algorithm is based on a uniform grid, and the edges may need particular chord discretization as in the SD7003 airfoil case, once the camber-line slope presents high gradients at the leading edge.
As illustrated, the linear approximation underestimates the camber slope coefficients when compared to the FDM-B algorithm.The difference between the coefficients obtained from both methods increases from  = 0 to  = 5, decreasing to nearly zero for  = 16, affecting especially the zero-lift angle of attack and the quarter-chord moment coefficient (Eqs.( 50) and ( 51)), as presented in Table 3.The geometric characteristics such as maximum camber and thickness obtained by both methods are also in Table 3 and show excellent agreement with Selig [44].
Airfoil data Bézier fitting FDM-B and Eq. ( 1) Fig. 13 Fitted surfaces of the SD7003 airfoil.Concerning the zero-lift angle of attack and quarter-chord moment coefficient, no detailed experimental data is found in [44] for the SD7003 airfoil.In this way, the XFOIL code [45] is employed for the inviscid aerodynamic analysis.Considering the pressure distribution coefficient as the composition of Eqs. ( 11) and ( 19), Fig. 14 compares the results from both methods against the   obtained by the XFOIL.Both methods demonstrate a similar performance in evaluating the pressure distribution along the chord, where the leading-edge pressure tends to be infinity -as expected in the thin airfoil theory.Yet, the results at the lower and the aft part of the upper surfaces are in excellent agreement with the XFOIL outcomes.As no analytic expressions are available for the SD7003 geometric composition, the FDM-B algorithm is a more suitable alternative than the linear approximation to obtain the main aerodynamic features based on the thin airfoil theory.On this wise, a second-order subsonic airfoil theory [46] may be applied to investigate the accuracy of the proposed method at the leading edge.Furthermore, the zero-lift angle and moment coefficient obtained from the XFOIL code are −1.737• and −0.0413, respectively.Thus, the relative errors of the FDM-B algorithm and linear approximation are 1.2% and 1.6% for the zero-lift angle and 3.6% and 4.6% for the moment coefficient, respectively.However, one should point out that the XFOIL code is founded on the panel method, whereas the proposed method uses the thin airfoil theory approximations.

IV. Conclusion
A straightforward method was implemented to obtain the camber line and the thickness distribution of general airfoils.The model is based on a finite differences method coupled with Bézier curve fitting without resorting to small-angle assumptions or additional geometric parameters.Thus, integrals in the aerodynamic thin airfoil theory involving the Bernstein basis are performed, giving rise to series coefficients of Gegenbauer polynomials.The FDM-B

Fig. 2
Fig. 2 Induced velocity   at a point (, ) by the vortex sheet.

Fig. 3
Fig. 3 Induced velocity   -and its Cartesian components -at the point (, ) by the source sheet.

Fig. 5
Fig. 5 Flowchart of the nonlinear geometric decomposition method.

Fig. 7 Fig. 8
Fig. 7 Comparison of the camber line obtained by the linear approximation and the FDM-B algorithm.

Fig. 9
Fig. 9 Comparison of the thickness distribution obtained by the linear approximation and the FDM-B algorithm.

Fig. 11 Fig. 12
Fig. 11 Comparison of the camber line and the thickness distribution obtained by the linear approximation and the FDM-B algorithm for the SD7003 airfoil.

Fig. 14
Fig. 14 Comparison of the pressure distribution obtained by the linear approximation, the FDM-B algorithm, and the XFOIL code.