Aerodynamic Shape Optimization with Grassmannian Shape Parameterization Method

Abstract

The conventional method of optimizing the aerodynamic performance of an airfoil heavily depends on the confines of the design space. The design variables create a non-normalized space that is fragmented into several different clusters of airfoils. An approach that is data-driven and deforms airfoils over a Grassmannian submanifold is utilized in the work that is being presented here. The affine deformation, which includes camber and thickness, can be uncoupled from the method that is currently in use, and the operations that are performed on the airfoil shape can be made smooth enough to prevent unreasonable shapes from being produced. The CST method is also a part of the current study so that a comparison can be made between the two. A new method to describe the airfoil geometries over the Grassmannian space was generated using a dataset that contained 7007 different shapes of airfoils. These two methods are used to parameterize the subsonic (NACA0012) and transonic (RAE2822) airfoils, and the new method cuts the number of design variables from twelve to six, resulting in a reduction in overall complexity. The findings demonstrate that the new method maintains a high degree of consistency regardless of the flow conditions.

1. Introduction

As the fundamental element of an aerodynamic device, the airfoil performs a crucial role. Many techniques, such as inverse design [1] and optimization, have been utilized to enhance the performance of aircraft and wind turbines, which is mostly determined by the shape associated with the pressure distribution. The first method necessitates a great deal of skill on the part of the designer, as the pressure distribution must correspond to the intended lift, drag, and moment. Due to the limitations of the designers’ cognition and the stringent numerical constraints required for a convergent and smooth airfoil shape [2], it is easy to disregard the new arrangement. Thus, attempts at airfoil design employing the boundary-layer methodology such as XFOIL [3] or Reynolds-averaged Navier–Stokes such as CFL3D [4], paired with gradient-based [5,6,7,8,9,10] or gradient-free [11,12,13,14,15,16] searching approaches, have emerged as a potent tool for facilitating an unskilled and inexperienced design process. In spite of the fact that the primary difference between gradient-based and gradient-free approaches is whether the optimal global position can be reached, the two approaches share a parameterization strategy. The numerical-simulation-based optimization typically consists of four steps: a specific mathematical function to parameterize the airfoil geometry; an adaptive mesh transformation to export an acceptable domain discretization; a low- or high-fidelity solver for numerical flow solution; and an optimum searching method. In most cases, the aerodynamic property of an airfoil can be affected by the shape parameterization, which will be the subject of this work.

The optimal solver typically considers the airfoil to be a formula with constraint conditions. A change in the input variables, also called design variables, can be used to change the airfoil profile and deform sensitive sections such as the leading edge, chamber, and thickness. Clearly, the design space’s range is proportional to the number of design variables. Two factors determine the precision of the parameterization method: minimizing the error between the parameterized airfoil and the original airfoil, accurately expressing the airfoil form with a limited number of design variables, and preserving the smoothness of the profile curves. Many methods for representing an airfoil have been developed. The parametric section (PARSEC) [17] uses a sixth-order polynomial to approximate each section of the airfoil, and the good news is that each of the twelve supplied parameters corresponds exactly to the geometry. The downside is that faithfulness cannot be supplemented with variety or adaptability [18]. Directly disturbing the mesh point is an alternative method that eliminates the need for mathematical formulas [5]. This, however, may result in a rough profile for mesh creation and, as a result, a failed solution. Hicks and Henne’s concept is to create a new airfoil profile by implementing the collection of single-signed sine functions to add a bump effect to a basic airfoil [19]. Using this strategy, numerous researchers have successfully optimized aerodynamic performance [20,21,22,23,24]. However, the sine functions at the locations where the leading and trailing edges are perturbed are weak. Kulfan proposes the class-shape function transformation (CST), which is defined by the Bernstein polynomial, to represent the airfoil [25]. The CST method is widely used in the design and application of aerodynamic devices [26,27,28,29].

According to Master et al. [18], who conducted a research case on transonic airfoil, a one-drag or lift count of aerodynamic convergence requires a number of design variables that is between 38 and 66. Meanwhile, the design space, which consists of the set of design variables, is very large, and most of the exploration by the optimal solver is wasteful and time-consuming due to numerical evaluation iteration. Unfortunately, the majority of parameterization methods pay less attention to this issue than they do to the method’s precision, and the ’curse of dimensionality’ will occur in the global optimization approaches for aerodynamic design. The optimized process will be hindered if the large design space and noisy landscapes are not reduced and filtered, respectively, [30].

Dimension reduction (DR) is a highly effective technique for condensing the design space, such as through proper orthogonal decomposition (POD). According to Sripawadku et al. [31], each set of design variables should correspond to a single airfoil shape, but the non-orthogonal design space may result in inadequate design space coverage. Using POD, specifically singular value decomposition (SVD), to rotate and scale the design space and derive the smallest number of shape degrees of freedom can be a straightforward method for identifying the inherent property of geometry data. Poole et al. [32,33] analyzed a several airfoil data with POD in order to derive generic airfoil design variables. Master et al. [18] presented a comparison between seven types of parameterized method and concluded that POD can give better results. Grey et al. [34] discovered the active subspace method (ASM), in which the internal method is still the POD approach, for reducing the design space. Li et al. [35] presented an efficient approach to applying ASM in high-dimensional aerodynamic shape optimization problems and validate this optimization approach in an ONERA M6 wing optimization case with 220 shape variables. Machine learning is shown to be powerful in solving the issues of conventional parameterization, such as Modal parameterization based on geometric anomalies, and almost avoids the curse of a dimensionality issue in wing shape optimization, since it reduces the number of design variables from hundreds to tens with little sacrifice in design effectiveness [36,37].

Another issue to be aware of is the smoothness of the new generated airfoil shape, as seen in Li et al. [35] who used the Laplacian smoothing algorithm to provide a reasonable pressure distribution. The work provided in this paper expands the use of POD to the aerodynamic optimization of airfoil by analyzing the parameterized technique from the perspective of the flow manifold. Doronina et al. [38] provided a novel representation of forms that is specified in a submanifold of the Grassmannian. The research demonstrates a number of benefits, including an enhanced low-dimensional parameter domain for inferential statistics informing design and manufacturing. The objective of this study was to apply this parameterization to a variety of real aerodynamic applications.

2. Shape Deformation on Discrete Section of Airfoil

As a subset of a two-dimensional Euclidean space, the airfoil shape $\mathcal{F}$ can be expressed as the boundary of the subset of a two-dimension real number system ${\mathbb{R}}^2$, which is to say that $\partial \mathcal{F}\subset \mathbb{R}$. The lower and upper curvilinear equations of the airfoil, $s_U$ and $s_L$, are mapping to a normalized spatial range $\mathit{l}\subset [0,1]$ from the leading to trailing edge, commonly. The airfoil shape $\partial \mathcal{F}$ can be expressed with a combination of the basis function, ${\varphi }_j:\mathbb{R}\to \mathbb{R}$, such that, for s as one of $s_U$ or $s_L$,

$$s\left(\mathit{l};r_1,r_2,\mathit{x}\right)=c_f\left(\mathit{l};r_1,r_2\right)s_f\left(\mathit{l};\mathit{x}\right)$$

(1)

where $r_1$ and $r_2$ are the class parameters with values of 0.5 and 1.0, respectively, for a round-nose and aft-end airfoil, and the vector $\mathit{x}\subset {\mathbb{R}}^m$ includes the polynomial coefficients, which corresponds to $\mathit{m}$ design variables(parameterized variables). The class function is given by:

$$c_f\left(\mathit{l},\mathit{r}\right)={\mathit{l}}^{r_1}{\left(1-\mathit{l}\right)}^{r_2}$$

(2)

Meanwhile, the shape function is a polynomial of degree $\mathit{m}-1$ in $\mathit{l}$:

$$s_f\left(\mathit{l},\mathit{x}\right)=\sum _{j=0}^{m-1}{x}_j\varphi {\left(\mathit{l}\right)}^j,\mathit{x}={\left[x_0,\cdots ,x_{m-1}\right]}^T$$

(3)

where ${\varphi }_j\left(\mathit{l}\right)$ is the component shape function and is represented by a Bernstein polynomial. The component shape function is written as:

$${\varphi }_j\left(\mathit{l}\right)=K_j^m{\mathit{l}}^m{(1-\mathit{l})}^{m-j}$$

(4)

where $\mathit{K}$ is the binomial coefficient, and defined as:

$$K_j^m=\frac{m!}{j!-(m-j)!}$$

(5)

Assuming that $t=\sqrt{\mathit{l}}$, the airfoil shape parameterized with the CST method can be written as:

$$s(t;\mathit{x})=x_{0}t-x_{m-1}{t}^{2m-1}+\sum _{j=0}^{m-1}(x_{j+1}-x_j)t^{2j+1}$$

(6)

From Equation (3), the control variables which are listed in vector $\mathit{x}$ do not have physical meanings on airfoil geometry. On the other hand, the polynomial in t is what primarily disturbs the control variables, causing the airfoil to change shape.

Commonly, the design of an airfoil starts from an existing base, such as SD7003, a low-Reynolds-number airfoil. The optimized method searches the optimal solution within the independent variables ranges. In the other words, the potential assumption is that the design space included the min–max value of the quality of interest (QOI). Several airfoils were described by the CST method, and the components of input vector $\mathit{x}$ were listed in

Table 1. The components of vector $\mathit{x}$ corresponding to the specific airfoils.

Airfoil	${\mathit{x}}_1$	${\mathit{x}}_2$	${\mathit{x}}_3$	${\mathit{x}}_4$	${\mathit{x}}_5$
	Lower Surface
SD7003	−0.1113	−0.0605	−0.0611	−0.0561	0.0217
NLF0215F	−0.1147	−0.0424	−0.1336	−0.1815	0.3003
RAE2822	−0.1305	−0.1309	−0.2340	−0.0665	0.0325
NLR7301	−0.2535	−0.0985	−0.3397	−0.1388	0.1041
SC1095	−0.1444	−0.0722	−0.1158	−0.0858	−0.0929
NPL9626	−0.1405	−0.1388	−0.1419	−0.1111	−0.1558
OA212	−0.1678	−0.0271	−0.1266	−0.1237	−0.1613
Airfoil	${\mathit{x}}_{\mathbf{6}}$	${\mathit{x}}_{\mathbf{7}}$	${\mathit{x}}_{\mathbf{8}}$	${\mathit{x}}_{\mathbf{9}}$	${\mathit{x}}_{\mathbf{10}}$
	Upper Surface
SD7003	0.1585	0.1642	0.1203	0.1240	0.1092
NLF0215F	0.2278	0.2581	0.4142	0.1965	0.3958
RAE2822	0.1277	0.1404	0.1877	0.1968	0.1990
NLR7301	0.3120	0.1051	0.3163	0.2542	0.2345
SC1095	0.1797	0.1334	0.1318	0.1417	0.1245
NPL9626	0.1708	0.1849	0.1763	0.1268	0.1840
OA212	0.1890	0.2806	0.1116	0.3548	0.0038

, and the airfoil profiles are plotted in Figure 1a. The number of the components is 10.

Obtaining some rules from the distribution and modification of the components for a better QOI with a general rule is extremely difficult. In order to demonstrate a generalized condition for the design space given by the CST method, Figure 1b depicts the variation of two of ten CST parameterized variables for seven airfoils listed in Table 1. Each colored square represents a region where the design variables are feasible. Clearly, the slices of different airfoils cannot exactly overlap, and a highly fragmented space for design variables is demonstrated.

For example, suppose the design baseline is an SD7003 airfoil and a 50% perturbation is performed on the baseline values, and the design target is a supercritical airfoil, such as RAE2822, which requires a flat upper surface and a bent lower surface to compensate for the lift loss caused by the nose’s low curvature. Figure 2 shows that the majority of RAE2822 lies outside of the design space, in the blush region.

Therefore, the issue raised above makes us aware that the target could be missed if the design variables are confined inside inappropriate ranges. One would believe that if the bound values were large enough, the issue could be resolved. Although expanding the bound values or adding more design variables would undoubtedly provide a closer approximation to the theoretical ideal values, the cost of searching, which is calculated as the time required to traverse the whole design space, would increase exponentially. Doronina et al. [38] treated the airfoil parameterization as a data-driven problem, and they separated the variables which control the deformation of airfoil. There are two distinct types of airfoil deformation: (1) The affine deformation, which has four structurally significant deformations to the shape of the airfoil, including changes in thickness, camber, chord length, and twist; (2) high-order deformation, which is a common operation provided by the CST and Hicks–Henne methods and includes undulations on the precise positions. As a linear transformation, the affine deformation can be written as:

$$\tilde{s}=M^{T}s\left(\mathit{x}\right)+\mathit{b}$$

(7)

where $M\in G{L}_2$, and $\mathit{b}\in {\mathbb{R}}^2$. M stands for an element from the set of all invertible $2\times 2$ matrices. Commonly, the airfoil shapes are represented with a discrete form, which can be regarded as a matrix $A\in {\mathbb{R}}^{n\times 2}$. That is, a real airfoil shape can be transformed with $AM+1_{n,2}diag\left(\mathit{b}\right)$, where $1_{n,2}$ stands for a $n\times 2$ matrix filled with one. The matrix $M\left(\eta \right)$ for $\eta \in (0,1)$ can be represented as

$$\begin{array}{cc}\hfill & \left(1\right)M\left(\eta \right)\triangleq \left[\begin{array}{cc}1& 0\\ 0& \eta \end{array}\right],\left(2\right)M\left(\eta \right)\triangleq \left[\begin{array}{cc}(1-\eta )& 0\\ 0& \eta \end{array}\right],\hfill \\ \hfill & \left(3\right)M\left(\eta \right)\triangleq \left[\begin{array}{cc}\eta & 0\\ 0& 1\end{array}\right],\left(4\right)M\left(\eta \right)\triangleq \left[\begin{array}{cc}cos\left(\frac{\eta \pi }{2}\right)& -sin\left(\frac{\eta \pi }{2}\right)\\ sin\left(\frac{\eta \pi }{2}\right)& cos\left(\frac{\eta \pi }{2}\right)\end{array}\right],\hfill \end{array}$$

(8)

respectively.

Obviously, the affine deformation can influence the aerodynamic performance of an airfoil. Due to structural constraints, however, these changes are typically not taken into account. The perturbation on the airfoil induced by the mathematical function, such as the CST and Hicks–Henne methods, has become an very popular means of improving the airfoil aerodynamic capability under some strong constraints, such as twist angle. These adjustments frequently result in sags and crests on the airfoil, which can cause performance variations. The Grassmannian space, denoted as $\mathcal{G}(n,q)$, is the subspace of all q-dimensional subspaces of ${\mathbb{R}}^n$, and it is invariant under $G{L}_q$ transformations. In other words, if $\tilde{A}\in {\mathbb{R}}_{*}^{n\times q}$ stand for a full-rank representative element of the equivalence class $[\tilde{A}]\in \mathcal{G}(n,q)$ of all matrices with equivalent column span, the ever element of Grassmannian is a full-rank matrix modulo $G{L}_q$ deformations. Therefore, the affine deformations cannot affect the deformation over $\mathcal{G}(n,q)$.

As is well-known, the Grassmannian can be viewed as a quotient topology of orthogonal subgroups, meaning that every element in Grassmannian space has sample covariance proportional to identity matrix [39]. Subsequently, the way given by Doronina et al. [38] to map the physical airfoil shapes, A, to the element in Grassmannian space, $\tilde{A}$, is landmark-affine standardization (LAS) [40]. The LAS method normalizes the airfoil shapes (profiles) in order to obtain a zero mean and a sample covariance proportional to identical matrix. For example, an airfoil shape $A\in {\mathbb{R}}_{*}^{n\times 2}$ has n landmarks, and a thin singular SVD is carried on the shape vector as

$${(A-B\left(A\right))}^{\top }=U\Sigma V^{\top },$$

(9)

where $B\left(A\right)=1_{n,2}diag\left(\right(1/n\left)A^{\top }{1}_{n,1}\right)$ is the discrete center of mass, which is corresponding to every row of A. Let $\tilde{M}=\Sigma U^{\top }$ denote a $2\times 2$ invertible matrix. The mapping between physical shape A and its corresponding part $\tilde{A}(=V)$ over the Grassmannian is the affine transformation, that is,

$$A-B\left(A\right)=\tilde{A}\tilde{M}.$$

(10)

It is very easy to find that ${\tilde{A}}^{\top }\tilde{A}=I_2$, which has been described as a basic feature of Grassmannian space, and we have $[\tilde{A}]\in \mathcal{G}(n,2)$. Additionally, the physical shape can be recovered with a linear operation via Equation (10). The advantage of the airfoil design space being sampled into the Grassmannian space is demonstrated by Figure 3.

The shadow in Figure 3 consists of the lower and upper bounds of 7007 airfoils, and this airfoil cluster comes from a thousand samples with a $\pm 20\%$ perturbation over the seven baseline airfoils. As shown in Figure 3a, physical airfoils, which are generally specified in Euclidean space, have some discontinuous sections, particularly around the trailing point. This corresponds to the disjoint regions shown in Figure 1. Moreover, the mean shape from the airfoil cluster also departs from the center line of the shadow region. After being aggregated with the LAS method, the scale and translation invariance are brought into the space. It is abundantly clear that the airfoil cluster over the Grassmannian space has normal coordinates, i.e., the maximum and minimum values are similar in different directions. Another advantage is that the varying space has been compacted around the mean shape. A very important step for the LAS method is the Procrustes Analysis, with which can optimally superimpose a set of objects, instead of superimposing them to an arbitrarily selected shape.

The next problem is the principal component analysis (PCA) over the design space to separate out the main elements which can have a large effect on the deformation of airfoils. Because Grassmannian is a curved space, which belongs to Riemannian space, the PCA cannot be directly used on the new dataset. Fletcher et al. [41] presented an idea to apply PCA on manifolds, such as the Grassmannian manifold. The first step is to obtain the intrinsic mean, which is defined as:

$$\mu =\underset{\tilde{A}\in \mathcal{G}(n,2)}{argmin}\sum _{i=1}^{n}d{(\tilde{A},{\tilde{A}}_i)}^2,$$

(11)

where $\mathcal{G}(n,2)$ is a Grassmannian manifold, $\tilde{A}$ stands for an arbitrary element in a collection of points ${\tilde{A}}_1,\cdots ,{\tilde{A}}_n\in \mathcal{G}(n,2)$ and $d(·,·)$ denotes the Riemannian distance on $\mathcal{G}(n,2)$. The way we obtained the intrinsic mean is given by Karcher [42], who treats the computation of the intrinsic mean as solving the minimization problem with the gradient descent algorithm. Defining the cost function as

$$f(\tilde{A})=\frac{1}{2n}\sum _{i-1}^{n}d{(\tilde{A},{\tilde{A}}_i)}^2,$$

(12)

and the target for us is to find a point a to satisfy $f(\tilde{A})$ which has a minimum value, i.e., zero. With a size step $\tau$ and an previous estimation intrinsic mean ${\mu }_j$, the next estimation is

$${\mu }_{j+1}={\mathbf{Exp}}_{{\mu }_j}\left(\frac{\tau }{N}\sum _{i=1}^n{\mathbf{Log}}_{{\mu }_j}({\tilde{A}}_i)\right),$$

(13)

where $\mathbf{Exp}$ and $\mathbf{Log}$ are the Riemannian exponential map and logarithmic map, respectively. The principal geodesic analysis (PGA) method is to find a sequence of nested geodesic sub-manifolds that have the maximum projected variance of the data, and it is necessary to find the tangent space of $\mathcal{G}(n,2)$, which is always denoted by ${\mathcal{T}}_{\mu }\mathcal{G}(n,2)$ with an intrinsic mean $\mu$. Then, constructing a matrix as

$$\mathbf{H}=\frac{1}{\sqrt{n-1}}\left({\mathbf{Log}}_{\mu }({\tilde{A}}_1),\cdots ,{\mathbf{Log}}_{\mu }({\tilde{A}}_N)\right),$$

(14)

and computing the reduced SVD on $\mathbf{H}$ up to $r\le dim\left(\mathcal{G}\right)$. The parameter r is the number of design variables corresponds to the first r modes which needs to be held to capture the most variability, e.g., 95%. The new basis for a subspace of the tangent space can be provided by Equation (14). Meanwhile, a neighborhood of zero is defined as ${\mathbf{U}}_r\subseteq {\mathcal{T}}_{\mu }\mathcal{G}(n,2)$, and hence the projection is well-defined for all geodesic submanifolds. As in the expression given by Doronina et al. [38], a compact space $\mathcal{U}$ which has the relation $\mathcal{U}\subset {\mathbb{R}}^r$ can be used to find a normal coordinates $\mathbf{t}\in \mathcal{U}$. Then, a set of all the linear assemblages of the principle components ${\mathbf{U}}_r\mathbf{t}$ defines a reduced dimensional domain (r-dimension) over the tangent space ${\mathcal{T}}_{\mu }\mathcal{G}(n,2)$ with the image of an Riemannian exponential map. Correspondingly, the change over $\mathbf{t}$ will give perturbations on the airfoil, which is, obviously, independent of the affine deformations.

The current dataset that contains 7007 airfoil shapes can be used to give an example. The Karcher mean is plotted in Figure 3b, and it is clear that the value is centered among the elements. Figure 4 is a plot of the first 10 modes of the $\mathbf{H}$ over the current dataset, and the first 6 modes capture more than 95% of the total variability, which is a significant reduction from the original 10 dimensions of the design variables. Moreover, over 85 % of the total variability is controlled with the first 4 modes, which must be very beneficial for a quick understanding of the large design space.

3. Airfoil Shape Optimization

In the previous section, the shape parameterization method was described, and in this section, two airfoil optimizations will be performed. Specifically, the optimization frame is comprised of an optimal search module implemented by the stimulated annealing toolbox in MATLAB${}^{\circledR }$, a computational fluid dynamics (CFD) module implemented by the open source code CFL3D [4], and a mesh deformation module implemented by Pointwise${}^{\circledR }$. During the search phase, all parameters of the optimization module are frozen. In addition, the turbulence model employed here is Menter’s $k-\omega$ shear stress transport model [43], and the only variable between cases is the freestream condition. The CST approach utilized in the present work employs 12 parameters to provide a fair description of the airfoil shape, and the upper and lower limits of each design parameter are limited to a 20% increase or decrease in the baseline value. As a data-driven parameterized method, it is impossible to arbitrarily specify the range of design variables. Therefore, the design variables’ limits are derived from the maximum and minimum values of the dataset. As shown in Figure 4, six design variables are utilized to complete the optimization using the Grassmannian shape parameterization (GSP) method.

3.1. Subsonic Airfoil Optimization

Maximizing the lift-to-drag ratio is a common design target for aircraft under subsonic flow conditions, and hence an takeoff freestream velocity, $Ma=0.15$, is used in this case. The Reynolds number based on the chord length is 6 million. The bound of farfield is located at 100c and 72,000 grid points is used to discretize the entire computational domain. The distance of the first grid point away from the airfoil is $1\times {10}^{-6c}$ to guarantee a normalized wall distance $y^{+}<1$.

The aerodynamic performance of the airfoils given by optimization was shown in Figure 5a,b. Obviously, the whole lift line is lifted up a lot by the airfoil given by the GSP method, but a little bit by the airfoil given by the CST method. As expressed in the previous section, the GSP method separates the affine deformation from the shape deformation if the matrix $\tilde{M}$ has no change, which means the scale and translation invariance during the shape deformation. From the view of low-speed airfoil design, a consensus is that the curvature of the upper surface affects the lift-to-drag ratio pretty much if the thickness cannot be changed due to the structural constant. In Figure 6a, the optimized airfoil given by GSP has a larger curvature on the upper surface line and camber on the tail of the lower surface line, and both of these two characters are beneficial to aerodynamic performance. The CST method provides an airfoil with a 1% camber increase for obtaining improved lift force. The slopes of the lift line plotted in Figure 5a from different airfoils are close to each other, which means a similar camber of airfoils. The experimental data come from Ladson [44]. The lift-to-drag ratio, K, demonstrated in Figure 5b, the maximum K is improved from 83 to 110 with the GSP method, and 87 with CST method, respectively. The distribution of pressure coefficient $C_p$ of the airfoil at $\alpha =10deg$ is demonstrated in Figure 6b, the suction peak is increased and the total area of $C_p$ is enlarged from the airfoil given by the GSP method.

The streamline at a large angle of attack, $\alpha =16deg$, is shown in Figure 7. The velocity contour is very high on the nose of the optimized airfoil from the GSP method. The positive velocity region is enlarged from the baseline to the optimized airfoils.

3.2. Transonic Airfoil Optimization

The most typical transonic airfoil which is used to validate the CFD code is RAE2822. An experiment is performed at $Ma=0.729$ with $\alpha =2.31deg$ [45]. The Reynolds number based on the chord length is 6.5 million.

The minimized drag coefficient $C_d$ with a constant lift coefficient $C_l$ is very common in the transonic airfoil design process. The supercritical airfoil proposed by Whitcomb et al. [46] was proved to be very good at such flow conditions. The shock wave often appears at the latter part due to the little curvature of the front part of the upper surface. The main principle of optimizing transonic airfoil is smoothing the adverse pressure gradient of the upper surface and increasing the camber of tail to recover the lift.

In Figure 8a, a small deformation in the airfoil shape can induce the pressure distribution plotted in Figure 8b, which has very large change, especially at the location with shock wave. The pressure distribution comes from simulations with lift coefficient $C_l=0.7429$, and the drag coefficients $C_d$ for baseline, CST and GSP models are 0.0123, 0.0101 and 0.099, respectively. Judging from the results, the new method brings very little improvement, however, it must be mentioned here that the shock wave is very sensitive to small-shape deformation. That is, a very little perturbation in the design variables could be enough and the drawback demonstrated in Figure 2 could be covered up. Hence, a general design must be affected by the bound of the design variables.

In the new case, the same flow condition is maintained. SD7003, a low-Reynolds number airfoil, which serves as the baseline. Figure 9a depicts the optimized airfoils, and the CST method produces little change on the baseline airfoil. The GSP method, on the other hand, produces an airfoil shape with the characteristics of supercritical airfoils. Meanwhile, the suction peak of a low Reynolds number airfoil is significantly higher than that of a supercritical airfoil, as shown in Figure 9b. Furthermore, the $Cp$ distribution produced by the GSP method has a very similar pattern to that shown in Figure 8b, demonstrating that it has a very stable performance even when the baseline airfoil is changed.

The gradient of density is used to demonstrate the strength of the shock wave on the airfoil. The shock wave in Figure 10 shows that the high suction peak induces strong adverse-pressure gradient, which corresponds to a strong shock wave. As the optimized module flattens upper surface, the peak drops and pressure gradient is released.

4. Conclusions

In this study, a method for shape parameterization was utilized to provide a precise and generalized description of airfoil. As a data-driven approach, a wide variety of airfoil types can be included prior to optimization, and this advantage can give the designer a sense of what the types of the optimized airfoil will be. However, the disjointed design space, which is caused by the non-normalized perturbation of design variables, could hinder the optimized method from achieving the best value. The benefit of transforming an airfoil from Euclidean space to Grassmannian space is that, as a result of the landmark-affine deformation, the design space can be closely matched. The dimensionality of the design space can also be reduced to a single digit using the PGA method. The Karcher-centered domain of normal coordinates provides an advantage in that the new method has the shortest routine to any element in the space, which, presumably, requires the least amount of time to locate the optimal shape in the design space. In the meantime, the affine deformation, which cannot be prevented with traditional methods such as CST, can be separated with the new method, and the design must be enhanced by the regularization against non-physical deformations.

We demonstrated that the new method is very robust at both the subsonic and transonic levels through the use of numerical processes of a variety of cases that were optimized. The disadvantage of the traditional parameterization method, which is that it heavily relies on the bounds of design variables, can now be overcome.