Cooperation of Thin-Airfoil Theory and Deep Learning for a Compact Airfoil Shape Parameterization

Abstract

An airfoil shape parameterization that can generate a compact design space is highly desirable in practice. In this paper, a compact airfoil parameterization is proposed by incorporating deep learning into the PAERO parameterization method based on the thin-airfoil theory. Following the PAERO parameterization, the mean camber line is represented by a number of aerodynamic performance parameters, which can be used to narrow down the design space according to the thin-airfoil theory. In order to further reduce the design space, the airfoil thickness distribution is represented by data-driven generative models, which are trained by the thickness distributions of existing airfoils. The trained models can automatically filter out the physically unreasonable airfoil shapes, resulting in a highly compact design space. The test results show that the proposed method is significantly more efficient and more robust than the widely used CST parameterization method for airfoil optimization.

1. Introduction

Airfoils play a critical role in the aerodynamic design of aircraft, wind turbines and gas turbines. Airfoil optimization design includes three components, namely airfoil parameterization, Computational Fluid Dynamics (CFD) simulation and optimization search. It is well known that during optimization, the number of CFD evaluations increases rapidly with the increase in design space, which is known as the “curse of dimensionality” [1,2]. Although the gradient-based methods such as the adjoint method [3,4] can make the number of CFD evaluations almost independent of design variables, the gradients are usually either not available or difficult to calculate. For the surrogate-based methods [5] and the gradient-free methods [6], an efficient parameterization method that can generate a compact design space is still highly desirable [7].

Traditionally, the airfoil parameterization methods can be classified into two categories. The first category consists of the direct construction methods such as parametric section (PARSEC) [8], class function/shape function transformation (CST) [9] and B-spline curves [10], and the second category consists of the perturbation deformation methods such as Hicks–Henne bump functions [11] and the free-form deformation (FFD) method [12]. In recent years, the data-driven methods have gradually been proved to be able to reduce the number of design variables in aerodynamic shape optimization. Masters et al. [13] studied a number of parameterization methods, concluding that data-driven singular value decomposition (SVD) parameterization can minimize the number of design variables representing airfoils. The idea of SVD is that a number of orthogonal modes can be extracted from the existing airfoil geometries, and then, the linear combination of these modes can be used to construct new airfoils. Li et al. [14] extended this idea by extracting the orthogonal modes from the airfoil thickness and the airfoil camber, respectively. More recently, Chen et al. [15] proposed Bézier Generative Adversarial Networks (Bézier-GANs) for airfoil shape parameterization. By learning from the existing UIUC [16] airfoil database, a Bézier-GAN can ensure a physically reasonable airfoil geometry, resulting in a compact design space for airfoil optimization. Du et al. [17,18] incorporated a surrogate model into a B-spline-based GAN to further speed up the airfoil design. In addition, Li et al. [19,20] adopted the UIUC airfoil database and a number of randomly deformed airfoils to train a discriminative model for unrealistic airfoils, resulting in a further speedup.

Although the above methods are proven to be useful in some applications, their design variables are usually independent of aerodynamic performances. In optimization design, aerodynamic performances often appear as goals or constraints. Therefore, incorporating aerodynamic performances into shape parameterization can provide extra guidance for optimal design. Recently, a number of studies introduced aerodynamic performances into shape parameterization by using deep learning technology. Achour et al. [21] trained a Conditional Generative Adversarial Network (CGAN) that can generate an airfoil geometry directly by a given lift–drag ratio and a given shape area. The results showed that about 75% of the generated airfoils maintained the given aerodynamic performances. Yilmaz et al. [22] adopted the airfoil stall performance and the drag pole curve as given design parameters to train a CGAN. Sekar et al. [23] proposed a deep learning model that can map a given pressure distribution to an airfoil geometry. Wang et al. [24] proposed an inverse design method for supercritical airfoils based on the Conditional Variational AutoEncoder (CVAE) and Wasserstein Generative Adversarial Networks (WGANs). This method first generates a Mach number distribution with given features, and then, the geometric coordinates are obtained from this Mach number distribution by a neural network model built for nonlinear mapping. In addition, Lei et al. [25] combined the generative model and the optimization algorithm to find a pressure coefficient distribution with given features and then mapped it to the corresponding airfoil geometry.

It is worth noting that the above trained models are often limited by the given training data, and they cannot provide a general shape parameterization method. Inspired by the classical thin-airfoil theory [26], Deng et al. [27] proposed a novel general parameterization method, named PAERO, which uses four aerodynamic performance parameters, namely the Mach number, the angle of attack, the lift coefficient and the pitching moment coefficient, and a number of optional parameters to represent the mean camber line of the airfoil. In the original method, the airfoil thickness distribution is represented by either the CST method or B-spline curves. In addition to being a general flexible shape parameterization, this method provides four aerodynamic performance parameters to manipulate the design space directly based on the thin-airfoil theory.

Following the previous work [26], this paper aims to extend the PAERO method by absorbing the advantages of deep learning technology. The airfoil thickness distribution will be represented by a data-driven generative model, which is built from the standard dataset of UIUC. The trained model can ensure a physically reasonable geometry due to its learning ability. The proposed parameterization method is expected to have the following advantages: (1) it can generate airfoil shapes with given aerodynamic performance parameters; (2) the aerodynamic performance parameters can be used as design variables, providing an efficient way to narrow down the design space; (3) the compact thickness representation can avoid a physically unreasonable geometry and reduce the design space as well.

2. Methodology

2.1. Problem Definition

The airfoil shape and coordinate system of interest are shown in Figure 1. It has a round nose with a radius equal to $r_{LE}$ and a trailing-edge thickness equal to $\Delta z_{TE}^{\left(t\right)}$. The coordinate system origin lies on its leading edge. Its chord line coincides with the x axis, with the chord length being c. In order to parameterize this airfoil, its shape was decomposed into two basic components, namely the mean camber line and the thickness distribution, denoted as $z^{\left(s\right)}\left(x\right)$ and $z^{\left(t\right)}\left(x\right)$, respectively. For brevity, the coordinates are normalized by the chord length, introduced as $X=\frac{x}{c}$, $Z^{\left(s\right)}\left(X\right)=\frac{z^{\left(s\right)}\left(X\right)}{c}$ and $Z^{\left(t\right)}\left(X\right)=\frac{z^{\left(t\right)}\left(X\right)}{c}$. Given the mean camber line and the thickness distribution, the upper surface of this airfoil can be written as

$$\frac{z_u}{c}=Z^{\left(s\right)}\left(X\right)+\frac{1}{2}{Z}^{\left(t\right)}\left(X\right),$$

(1)

and the lower surface of this airfoil can be written as

$$\frac{z_l}{c}=Z^{\left(s\right)}\left(X\right)-\frac{1}{2}{Z}^{\left(t\right)}\left(X\right).$$

(2)

The following sections will complete the defining equations for the mean camber line $Z^{\left(s\right)}\left(X\right)$ and the thickness distribution $Z^{\left(t\right)}\left(X\right)$.

2.2. Representation of Mean Camber Line

According to the thin-airfoil theory, the mean camber line of an airfoil is related to its aerodynamic performance parameters, including the lift coefficient and the pitching moment. The idea of the PAERO method [26] is to use these aerodynamic performance parameters as design variables. The mean camber line, $Z^{\left(s\right)}\left(X\right)$, can be expressed as follows:

$$Z^{\left(s\right)}\left(X\right)=\left(\alpha -A_0\right)X-\frac{1}{\pi }{{\displaystyle \int }}_0^1\left[{{\displaystyle \sum }}_{n=1}^N{A}_n\sin \left(n\times \mathrm{acos}\left(2X^{\prime }-1\right)\right)\right]ln\left|\frac{X-X^{\prime }}{X}\right|d{X}^{\prime },$$

(3)

where the first three coefficients, namely $A_0$, $A_1$ and $A_2$, can be calculated as follows:

$$A_0=3\left[\alpha -\frac{\Delta z_{TE}^{\left(s\right)}}{c}-\frac{1}{3\pi }{C}_L-\frac{4}{3\pi }{C}_M+{{\displaystyle \sum }}_{i=2}^{N/2}\frac{A_{2i}}{\left(2i-1\right)\left(2i+1\right)}\right],$$

(4)

$$A_1=\frac{1}{\pi }{C}_L-2A_0,$$

(5)

$$A_2=-\frac{1}{\pi }{C}_L-\frac{4}{\pi }{C}_M+2A_0,$$

(6)

where the angle of attack, $\alpha$, the lift coefficient, $C_L$, and the pitching moment coefficient around the quarter-chord point, $C_M$, are design variables. In addition, $\Delta z_{TE}^{\left(s\right)}$ is the height of the camber line on the trailing edge, and $c$ is the chord length.

When the freestream Mach number, denoted as $Ma$, is less than one, the Prandtl Glauert rule can be used to take the compressibility effect into account, and the effective aerodynamic coefficients are written as follows:

$$C_L=C_{L,Ma}\times \sqrt{1-M{a}^2},$$

(7)

$$C_M=C_{M,Ma}\times \sqrt{1-M{a}^2}.$$

(8)

Finally, the design variables for the mean camber line can be expressed as follows:

$$\mathit{D}{\mathit{V}}^{\left(s\right)}=\left(M{a}_{DV},{\alpha }_{DV},C_{L,DV},C_{M,DV},A_3,\dots \dots A_N,\Delta z_{TE}^{\left(s\right)}\right),$$

(9)

where the subscripts “DV” are added to avoid confusion with the actual aerodynamic performance parameters. The design Mach number, $M{a}_{DV}$, can be given according to the design requirements. Because the design lift coefficient, $C_{L,DV}$, depends on both the angle of attack, ${\alpha }_{DV}$, and the mean camber line, either $C_{L,DV}$ or ${\alpha }_{DV}$ can be fixed in the optimization design.

2.3. Representation of Airfoil Thickness Distribution

2.3.1. CST Method

The CST method is a widely used parameterization method. It combines a class function and a shape function. The surface coordinates of an airfoil with a trailing-edge thickness of $\Delta z_{TE}^{\left(t\right)}$ can be expressed as

$$Z^{\left(t\right)}\left(X\right)=2\times \left[\left(X^{0.5}\times \left(1-X\right)\right)\times S\left(X\right)+\frac{1}{2}\times X\times \frac{\Delta z_{TE}^{\left(t\right)}}{c}\right],$$

(10)

where $S\left(X\right)$ is the shape function defined as the linear combination of Bernstein polynomials, written as follows:

$$S\left(X\right)={{\displaystyle \sum }}_{i=0}^n{a}_i{B}_{i,n}\left(X\right),$$

(11)

$$B_{i,n}\left(X\right)=\frac{n!}{i!\left(n-i\right)!}\times X^i\times {\left(1-X\right)}^{n-i},$$

(12)

where $B_{i,n}\left(X\right)$ is a Bernstein polynomial of degree n.

2.3.2. SVD Method

The idea of SVD is to extract the orthogonal basic modes through a proper orthogonal decomposition of the chosen dataset. A matrix is first assembled by using the Z coordinates of the airfoil thickness distribution, while the corresponding X coordinates are kept the same for all airfoils. The matrix, H, is written as follows:

$$\mathit{H}=\left[\begin{array}{ccc}{Z}_{11}^{\left(t\right)}& \cdots & Z_{m1}^{\left(t\right)}\\ ⋮& \ddots & ⋮\\ Z_{1n}^{\left(t\right)}& \cdots & Z_{m1}^{\left(t\right)}\end{array}\right],$$

(13)

where $n$ is the number of coordinate points and $m$ is the number of samples. Then, this matrix is further decomposed by singular value decomposition as follows:

$$\mathit{H}=\mathit{U}\Sigma {\mathit{V}}^T,$$

(14)

where the columns of $\mathit{U}$, namely $\left(U_1,U_2,\cdots U_{\min \left(m,n\right)}\right)$, represent the orthogonal thickness modes. A new thickness distribution can be expressed as a linear combination of these modes, written as follows:

$$Z^{\left(t\right)}\left(X\right)={{\displaystyle \sum }}_{i=1}^k{a}_i{U}_i\left(X\right)+X\times \frac{\Delta z_{TE}^{\left(t\right)}}{c}.$$

(15)

Figure 2 shows the first ten thickness modes extracted from the UIUC dataset. Note that the number of samples, $m$, is equal to 1433 and the number of coordinate points, $n$, is equal to 121. The coordinate distribution is expressed as follows:

$$X_i=\frac{1}{2}\left(\cos \frac{2\pi \left(i+119\right)}{240}+1\right),i=1,2,\cdots ,121.$$

(16)

2.3.3. Deep Learning Method

The Generative Adversarial Network (GAN) [28] is a widely used model that can learn the main features of data and generate a dataset similar to the original data. The architecture of the standard GAN is shown in Figure 3a, where S is the training dataset.

It is composed of two players, namely the generator, G, and the discriminator, D. G obtains the random noise, $\mathit{w}$, from the random distribution, P_w, and generates the data, $G\left(\mathit{w}\right)$. D judges whether the input data come from the real distribution, P_data, or are generated by G. G constantly adjusts parameters to cheat D, making it unable to determine the source. The two players play games with each other and finally reach the Nash equilibrium, so that the generative distribution, P_G, of G is the same as the distribution of the real data, P_data. Generally, the training loss function of GAN can be expressed as follows:

$${\min }_G{\max }_D\mathrm{V}\left(D,G\right)={\mathbb{E}}_{s~P_{data}}\left[\log D\left(s\right)\right]+{\mathbb{E}}_{w~P_w}\left[\log \left(1-D\left(G\left(\mathit{w}\right)\right)\right)\right].$$

(17)

In fact, GAN is very easy to collapse during training, with the result that the generator, G, can only generate a small part of the training data to completely cheat the discriminator, D. Therefore, the WGAN [29] model was proposed to solve this problem by improving the activation function and the training process. The loss function of WGAN can be expressed as follows:

$${\min }_G{\max }_D\mathrm{V}\left(D,G\right)={\mathbb{E}}_{s~P_{data}}\left[D\left(s\right)\right]-{\mathbb{E}}_{w~P_w}\left[D\left(G\left(\mathit{w}\right)\right)\right],$$

(18)

$$\parallel \partial D/\partial \left(v^D,b^D\right)\parallel \le \gamma ,$$

(19)

where $v^D$ and $b^D$ are the training parameters of D and $\gamma$ is the weight gradient clipping boundary used to ensure the stability of model training.

Because the noise variables, $\mathit{w}$, of the standard GAN are uninterpretable, the InfoGAN [30] model was proposed to add a number of latent variables that are interpretable. The architecture of the InfoGAN model is shown in Figure 3b. Like GAN, this model also has a G and a D that play games with each other and ultimately enable the G to synthesize data similar to the training data. Unlike GAN, the input variables of G include the latent variables, $\mathit{u}$. InfoGAN adds the auxiliary distribution, Q, to restore the $\mathit{u}$ inputted to the G. Both Q and D can share the convolution layers. InfoGAN maximizes the mutual information, $I\left(\mathit{u};G\left(\mathit{w},\mathit{u}\right)\right)$, of the latent variables,$\mathit{u}$, and the generated samples, $G\left(\mathit{w},\mathit{u}\right)$. By training, the main features of the data can be extracted into $\mathit{u}$.

In practice, it is difficult to calculate the mutual information directly, so the lower bound, $L_1$, of the mutual information is generally calculated as follows:

$$L_1\left(G,D\right)={\mathbb{E}}_{s~G\left(w,u\right)}\left[{\mathbb{E}}_{u^{\prime }~P\left(u|s\right)}\left[\log Q\left({\mathit{u}}^{\prime }|s\right)\right]\right]+H\left(\mathit{u}\right)\le I\left(\mathit{u};G\left(\mathit{w},\mathit{u}\right)\right),$$

(20)

where $H\left(\mathit{u}\right)$ is the entropy of $\mathit{u}$. If the distribution of $\mathit{u}$ is fixed, $H\left(\mathit{u}\right)$ becomes a constant. Combined with the mutual information, the loss function for the InfoGAN training can be expressed as follows:

$${\min }_{G,Q}{\max }_D{V}_I\left(D,G,Q\right)={\mathbb{E}}_{s~P_{data}}\left[\log D\left(s\right)\right]+{\mathbb{E}}_{\mathit{u}~P_{\mathit{u}},\mathit{w}~P_{\mathit{w}}}\left[\log \left(1-D\left(G\left(\mathit{w},\mathit{u}\right)\right)\right)\right]-\lambda L_1\left(G,Q\right),$$

(21)

where $\lambda$ is introduced to balance the loss function and $P_{\mathit{u}}$ is the prior probability distribution of the latent variables.

The dataset S used for training the generative model is from the UIUC airfoil database. The thickness distribution of an airfoil can be regarded as showing a symmetric airfoil. The generator, G, obtains the latent variables, $\mathit{u}$, and the noise variables, $\mathit{w}$, from the random distribution and then synthesizes the airfoils. The discriminator, D, determines whether the airfoil is correct. Q attempts to restore the latent variables, ${\mathit{u}}^{\prime }$, by extracting the features. Eventually, G can synthesize a reasonable symmetric airfoil through the latent variables and the noise variables, and almost all geometric information is extracted into the latent variables. Therefore, the latent variables, $\mathit{u}$, of InfoGAN can be used for the airfoil thickness representation.

For a symmetric airfoil with a trailing-edge thickness of $\Delta z_{TE}^{\left(t\right)}$, its coordinates, $Z^{\left(t\right)}\left(X\right)$, can be expressed as follows:

$$Z^{\left(t\right)}\left(X\right)=G\left(X;\mathit{u}\right)+X\times \frac{\Delta z_{TE}^{\left(t\right)}}{c}.$$

(22)

Finally, the design variables for the airfoil thickness distribution can be expressed as follows:

$$\mathit{D}{\mathit{V}}^{\left(t\right)}=\left(\mathit{u},\Delta z_{TE}^{\left(t\right)}\right)=\left(u_0,u_1,u_2,\cdots \cdots u_N,\Delta z_{TE}^{\left(t\right)}\right).$$

(23)

2.4. Summary

The schematic diagram of the proposed airfoil parameterization is shown in Figure 4. The PAERO method is used to represent the mean camber line, and the generative model is used to represent the airfoil thickness distribution. The parameters used for the camber representation include two types, namely the aerodynamic performance parameters and the optional parameters. The parameters used for the thickness representation are the latent variables of the generator, G.

For an airfoil with a trailing-edge height of $\Delta z_{TE}^{\left(s\right)}$ and a trailing-edge thickness of $\Delta z_{TE}^{\left(t\right)}$, the design variables can be expressed as follows:

$$\mathit{D}\mathit{V}=\left(\mathit{D}{\mathit{V}}^{\left(s\right)},\mathit{D}{\mathit{V}}^{\left(t\right)}\right)=\left(M{a}_{DV},{\alpha }_{DV},C_{L,DV},C_{M,DV},A_3,\dots \dots A_N,\Delta z_{TE}^{\left(s\right)},u_0,u_1,u_2,\cdots \cdots u_N,\Delta z_{TE}^{\left(t\right)}\right).$$

(24)

Note that when the method is used to generate airfoil families, four aerodynamic performance parameters need to be specified; when it is used for aerodynamic shape optimization, only two aerodynamic performance parameters are needed, due to the fact that the lift coefficient is dependent on the angle of attack. The optional parameters for the mean camber line can be selected according to the design requirements. The number of latent variables needed to generate the airfoil thickness distribution will be discussed later.

In summary, the coordinates of the upper surface, $z_u$, and the coordinates of the lower surface, $z_l$, for an airfoil can be expressed as follows:

$$\frac{z_u}{c}=Z^{\left(s\right)}\left(X;\mathit{D}{\mathit{V}}^{\left(s\right)}\right)+\frac{1}{2}{Z}^{\left(t\right)}\left(X;\mathit{D}{\mathit{V}}^{\left(t\right)}\right),$$

(25)

$$\frac{z_l}{c}=Z^{\left(s\right)}\left(X;\mathit{D}{\mathit{V}}^{\left(s\right)}\right)-\frac{1}{2}{Z}^{\left(t\right)}\left(X;\mathit{D}{\mathit{V}}^{\left(t\right)}\right).$$

(26)

It should be pointed out that the chosen formulation does not exactly follow the theory of camber and thickness, as it neglects the X component of the normal to the camber line when adding the thickness. We adopt this common simplification to avoid negative X coordinates at the leading edge.

3. Results and Discussion

3.1. Model Training and Airfoil Thickness Synthesis

3.1.1. Data Preparation

The UIUC airfoil database includes most airfoils used in a wide range of scenarios, including aircraft, gas turbines, wind turbines, etc. Many studies [15,17,19,23] have adopted this database as the training basis for a neural network model. However, some of the airfoil data in the database are incomplete or impractical. Therefore, the airfoil database is filtered and normalized by using the method suggested by Li et al. [31]. The process is divided into the following four steps:

(1)

Remove the airfoils with incomplete data, and order the airfoil coordinates in a counterclockwise direction starting from the trailing edge, as shown in Figure 5.

(2)

Convert a blunt trailing edge into a sharp trailing edge, adjust the angle of attack to 0 degree, and synchronize the x coordinates.

(3)

Remove the airfoils with more than one local maximum thickness.

(4)

Remove the airfoil mean cambers.

Finally, 1433 symmetrical airfoils were selected from the UIUC database with a total of 1550 existing airfoils. Typical samples obtained from the database are shown in Figure 6. The database only stores the normalized Z coordinates of the airfoils, and the corresponding normalized X coordinate distribution is as follows:

$$X_i=\frac{1}{2}\left(\cos \frac{2\pi \left(i-1\right)}{240}+1\right),i=1,2,\cdots ,241.$$

(27)

3.1.2. Model Training

A symmetrical airfoil can be defined by a number of key geometric features, such as the leading edge radius, the airfoil thickness, and the location of the maximum thickness. An InfoGAN model is established to extract these interpretable features as design variables for airfoil parameterization. The schematic diagram of the InfoGAN model is shown in Figure 7.

The training set contains a total of 1433 symmetrical airfoils, each of which is represented as a vector in the format of 241 × 1. The latent variables, $\mathit{u}$, and the noise variables, $\mathit{w}$, are used as the input of the generator, G, where $\mathit{u}$ follows a uniform distribution in [0, 1], and $\mathit{w}$ follows a normal distribution with a mean of 0 and a variance of 1. G includes the multi-layer full connected layers and a mode layer. The modes come from the SVD decomposition of symmetric airfoils. The research by Li et al. [14] showed that the first 20 modes were already sufficient for fitting the airfoils. In addition, using the modes as the last layer of G is helpful for ensuring a smooth surface and a symmetrical geometry. The discriminator, D, is composed of a number of convolution layers and a full connected layer. The convolution layers extract the features of the airfoil, and the full connected layer integrates them for judgment. An independent full connected layer, Q, is set behind the convolution layers of D to extract the features of the generated airfoil and tries to restore the input information, $\mathit{u}$.

The parameters for the neural network were chosen by experience [32]. The final network parameters of the generator and the discriminator are shown in

Table 1. Neural network parameters for the generator and the discriminator.

Generator			Discriminator
Layer	Neurons	AF	Layer	Filters	AF
FC layer	128	ReLU	C layer	10	LeakyReLU
FC layer	128	ReLU	C layer	20	LeakyReLU
FC layer	64	ReLU	C layer	40	LeakyReLU
FC layer	20	ReLU	C layer	80	LeakyReLU
SVD mode layer			FC layer	1
			FC layer(Q)	3	Sigmoid

, where FC layer means a full connected layer and C layer means a convolution layer. The models were implemented in the open-source software PyTorch [33], and the final model parameters were determined by trial and error [32]. The parameters related to D’s convolution kernel are set as follows: kernel size = 7, step size = 2 and padding = 3. The training parameters are set as follows: batch size = 64, learning rate = 0.0001 and training round epochs = 500.

3.1.3. Symmetrical Airfoil Synthesis

The input variables of symmetrical airfoil generator include the latent variables, $\mathit{u}$, and the noise variables, $\mathit{w}$. Here, we chose $d_{latent}=3$ and $d_{noise}=10$ as an example to study the generated samples, as shown in Figure 8.

Figure 8a shows the results of airfoil synthesis by changing each latent variable, while the noise variables, $\mathit{w}$, remain 0. It can be seen that each component of $\mathit{u}$ is related to one of the geometric features. $u_1$ is related to the leading edge radius, $u_2$ is related to the airfoil thickness, and $u_3$ is related to the location of the maximum thickness. Figure 8b shows the results of airfoil synthesis by changing the first three components of $\mathit{w}$, while the latent variables, $\mathit{u}$, remain 0.5. It can be seen that there is no significant difference between the generated airfoils, indicating that the features of symmetric airfoils have been mainly extracted into the latent variables, $\mathit{u}$.

Figure 9 shows the symmetric airfoils randomly generated by the InfoGAN method and the traditional parameterization methods. The InfoGAN method adopts a generator with six latent variables. The CST method adopts a 5th-order polynomial, so that there are six design variables. The SVD method adopts the first six modes. The variable range of the generative model is [0, 1], and the variable range of CST and SVD comes from the fitting results of the dataset. It can be seen that all of the symmetric airfoils generated by the InfoGAN generator maintain the features of the existing airfoils, and a number of the symmetric airfoils generated by the CST and SVD methods are actually physically unreasonable.

In order to test the accuracy of the trained model, the existing airfoils in the database have been fitted. Fitting an existing airfoil can be turned into an unconstrained optimization problem. The latent variables of the generator are taken as the design variables, and Powell’s BOBYQA [6] algorithm is used to minimize the root mean square error (RMSE), which is expressed as follows:

$$RMSE=\sqrt{\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{\left(Z_{i,G}\left(X\right)-Z_{i,data}\left(X\right)\right)}^2},$$

(28)

where $n$, equal to 241, is the total number of the coordinate points, and $Z_{i,G}\left(x\right)$ and $Z_{i,data}\left(x\right)$ are the normalized Z coordinates of the generated airfoil and the airfoil to be fitted, respectively.

While keeping the noise variables, $\mathit{w}$, at 0, the optimization algorithm searches the design space spanned by the latent variables, $\mathit{u}$, to minimize the RMSE. The average fitting errors of all the symmetric airfoils in the database are shown in

Table 2. Errors of fitting symmetric airfoil database.

dlatent	$\mathbf{RMSE} (\times {10}^{-4})$
2	33.11
3	12.01
4	7.253
6	6.234
8 10	5.067 3.367

. It can be seen that when the number of latent variables is greater than 8, the RMSE is less than $5\times {10}^{-4}$, which is the tolerance for wind tunnel geometries [34].

Note that in some cases, the upper surface and the lower surface will cross each other near the trailing edge, as shown in Figure 10. This phenomenon was also observed by Chen et al. [15], who suggested adding a regularization term to the loss function to avoid this risk. We find that this problem can be easily solved by fitting the airfoil to a high-order CST polynomial, as shown in Figure 10.

3.2. Verification of Aerodynamic Performance Parameters

This section will test the accuracy of aerodynamic performance parameters. CFL3D [35], an open-source CFD solver based on multi-block structured grids, was used for high-fidelity aerodynamic performance calculations. The CFD solver was first validated by the wind tunnel test of the airfoil E387 [36]. The flow conditions are Ma = 0.1, α = 0° and Re = 0.2 × 10⁶. Note that the wind tunnel test data were measured by allowing a free transition. The Spalart–Allmaras one-equation turbulence model coupled with an e^N empirical method for the transition prediction was adopted. Three types of grids were used to calculate the aerodynamic performances. The farfield is 15 times the chord length, and the height of the first layer is set to 1 × 10⁻⁵ to ensure y+ < 1.0. As shown in

Table 3. Grid sensitivity study for the E387 airfoil (Ma = 0.1, α = 0° and Re = 0.2 × 10⁶).

No.	Grid	CL	CD	CM
1	$143\times 49$	0.3773	0.01012	−0.07740
2	$257\times 97$	0.3817	0.009681	−0.07792
3	$451\times 97$	0.3876	0.009594	−0.07860

, the results of grid No. 2 are close to those of the densest grid, No. 3. Considering the computational cost, grid No. 2, shown in Figure 11, was adopted in the simulations. The comparison of the numerical simulation and the wind tunnel test is shown in Figure 12, indicating a satisfying accuracy.

The airfoil performance parameters at a low Mach number were first tested. The aerodynamic performance parameters were set as follows: $M{a}_{DV}=0.2$, ${\mathsf{\alpha }}_{DV}=2.0$, $C_{L,DV}=0.40$ and $C_{M,DV}=-0.05$. The optional parameters are not considered here. The generator, G, synthesized the thickness distributions based on the latent variables, $\mathit{u}$, of a random distribution. The dimension of $\mathit{u}$ is three, and the value range is [0, 1]. High-fidelity CFD simulations were used to verify the accuracy of the aerodynamic performance parameters. The simulation conditions were $Ma=0.2$, $\mathsf{\alpha }=2.0$ and $Re=4\times {10}^6$. The simulations were performed assuming a fully turbulent flow. A total of 50 airfoils were randomly generated, and the test results are shown in Figure 13. The given aerodynamic performance parameters are represented by a star, and the airfoils are displayed on their corresponding aerodynamic performance coordinates. It can be seen that due to the viscous effect, $C_L<C_{L,DV}$ and $C_M>C_{M,DV}$. However, the aerodynamic performances of most airfoils are close to the given parameters.

Then, the airfoil performance parameters at a high Mach number were tested. The aerodynamic performance parameters were set as follows: $M{a}_{DV}=0.6$, ${\mathsf{\alpha }}_{DV}=2.0$, $C_{L,DV}=0.50$ and $C_{M,DV}=-0.07$. In addition, the optional parameters, A₃ and A₄, were considered as well. Following the work of [26], the recommended range of A₃ and A₄ is [−0.02, 0.02]. The flow simulation conditions were $Ma=0.6$, $\mathsf{\alpha }=2.0$ and $Re=6.5\times {10}^6$. The simulation results of 50 randomly generated airfoils are shown in Figure 14. It can be seen that for most of the generated airfoils, $C_L<C_{L,DV}$, $C_M>C_{M,DV}$, and their aerodynamic performances are still close to the given parameters. In summary, it can be concluded that the aerodynamic performance parameters can significantly determine the aerodynamic performances of the airfoil family even at a high Mach number. In contrast, the optional parameter A_i does not significantly affect the aerodynamic performances.

3.3. Airfoil Shape Optimization

In this section, the capabilities of the proposed parameterization method will be tested in the context of airfoil shape optimization. Three parameterization methods are adopted for comparison, namely the combination of PAERO and InfoGAN, the combination of PAERO and CST, and CST for the entire airfoil surface. Note that the PAERO method is only adopted for the camber representation, the InfoGAN generative model is only adopted for the thickness representation, and the CST method is a general parameterization method, which can be used for representing either the airfoil thickness distribution or the airfoil surface.

The problem of airfoil drag minimization is considered here. The objective function, the design variables and the constraints are described in

Table 4. Statement of airfoil optimization for drag minimization.

	Functions	Description	Quantity
Minimize	CD	Drag coefficient	1
With respect to	DV α	Design variables Angle of attack	12 1
Subject to	CL = 0.5 CM ≥ −0.07 area ≥ 0.07	Lift coefficient Quarter-chord pitching moment Area constraint	1 1 1

. The design conditions are $Ma=0.6$, $C_L=0.5$ and $Re=6.5\times {10}^6$. The chosen optimization algorithm is Powell’s BOBYQA [6], a gradient-free algorithm based on quadratic approximation. The convergence is thought to be reached when the change in the objective function is less than a tolerance of $1\times {10}^{-5}$. The maximum number of iterations is set to 500. The constraints are as follows: the airfoil area is not less than 0.07, and the pitching moment coefficient around the quarter-chord point is not less than −0.07. The aerodynamic constraints and the geometric constraints are added to the objective function as the penalty functions, so the objective function is expressed as follows:

$$\mathrm{Obj}={\lambda }_0{C}_D+{\lambda }_{1}max\left(1+\frac{0.07}{C_M},0\right)+{\lambda }_{2}max\left(1-\frac{aera}{0.07},0\right),{\lambda }_0=100,{\lambda }_{1,2}=10,$$

(29)

where ${\lambda }_0$, ${\lambda }_1$ and ${\lambda }_2$ are the Lagrange multipliers.

The value range of design variables is shown in

Table 5. Design variables and their value ranges.

Variable Type	Symbol	Value Range
Camber representation	$M{a}_{DV}$ $C_{L,DV}$ ${\mathsf{\alpha }}_{DV}$ $C_{M,DV}$ $A_3,A_4,A_5,A_6$	0.6 0.5 [1.8, 2.2] [−0.06, −0.08] [−0.02, 0.02]
Thickness representation (InfoGAN)	$u_0,u_1,\cdots ,u_n$	[0, 1]
Thickness representation (CST)	$a_0,a_1,\cdots ,a_n$	[0.01, 0.4]

. $M{a}_{DV}$ and $C_{L,DV}$ are fixed, and ${\mathsf{\alpha }}_{DV}$ and $C_{M,DV}$ are allowed to be changed. The optional parameters, namely $A_3$, $A_4$, $A_5$ and $A_6$, are added as well. The thickness representation maintains the same degree of freedom as the camber representation. The thickness distribution is generated by either the generative model controlled by six latent variables or the CST method with a 5th-order polynomial. The value range for the CST method is derived from the fitting results of the UIUC symmetric airfoil database, and all design variables are assumed to be positive. For representing the entire airfoil surface, the CST method adopts two 5th-order polynomials for the upper surface and the lower surface, respectively, and the value ranges come from the fitting results of the UIUC airfoil database as well. The value range for the upper surface is [0.01, 0.4], and the value range for the lower surface is [−0.4, −0.01].

The optimization tests were repeated 10 times for each method by using random initial design variables. The convergence histories are shown in Figure 15. The shaded regions in the figure show the convergence histories of different tests, and the solid lines are the average convergence histories. It can be seen from the convergence histories that the method of using the InfoGAN model for thickness representation has the best performance. It converged near the minimum point after about 30 iterations, far fewer than the other two methods. Note that the large initial objective function shown in the figure is due to the fact that the area of the initial shape does not satisfy the design requirements and a large penalty function is imposed. The final designs of airfoils are shown in Figure 16. The optimized results of the InfoGAN model maintain the features of the existing airfoils and are all physically reasonable. On the contrary, the optimized results of the CST method are divergent, indicating a less compact and bumpier design space. Thus, it can be concluded that the proposed method is significantly more efficient and more robust than the traditional CST method.

4. Conclusions

In this paper, a compact airfoil shape parameterization is proposed by adopting the recent PAERO method for the camber representation and the data-driven InfoGAN model for the thickness representation. The InfoGAN model shows an excellent ability to learn the thickness distributions of existing airfoils. Its generator can compactly represent the thickness distribution and relate some of the key geometric features to the latent variables. Moreover, it can filter the unrealistic thickness distribution, resulting in a highly compact design space for airfoil design. In addition, it is confirmed that the aerodynamic performance parameters can significantly determine the aerodynamic performances of the airfoil family at subsonic flows and the optional parameters do not significantly affect the aerodynamic performances. In the test of airfoil shape optimization, it is found that the proposed method is significantly more efficient and more robust than the traditional CST method.