Adaptive Free-Form Deformation Parameterization Based on Spring Analogy Method for Aerodynamic Shape Optimization

Abstract

An adaptive Free-Form Deformation parameterization method based on a spring analogy is presented for aerodynamic shape optimization problems. The proposed method effectively incorporates the gradients of the objective and constraint functions, achieving automatic control point adjustment based on variances in design variable components. To evaluate the performance of the adaptive FFD parameterization method, two 2D airfoil optimization design problems are examined. The optimization of the RAE2822 airfoil with 12, 18 and 24 design variables demonstrates superior results for the adaptive method compared to uniform parameterization. The adaptive method requires fewer iterations and achieves lower objective function values. Additionally, the optimization design from NACA0012 to RAE2822 airfoil with 18 design variables shows that the adaptive parameterization method achieves a lower drag coefficient while satisfying the optimization objective. This validates the method’s capability to finely adjust airfoil shapes and capture more optimal design points by exerting stronger control over local shapes. The proposed adaptive FFD parameterization method proves highly effective for optimizing aerodynamic shapes, offering stability and efficiency in the early stages of optimization, even with a limited number of design variables.

1. Introduction

Aerodynamic shape optimization design is fundamentally a mathematical optimization problem, comprising essential elements such as the objective function, design variables, constraint functions, and variable limits [1]. For aircraft-related aerodynamic shape optimization, objective function values or constraints are often represented by aerodynamic forces, such as drag and the drag-to-lift ratio. In CFD-based optimization, evaluating these objective function values or constraints is time-consuming, and multiple optimization iterations impose high demands on computational resources.

Optimization methods can be classified into non-gradient-based and gradient-based approaches, depending on whether gradient information is used [2,3,4]. Non-gradient-based optimization methods offer global optimization capabilities, theoretically identifying optimal solutions within the design space [5,6]. The choice of design variables, however, directly influences the problem’s global optimization, creating different optimization challenges. Non-gradient-based methods can locate an optimal solution for a given problem but cannot guarantee the best aerodynamic shape for the design physics.

The efficiency of non-gradient-based optimization methods is directly related to the number of design variables. In three-dimensional aerodynamic optimization, hundreds of design variables are often required, resulting in substantial computational loads. Surrogate models can partially reduce the need for objective function evaluations, though the accuracy of these models remains a concern [7,8].

In comparison, gradient-based methods offer significant efficiency benefits. These methods require additional gradient calculation steps, and currently, adjoint methods have become the mainstream approach for gradient calculation [9,10]. The computational cost of gradient calculation based on adjoint methods is almost independent of the number of design variables. However, projecting grid sensitivities onto design variables, along with calculating geometric quantity gradients via finite differences methods, can lead to significant computational costs with many design variables [11].

Moreover, excessive design variables in certain optimization problems can complicate finding optimal aerodynamic shapes, potentially requiring additional optimization iterations [12].

Both non-gradient-based and gradient-based methods rely on the same design variable definition approach. Numerical parameterization methods are used to convert many grid points into a manageable parameter set. Common parameterization methods include Hicks–Henne, B-spline, Class Shape Transformation (CST), and Free-Form Deformation (FFD), with FFD widely applied due to its adaptability in 2D and 3D aerodynamic design optimization [13,14].

In FFD, control point displacements are the design variables. By mapping control points to grid points, changes in aerodynamic shape are linked to FFD control point displacements. However, there is a trade-off in the aerodynamic shape optimization process. Using a larger number of control points (design variables) can lead to better optimization results but increases computational complexity. Conversely, using a smaller number of control points can improve the convergence speed of optimization but may not achieve the desired optimization performance. Thus, there is still room for research and consideration in defining the number and positions of control points in the FFD parameterization method.

Design variables play a crucial role in formulating optimization problems and are closely related to the efficiency and quality of optimization results. In order to improve the effectiveness of aerodynamic shape optimization, some adaptive parameterization methods have been adopted, such as progressive parameterization method [15], nested and self-adaptive Bezier parameterization approach [16], adaptive multilevel algorithm [17], B-spline knot insertion method [18] and so on. In adaptive FFD parameterization research, the adjustment or insertion of control points depends on the gradient information of the corresponding design variables with respect to the objective function, as mentioned earlier. However, in optimization algorithms based on quadratic programming, the solution for design variables is based on the Karush–Kuhn–Tucker (KKT) conditions, which not only include the gradient of the objective function but also the gradient of constraint functions [19]. In other words, changes in the aerodynamic shape are related to both the gradients of the objective function and constraint functions. Therefore, both gradients should be considered when adjusting control point density.

The updated design variables after each optimization sub-iteration reflect the combined effects of the gradients. Hence, using the updated design variable solution as the basis for the adaptive positioning of FFD control points is more reasonable. The objective of this study is to apply the FFD parameterization method with a fixed number of control points. Based on the updated design variable information, the spring analogy method is used to automatically adjust the chord-wise positions of control points. This adjustment enables control points with a significant impact on the objective function and constraint functions to gather closely, while control points with a lesser impact are more widely spaced. The newly positioned control points are then re-parameterized, and the performance of this parameterization strategy in aerodynamic shape optimization design for airfoils is analyzed.

The organization of this paper is as follows: Section 2 introduces the SU2 optimization framework used in this study, the principles of the adaptive FFD parameterization method based on spring analogy, and the validation of adaptive control point adjustment for spline curve fitting capability. Section 3 discusses the results of the 2D airfoil optimization design based on adaptive FFD parameterization, including optimization of the Rae2822 airfoil and the design transformation from NACA0012 to Rae2822.

2. Methodology

The open–source suite SU2 has been used as the optimization framework in this study and has been applied in numerous aerodynamic optimization design researches [20]. This section provides a concise explanation of the different functional modules incorporated within the SU2 optimization framework, as well as the fundamental principle underlying the adaptive FFD parameterization method, which is based on the spring analogy approach.

2.1. Optimization Frame

In the SU2 optimization framework, the CFD module serves as an unstructured flow solver. The flow is governed by the compressible RANS equations and employs the Spalart–Allmaras (SA) turbulence model [21]. To handle the flow’s inviscid terms, the Jameson–Schmidt–Turkel scheme (second order in space) is utilized, striking a favorable balance between accuracy and efficiency for the optimization problem addressed in this paper. For the turbulence model, a scalar upwind solver is employed. The Green–Gauss method [22] is applied to estimate the spatial gradients of the flow at the cell faces, which are crucial for determining the viscous fluxes. The time integration of the flow equations and turbulence model employs the Euler implicit scheme. To solve the linear system, the Flexible Generalized Minimum Residual (FGMRES) method [23] is employed with ILU preconditioning. The function sensitivity of the flow variables is obtained through discrete adjoint equations using Algorithmic Differentiation.

Mesh deformation operation is applied to the design variables, which are the mesh coordinates parameterized by FFD method. The linear elasticity equation is used to deform the volume mesh surrounded FFD control points. Cell stiffness is scaled by distance to nearest solid surface.

FFD based on B-spline basis function is adopted as the parameterization method of optimization problem in this work and the relationship between global and parametric mesh coordinates is described by the following equation (Equation (1)):

(1)$X(u,v,w)=\sum _{i=0}^{l-1}\sum _{j=0}^{m-1}\sum _{k=0}^{n-1}{N}_{i,p}\left(u\right)N_{j,q}\left(v\right)N_{k,r}\left(w\right)B_{i,j,k}$

where N is the open uniform B-spline basis function and B denotes the control point of FFD (namely, design variable). When a parametric coordinate (u,v,w) is given, a global mesh coordinate (X(u,v,w)) is determined and basis functions define the projection relationship of the two coordinates. The parameter coordinates are iteratively solved using Newton’s method.

The optimization method employed in this study is the Sequential Least Squares Programming (SLSQP) algorithm [24] implemented in Scipy library. The optimizer uses a quasi-Newton Hessian approximation and an L1-test function in the line search algorithm.

2.2. The Adaptive FFD Parameterization Strategy Based on the Spring Analogy Method

2.2.1. Spring Analogy Method

Achieving adaptive positioning of parameterized control points, Gnoffo introduced the spring analogy method [25]. In this method, the parameterized control points are modeled as nodes connected by springs of varying stiffness. Specifically, considering a 1-D spring system where the nodes move solely in the x-direction, the equilibrium state of the system results in a balance equation involving the elastic forces between the $i-1$, i and $i+1$ spring nodes, given by the equation

(2)$K_{i+1/2}(x_{i+1}-x_i)=K_{i-1/2}(x_i-x_{i-1})$

where $x_i$ is the 1-D coordinate of the spring node, $K_{i+1/2}$ is the stiffness of spring between the ith and $(i+1)$th spring knot and is modeled as

(3)$K_{i+1/2}=1+(A-1)F\left({\displaystyle \frac{f_i+f_{i-1}}{2}}\right)$

In Equation (3), parameter A represents the ratio between the maximum and minimum node spacing in the spring system, reflecting the density of the adapted parameterized control points. $f_i=(a_i-a_{min})/(a_{max}-a_{min})$, where $a_i$ and $a_{i-1}$ determine the relative magnitude of stiffness $K_{i-1/2}$. For general aerodynamic shape optimization problems, they represent the sensitivities of the objective function (or constraint) with respect to the design variables. $a_{min}$ and $a_{max}$ are the maximum and minimum values of sensitivity among all design variables, respectively. After normalization, $1/2(f_i+f_{i-1})$ is within the range of $[0,1]$, ensuring that $F\left(f\right)$ also falls within $[0,1]$. Nakahashi [26] suggested $F\left(f\right)=fB$, where B is a positive constant. In this study, $F\left(f\right)=min\{1.0,2sin(f\left)\right\}$ so that the ratio of the maximum and minimum node spacing in the spring system depends solely on A. Additionally, this setting allows for the maximum node spacing ratio and a wider adaptive range for parameterized control points. Similarly, for the 1-D problem, the unknown positions of the nodes, to be solved using the spring analogy method, form a linear equation system. To solve this system, $K_{i+1/2}$ is first calculated and stored, followed by solving a tri-diagonal linear equation system concerning $x_i$, as shown below

(4)$\left[\begin{array}{ccccc}{K}_{1/2}+K_{3/2}& -K_{3/2}& & & \\ -K_{3/2}& K_{3/2}+K_{5/2}& -K_{5/2}& & \\ & \ddots & \ddots & \ddots & \\ & & -K_{N-5/2}& K_{N-5/2}+K_{N-3/2}& -K_{N-3/2}\\ & & & -K_{N-3/2}& K_{N-1/2}+K_{N-3/2}\end{array}\right]\left(\begin{array}{c}{x}_2\\ x_3\\ ⋮\\ x_{N-1}\end{array}\right)=\left(\begin{array}{c}{K}_{1/2}{x}_1\\ ⋮\\ K_{N-1/2}{x}_N\end{array}\right)$

In Equation (4), $x_1$ and $x_N$ are the two endpoints of the spring system. In general, these two points are kept fixed.

2.2.2. Verification of Control Points Adaptive Spline Curve Fitting Capability

The FFD parameterization method uses spline functions as basis functions, so the influence of control point positions on the basis functions can be examined through their impact on B-spline curves.

• B-spline Curve fitting

A B-spline curve is given by $P\left(t\right)={\sum }_{i=1}^{n+1}{B}_i{N}_{i,k}\left(t\right)$, where $N_{i,k}\left(t\right)$ denotes the basis function, $B_i$ represents the control points, and $D_i$ is the fitting target point. Their relationship is shown schematically in Figure 1. Utilizing the gradient descent algorithm, we compute a set of t values that minimize the root mean square of the difference between the y-coordinates of the given points $D_i$ and the points on the B-spline curve (with $B_i$ as the control points) with the same y-coordinates. This implies that, even when the fitting points remain unaltered, different control points yield different parameter values t for the B-spline curve. Consequently, in order to generate points on the B-spline curve with x-coordinates matching the fitting points (within a tolerance of 0.01), adjustments in the control points become necessary. The number of control points $B_i$ significantly influences the ability of the B-spline to accurately describe the shape. In cases where the number of control points $B_i$ is fixed, modifying their positions allows for control over the shape of the spline curve. Moreover, the density or sparsity of the control points plays a vital role in the local shape representation of the B-spline. The adaptive control point B-spline curve fitting approach presented in this paper is based on these principles.

• b.. Verification test

In this study, the method is tested by fitting a segmented function defined as follows:

(5)$\left\{\begin{array}{cc}-{(x-0.5)}^2+0.25& 0\le x<0.4\cup 0.6\le x\le 1\\ -8\ast {(x-0.5)}^2+0.32& 0.4\le x<0.6\end{array}\right.$

The above function consists of three segments of parabolic curves, where the curvature varies as a function of x. The locations of high curvature correspond to the breakpoints of the segmented function and the extremal points of the function. These positions exhibit rapid changes in the function gradient, making it challenging to fit this curve using a B-spline curve with a small number of control points. The curve and its curvature are illustrated in Figure 2.

Initially, a uniform open B-spline curve of degree 3 with 10 control points is used to fit this curve. A total of 50 points are uniformly sampled on the curve, and the proposed method is employed to generate 50 points on the desired B-spline curve. These points are adjusted to match the x-coordinates of the given points (within a certain tolerance). Subsequently, the y-coordinates of the 50 points on the B-spline curve are calculated to approximate the given points. This process leads to an optimization problem, where the design variables are the y-coordinates of the B-spline curve, and the objective function f represents the root mean square difference between the y-coordinates of the 50 points on the B-spline curve and the given points. Two optimization algorithms, namely the gradient-based method and the non-gradient-based method, are used to solve this optimization problem. It is verified that both methods converge to the same minimum value of the objective function, indicating that the objective function for this type of problem is unimodal.

The fitting result of the B-spline curve with 10 control points, uniform open, and degree 3 is shown in Figure 3. The root mean square error between the B-spline fitting points and the data points is 0.0074. It can be observed that at positions where the curvature of the data points’ curve changes slightly, the difference between the B-spline fitting and the data points is minimal. Larger fitting errors occur at positions with significant curvature variations. Next, based on the spring analogy method, the control points of the B-spline curve are adaptively adjusted to fit the data points. The adaptivity is based on the curvature of the data points’ curve. Since the data points do not coincide with the control points of the B-spline curve, linear interpolation is used to map the curvature values from the data points to the control points. The different curvature values of the control points are reflected in the stiffness of the spring system, where higher curvature values correspond to greater stiffness, resulting in shorter distances between the two nodes of the spring after system equilibrium. Figure 4 presents the fitting result using adaptive control point B-spline curve fitting with 10 control points, uniform open, and degree 3. The root mean square error between the B-spline fitting points and the data points is 0.0027. The fitting error is significantly reduced compared to the non-adaptive approach, indicating an improved fitting accuracy.

2.2.3. Optimization Process for Adaptive Control Point Parameterization Based on FFD

Adaptive control point parameterization requires an indicator. In this study, the indicator is chosen as the updated design variables corresponding to each control point after the main iteration of the quadratic programming optimizer. Therefore, the updated position of control points comprehensively reflects the impact of the objective function and the gradients of active constraints. The control point adaptation operation occurs in the optimization framework after the main iteration step. The updated design variables are introduced into the spring system to calculate the new positions of the control points. Since FFD control point positions have changed, the aerodynamic surface grid points need to be re-parameterized, and the gradients of the objective function and constraints are calculated. It should be noted that based on ${\displaystyle \frac{\partial J}{\partial X}}={\displaystyle \frac{\partial J}{\partial S}}{\displaystyle \frac{\partial S}{\partial X}}$, where ${\displaystyle \frac{\partial J}{\partial S}}$ is the adjoint derivative, it is independent of the parameterization method and does not require resolving. The FFD control point adaptive process does not incur significant computational costs. The detailed optimization process is as follows:

Step 1: Evaluate the objective function, constraint function, and their gradients with respect to the design variables.

Step 2: Determine convergence based on the KKT conditions.

Step 3: If convergence is achieved, stop the optimization; otherwise, update the design variables using the SLSQP algorithm.

Step 4: When the adaptive condition for control points is satisfied, use the spring analogy method with the current design variables to adaptively adjust the control point positions as shown in Figure 5, and then re-parameterize and re-project the gradients.

Step 5: Solve for new design variables based on the updated gradient information.

Step 6: Update the grid based on the new design variables and repeat steps 1–6.

A schematic diagram of the optimization framework is shown in Figure 6.

3. Optimization Results

To evaluate the performance of adaptive FFD parameterization in aerodynamic shape optimization, two optimization design cases were implemented: the RAE2822 airfoil under cruise conditions and the NACA0012 airfoil, with the RAE2822 airfoil under cruise conditions as the design point.

3.1. Optimization Design of RAE2822 Airfoil Based on Adaptive FFD Parameterization Method

In this study, the optimization objective for all cases is the RAE2822 airfoil in viscous, transonic flow, framed as an optimization problem. The flow conditions are defined by a freestream Mach number of 0.734 and Reynolds number of 6.5 $\times {10}^6$. To achieve the optimal design, several constraints are imposed, including a maximum lift coefficient of 0.824, a minimum pitching moment coefficient (evaluated at the quarter-chord) of −0.092, and a requirement to maintain the airfoil area during optimization. The optimization problem is mathematically summarized in Table 1.

The computational grid adopts a hybrid form, with a quadrilateral grid near the airfoil and gradually transitioning to a fully triangular grid in the far-field direction. The grid resolution is adjusted by controlling the number of grid points on the airfoil and the stretching ratio of the viscous grid. In this study, five sets of grids are generated to verify grid convergence. Table 2 shows the number of elements and the coefficient results for drag and moment of the four grids. Considering both computational effort and accuracy, all cases in this study use Level 3 as the initial grid for the optimization design as shown in Figure 7.

In the distribution of FFD control points, the upper and lower control points have the same vertical position, forming a rectangular frame. This point distribution method, compared to placing points along the upper and lower surfaces of the airfoil, provides similar control capabilities for the grid points, and the optimization results are essentially the same. The design variables in this study are the displacement of FFD control points in the y-direction. The control points at the leading and trailing edges of the airfoil are thickness variables, moving in opposite directions to ensure that the positions of the leading and trailing edge points of the airfoil remain essentially unchanged. In addition to the design variables related to FFD control points, there is also an angle of attack variable. The angle of attack can vary in each optimization iteration step to achieve the target lift coefficient.

In Figure 8, FFD control points enclosing the RAE2822 airfoil shows the distribution of FFD control points enclosing the RAE2822 airfoil. The design variables generated involve displacements in the z-direction for all control points. To prevent large curvatures at the leading and trailing edges of the airfoil, the displacements of the upper and lower control points at the leading and trailing edges are set to be opposite in sign.

3.1.1. Uniform FFD Control Point

Before evaluating the optimization effectiveness of the adaptive FFD parameterization method, the influence of different numbers of FFD control points on the optimization results is analyzed. In total, three different control point quantities, namely 14, 20, and 26, are selected. Based on the previously defined design variable formulation, the number of design variables for each case is 12, 18, and 24, respectively. The convergence history of the drag coefficient for these three optimization cases is shown in Figure 9. It is observed that as the number of design variables increases, the converged drag coefficient becomes smaller, but a larger number of iterations are required for convergence. The peaks in the curve are due to the SLSQP line search encountering infeasible points. Regarding the pressure distribution, all three design variable cases eliminate the shockwave at the 0.6c location on the upper surface of the airfoil, exhibiting similar upper surface pressure distribution and airfoil shape as shown in Figure 10. The major difference lies in the airfoil’s leading edge, where the optimized shapes do not exhibit a strong suction peak in the pressure distribution compared to the initial shape. Moreover, increasing the number of design variables reduces the suction at that location. Significant differences are observed in the pressure distribution on the lower surface of the airfoil among the different design variable optimization cases. All three results exhibit increased leading edge loading and reduced trailing edge loading compared to the initial shape, indicating that the optimization process tends to raise the airfoil nose to reduce drag. Consequently, the optimal values appear at the boundary of the feasible domain in the optimization problem. Analyzing the FFD control point displacements corresponding to the optimal values, it is observed that the displacements of the control points on the upper surface of the airfoil are relatively small, as shown in Figure 11. However, these small displacements have a significant impact on reducing drag by weakening the shockwave on the airfoil. In this case, complex changes to the airfoil shape are not necessary, and the three cases result in almost identical upper surface curves. There are no significant differences observed among the different design variables. As for the control points on the lower surface, the FFD control points corresponding to the 24 design variables are more densely distributed, providing stronger control capability over the leading and trailing edge curves and resulting in more complex curvature changes. The variations in the displacements of the FFD control points along the chordwise direction in these cases indicate that the optimization trend for the aerodynamic shape of the airfoil may not result in significant curvature changes across the entire chordwise span. Intense changes in the shape curvature may only occur locally (typically at the leading and trailing edges) for two-dimensional problems. Therefore, to achieve effective control over local shapes and enhance drag reduction effects, positions with large curvature changes should have more control points, while fewer control points can meet the shape control requirements at other positions.

3.1.2. Adaptive FFD Control Point

The purpose of adaptive FFD parameterization is to automatically adjust the distribution of FFD control points by aggregating and dispersing them based on a relatively small number of design variables. The objective is to have more control points at locations with significant curvature changes and fewer control points at locations with minor curvature changes, while keeping the control points unchanged. In the implementation of the spring analogy method described in this paper, the local density of control points is adjusted by tuning parameter A. Based on the analysis presented earlier, certain local positions require a sufficient number of control points to achieve more precise shape control. In the case with 14 design variables, the number of control points is too small. Computational experiments have shown that even after local aggregation, the control points still do not reach the desired density. Increasing the parameter A to further increase the local control point density is not feasible because there is an upper limit to the ratio between the maximum and minimum distances between control points. When this limit is reached, increasing parameter A cannot further increase the local control point density. Taking these factors into consideration, this paper selects the case with 20 control points (18 design variables) as the evaluation case for adaptive FFD parameterization. The goal is to adaptively aggregate a certain number of control points locally by adjusting parameter A, so that the minimum distance between control points in this case is similar to the distance between the 26 uniformly distributed control points. The aim is to examine the improvement in drag reduction achieved by this method.

In addition to the issue of control point density, there is also the question of when to initiate the adaptive process. If the adaptive process starts too late, even close to convergence, although further improvement in the objective function can be achieved through control point adjustments, the optimization efficiency is relatively low, and more iterations are still required. On the other hand, if the adaptive process starts too early, when the aerodynamic shape is still in the preliminary adjustment phase, adapting the control points is unnecessary. In the context of this paper’s case study, during the initial stages of optimization, the most sensitive region for drag reduction is the shock wave area on the upper surface of the airfoil. The optimization process at this stage aims to eliminate the shock wave, and the upper surface undergoes less pronounced curvature changes. A coarser distribution of control points can effectively control the shape, and a few iterations are sufficient to achieve shockfree optimization.

To determine the timing for parameterized adaptation, a common criterion is to use the slope of the objective function convergence history. In the early stages of optimization, the objective function decreases rapidly, resulting in a large slope, while in the later stages, the objective function changes less, with the slope approaching zero. The slope criterion for initiating adaptation should lie between these two extremes. In this paper, the ratio between the current iteration step’s slope and the previous iteration step’s slope is selected as the criterion. The slope is calculated using forward difference. When this ratio, $S_i/S_{i-1}<0.1$, the adaptive control point positioning begins. The threshold of 0.1 indicates that the current step’s reduction in the objective function is 0.1 times that of the previous step. This ratio reflects the changing trend of the objective function for three adjacent iterations. In the early stages of optimization, this ratio is close to 1, indicating a continuous decrease in the objective function. When this ratio suddenly decreases, it indicates that the design variables have started to undergo minor updates, and the optimization process has entered the fine-tuning phase of the shape, which is a suitable time to initiate adaptation.

The following discussion focuses on the performance of aerodynamic shape optimization using the slope ratio as the criterion for FFD parameterized control point adaptation in the case with 18 design variables. The convergence history of the objective function for three optimization processes is shown in Figure 12: the case with 18 design variables and uniform control points, the case with 18 design variables and adaptive control points, and the case with 24 design variables and uniform control points. In the case of adaptive control points, the automatic adjustment of chordwise position of the control points was performed only once, occurring at the 33th iteration step, where the slope ratio, $S_i/S_{i-1}=0.1$, satisfying the condition for starting the adaptation process. The specific operation of the adaptation process involves inputting the updated design variable solution from the current step into the spring system, solving the system to generate new chordwise positions for the control points, performing FFD parameterization based on the new control points and the aerodynamic shape, and then starting the optimization process to find the optimal solution for the new problem.

In this paper, the parameter A is set to 20, and the adaptive control points tend to cluster around two regions: the leading edge and approximately 80% of the chord length. This clustering results in local minimum chordwise normalized distances of 0.061 and 0.055, which are smaller than the uniform control point spacing of 0.083 in the case with 24 design variables, as shown in Figure 13. From the perspective of optimization convergence history, the objective function of the adaptive case continues to decrease after the control point adjustment, and the objective function of subsequent iteration steps is generally lower than that of the case with 24 design variables. The final convergence result of the objective function is better for the adaptive case than for the case with 24 design variables. This indicates that the proposed adaptive method for control point chordwise positions not only maintains the characteristics of rapid objective function decrease with fewer design variables but also achieves lower objective function convergence results with more design variables.

The pressure distribution and airfoil shape in Figure 14 reveal that the forward loading is further enhanced in the adaptive case. The positive pressure on the lower surface near the leading edge is higher compared to the cases with 18 design variables and 24 design variables. However, in general, the chordwise scale of the positive pressure region is smaller in the adaptive case compared to the other two cases, indicating that the generated moment in this region may not change significantly. At approximately 80% of the chord length, where the control point spacing is at a minimum after adaptation, there is some fluctuation in the pressure distribution on the lower surface, which can also be observed on the airfoil shape. The possible reasons for this phenomenon are that the control point spacing is small, leading to larger differences in spacing values. Additionally, the spline basis function used in the FFD parameterization in this paper has an order of 4, which may exacerbate this issue. Therefore, reducing the value of A and increasing the order of the basis function may help alleviate this problem.

3.2. From NACA0012 to RAE2822

The purpose of this section is to demonstrate the aerodynamic optimization design process from an arbitrary airfoil to an optimized airfoil that meets specific design requirements using the adaptive FFD parameterization framework. In this particular case, starting with the NACA0012 airfoil, the objective and constraint conditions for the RAE2822 airfoil optimization problem mentioned earlier are given. Through the optimization design process based on adaptive FFD parameterization, the optimized RAE2822 airfoil is obtained. The formulation of this optimization design problem is similar to that in Section 3.1, with 20 FFD control points corresponding to 18 design variables.

This optimization design process raises the following questions: Given the significant differences between the NACA0012 airfoil and the RAE2822 airfoil, can the optimization algorithm find the global minimum of the objective function while satisfying the design conditions? Can the adaptive FFD parameterization further reduce the extremum of the objective function, and is this extremum the global minimum in the design space?

As mentioned earlier, the design variables in the optimization process in this paper include both the displacement of the FFD control points and the angle of attack. In each iteration step, the angle of attack changes to satisfy the lift coefficient constraint, while other constraints and the objective function form the Lagrangian function for solving the optimization problem. The optimization problem-solving approach used in this study may lead to the following situations: First, if the target lift coefficient is too large and exceeds the maximum lift coefficient of the initial airfoil, it may not be achievable by adjusting the angle of attack of the initial airfoil. Second, even if the target lift coefficient is satisfied, the required angle of attack may be very large. This was the case in the optimization from NACA0012 to RAE2822 in this section. For the flow solver and numerical methods used in this research, the angle of attack for the NACA0012 airfoil at a lift coefficient of 0.824 is 21.1°, while the optimization results in the previous section indicate that the angle of attack for the optimized RAE2822 airfoil is between 2 and 4°. The significant difference between these values necessitates a large number of iterations and may even result in a suboptimal solution.

To address this issue, the optimization process in this study is divided into two stages: In the first stage, while keeping the constraint conditions unchanged, a relatively small target lift coefficient is set to form an optimization problem, and the optimal solution for this problem is obtained. In the second stage, using the optimal shape from the first stage as the starting point, the final optimization process of the airfoil is performed with the target lift coefficient as a constraint.

For the optimization from NACA0012 to RAE2822, a target lift coefficient of 0.6 was used for the first stage. The corresponding angle of attack for NACA0012 to achieve this lift coefficient is 14°. Figure 15 (left) shows the iteration history of the drag and moment coefficients for the first optimization stage, representing the objective function and constraints, respectively. At 82nd iteration, the KKT conditions are satisfied, with a drag coefficient of 0.01095 and an angle of attack of 1.93°. The moment coefficient gradually approaches the given constraint value of −0.092. At the end of the iteration, the moment coefficient is smaller than the constraint value, indicating that the optimization result is not on the boundary of the feasible domain, and the moment constraint did not take effect. The first stage is an intermediate process in the overall optimization, and finding a local minimum that is not on the boundary of the feasible domain is acceptable as it does not affect the final optimization result in the second stage. Figure 15 (right) shows the iteration history of the drag and moment coefficients of the second optimization stage. In this stage, the target lift coefficient is 0.824, and the initial iteration starts with an angle of attack of 3.22°. At the 157th iteration, the KKT conditions are satisfied, with an angle of attack of 3.01°, a drag coefficient of 0.01329, and a converged moment coefficient of 0.0088, which is smaller than the given value. This indicates that this extremum point is not on the boundary of the feasible domain. Figure 16 also shows the evolution of the airfoil in the two optimization processes.

There are some abrupt iteration points in the convergence history of the drag and moment coefficients. This is due to the large step size in the online search process, which leads to significant changes in the airfoil shape and results in large fluctuations in the coefficient values. When this occurs, the optimizer automatically reduces the step size until the objective function decreases compared to the previous iteration step.

Regarding the second question raised at the beginning of this section—whether the adaptive FFD parameterization method can find design points with lower drag coefficients compared to uniform control points while satisfying the optimization constraints—this section explores this using the adaptive criterion introduced in Section 3.1. The adaptive parameterization is only implemented in the second stage of the optimization design. Based on the criterion that the slope ratio of the objective function with respect to the iteration step number is $S_i/S_{i-1}=0.1$, the adaptive parameterization is performed at the 44th iteration. Figure 17 shows the convergence history of the drag coefficient for the adaptive parameterization and uniform parameterization optimization design cases. It can be observed that the drag coefficient continues to decrease after the adaptive parameterization is performed, and it satisfies the convergence condition at 438th iteration step with a converged drag coefficient of 0.01320, which is further reduced compared to the optimized drag coefficient of the uniform parameterization case. The optimized airfoils and pressure distributions for the two optimization cases are shown in Figure 18. From the optimized airfoil perspective, the upper surface shapes of the two cases are almost identical, while there are significant differences in the lower surface shape, especially in the leading edge. The adaptive parameterization case shows greater variation in the curvature of the airfoil shape, and the concavity of the leading edge is more pronounced. The pressure distribution also reflects the differences in the optimization results. The pressure distributions on the upper surfaces of the two cases are generally similar. For the adaptive parameterization case, the pressure is higher at the location corresponding to the concave region of the lower surface, indicating a certain degree of leading edge suction. These optimization results demonstrate that adaptive parameterization, with its enhanced control over local shapes, enables finer adjustments to the airfoil, achieving improved design points compared to uniform parameterization.

4. Conclusions

This paper proposes an adaptive Free-Form Deformation parameterization method, based on a spring analogy, to address the balance between design variable count, gradient calculation cost, and optimization algorithm stability in gradient-based methods for aerodynamic shape optimization. In this approach, FFD control points are treated as spring nodes, with stiffness values that adjust according to the updated design variables during optimization. The resulting spring system forms a tridiagonal linear equation set, and its solution provides updated FFD control point positions. By using design variable solutions as the basis for spring stiffness calculations, this method comprehensively reflects the influence of objective and constraint function gradients on control point adjustment, allowing automatic density adjustments of control points based on variations in design variable values.

To evaluate the performance of this adaptive FFD parameterization method, two 2D airfoil optimization cases were studied and analyzed. First, the RAE2822 airfoil was optimized with 12, 18, and 24 design variables. Starting conditions for adaptation were set with a slope ratio of 0.09 for the objective function and a control point distance ratio of 20. Convergence history showed a continuous decrease in the objective function post-adaptation, achieving lower values than the 24-variable uniform case. The comparison of the convergence results of the objective function shows that the adaptive parameterization method is superior to the uniform parameterization method and requires fewer iterations. This demonstrates the proposed adaptive method’s effectiveness in maintaining rapid objective function reduction with fewer design variables and achieving lower convergence results with more variables, surpassing the uniform parameterization approach.

In the second case, optimization from NACA0012 to RAE2822 air-foil was performed with 18 design variables, applying the same adaptive criterion as in the first case. The optimization results show that the case using the adaptive parameterization method achieves a lower drag coefficient for the optimized airfoil compared to the uniform parameterization case while satisfying the optimization objective. This also illustrates that the adaptive parameterization method, through its stronger control capability over local shapes, can fine-tune the airfoil and capture more optimal design points. In conclusion, the proposed adaptive FFD parameterization method demonstrates strong adaptability. It provides stability and efficiency in early optimization stages with fewer design variables, while its adaptive adjustments in later stages facilitate the search for more optimal design points.

Author Contributions

Conceptualization, J.Z.; Methodology, J.Z.; Software, J.Z.; Validation, J.Z.; Resources, X.W.; Writing—original draft, J.Z.; Writing—review & editing, H.J.; Supervision, J.Y. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. B-spline curve fitting.

Figure 2. The curve and its curvature of the segmented function defined by Equation (5).

Figure 3. Fitting result of the B-spline curve with 10 control points, uniform open, and degree 3.

Figure 4. Fitting result using adaptive control point B-spline curve fitting with 10 control points, uniform open, and degree 3.

Figure 5. Diagram of adaptive FFD control points based on spring analogy method.

Figure 6. Optimization frame of Control Point Adaptive Parameterization based on FFD.

Figure 7. Original mesh for RAE2822 optimization.

Figure 8. FFD control points enclosing the RAE2822 airfoil.

Figure 9. Convergence history of the drag coefficient for these three optimization cases with uniform FFD control points.

Figure 10. Pressure distribution and airfoil of original and optimized case with uniform FFD control points ((Left): Pressure distribution, (Right): Airfoil).

Figure 11. Control points position of three optimized cases.

Figure 12. Convergence history of the drag coefficient for these three optimization cases with uniform FFD control points.

Figure 13. Adaptive FFD control points and optimized result.

Figure 14. Pressure distribution and airfoil of original and optimized case with adaptive FFD control points ((Left): Pressure distribution, (Right): Airfoil).

Figure 15. Iteration history of drag coefficient and moment coefficient for optimization from NACA0012 to RAE2822 airfoil ((Left): First optimization stage, (Right):Second optimization stage).

Figure 16. Airfoil evolution in the two optimization design processes.

Figure 17. Convergence history of the drag coefficient for the adaptive parameterization and uniform parameterization optimization design cases.

Figure 18. Optimized airfoils and pressure distributions for the two optimization cases ((Left): Pressure distribution, (Right): Airfoil).

Table 1

Optimization Problem.

	Function/Variable Description		Quantity
Minimize	$C_D$	Drag coefficient	-
With respect to	$\alpha$	Angle of attack	1/1/1
	z	FFD control point z-coordinates	6/12/24
Subject to	$C_L$ = 0.824	Lift coefficient constraint	1
	$C_m>-0.092$	Moment coefficient constraint	1
	$A\ge A_{base}$	Area constraint	1

Table 2

Mesh convergence.

Level	Mesh Cells Quantity	$\mathit{\alpha }{/}^{\circ }({\mathit{C}}_{\mathit{L}}=0.824)$	${\mathit{C}}_{\mathit{D}}$	${\mathit{C}}_{\mathit{m}}$
1	3205	3.01	0.022083	−0.093304
2	7698	2.94	0.021246	−0.094703
3	12,537	3.24	0.023112	−0.085896
4	23,426	3.22	0.022827	−0.086603
5	29,812	3.20	0.022580	−0.086949

References

1. Martins, J.R. Aerodynamic design optimization: Challenges and perspectives. Comput. Fluids; 2022; 239, 105391. [DOI: https://dx.doi.org/10.1016/j.compfluid.2022.105391]

2. Kim, J.E.; Rao, V.N.; Koomullil, R.P.; Ross, D.H.; Soni, B.K.; Shih, A.M. Development of an efficient aerodynamic shape optimization framework. Math. Comput. Simul.; 2009; 79, pp. 2373-2384. [DOI: https://dx.doi.org/10.1016/j.matcom.2009.01.012]

3. Nejat, A.; Mirzabeygi, P.; Shariat Panahi, M. Airfoil shape optimization using improved Multiobjective Territorial Particle Swarm algorithm with the objective of improving stall characteristics. Struct. Multidiscip. Optim.; 2014; 49, pp. 953-967. [DOI: https://dx.doi.org/10.1007/s00158-013-1025-3]

4. Yiğit, S.; Abuhanieh, S.; Biçer, B. An Open-Source Aerodynamic Shape Optimization Application for an Unmanned Aerial Vehicle (UAV) Propeller: An open-source aerodynamic shape optimization application. J. Aeronaut. Space Technol.; 2022; 15, pp. 1-12.

5. Li, J.; Du, X.; Martins, J.R. Machine learning in aerodynamic shape optimization. Prog. Aerosp. Sci.; 2022; 134, 100849. [DOI: https://dx.doi.org/10.1016/j.paerosci.2022.100849]

6. Raul, V.; Leifsson, L. Surrogate-based aerodynamic shape optimization for delaying airfoil dynamic stall using Kriging regression and infill criteria. Aerosp. Sci. Technol.; 2021; 111, 106555. [DOI: https://dx.doi.org/10.1016/j.ast.2021.106555]

7. Han, Z.H.; Abu-Zurayk, M.; Görtz, S.; Ilic, C. Surrogate-based aerodynamic shape optimization of a wing-body transport aircraft configuration. Proceedings of the AeroStruct: Enable and Learn How to Integrate Flexibility in Design: Contributions to the Closing Symposium of the German Research Initiative AeroStruct; Braunschweig, Germany, 13–14 October 2015; Springer: Berlin/Heidelberg, Germany, 2018; pp. 257-282.

8. Skinner, S.N.; Zare-Behtash, H. State-of-the-art in aerodynamic shape optimisation methods. Appl. Soft Comput.; 2018; 62, pp. 933-962. [DOI: https://dx.doi.org/10.1016/j.asoc.2017.09.030]

9. Kenway, G.K.; Martins, J.R. Multipoint aerodynamic shape optimization investigations of the common research model wing. AIAA J.; 2016; 54, pp. 113-128. [DOI: https://dx.doi.org/10.2514/1.J054154]

10. Lyu, Z.; Martins, J.R. Aerodynamic shape optimization of an adaptive morphing trailing-edge wing. J. Aircr.; 2015; 52, pp. 1951-1970. [DOI: https://dx.doi.org/10.2514/1.C033116]

11. Jameson, A.; Kim, S. Reduction of the adjoint gradient formula for aerodynamic shape optimization problems. AIAA J.; 2003; 41, pp. 2114-2129. [DOI: https://dx.doi.org/10.2514/2.6830]

12. He, X.; Li, J.; Mader, C.A.; Yildirim, A.; Martins, J.R. Robust aerodynamic shape optimization—from a circle to an airfoil. Aerosp. Sci. Technol.; 2019; 87, pp. 48-61. [DOI: https://dx.doi.org/10.1016/j.ast.2019.01.051]

13. Masters, D.A.; Taylor, N.J.; Rendall, T.; Allen, C.B.; Poole, D.J. Geometric comparison of aerofoil shape parameterization methods. AIAA J.; 2017; 55, pp. 1575-1589. [DOI: https://dx.doi.org/10.2514/1.J054943]

14. Wu, X.; Zhang, W.; Peng, X.; Wang, Z. Benchmark aerodynamic shape optimization with the POD-based CST airfoil parametric method. Aerosp. Sci. Technol.; 2019; 84, pp. 632-640. [DOI: https://dx.doi.org/10.1016/j.ast.2018.08.005]

15. Anderson, G.R.; Aftosmis, M.J. Adaptive Shape Parameterization for Aerodynamic Design; NASA Technical Report NAS-2015-02; National Aeronautics and Space Administration, Ames Research Center: Moffett Field, CA, USA, 2015.

16. Désidéri, J.A.; Abou El Majd, B.; Janka, A. Nested and self-adaptive Bézier parameterizations for shape optimization. J. Comput. Phys.; 2007; 224, pp. 117-131. [DOI: https://dx.doi.org/10.1016/j.jcp.2006.12.016]

17. Abou El Majd, B.; Duvigneau, R.; Désidéri, J.A. Aerodynamic shape optimization using a full and adaptive multilevel algorithm. Proceedings of the ERCOFTAC Conference Design Optimization: Methods and Applications; Canary Island, Spain, 31 March–2 April 2004.

18. Han, X.; Zingg, D.W. An adaptive geometry parametrization for aerodynamic shape optimization. Optim. Eng.; 2014; 15, pp. 69-91. [DOI: https://dx.doi.org/10.1007/s11081-013-9213-y]

19. Yu, Y.; Lyu, Z.; Xu, Z.; Martins, J.R. On the influence of optimization algorithm and initial design on wing aerodynamic shape optimization. Aerosp. Sci. Technol.; 2018; 75, pp. 183-199. [DOI: https://dx.doi.org/10.1016/j.ast.2018.01.016]

20. Yang, G.; Da Ronch, A.; Drofelnik, J.; Xie, Z.T. Sensitivity assessment of optimal solution in aerodynamic design optimisation using SU2. Aerosp. Sci. Technol.; 2018; 81, pp. 362-374. [DOI: https://dx.doi.org/10.1016/j.ast.2018.08.012]

21. Allmaras, S.R.; Johnson, F.T. Modifications and clarifications for the implementation of the Spalart–Allmaras turbulence model. Proceedings of the Seventh International Conference on Computational Fluid Dynamics (ICCFD7); Big Island, HI, USA, 9–13 July 2012; Volume 1902, pp. 1-11.

22. Shima, E.; Kitamura, K.; Haga, T. Green–gauss/weighted-least-squares hybrid gradient reconstruction for arbitrary polyhedra unstructured grids. AIAA J.; 2013; 51, pp. 2740-2747. [DOI: https://dx.doi.org/10.2514/1.J052095]

23. Saad, Y. A flexible inner-outer preconditioned GMRES algorithm. SIAM J. Sci. Comput.; 1993; 14, pp. 461-469. [DOI: https://dx.doi.org/10.1137/0914028]

24. Lyu, Z.; Xu, Z.; Martins, J. Benchmarking optimization algorithms for wing aerodynamic design optimization. Proceedings of the 8th International Conference on Computational Fluid Dynamics; Chengdu, China, 11–13 November 2014; Volume 11, 585.

25. Gnoffo, P. A finite-volume, adaptive grid algorithm applied to planetary entry flowfields. AIAA J.; 1983; 21, pp. 1249-1254. [DOI: https://dx.doi.org/10.2514/3.8236]

26. Nakahashi, K.; Deiwert, G.S. Self-adaptive-grid method with application to airfoil flow. AIAA J.; 1987; 25, pp. 513-520. [DOI: https://dx.doi.org/10.2514/3.9655]