A New Method of Sectional Expression of Wind Turbine Airfoil Function

Abstract

To expand the airfoil database and realize the reverse digital airfoil design, this paper proposes a new sectional expression function of wind turbine airfoil based on the Joukowski transformation and derives the function equation for the novel airfoil. Compared with the existing airfoil function, the new airfoil function can adjust the parameter values to control the relative thickness, relative thickness position, relative camber, and relative camber position of the airfoil. Taking the common airfoil (NACA63-212) as an example, the existing airfoil is fitted according to the expected direction. The maximum fitting deviation is $1.22135\times {10}^{-3}$, better than the four existing methods for airfoil shape fitting, which demonstrates well the practicability and accuracy of this method. It also shows that changing the parameter values can generate a new airfoil to enrich the airfoil database. Therefore, the proposed method broadens the ideas for airfoil design and method research.

1. Introduction

Energy is an important material basis for the survival and development of human society, which is related to economic development, social progress, and national security, and also occupies a major position in the development strategies of various countries [1]. Fossil fuels have played an extremely important role in the process of global economic and social development. However, due to the declining reserves, alternatives need to be sought, the most important of which is renewable energy. Among the many types of renewable energy, wind energy has unique advantages. On the one hand, wind energy has large reserves and wide distribution with broad development prospects; on the other hand, wind energy has high energy conversion efficiency which is suitable for large-scale industrial development [2].

Wind power generation is the main form of wind energy utilization [3]. Wind turbines are rotating mechanical devices that convert wind energy into electrical energy. For wind turbines, the airfoil is the basis for the blade shape and structure design. Moreover, the airfoil shape directly affects the aerodynamic and structural characteristics of the blade and determines the pressure distribution of the entire flow field, which is the basis and key of flow analysis. Therefore, the expression method of the airfoil shape is very important, which affects the efficiency of the subsequent calculation process and the comprehensive performance of the airfoil [4].

The parameterization of airfoil means that the airfoil shape is represented by a series of parameters, and the airfoil shape can be changed by adjusting each parameter value. The existing parametric expression methods of wind turbine airfoil include the spline curve method, Hicks-Henne method, PARSEC method, mesh deformation method, class/shape function transformation method (CST), etc.

The spline curve method refers to the method of fitting the airfoil shape by linearly weighting basis functions, such as B-spline [5,6], non-uniform rational B-spline (NURBS) [7,8], and Bernstein polynomial. SF Hosseini [9] and Pérez-Arribas [10] proposed a method for HAWT blade shape design based on B-spline surfaces, and provided detailed graphics, geometric parameters, and the descriptions and their influence on airfoil shape. Their main objective was to generate smooth blade surfaces with minimal curvature variation. Wei [11] proposed a high-precision airfoil parameterization and multi-objective optimization program. The coordinates of eight control points are selected as the design variables by the Bezier curve in the airfoil parameterization to ensure the fitting accuracy. Hao [12] held that the advantages of the Bezier curve can reduce the control points number required to compose the curve, which can promote neural network training to optimize the airfoil. Fitting the airfoil S809 with a Bezier curve, a new airfoil was generated with better performance than S809. Zhi [13] used regularization technology to improve the Bezier curve by comparing the different Bezier parameters of the airfoil and concluded that the fitting quality of the Bezier curve depends on the way the abscissa of the control point is specified.

The Hicks-Henne method was proposed by R.M. Hicks and P.A. Henne [14]. Its principle was to use several different functions to describe the airfoil and obtain the parametric representation of the airfoil. This method superimposes the disturbance function described by the parameters on the original airfoil to form a new airfoil. The trailing edge shape of the airfoil affects the fluid flow on the airfoil surface, which in turn affects the aerodynamic performance of the airfoil. To change the angle of the profile line at the trailing edge of the airfoil, Wang [15] added a new basis function at the trailing edge. The improved basis function had a larger amplitude at the trailing edge of the airfoil, which ensures the new airfoil shape changes at a certain amount at the trailing edge and improves the overlap of the new airfoil shape at the trailing edge. Using the improved Hicks-Henne function method, Li [16] parameterized the airfoil shape, keeping the lift-drag ratio of the airfoil as the optimization objective. Moreover, the low-speed wind turbine airfoil is optimized by considering the basis function coefficients of the Hicks-Henne function as the design variable.

The PARSEC method was proposed by Sobieczky [17]. The airfoil shape was described by the eleven geometric parameters, such as the leading edge radius, trailing edge longitudinal ordinate, trailing edge thickness, maximum thickness coordinates of upper and lower surfaces, curvature, trailing edge angle, and trailing edge direction angle of the airfoil. Vecchia [18] proposed a new airfoil shape design procedure based on PARSEC parameterization and genetic algorithm optimization, which studied new airfoil function expressions with good aerodynamic characteristics. The mesh deformation method [19] is based on the finite element mesh model and selects the control points of the airfoil by changing the coordinates of the control points to change the airfoil shape. However, the adjacent relationship between the grid nodes does not change, and there is no need to re-mesh. The method is intuitive, simple to operate, and simplifies the optimization process without repeated meshing. However, the selection of the control points requires a certain design experience and aerodynamic theory basis.

The CST method was proposed by Kulfan and Bussoletti [20] in 2006 and is widely applied because of advantages such as better parameter adjustment, fewer design variables, higher shape fitting accuracy, and wide design range. Wu et al. [21] converted the basis function of the CST parameter method into a simplified orthogonal POD basis function to describe the aerodynamic shapes changes. The results show that, compared with the original high-dimensional CST method, the developed parametric method fits the airfoil well with fewer parameters and maintains almost the same ability, thus overcoming the contradiction of high-dimensional design parameters. Bu et al. [22] proposed the CST method of the superimposed thickness distribution of the mid-arc line. A new airfoil was obtained by using CST equation to express the superimposed thickness distribution of the mid-arc line of the airfoil. The feature of this method is that due to the excessive range of design variables, it can prevent the crossover of the airfoil geometry (especially near the trailing edge). The research shows that this method has a higher efficiency than the direct CST method.

In addition, there is also a compound mapping method. Chen [23] proposed a general integral expression of the airfoil shape based on the generalized functional and the Joukowski conformal transformation. The airfoil profile can be determined by changing the coefficients of its shape function, but the relative thickness position and the relative camber position cannot be transformed.

In general, the spline curve method can locally control and smooth the curve shape, but if the control points are not selected properly, wavy airfoils will easily appear. The Bézier curve method has the advantages of intuitiveness and convexity, but the Bézier curve cannot be modified locally, and changing any control point position will affect the whole fitted curve. The Hicks-Henne method has good control over the airfoil shape, but the variable range should not be too large; otherwise, the airfoil will appear unsmooth, which is suitable for the small adjustments of the airfoil. PARSEC method has good robustness and a wave shape is unlikely to appear, but the ability to describe the details is poor, so it is suitable for the rough airfoil design. Furthermore, the CST method can be used to describe a larger design space, but this method is less robust. To this end, this research proposes a new method of segmental expression of the wind turbine airfoil functions, which can both fit existing airfoils and generate new airfoils.

This paper is organized as follows: Section 2 introduces the Joukowski Transformation, Section 3 derives the function equations of the airfoil function segmental expression, Section 4 verifies the feasibility and effectiveness of the segmental expression model of the airfoil function, and Section 5 provides the conclusion.

2. Joukowski Transformation

The Joukowski transformation [24] is a conformal mapping. In a complex variable function, one plane is transformed into another plane through a functional relationship. Such a transformation is a conformal transformation where the derivatives of all analytic functions are not zero. The Joukowski transformation is a very important transformation function in conformal transformation. The analytic function is $z=\zeta +C^2/\zeta$, where $C$ is a positive real number.

As shown in Figure 1, the circle with radius $R$ and center in the second quadrant is transformed by Joukowski transformation on plane $\zeta$ to obtain an airfoil with certain relative thickness and relative camber on plane $z$. The distance between the circle center and the origin of the coordinates is $p$, $p=\sqrt{m^2+n^2}$, $m$ is the abscissa of the circle center in the second quadrant, and $n$ is the ordinate of the circle center. On the plane $\zeta$, the points on the contour are represented by the polar coordinates: $\zeta =R{e}^{ia}$, and $R$ can be expressed as a function of $\alpha$. Substitute $\zeta$ into the analytic function of the Joukowski transformation:

$$z=R{\mathrm{e}}^{ia}+\frac{{(\frac{C}{2})}^2}{R}{e}^{-ia}=\left[R+\frac{{(\frac{C}{2})}^2}{R}\right]\cos \alpha +\mathrm{i}\left[R-\frac{{(\frac{C}{2})}^2}{R}\right]\sin \alpha$$

(1)

Thus, the parameter equation of the airfoil line on the plane $z$ with $\alpha$ is as follows:

$$\left\{\begin{array}{l}X=\left[R+\frac{{(\frac{C}{2})}^2}{R}\right]\cos \alpha \\ Y=\left[R-\frac{{(\frac{C}{2})}^2}{R}\right]\sin \alpha \end{array}\right.\alpha \in [0,2\mathit{\pi }]$$

(2)

which is the parameter equation of the upper airfoil line when $\alpha \in [0,\mathit{\pi }]$ and it is the lower airfoil line parameter equation when $\alpha \in [\mathit{\pi },2\mathit{\pi }]$. $(X,Y)$ is the airfoil coordinate, $R$ is the vector length of the airfoil in the complex plane $\zeta$, $\alpha$ is the argument, and $C$ is a positive real number.

The flow pattern research of the airfoil is based on the Joukowski transformation of conformal theory, which can transform the graph on certain plane to another plane and can obtain another graph. The transformed graph not only depends on the specific transformation formula, but also depends on the position of the original graph. The parameter of $m$ controls the relative thickness of the airfoil, and the parameter of $n$ controls the relative camber value of the airfoil. As shown in Figure 2a, the same value of $m$ means that the relative thickness of the airfoil is the same. As the value of $n$ increases, the relative camber of the airfoil becomes larger; in Figure 2b, the same value of $n$ means that the relative camber of the airfoil is the same. As the value of $m$ increases, the relative thickness of the airfoil becomes larger; in Figure 2c, under different values of $m$ and $n$, the relative thickness and relative camber of the airfoil are obtained after the Joukowski transformation.

However, whatever the parameter values are changed to, the relative thickness position of the obtained airfoil is at 25% of the chord and the relative camber position is at 50% of the chord. The relative thickness position and relative camber position are unchanged, and the relative thickness value and relative camber value of the airfoil can only be changed at a specific position.

3. Sectional Expression Model of Wind Turbine Airfoil Function

As shown in Figure 3, to solve the situation that the relative thickness position and the relative camber position are fixed, the circle in the second quadrant is modified to a rotating elliptical circle. Then, the Joukowski transformation is applied to obtain the airfoil with variable relative camber, relative thickness, relative camber position, and relative thickness position on the plane $z$.

The ellipse is moved a certain distance along the real axis, and the airfoil with a variable relative thickness position is obtained by the Joukowski transformation. The relationship between the major axis and the minor axis of the ellipse determines whether the relative thickness position is close to the leading edge or the trailing edge. Thus, a symmetric airfoil with different relative thickness positions can be obtained by moving the rotating ellipse for a certain distance along the imaginary axis, and the airfoil with a variable relative camber position can be obtained through the Joukowski transformation. The size of the major axis determines the relative camber position of the airfoil. The transformed airfoil parameter equations are added to obtain the variable relative camber position and the airfoil parameter equation with the variable relative thickness position.

3.1. Determination of Relative Thickness Position

As shown in Figure 4a, the parametric equation of the ellipse at the coordinate origin ($A$) on the plane $\zeta$ is as follows:

$$\left\{\begin{array}{l}x=a\cos \theta \\ y=b\sin \theta \end{array}\right.\theta \in [0,2\mathit{\pi }]$$

(3)

where $a$ is the major axis of the ellipse, $b$ is the minor axis of the ellipse, and $\theta$ is the parameter.

If the ellipse is translated the distance of $m$ along the real axis ($m$ is a positive real number), then the ellipse parameter equation is as follows:

$$\left\{\begin{array}{l}X=a\cos \theta -m\\ Y=b\sin \theta \end{array}\right.\theta \in [0,2\mathit{\pi }]$$

(4)

Among them, $m$ is the abscissa of the center of the ellipse in the second quadrant. After translating the value of $m$, the value of the major axis of the ellipse is $a=m+\frac{C}{2}$.

$R$ is the distance between the center of circle $A$ and any point on the circumference of the ellipse:

$$R^2=X^2+Y^2$$

(5)

The distance between the center of the circle $B$ and any point on the circumference of the ellipse is $D$:

$$D^2={(X+m)}^2+Y^2$$

(6)

From the three-sided relationship of the triangle, the angle $\angle BAC$ between any point $C$ on the circumference of the ellipse and the center $A$ and $B$ is as follows:

$$\angle BAC=\arccos (\frac{m^2+R^2-D^2}{2mR})$$

(7)

Then, the angles are different corresponding to different points on the circumference of the ellipse, so the angle $\alpha$ between the line of $AC$ and the real axis is as follows:

$$\alpha =\left\{\begin{array}{l}\mathit{\pi }-\angle BAC,\theta \in [0,\mathit{\pi }]\\ \mathit{\pi }+\angle BAC,\theta \in [\mathit{\pi },2\mathit{\pi }]\end{array}\right.$$

(8)

Combining the above Formulas (7) and (8), the argument $\alpha$ of the circumference of the ellipse is as follows:

$${\alpha }^{\delta }=\left\{\begin{array}{l}\mathit{\pi }-\arccos (\frac{m-a\cos \theta }{\sqrt{{(a\cos \theta -m)}^2+{(b\sin \theta )}^2}}),\theta \in [0,\mathit{\pi }]\\ \mathit{\pi }+\arccos (\frac{m-a\cos \theta }{\sqrt{{(a\cos \theta -m)}^2+{(b\sin \theta )}^2}}),\theta \in [\mathit{\pi },2\mathit{\pi }]\end{array}\right.$$

(9)

Inputting Formulas (5) and (9) into Formula (2), the parametric equation of the airfoil is as follows:

$$\begin{array}{l}\left\{\begin{array}{l}{X}_u^{\delta }=\left[\sqrt{{(acos\theta -m)}^2+{(b\sin \theta )}^2}+\frac{{(\frac{C}{2})}^2}{\sqrt{{(acos\theta -m)}^2+{(b\sin \theta )}^2}}\right]\cos {\alpha }_u^{\delta }\\ Y_u^{\delta }=\left[\sqrt{{(acos\theta -m)}^2+{(b\sin \theta )}^2}-\frac{{(\frac{C}{2})}^2}{\sqrt{{(acos\theta -m)}^2+{(b\sin \theta )}^2}}\right]\sin {\alpha }_u^{\delta }\end{array}\right.{\alpha }_u^{\delta }\in [0,\mathit{\pi }]\\ \left\{\begin{array}{l}{X}_l^{\delta }=\left[\sqrt{{(acos\theta -m)}^2+{(b\sin \theta )}^2}+\frac{{(\frac{C}{2})}^2}{\sqrt{{(acos\theta -m)}^2+{(b\sin \theta )}^2}}\right]\cos {\alpha }_l^{\delta }\\ Y_l^{\delta }=\left[\sqrt{{(acos\theta -m)}^2+{(b\sin \theta )}^2}-\frac{{(\frac{C}{2})}^2}{\sqrt{{(acos\theta -m)}^2+{(b\sin \theta )}^2}}\right]\sin {\alpha }_l^{\delta }\end{array}\right.{\alpha }_l^{\delta }\in (\mathit{\pi },2\mathit{\pi }]\end{array}$$

(10)

Among them is the upper airfoil line parameter equation with a variable relative thickness position when ${\alpha }_u^{\delta }\in [0,\mathit{\pi }]$, which is the lower airfoil line parameter equation with a variable relative thickness position when ${\alpha }_l^{\delta }\in [0,\mathit{\pi }]$.

In the same way, as shown in Figure 4b, after the ellipse translating the value of $m$, the minor axis of the ellipse is $b=m+\frac{C}{2}$, and the sectional function equation of the upper and lower airfoils with variable relative thickness positions can also be obtained.

Take the same value of $m$ and different minor axis values of $b$ to get the airfoils with different relative thicknesses and different relative thickness positions without the values of the camber. As shown in Figure 5a, when the value of $b$ is less than the value of $a$ and closer to the value of $a$, the relative thickness of the airfoil is closer to 25% of the chord, and when the value of $a$ is smaller, the relative thickness of the airfoil is closer to the leading edge. In Figure 5b, when the value of $b$ is greater than the value of $a$ and closer to the value of $a$, the relative thickness of the airfoil is closer to 25% of the chord, and if the value of $a$ is greater, the relative thickness of the airfoil is closer to the trailing edge.

3.2. Determination of Relative Camber Position

As shown in Figure 6a, the ellipse is translated along the imaginary axis by $n$ ($n$ is a positive real number); the ellipse parameter equation is as follows:

$$\left\{\begin{array}{l}{x}^{\text{'}}=c\cos \theta \\ y^{\text{'}}=d\sin \theta +n\end{array}\right.\theta \in [0,2\mathit{\pi }]$$

(11)

where $c$ is the major axis of the ellipse, $d$ is the minor axis of the ellipse, and $\theta$ is the parameter.

Define $d^2=n^2+{(\frac{C}{2})}^2$, then the counterclockwise rotation angle $\delta$ of the ellipse parameter equation is as follows:

$$\delta =\arctan \frac{C}{2n}$$

(12)

The parameter equation of the ellipse after rotation is as follows:

$$\left[\begin{array}{l}X\\ Y\end{array}\right]=\left[\begin{array}{l}\cos \delta -\sin \delta \\ \sin \delta \cos \delta \end{array}\right]\left[\begin{array}{l}{x}^{\text{'}}\\ y^{\text{'}}\end{array}\right]$$

(13)

The distance $R$ between the center $A$ and any point on the circumference of the rotating ellipse is as follows:

$$R^2=X^2+Y^2$$

(14)

The distance $D$ between any point on the circumference of the ellipse after rotation and the center $B$ is as follows:

$$D^2=X^2+{(Y-n)}^2$$

(15)

From the three-sided relationship of the triangle, the angle $\angle BAC$ between any point $C$ on the circumference of the ellipse and the center $A$ and $B$ is as follows:

$$\angle BAC=\arccos (\frac{n^2+R^2-D^2}{2nR})$$

(16)

The argument is different corresponding to different points on the circumference, so the angle $\alpha$ between $AC$ and the real axis is as follows:

$$\alpha =\left\{\begin{array}{l}\frac{\mathit{\pi }}{2}-\angle BAC,\theta \in [0,\frac{\mathit{\pi }}{2}]\\ \frac{\mathit{\pi }}{2}+\angle BAC,\theta \in (\frac{\mathit{\pi }}{2},\mathit{\pi }]\\ \frac{\mathit{\pi }}{2}+\angle BAC,\theta \in (\mathit{\pi },\frac{3\mathit{\pi }}{2}]\\ \frac{5\mathit{\pi }}{2}-\angle BAC,\theta \in (\frac{3\mathit{\pi }}{2},2\mathit{\pi }]\end{array}\right.$$

(17)

The argument $\alpha$ of the circle of the ellipse after rotation is as follows:

$${\alpha }^f=\left\{\begin{array}{l}\frac{\mathit{\pi }}{2}-\arccos (\frac{c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta }{\sqrt{{(c\cos \theta \cos \delta -d\sin \theta \sin \delta -n\sin \delta )}^2+{(c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta )}^2}}),\theta \in [0,\frac{\mathit{\pi }}{2}]\\ \frac{\mathit{\pi }}{2}+\arccos (\frac{c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta }{\sqrt{{(c\cos \theta \cos \delta -d\sin \theta \sin \delta -n\sin \delta )}^2+{(c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta )}^2}}),\theta \in (\frac{\mathit{\pi }}{2},\mathit{\pi }]\\ \frac{\mathit{\pi }}{2}+\arccos (\frac{c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta }{\sqrt{{(c\cos \theta \cos \delta -d\sin \theta \sin \delta -n\sin \delta )}^2+{(c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta )}^2}}),\theta \in (\mathit{\pi },\frac{3\mathit{\pi }}{2}]\\ \frac{5\mathit{\pi }}{2}-\arccos (\frac{c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta }{\sqrt{{(c\cos \theta \cos \delta -d\sin \theta \sin \delta -n\sin \delta )}^2+{(c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta )}^2}}),\theta \in (\frac{3\mathit{\pi }}{2},2\mathit{\pi }]\end{array}\right.$$

(18)

Inputting Formulas (14) and (17) into Formula (2), the parametric equation of the airfoil is as follows:

$$\begin{array}{l}\left\{\begin{array}{l}{X}_u^f=\left[\sqrt{{(c\cos \theta \cos \delta -d\sin \theta \sin \delta -n\sin \delta )}^2+{(c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta )}^2}+\frac{{(\frac{C}{2})}^2}{\sqrt{{(c\cos \theta \cos \delta -d\sin \theta \sin \delta -n\sin \delta )}^2+{(c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta )}^2}}\right]\cos {\alpha }_u^f\\ Y_u^f=\left[\sqrt{{(c\cos \theta \cos \delta -d\sin \theta \sin \delta -n\sin \delta )}^2+{(c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta )}^2}-\frac{{(\frac{C}{2})}^2}{\sqrt{{(c\cos \theta \cos \delta -d\sin \theta \sin \delta -n\sin \delta )}^2+{(c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta )}^2}}\right]\sin {\alpha }_u^f\end{array}\right.{\alpha }_u^f\in [0,\mathit{\pi }]\\ \left\{\begin{array}{l}{X}_l^f=\left[\sqrt{{(c\cos \theta \cos \delta -d\sin \theta \sin \delta -n\sin \delta )}^2+{(c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta )}^2}+\frac{{(\frac{C}{2})}^2}{\sqrt{{(c\cos \theta \cos \delta -d\sin \theta \sin \delta -n\sin \delta )}^2+{(c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta )}^2}}\right]\cos {\alpha }_l^f\\ Y_l^f=\left[\sqrt{{(c\cos \theta \cos \delta -d\sin \theta \sin \delta -n\sin \delta )}^2+{(c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta )}^2}-\frac{{(\frac{C}{2})}^2}{\sqrt{{(c\cos \theta \cos \delta -d\sin \theta \sin \delta -n\sin \delta )}^2+{(c\cos \theta \sin \delta +d\sin \theta \cos \delta +n\cos \delta )}^2}}\right]\sin {\alpha }_l^f\end{array}\right.{\alpha }_l^f\in (\mathit{\pi },2\mathit{\pi }]\end{array}$$

(19)

Among them is the upper airfoil line parameter equation with a variable relative thickness position when ${\alpha }_u^f\in [0,\mathit{\pi }]$, while is the lower airfoil line parameter equation with a variable relative thickness position when ${\alpha }_l^f\in [0,\mathit{\pi }]$.

Similarly, as shown in Figure 6b, after the ellipse shifting the value of $n$, the major axis of the ellipse is $c^2=n^2+{(\frac{C}{2})}^2$, and the sectional function equation of the upper and lower airfoil lines with variable relative camber positions can also be obtained.

As shown in Figure 7a, when taking the same value of $n=0.04$, the value of $c$ is greater than the value of $d$. When the value of $c$ is larger, the maximum camber position is closer to the trailing edge, and when the value of $c$ is smaller, the maximum camber position is closer to the 50% chord. In Figure 7b, when the value of $c$ is less than the value of $d$ and the value of $c$ is smaller, the maximum camber position is closer to the leading edge, if the value of $c$ is larger, the maximum camber position is closer to the 50% chord.

3.3. Influence of the Transformation Parameters on the Airfoil Shape

When the ellipse translates the same value of $m$ and $n$ along the real axis and the imaginary axis, the relative thickness position of the airfoil is forward and the relative camber position is backward, which can be seen in Figure 8a; in Figure 8b, the relative thickness of the airfoils at the back, and the relative camber is at the back; the relative thickness of the airfoils at the front and the relative camber is at the front in Figure 8c ; the relative thickness position of the airfoil is rearward, and the relative camber position is forward in Figure 8d. The function equations of the lines of different upper and lower airfoils can be combined to produce a richer airfoil structure, as shown in Figure 8e,f (that is, the upper and lower profiles can adopt different parameters).

4. Case Verification and Analysis

4.1. Fitting Existing Airfoils

The NACA63-212 airfoil is adopted to be fitted. The parameters of the geometric characteristics are the maximum thickness of 12%, the maximum thickness position of 35.1%, the maximum camber of 1.10%, and the maximum camber position of 50%. The original airfoil coordinates are connected sequentially with a smooth curve to form a lattice image, and the transformation parameters of the airfoil transformation sectional function are used to approximate the original airfoil lattice. According to the difference between the two graphics, the values of the transformation parameters can be adjusted to gradually approach the original airfoil lattice according to the change law of the influence of the transformation parameter on the airfoil shape. As shown in Figure 9a, the upper airfoil transformation parameters are m = 0.035, n = 0.02, a = 0.285, b = 0.275, c = 0.27, and d = 0.2508, and the lower airfoil transformation parameters are m = 0.042, n = 0.022, a = 0.2920, b = 0.3, c = 0.275, and d = 0.2510. The deviation value is used to judge the geometric convergence of the discrete coordinate points of the fitted airfoil and the original airfoil, as shown in Figure 9b. The maximum deviation of the fit is $1.22135\times {10}^{-3}$. Compared with the literature [25], as shown in

Table 1. Comparative analysis of the fitted airfoils.

Methods	Hicks-Henne	Bezier	Spline	B-Spline	This Paper
Max Deviation (10-3)	1.5652	9.2297	1.7590	6.4237	1.22135

, the new method is better than the other four methods, thus verifying the effectiveness of the model.

The XFOIL software is a commercial software used especially for wind turbine airfoil design and aerodynamic characteristics analysis. This paper used the XFOIL software to analyze the airfoil and test the aerodynamic characteristics of the fitted airfoil. Figure 10 shows the changes in the lift coefficient $Cl$ and drag coefficient $Cd$ with the angle of attack $\alpha$. At the same time, comparing the calculation results with the original airfoil under the same working condition (Re = 500,000), the fitting deviation of the lift coefficient is less than 0.03, and the fitting deviation of the drag coefficient is less than 0.0015. The aerodynamic characteristics of the fitted airfoil under the same working conditions show very little deviation from the original airfoil data, and the lift and drag characteristics of the fitted airfoil can be well-matched with the original airfoil. As shown in Figure 11, the pressure distribution of the original airfoil and the fitted airfoil is analyzed under the same working conditions, there is a deviation at the leading edge of the airfoil. The reason for the error is that there is a superposition process during the transformation process, and the imaginary part needs to be omitted. Moreover, the pressure distribution at the other positions is consistent. Therefore, the fitted airfoil can be used to replace the original airfoil for wind turbine blade design and analysis.

4.2. Generating a New Airfoil

By setting different transformation parameters, the geometric characteristic parameters of the airfoil can be fine-tuned to design a new airfoil. According the transformation parameters in

Table 2. Parameter values for generating new airfoil.

		m	n	a	b	c	d
Airfoil1	Upper1	0.03	0.04	0.28	0.295	0.28	0.2532
	Lower1
Airfoil2	Upper2	0.03	0.04	0.28	0.265	0.25	0.2532
	Lower2	0.03	0.04	0.28	0.27	0.24	0.2532

, two new airfoil shapes are generated in Figure 12. The geometric characteristics of airfoil 1 are as follows: the maximum thickness is 14.07%, the maximum thickness position is 39.4%, the maximum camber is 4.25%, and the maximum camber position is 49.8%. The geometric characteristics of airfoil 2 are as follows: the maximum thickness is 5.8%, the maximum thickness position is 29.8%, the maximum camber is 2.74%, and the maximum camber position is 52.4%.

The XFOIL software is also used to calculate the lift-drag characteristics and pressure distribution of the newly generated two airfoils at a Reynolds number of 500,000. As shown in Figure 13, the maximum lift coefficient of airfoil 1 is 1.17813 and the minimum drag coefficient is 0.00968; the maximum lift coefficient of airfoil 2 is 1.06386 and the minimum drag coefficient is 0.00487. Figure 14 shows the pressure distribution of the two newly generated airfoils at 3° and 8° angles of attack. At an angle of attack of 3°, the pressure distribution of airfoil 1 fluctuates at the chord length of 0.65–0.7, and airfoil 1 can also be further optimized. The pressure distribution of airfoil 2 is relatively smooth, and it can be directly selected according to the actual situation. At an angle of attack of 8°, the pressure distribution of airfoil 1 and airfoil 2 is normal, and there is no major fluctuation. It can be directly selected into the airfoil library for subsequent use, or it can be used as a basic airfoil to further optimize the required airfoil performance, and it can then be considered whether it should be selected.

5. Conclusions

(1) This paper proposes a new segmental expression function of wind turbine airfoil based on the Joukowski transformation, and the functional equation for generating airfoil is derived. It shows that the translational parameter values along the real and imaginary axes determine the relative thickness and camber of the airfoil, and the size of the major and minor axis determine the relative thickness and camber position of the airfoil.

(2) Taking the commonly used NACA63-212 airfoil as an example, the geometric and lift-drag characteristics and pressure distribution of the fitted and original airfoil are studied. The results show that the geometric deviation of the fitted airfoil is $1.22135\times {10}^{-3}$, the maximum deviation of the lift coefficient is 0.03, the maximum deviation of the drag coefficient is 0.00117, and the pressure distribution has slight fluctuations at the leading edge, which demonstrate the validity of the model.

(3) Two new airfoils are generated by setting different parameter values that can be analyzed, namely, the geometric characteristics, lift-drag characteristics of the new airfoils, and the pressure distribution at 3° and 8° angles of attack. The newly generated airfoil may have pressure distribution fluctuations, which can be used as a basic airfoil to be further optimized so as to achieve the desired performance of the airfoil shape. The airfoil with smooth pressure distribution can be directly added into the airfoil library to provide a reference for subsequent airfoil selection.

Future work: In this paper, the obtained airfoil shape is quite different when the setting parameter values are relatively similar. Hence, the next step will be to study the changing regular between the parameter value and the airfoil shape. The proposed segmented airfoil function will be further applied to the actual blade design to validate the effectiveness of the new methodology.