Airfoil Analysis and Optimization Using a Petrov–Galerkin Finite Element and Machine Learning

Abstract

For the analysis of low-speed incompressible fluid dynamics with turbulence around airfoils, we developed a finite element formulation based on a stabilized pressure and velocity formulation. To shape the optimization of bidimensional airfoils, this formulation is applied using machine learning (TensorFlow) and public domain global optimization algorithms. The goal is to maximize the lift-over-drag ratio by using the class-shape function transformation (CST) parameterization technique and machine learning. Specifically, we propose equal-order stabilized three-node triangles for the flow problem, standard three-node triangles for the approximate distance function (ADF) required in the turbulence stage, and stabilized three-node triangles for the Spalart–Allmaras turbulence model. The backward Euler time integration was employed. An implicit time-integration algorithm was adopted, and a solution was obtained using the Newton–Raphson method. This was made possible in the symbolic form via Mathematica with the AceGen package. Three benchmarks are presented, with Reynolds numbers up to $1\times {10}^7$, demonstrating remarkable robustness. After the assessment of the new finite element, we used machine learning and global optimization for four angles of attack to calculate airfoil designs that maximized $C_L/C_D$.

1. Introduction

Research on airfoil shape optimization with classical gradient-based methods has been published for nearly five decades [1]. With the advent of automatic differentiation algorithms for machine learning (see the review article [2]) and backpropagation, it became advantageous to use machine learning algorithms for problems where sensitivity analyses were intricate, see, e.g., [3].

For airfoil optimization, it is customary to use finite-volume software to perform the CFD analysis [4], but PINN has also been successfully adopted [5]. The finite element (FE) choice in this work is supported by the need to further increase the application range of the simulation to a fluid–structure interaction (FSI).

We focus on three main components in the approach: (i) airfoil parameterization, (ii) CFD discretization and time integration, (iii) and machine learning/optimization. In [6], a TLBO (teaching–learning-based) optimization algorithm was proposed for shape optimization. The authors used class-shape function transformation (CST) parameterization for 3D numerical tests. For CFD analysis, a finite-volume method was adopted. In the work by Deng and Yi [7], a machine learning algorithm was trained to generate pressure distributions. Their objective was analysis efficiency. Latin hypercubes and CST parameterization were adopted for 2D numerical tests. In the work by Du et al. [8], a convolutional neural network algorithm for airfoil design and performance prediction (DPCNN) was established, A Bézier representation of the airfoil geometry was adopted, and the UIUC database was adopted in the training process [9]. In [8], ANSYS Fluent software was used to analyze the aerodynamic performance of the airfoil. In that paper, optimization was performed in the machine learning framework to maximize the function denoted by the ratio of the lift coefficient to the drag coefficient $C_L/C_D$. In Hui’s study [10], reinforced learning was combined with multiobjective optimization. A constraint on the airfoil thickness was enforced to prevent the thinning observed in the current paper. A free-form deformation (FFD) algorithm was adopted for the shape motion and the CFL3D CFD code was employed for the simulations. In Karali et al. [11], a lifting line method was introduced for modeling and optimizing small unmanned vehicles. In a similar approach to that of this work, Karali et al. tested a range of angles of attack and provided a black box function for both analysis and optimization. TensorFlow functions and parameters are similar to other contributions, but only 250 epochs were used. In [12], a deep convolutional generative adversarial network was adopted, which produced more realistic airfoils than normal, and a detector of geometric abnormalities was introduced. FFD parametrization was adopted in [12], with the optimizer being a variant of the sequential quadratic programming (SQP) algorithm. Tyan et al. [13] proposed an inverse design deep neural network (IDNN), dimensional reduction, and the NACA four-series airfoil geometry representation via 2D examples. In an alternative approach to optimization, Xu [14] proposed machine learning (ML) for the adjoint-variable method. In a parallel line of developments, it is worth mentioning the work by Gutierrez et al. [15], which focused on the optimization of propellers for high-altitude pseudo-satellites. More sophisticated modeling approaches, such as physics-informed neural networks (PINNs) [16], have also been adopted. A review of the state of the art was provided by Li, Du, and Martins [17], who also discussed filtering and dimensional reduction techniques.

Here, we propose distinct approaches: the panel solver XFOIL [18] and the mixed-Finite Element software SIMPLAS [19] (developed by the first author). CST parameterization with 10 parameters was adopted for each airfoil, which was read from the UIUC database [9], and the control volume was meshed using Triangle [20]; both SIMPLAS and XFOIL produced the $C_L/C_D$ values. Four values of the attack angle were tested for each airfoil: $\alpha =0^{\circ }$, $1.{25}^{\circ }$, $2.5^{\circ }$ and $5^{\circ }$. After this, the TensorFlow machine learning library was employed to produce a function that returned $C_L/C_D$ from the 10 CST parameters. A dual annealing algorithm was then employed to estimate the optimized airfoil.

This paper, which is aimed at optimizing recreational aircraft airfoils, is organized as follows: Section 2 discusses the CST parameterization. Section 3 focuses on a RANS finite element formulation based on equal-order velocity/pressure interpolation and stabilization terms. This is a variant of the Petrov–Galerkin formulation, which was created by the first author. In Section 4, the Spalart–Allmaras transport equation for the eddy viscosity is described and discretized. The machine learning algorithm for the solution and optimization is described in Section 5. The results are presented in Section 6.

2. Parameterization

The parameterization of airfoil geometry, although not strictly required from the simulation perspective, reduces the required information needed to feed the machine learning algorithm [21]. Among the established parameterization algorithms, the cubic version of the CST [22] (class function/shape function transformation method) algorithm has been profusely adopted. Although CST is not suitable for all profiles in the repositories, the greatest advantages of CST are as follows: (i) it is applicable to 3D airfoils and (ii) it only requires a small number of parameters. The latter is of special importance, as only a small number of profiles are available in the databases [9]. For traditional optimization algorithms, see, e.g., [23], the CST parameterization is also commonly adopted due to its convenience. Figure 1 shows the essential nomenclature of an airfoil.

We directly implemented the CST functions proposed by Kulfan [22], which are capable of analytically defining a variety of airfoils with a reduced number of parameters. CST equations are based on a Bézier curve with an added class function term. Lower and upper surfaces are defined, respectively, by:

$$\begin{array}{cc}\hfill {\zeta }_l\left(\chi \right)& ={\chi }^{0.5}\left(1-\chi \right)S_l\left(\chi \right)+\chi \mathrm{\Delta }{\zeta }_l\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\zeta }_u\left(\chi \right)& ={\chi }^{0.5}\left(1-\chi \right)S_u\left(\chi \right)+\chi \mathrm{\Delta }{\zeta }_u\hfill \end{array}$$

(2)

where ${\zeta }_l$ is the upper surface, ${\zeta }_u$ is the lower surface, $\mathrm{\Delta }{\zeta }_u$ and $\mathrm{\Delta }{\zeta }_L$ define the thicknesses of the trailing edge, and $S_u\left(\chi \right)$ and $S_l\left(\chi \right)$ are the shape functions obtained from the Bernstein polynomials, ${\mathit{N}}_u\left(\chi \right)$ ($N_u+1-$dimensional) and ${\mathit{N}}_l\left(\chi \right)$ ($N_l+1-$dimensional), for the upper and lower surfaces, respectively:

$$\begin{array}{cc}\hfill {\mathit{N}}_l\left(\chi \right)& =\left\{\cdots ,\frac{N_l!}{i!\left(N_l-i\right)!}{\chi }^i{\left(1-\chi \right)}^{N_l-i},\cdots \right\},i=0,\dots ,N_l\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\mathit{N}}_u\left(\chi \right)& =\left\{\cdots ,\frac{N_u!}{i!\left(N_u-i\right)!}{\chi }^i{\left(1-\chi \right)}^{N_u-i},\cdots \right\},i=0,\dots ,N_u\hfill \end{array}$$

(4)

The dimensionless coordinates are defined by:

$$\left\{\begin{array}{c}\chi \\ \zeta \end{array}\right\}=\frac{1}{c}\left\{\begin{array}{c}x\\ z\end{array}\right\}$$

(5)

where c represents the chord length, x represents the coordinate along the chord length, and z represents the coordinate along the profile thickness. We now use all nodal images of ${\zeta }_l$ as ${\mathit{\zeta }}_{lN}$ and of ${\zeta }_u$ as ${\mathit{\zeta }}_{uN}$. Solving the NU and NL equations and applying the method of least squares, we obtain

$$\begin{array}{cc}\hfill {\mathit{A}}_l& ={\left({\mathit{C}}_l^T·{\mathit{C}}_l\right)}^{-1}·{\mathit{C}}_l^T·\left({\mathit{\zeta }}_{lN}-{\mathbf{\chi }}_{lN}\mathrm{\Delta }{\zeta }_l\right)\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\mathit{A}}_u& ={\left({\mathit{C}}_u^T·{\mathit{C}}_u\right)}^{-1}·{\mathit{C}}_u^T·\left({\mathit{\zeta }}_{uN}-{\mathbf{\chi }}_{uN}\mathrm{\Delta }{\zeta }_u\right)\hfill \end{array}$$

(7)

So now we have everything to obtain the following:

$$\begin{array}{cc}\hfill S_l\left(\chi \right)& ={\mathit{N}}_l\left(\chi \right)·{\mathit{A}}_l\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill S_u\left(\chi \right)& ={\mathit{N}}_u\left(\chi \right)·{\mathit{A}}_u\hfill \end{array}$$

(9)

Finally, we obtain the parameters related to the thicknesses of the lower and upper surfaces, respectively:

$$\begin{array}{cc}\hfill z_l\left(\chi \right)& =c{\chi }^{0.5}\left(1-\chi \right){\mathit{N}}_l\left(\chi \right)·{\mathit{A}}_l+c\chi \mathrm{\Delta }{\zeta }_l\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill z_u\left(\chi \right)& =c{\chi }^{0.5}\left(1-\chi \right){\mathit{N}}_u\left(\chi \right)·{\mathit{A}}_u+c\chi \mathrm{\Delta }{\zeta }_u\hfill \end{array}$$

(11)

3. RANS: Consistent Stabilized Petrov–Galerkin Formulation with Equal-Order Interpolation

The recreational aircraft airfoil flow is clearly within the domain $Ma\in [0,0.3]$, which means it can be modeled with the following constitutive assumptions (see, e.g., [24]): (i) incompressible flow and (ii) turbulent behavior. SIMPLAS, developed by P. Areias, is used; see [19]. The specific elements and turbulence model are created in Wolfram Mathematica [25] with the AceGen add-on [26]. We use the following tools for the solution of the flow problem:

• Equal-order stabilized three-node triangles for the flow problem.
• Standard three-node triangles for the distance function (ADF) required in the turbulence stage.
• Stabilized three-node triangles for the Spalart–Allmaras turbulence model.
• Backward Euler time integration is employed.
• A fully-implicit algorithm is adopted, and the solution is obtained using the Newton–Raphson method.

We now consider an incompressible flow with Newtonian constitutive behavior and incorporate the effect of turbulence. Given a fixed control volume $\mathrm{\Omega }$, and for each point $\mathit{x}\in \mathrm{\Omega }$, we have a velocity field $\mathit{v}\left(\mathit{x}\right)$, which describes the fluid flow. A pressure field $p\left(\mathit{x}\right)$ is necessary, as it is a Lagrange multiplier for the incompressibility condition. Mass density is identified as usual by ${\rho }_0$ and the volume force by $f_i$. Omitting the dependence on $\mathit{x}$, Newtonian fluids follow the constitutive law:

$${\sigma }_{ij}=2\mu {\dot{\epsilon }}_{ij}-p{\delta }_{ij}$$

where ${\sigma }_{ij}$ is the Cauchy stress tensor, $\mu$ is the dynamic viscosity, $v_{i,j}=\partial v_i\left(\mathit{x}\right)/\partial x_j$, with ${\dot{\epsilon }}_{ij}=1/2\left(v_{i,j}+v_{j,i}\right)$ being the strain rate. The equations of motion and the incompressibility condition are given, respectively, by:

$$\begin{array}{cc}\hfill {\rho }_0{\mathring{v}}_i& +{\left(p{\delta }_{ij}-2\mu {\dot{\epsilon }}_{ij}\right)}_{,j}-{\rho }_0{f}_i=0\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill v_{i,i}& =0\hfill \end{array}$$

(13)

where

$${\mathring{v}}_i={\dot{v}}_i+v_{i,j}{v}_j.$$

(14)

In general, the convective time-derivative (14) corresponds to the following operator on a generic argument $\mathsf{x}$:

$$\mathring{\mathsf{x}}=\dot{\mathsf{x}}+\nabla \mathsf{x}·\mathit{v}$$

(15)

where ${\left[\nabla \mathsf{x}\right]}_i=\partial \mathsf{x}/\partial x_i$. Boundary conditions for the flow problem are essential in ${\mathrm{\Gamma }}_v$ and natural in ${\mathrm{\Gamma }}_{\sigma }$:

$$\begin{array}{cc}\hfill v_i& =f_i^v\left(t\right)\mathrm{on}{\mathrm{\Gamma }}_v\hfill \\ \hfill {\sigma }_{ij}\left(t\right)n_j& =f_i^{\tau }\left(t\right)\mathrm{on}{\mathrm{\Gamma }}_{\sigma }\hfill \end{array}$$

Using the test functions $q\left(\mathit{x}\right)$ and $w\left(\mathit{x}\right)$, we project (12) and (13) to obtain the weak form:

$$\begin{array}{cc}\hfill {\int }_{\mathrm{\Omega }}{\rho }_0{w}_i{\mathring{v}}_i\mathrm{d}V& +{\int }_{\mathrm{\Omega }}\left(2\mu w_{i,j}{\dot{\epsilon }}_{ij}-p{w}_{i,i}-{\rho }_0{w}_i{f}_i\right)\mathrm{d}V-{\int }_{\mathrm{\Gamma }}{w}_i{f}_i^{\tau }\mathrm{d}A=0\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\int }_{\mathrm{\Omega }}q{v}_{i,i}\mathrm{d}V& =0\hfill \end{array}$$

(17)

with $v_i=f_i^v\left(t\right)\mathrm{on}{\mathrm{\Gamma }}_v$. Finite element discretization using equal-order interpolation (the pressure and velocity share the shape functions) is performed as described in T.E. Tezduyar’s papers [27,28]. We consider an element with domain ${\mathrm{\Omega }}_e$ and boundary ${\mathrm{\Gamma }}_e$ with shape functions $N_K^v\left(\mathit{\xi }\right)$ for $v_i$ and $w_i$ and $N_K^p\left(\mathit{\xi }\right)$ for p and q:

$$\begin{array}{cc}\hfill & {\int }_{{\mathrm{\Omega }}_e}{w}_{Ki}\left[{\rho }_0{N}_K^v{\mathring{v}}_i+N_{K,j}^v\left(2\mu {\dot{\epsilon }}_{ij}-p{\delta }_{ij}\right)-{\rho }_0{N}_K^v{f}_i\right]\mathrm{d}{V}_e\hfill \end{array}$$

$$\begin{array}{cc}\hfill & +{\int }_{{\mathrm{\Gamma }}_e}{w}_{Ki}{N}_K^v{f}_i^{\tau }\mathrm{d}{A}_e=0\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & {\int }_{{\mathrm{\Omega }}_e}{q}_K{N}_K^p{v}_{i,i}\mathrm{d}{V}_e=0\hfill \end{array}$$

(19)

where

$$\begin{array}{cc}\hfill {\mathring{v}}_i& =N_L^v{\dot{v}}_{Li}+N_M^v{N}_{L,j}^v{v}_{Mj}{v}_{Li}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\dot{\epsilon }}_{ij}& =\frac{1}{2}\left(N_{L,j}^v{v}_{Li}+N_{L,i}^v{v}_{Lj}\right)\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill v_{i,i}& =N_{K,i}^v{v}_{Ki}\hfill \end{array}$$

(22)

in (18) and (19), $V_e$ is the volume of element e, which is assumed to have unit thickness. With a similar notation, $A_e$ is the area of the boundary where natural conditions are applied. In a steady-state analysis, which is used for the optimization stage, the following condition is enforced:

$${\mathring{v}}_i\cong v_{i,j}{v}_j$$

(23)

which reduces the acceleration term to be equal to the transport term. We note that stabilization is required for (18) and (19) and, therefore, a non-symmetric term is inserted, which is a penalty for the strong form of the equations of motion:

$$r_{\mathrm{gls}}={\int }_{{\mathrm{\Omega }}_e}\tau \left(\mathit{v}\right)\left({\rho }_0{\mathring{w}}_i+q_{,i}\right)\left({\rho }_0{\mathring{v}}_i+p_{,i}-{\rho }_0{f}_i\right)\mathrm{d}{V}_e$$

(24)

with $\tau \left(\mathit{v}\right)$ being the penalty term, depending on the velocity $\mathit{v}$:

$$\tau \left(\mathit{v}\right)=\frac{1}{{\rho }_0}{\left[{\left(\frac{2}{\mathrm{\Delta }t}\right)}^2+{\left(\frac{2\parallel \mathit{v}\parallel }{h}\right)}^2+{\left(\frac{4\nu }{h^2}\right)}^2\right]}^{-1/2}$$

(25)

An important aspect of the penalty term (25) is that, for high values of the Reynolds number, the iterative velocity $\mathit{v}$ must be used and not the converged value. In [29], a remark is made concerning the use of converged values, and we notice the loss of the Newton–Raphson convergence in the presence of high Reynolds numbers. The shortcoming of using iterative velocities in the stabilization parameter is the required derivative. Since we are performing these calculations with Mathematica [25] and AceGen [26], this additional task is easily performed here. This stabilization and several variants have been extensively studied and benchmarked by T.E. Tezduyar and co-workers (see [27,28,30]). It is now an established solution for the finite element analysis of fluid flow problems. This and alternative techniques are discussed in standard textbooks [31]. In (25), the kinematic viscosity is given by:

$${\displaystyle \nu =\frac{\mu }{{\rho }_0}}$$

(26)

This results in a Petrov–Galerkin formulation that allows equal-order interpolation, which otherwise would fail the Babuska–Brezzi condition. After introducing the stabilization term and performing the discretization with shape functions $N_K^v\left(\mathit{\xi }\right)$ for the velocities and $N_L^p\left(\mathit{\xi }\right)$ for the pressures (which coincide in this case, but we retain generality), the following discrete form is obtained:

$$\begin{array}{cc}\hfill & {\int }_{{\mathrm{\Omega }}_e}{w}_{Ki}{\rho }_0{N}_K^v{\mathring{v}}_i\mathrm{d}{V}_e\hfill \\ \hfill & +{\int }_{{\mathrm{\Omega }}_e}{w}_{Ki}\left[N_{K,j}^v\left(2\mu {\dot{\epsilon }}_{ij}-p{\delta }_{ij}\right)-{\rho }_0{N}_K^v{f}_i\right]\mathrm{d}{V}_e\hfill \\ \hfill & +{\int }_{{\mathrm{\Omega }}_e}{w}_{Ki}\tau \left(\mathit{v}\right){\rho }_0\left(\partial \dot{w}/\partial w{N}_K^v+N_{K,l}^v{v}_l\right)\left({\rho }_0{\mathring{v}}_i+p_{,i}-{\rho }_0{f}_i\right)\mathrm{d}{V}_e\hfill \\ \hfill & -{\int }_{{\mathrm{\Gamma }}_e}{w}_{Ki}{N}_K^v{f}_i^{\tau }\mathrm{d}{A}_e\hfill & \hfill =0\end{array}$$

(27)

$$\begin{array}{ccc}\hfill & {\int }_{{\mathrm{\Omega }}_e}{q}_K{N}_K^p{v}_{i,i}\mathrm{d}{V}_e+{\int }_{{\mathrm{\Omega }}_e}{q}_K{N}_{K,i}^p\tau \left(\mathit{v}\right)\left({\rho }_0{\mathring{v}}_i+p_{,i}-{\rho }_0{f}_i\right)\mathrm{d}{V}_e\hfill & \hfill =0\end{array}$$

(28)

From (27) and (28), we calculate the forces ${\mathit{f}}^w$ (for the velocity degree-of-freedom) and ${\mathit{f}}^q$ (for the pressure degree-of-freedom). Since the Newton–Raphson method is adopted to solve the discrete system, the Jacobian matrix is required. The Jacobian matrix, ${\mathit{K}}^e$, is obtained as follows:

$${\mathit{K}}^e=\left[\begin{array}{cc}{\left[{\mathit{K}}^e\right]}^{wv}& {\left[{\mathit{K}}^e\right]}^{wp}\\ {\left[{\mathit{K}}^e\right]}^{qv}& {\left[{\mathit{K}}^e\right]}^{qp}\end{array}\right]=\left[\begin{array}{cc}-\partial {\mathit{f}}^w/\partial \mathit{v}& -\partial {\mathit{f}}^w/\partial \mathit{p}\\ -\partial {\mathit{f}}^q/\partial \mathit{v}& -\partial {\mathit{f}}^q/\partial \mathit{p}\end{array}\right]$$

(29)

with

$$\begin{array}{cc}\hfill {\left[{\mathit{K}}^e\right]}_{KiLj}^{wv}& ={\int }_{{\mathrm{\Omega }}_e}{\rho }_0{N}_K^v\left[{\delta }_{ij}\partial \dot{v}/\partial v{N}_L^v+N_L^v{v}_{i,j}+{\delta }_{ij}{N}_{L,l}^v{v}_l\right]\mathrm{d}{V}_e\hfill \\ \hfill & +{\int }_{{\mathrm{\Omega }}_e}\mu \left(N_{K,i}^v{N}_{L,j}^v+N_{K,j}^v{N}_{L,i}^v\right)\mathrm{d}{V}_e+{\left[{\mathit{K}}_{\mathrm{stab}}^e\right]}_{KiLj}^{wv}\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\left[{\mathit{K}}^e\right]}_{KiL}^{wp}& =-{\int }_{{\mathrm{\Omega }}_e}{N}_{K,i}^v{N}_L^p\mathrm{d}{V}_e+{\left[{\mathit{K}}_{\mathrm{stab}}^e\right]}_{KiL}^{wp}\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\left[{\mathit{K}}^e\right]}_{KLi}^{qv}& =-{\int }_{{\mathrm{\Omega }}_e}{N}_K^p{N}_{L,i}^p\mathrm{d}{V}_e+{\left[{\mathit{K}}_{\mathrm{stab}}^e\right]}_{KLi}^{qv}\hfill \\ \hfill {\left[{\mathit{K}}^e\right]}_{KL}^{qp}& ={\left[{\mathit{K}}_{\mathrm{stab}}^e\right]}_{KL}^{qp}\hfill \end{array}$$

(32)

where the four stabilization matrices are calculated using Mathematica [25] with the AceGen add-on [26].

$$\begin{array}{cc}\hfill {\left[{\mathit{K}}_{\mathrm{stab}}^e\right]}_{KiLj}^{wv}& ={\int }_{{\mathrm{\Omega }}_e}\tau \left(\mathit{v}\right){\rho }_0{N}_{K,j}^v{N}_L^v\left[{\rho }_0\left({\dot{v}}_i+v_{i,k}{v}_k\right)+p_{,i}-f_i\right]\mathrm{d}{V}_e\hfill \\ \hfill & +{\int }_{{\mathrm{\Omega }}_e}\tau \left(\mathit{v}\right){\rho }_0\left(\partial \dot{w}/\partial w{N}_K^v+N_{K,k}^v{v}_k\right){\rho }_0\left(\partial \dot{v}/\partial v{N}_L^v{\delta }_{ij}+N_{L,k}^v{v}_k{\delta }_{ij}+v_{i,j}{N}_L^v\right)\mathrm{d}{V}_e\hfill \\ \hfill & +{\int }_{{\mathrm{\Omega }}_e}\partial \tau \left(\mathit{v}\right)/\partial v_j{N}_L^v{\rho }_0\left(\partial \dot{w}/\partial w{N}_K^v+N_{K,k}^v{v}_k\right)\left\{{\rho }_0\left[{\dot{v}}_i+\left(v_{i,k}{v}_k\right)\right]+p_{,i}-f_i\right\}\mathrm{d}{V}_e\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\left[{\mathit{K}}_{\mathrm{stab}}^e\right]}_{KiL}^{wp}& ={\int }_{{\mathrm{\Omega }}_e}\tau \left(\mathit{v}\right){\rho }_0\left(\partial \dot{w}/\partial w{N}_K^v+N_{K,k}^v{v}_k\right){\delta }_{ij}{N}_{L,i}\mathrm{d}{V}_e\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill {\left[{\mathit{K}}_{\mathrm{stab}}^e\right]}_{KLj}^{qv}& ={\int }_{{\mathrm{\Omega }}_e}{N}_{K,i}^p\tau \left(\mathit{v}\right)\left\{{\rho }_0\left[{\dot{v}}_i\left(v_i\right)+v_{i,k}{v}_k\right]+p_{,i}-f_i\right\}\mathrm{d}{V}_e\hfill \\ \hfill & +{\int }_{{\mathrm{\Omega }}_e}{N}_{K,i}^p\tau \left(\mathit{v}\right)\left(\partial \dot{v}/\partial v{N}_L^v{\delta }_{ij}+N_{L,k}^v{v}_k{\delta }_{ij}+v_{i,j}{N}_L^v\right)\mathrm{d}{V}_e\hfill \\ \hfill & +{\int }_{{\mathrm{\Omega }}_e}{N}_{K,i}^p\partial \tau \left(\mathit{v}\right)/\partial v_j{N}_L^v\left\{{\rho }_0\left[{\dot{v}}_i\left(v_i\right)+\left(v_{i,k}{v}_k\right)\right]+p_{,i}-f_i\right\}\mathrm{d}{V}_e\hfill \end{array}$$

(35)

and

$$\begin{array}{cc}\hfill {\left[{\mathit{K}}_{\mathrm{stab}}^e\right]}_{KL}^{qp}& ={\int }_{{\mathrm{\Omega }}_e}{N}_{K,i}^p\tau \left(\mathit{v}\right)N_{L,i}^p\mathrm{d}{V}_e\hfill \end{array}$$

(36)

The derivative of $\tau \left(\mathit{v}\right)$ with respect to a component $v_i$ is calculated as follows:

$$\partial \tau \left(\mathit{v}\right)/\partial v_i=-\frac{4{\tau }^3}{h^2}{v}_i$$

(37)

A simple backward Euler (implicit) method is used for the time integration. Therefore, for step n, we have:

$${\dot{\mathit{v}}}^n\cong \frac{1}{\mathrm{\Delta }t}\left({\mathit{v}}^n-{\mathit{v}}^{n-1}\right)$$

(38)

This results, component-wise, as follows:

$$\begin{array}{cc}\hfill {\dot{v}}_i\left(v_i\right)& =v_i/\mathrm{\Delta }t+C_i\hfill \\ \hfill \mathrm{with}& C_i=-v_i^{n-1}/\mathrm{\Delta }t\hfill \end{array}$$

(39)

and, therefore,

$$\partial {\dot{v}}_i/\partial v_i=\partial {\dot{w}}_i/\partial w_i=1/\mathrm{\Delta }t$$

(40)

The source code for these equations is available in GitHub [32].

4. Turbulence with the Spalart–Allmaras Transport Equation

Turbulence models make use of the distance of a given point in the control volume to the wall, or to the airfoil. It is well known that, using the Poisson equation, Varadhan [33] established that an ADF can be generated by the solution of the following problem:

$$\begin{array}{cc}\hfill c_L^2{\nabla }^2\varphi \left(\mathit{x}\right)-\varphi \left(\mathit{x}\right)& =0\mathrm{in} \mathrm{\Omega }\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill \varphi \left(\mathit{x}\right)& =1\mathrm{on} \mathrm{\Gamma }\hfill \end{array}$$

(42)

with $g\left(\mathit{x}\right)=-c_{L}log\left[\varphi \left(\mathit{x}\right)\right]$ and $d\left(\mathit{x}\right)={lim}_{c_L\to 0}g\left(\mathit{x}\right).$ Therefore, we can use this ADF to estimate $d\left(\mathit{x}\right)$ for use in turbulence models. Note that, in contrast with geometric algorithms, the precise point corresponding to the orthogonal projection on $\mathrm{\Gamma }$ is not obtained, only the distance.

The Spalart–Allmaras turbulence model is a one-PDE turbulence model specific to aerodynamic flows, developed by two staff members of the Boeing group, see [34,35]. The model is shown to work well with a variety of turbulent flows, including boundary layers, separated flows, and wakes. It has been successfully used for aeroacoustic simulations and conjugated heat transfer simulations. We follow the guidelines provided by the NASA Langley Research Center Turbulence Modeling Resource [36]. A single equation is solved for $\tilde{\nu }$, which is the mapped kinematic eddy viscosity:

$$\begin{array}{cc}\hfill \mathring{\tilde{\nu }}=& c_{b1}\tilde{S}\tilde{\nu }+\frac{1}{\sigma }\left\{\nabla ·\left[\left(\nu +\tilde{\nu }\right)\nabla \tilde{\nu }\right]+c_{b2}\nabla \tilde{\nu }·\nabla \tilde{\nu }\right\}-c_{w1}{f}_w{\left(\frac{\tilde{\nu }}{d}\right)}^2+f_{t1}\mathrm{\Delta }\mathit{v}·\mathrm{\Delta }\mathit{v}\mathrm{in}\mathrm{\Omega }\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill \tilde{\nu }=& 0\mathrm{on}\mathrm{\Gamma }\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill {\left.\tilde{\nu }\right|}_{t=0}=& \nu /10\mathrm{in}\mathrm{\Omega }\hfill \end{array}$$

(45)

the reason to use $\tilde{\nu }$ instead of ${\nu }_t$ is that the behavior of the latter near the fixed boundaries is quartic, and a large number of elements would be required near the boundaries. The significance of terms in Equation (43) is as follows:

• The left-hand-side is the convective derivative of $\tilde{\nu }$, which provides the rate of change of $\tilde{\nu }$.
• The first term on the right-hand-side is the shear-driven turbulence generation, which occurs due to the velocity profile near solid boundaries.
• The second term on the right-hand-side is the nonlinear diffusion term, which has a linear part ($\nabla ·\left[\left(\nu +\tilde{\nu }\right)\nabla \tilde{\nu }\right]/\sigma$) and a nonlinear part ($c_{b2}\nabla \tilde{\nu }·\nabla \tilde{\nu }/\sigma$). The nonlinear part was fine-tuned using the parameter $c_{b2}$ to correctly represent the spreading of turbulence.
• The third term (with a negative sign) is the destruction term, whose purpose is to destruct turbulence near a wall. It obviously depends on the distance to the wall d.
• The fourth term is the source, which depends on the square of the velocity difference $\mathrm{\Delta }v$ and the vorticity, by means of $f_{t1}$. Note that for fixed walls, $\mathrm{\Delta }\mathit{v}=\mathit{v}$. Without this term, for homogeneous initial and boundary conditions, the equation would produce a null value of $\tilde{\nu }$.

The solution of (43) allows the determination of the turbulent eddy viscosity as follows:

$$\begin{array}{cc}\hfill {\mu }_t& =\rho \tilde{\nu }{f}_{v1}\mathrm{with}\hfill \\ \hfill f_{v1}& =\frac{{\chi }^3}{{\chi }^3+c_{v1}^3}\mathrm{and}\hfill \\ \hfill \chi & =\tilde{\nu }/\nu \hfill \end{array}$$

(46)

Note that the cube of $\chi$, multiplied by $\tilde{\nu }$, will produce the intended quartic profile of $\nu$ near the walls. The total viscosity is simply given as the following sum:

$${\mu }_T=\mu +{\mu }_t$$

(47)

This value is a replacement for $\mu$ in the equations of the previous section. The remaining quantities are determined as follows:

$$\begin{array}{cc}\hfill f_{t1}& =c_{t1}{g}_{t}exp\left[-c_{t2}{\parallel \mathit{W}\parallel }^2/{\parallel \mathit{v}\parallel }^2\right]\left[d^2\left(1+g_t^2\right)\right]\hfill \\ \hfill f_{t2}& =c_{t3}exp\left[-c_{t4}{\chi }^2\right]\hfill \\ \hfill f_{v2}& =1-\frac{\chi }{1+\chi f_{v1}}\hfill \\ \hfill f_w& =g{\left(\frac{1+c_{w3}^6}{g^6+c_{w3}^6}\right)}^{1/6}\hfill \\ \hfill g& =r+c_{w2}\left(r^6-r\right)\hfill \\ \hfill g_t& =min\left[0.1,\parallel \mathit{v}\parallel /h\parallel \mathit{W}\parallel \right]\hfill \\ \hfill r& =min\left[\frac{\tilde{\nu }}{\tilde{S}{k}^2{d}^2},10\right]\hfill \\ \hfill \tilde{S}& =max\left[\sqrt{2\mathit{W}:\mathit{W}}+\frac{\tilde{\nu }}{k^2{d}^2}{f}_{v2},0.3\parallel \mathit{W}\parallel \right]\hfill \\ \hfill \mathit{W}& =\frac{1}{2}\left[\left(\nabla \mathit{v}\right)-{\left(\nabla \mathit{v}\right)}^T\right]\hfill \end{array}$$

(48)

where the constants are given in [36], with d being the distance to the wall, as obtained from the ADF equation. The weak form of (43) is obtained, after integrating by parts, as follows:

$$\begin{array}{cc}\hfill & {\int }_{{\mathrm{\Omega }}_e}w\left(\mathring{\tilde{\nu }}-c_{b1}\tilde{S}\tilde{\nu }-f_{t1}\mathrm{\Delta }\mathit{v}·\mathrm{\Delta }\mathit{v}\right)\mathrm{d}{V}_e+{\int }_{{\mathrm{\Omega }}_e}w{c}_{w1}{f}_w{\left(\frac{\tilde{\nu }}{d}\right)}^2\mathrm{d}{V}_e\hfill \\ \hfill & -{\int }_{{\mathrm{\Omega }}_e}\frac{1}{\sigma }\left[-\left(\nu +\tilde{\nu }\right)\nabla w·\nabla \tilde{\nu }+c_{b2}w\nabla \tilde{\nu }·\nabla \tilde{\nu }\right]\mathrm{d}{V}_e=0\hfill \end{array}$$

(49)

where $w\left(\mathit{x}\right)$ is a test function. The stabilization of (49) follows the same pattern as before:

$$\begin{array}{cc}\hfill r_{\mathrm{gls}}^{\nu }& ={\int }_{{\mathrm{\Omega }}_e}\tau \left(\mathit{v}\right)\mathring{w}\left[\mathring{\tilde{\nu }}-c_{b1}\tilde{S}\tilde{\nu }+c_{w1}{f}_w{\left(\frac{\tilde{\nu }}{d}\right)}^2-f_{t1}\mathrm{\Delta }\mathit{v}·\mathrm{\Delta }\mathit{v}-\frac{c_{b2}}{\sigma }\nabla \tilde{\nu }·\nabla \tilde{\nu }\right]\mathrm{d}{V}_e\hfill \end{array}$$

(50)

This results in the stabilized weak form:

$$\begin{array}{cc}\hfill & {\int }_{{\mathrm{\Omega }}_e}w\left(\mathring{\tilde{\nu }}-c_{b1}\tilde{S}\tilde{\nu }\right)\mathrm{d}{V}_e\hfill \\ \hfill & +{\int }_{{\mathrm{\Omega }}_e}w{c}_{w1}{f}_w{\left(\frac{\tilde{\nu }}{d}\right)}^2\mathrm{d}{V}_e\hfill \\ \hfill & -{\int }_{{\mathrm{\Omega }}_e}\frac{1}{\sigma }\left[-\left(\nu +\tilde{\nu }\right)\nabla w·\nabla \tilde{\nu }+c_{b2}w\nabla \tilde{\nu }·\nabla \tilde{\nu }\right]\mathrm{d}{V}_e\hfill \\ \hfill & +{\int }_{{\mathrm{\Omega }}_e}\tau \left(\mathit{v}\right)\mathring{w}\left[\mathring{\tilde{\nu }}-c_{b1}\tilde{S}\tilde{\nu }-f_{t1}\mathrm{\Delta }\mathit{v}·\mathrm{\Delta }\mathit{v}+c_{w1}{f}_w{\left(\frac{\tilde{\nu }}{d}\right)}^2-\frac{c_{b2}}{\sigma }\nabla \tilde{\nu }·\nabla \tilde{\nu }\right]\mathrm{d}{V}_e\hfill \\ \hfill & =0\hfill \end{array}$$

(51)

After discretization, the residual reads:

$$\begin{array}{cc}\hfill f_K^{\nu }=& {\int }_{{\mathrm{\Omega }}_e}{N}_K\left(\mathring{\tilde{\nu }}-c_{b1}\tilde{S}\tilde{\nu }\right)\mathrm{d}{V}_e\hfill \\ \hfill & +{\int }_{{\mathrm{\Omega }}_e}{N}_K{c}_{w1}{f}_w{\left(\frac{\tilde{\nu }}{d}\right)}^2\mathrm{d}{V}_e\hfill \\ \hfill & -{\int }_{{\mathrm{\Omega }}_e}\frac{1}{\sigma }\left[-\left(\nu +\tilde{\nu }\right)N_{Ki}\frac{\partial \tilde{\nu }}{\partial x_i}+c_{b2}{N}_K\nabla \tilde{\nu }·\nabla \tilde{\nu }\right]\mathrm{d}{V}_e\hfill \\ \hfill & +{\int }_{{\mathrm{\Omega }}_e}\tau \left(\mathit{v}\right)\left(\frac{1}{\mathrm{\Delta }t}{N}_K+v_i{N}_{Ki}\right)\left[c_{w1}{f}_w{\left(\frac{\tilde{\nu }}{d}\right)}^2-\frac{c_{b2}}{\sigma }\nabla \tilde{\nu }·\nabla \tilde{\nu }\right]\mathrm{d}{V}_e\hfill \\ \hfill & +{\int }_{{\mathrm{\Omega }}_e}\tau \left(\mathit{v}\right)\left(\frac{1}{\mathrm{\Delta }t}{N}_K+v_i{N}_{Ki}\right)\left[\mathring{\tilde{\nu }}-c_{b1}\tilde{S}\tilde{\nu }-f_{t1}\mathrm{\Delta }\mathit{v}·\mathrm{\Delta }\mathit{v}\right]\mathrm{d}{V}_e\hfill \end{array}$$

(52)

The solution of (49) makes use of the Newton method, with  a very intricate Jacobian. The specific source code is available in GitHub [32] as files turbulence.nb and turbulence.f90. For the driven cavity, see Figure 2; streamlines for the turbulent case are shown in Figure 3.

We adopted a control volume consisting of a rectangle, including the airfoil geometry, as shown in Figure 4. The chord line has a 1 m length. For $\alpha =5^{\circ }$ and the optimized airfoil, we present the results in Figure 5 for the optimized airfoil. The velocity contour plots for the conditions are presented in Figure 4 along with the velocity lines in the neighborhood of the airfoils. A video of the optimized case flow for $\mathsf{Re}=1\times {10}^7$ is available on the video repository [37].

The convergence of results for the eddy viscosity $\tilde{\mu }$ and pressure p using three distinct meshes (17,400, 34,737, and 57,211 nodes) shows that smoothness increases with mesh refinement. Figure 6 shows the convergence of these fields.

5. Machine Learning and Optimization

The optimized airfoil optimization procedure consists of two stages: in the first stage, the relation between the CST parameters and the lift/drag ratio is established by machine learning. In the second stage, this function is maximized. In the first stage, we resort to the machine learning algorithm to actually construct a regression relating $C_L/C_D$ with 10 parameters defining the upper and lower faces of the airfoil. The optimization processes in the second stage are based on the dual annealing algorithm. We make use of the libraries TensorFlow [38,39], SciPy, and NumPy. Optimization is performed with the function SciPy.optimize.dual_annealing. It uses a generalization of simulated annealing with a final gradient-based local search stage. $C_L/C_D$ gradients are available from the machine learning model. In practice, the machine learning system makes use of a large set of input and output data, relates it statistically, and establishes the rules for regression. More data and better data will produce a better regression [39]. We note that, although possible (see, e.g., [10,40], the intent here is not to use full-feature machine learning to replicate the finite element simulations, but rather to extract the $C_L/C_D$ ratio.

In summary, the machine learning algorithm involves three parts:

• Parameter inputs, which are the 10 parameters of the CST representation for the complete set of UIUC airfoils.
• Outputs, which are images under a function of the inputs, in this case, $C_L/C_D$.
• Assessment method: this is required to measure the evolution of the iterations toward the result.

The data collection was performed with the flowchart in Figure 7. A total of 1598 airfoils from the UIUC database were considered valid for this application. Specific features of the machine learning procedure are as follows (see also [41]):

• CST [22] parameterization with 10 parameters, which are both the inputs in the learning process and the design variables.
• Data source: UIUC database [9].
• Four angles of attack were tested, 0, 1.25, 2.50, and 5.00 degrees.
• Set partitioning is as follows: 80% for training and 20% for testing.
• The type is “Regression Classification”.
• Activation function is “ReLU”.
• Kernel initializer is a uniform distribution with HeUniform class.
• Four dense layers.
• 6000 epochs.
• Loss function is “MAE”.
• Adam optimizer for the stochastic gradient descent.

We make use of the deep learning framework, the Keras CNN dense neural network, which is adequate for regression. The network diagram is shown in Figure 8 for 8 inputs (omitting $\mathrm{\Delta }{\zeta }_l$ and $\mathrm{\Delta }{\zeta }_u$).

The reasons for these options are as follows. The UIUC database contains a total of 1621 airfoils but we found that 23 were unsuitable for parameterization by the CST algorithm. Since at least a single CST patch is needed for each surface, the minimum number of CST parameters is 10. Further increasing the number of parameters would create difficulties for the machine learning stage. The number and values of angles of attack are aligned with what is typical in recreational airplanes, ensuring that simulations preclude stalling. The type of procedure is “Regression Classification” since the purpose of simulation here is to be reproduced by the machine learning model. The “ReLU” activation function is adopted due to the vanishing gradient issue that occurred with other functions. Since an optimization algorithm (with a final gradient-based stage) is run on top of the machine learning model, this advantage is important. In addition, due to sparsity, good performance is achieved with the “ReLU” activation function. The layer and epoch numbers are the results of the experimentation.

A flowchart for the algorithm following the data collection depicted in Figure 7 is shown in Figure 9.

6. Results and Discussion of Airfoil Optimization

With a lightweight airplane, such as the MC12 Cri-Cri airplane shown in Figure 10, efficiency is paramount, even at the expense of lift. High values of $C_L/C_D$ can be observed in recreational airplanes. By complying with this requirement, we are fulfilling two requirements, one is to maximize efficiency and the other is to avoid stalling. In this type of analysis, attention is focused on the behavior of the boundary layer, known as the layer of fluid close to a surface, a layer that strongly depends on the Reynolds number, the shape of the airfoil, and its inclination. Evaluating the boundary layer has the objective of reducing resistance measured by the resistance coefficient $C_D$, composed of two effects, pressure resistance ${\delta }^{*}$ and frictional resistance $C_L$ [42].

Using Keras/TensorFlow [38,39] combined with global optimization, it is possible to train and then optimize the value of $C_L/C_D$ using (i) XFOIL [18] results or (ii) FEM RANS [19]. In terms of optimization, the toolkit from SciPy was adopted with the function optimize.dual_annealing. We show the evolution of loss and MAE (mean absolute error) in TensorFlow [38], as functions of the epochs, from which we have 6000, in Figure 11 for the four considered angles of attack: $\alpha =0^{\circ },1.{25}^{\circ },2.{50}^{\circ }$ and $5^{\circ }$. The precision is sufficient for obtaining airfoil profiles in parametric form, which are superior to those in the UIUC database [9].

Optimized airfoils, each exhibiting positive curvature, are depicted in Figure 12. When analyzing the optimized airfoil for $\alpha =0^{\circ }$, the positive curvature stands out. This configuration, despite being obtained for $\alpha$ = $0^{\circ }$, subjects the airfoil to higher pressure on the lower surface compared to the upper surface, producing a positive lift. There is a certain bounce on the leading edge, which, together with the curvature and its thickness, alter the behavior of the boundary layer, trying to keep it close to the airfoil, inhibiting separation. It should be noted that, as the angle of attack increases, the curvature tends to not be as pronounced; this happens due to the progressive increase in the contribution of the angle of attack to the increase in lift, leaving the airfoil to have such a pronounced curvature to compensate for this need. This observation is even clearer for the $5^{\circ }$ case. It is known [42] that, at small angles of attack, thinner airfoils are preferable; at higher angles of attack, thicker airfoils are preferable. Comparing with the CST-parameterized versions of the UIUC database (pristine versions can be significantly different), both FEM and XFOIL show better $C_L$/$C_D$ results, see

Table 1. Optimized $C_L/C_D$ and corresponding $C_L$ for the optimized airfoils.

(a) ${\mathit{C}}_{\mathit{L}}/{\mathit{C}}_{\mathit{D}}$ for $\mathtt{Re}=1.08333\times {10}^6$
	Optimized
$\mathit{\alpha }$	XFOIL CST${\mathit{C}}_{\mathbf{L}}=$	FEM CST${\mathit{C}}_{\mathbf{L}}=$	Best-in-Database ${}^{★}$ ${\mathit{C}}_{\mathbf{L}}=$
0	$\begin{array}{|c|}\hline 160\\ \hline\end{array}$	$\begin{array}{|c|}\hline 160\\ \hline\end{array}$	105
$1.25$	$\begin{array}{|c|}\hline 223\\ \hline\end{array}$	222	219
$2.50$	190	$\begin{array}{|c|}\hline 200\\ \hline\end{array}$	179
$5.00$	$\begin{array}{|c|}\hline 174\\ \hline\end{array}$	$\begin{array}{|c|}\hline 174\\ \hline\end{array}$	173
(b) ${\mathit{C}}_{\mathbf{L}}$ for $\mathtt{Re}=\mathbf{1}.\mathbf{08333}\times {\mathbf{10}}^{\mathbf{6}}$
$\mathit{\alpha }$	XFOIL CST	FEM CST
0	$\begin{array}{|c|}\hline 0.98\\ \hline\end{array}$	$0.61$
$1.25$	$\begin{array}{|c|}\hline 1.12\\ \hline\end{array}$	$1.08$
$2.50$	$0.93$	$\begin{array}{|c|}\hline 1.04\\ \hline\end{array}$
$5.00$	$\begin{array}{|c|}\hline 1.49\\ \hline\end{array}$	$\begin{array}{|c|}\hline 1.49\\ \hline\end{array}$

, with similar performances for this target. In terms of $C_L$ alone, this table also shows that significant differences exist for very similar airfoils. Specific parameters are shown in

Table 2. Optimal parameters for each $\alpha$ (8 out of 10, since two parameters are always close to zero for the UIUC database [9]).

Method	$\mathit{\alpha }$	p
XFOIL	0	$\{-0.06220,+0.17723,+0.01487,+0.35473,+0.13862,+0.26888,+0.21769,+0.48022\}$
	$1.25$	$\{-0.05168,+0.16828,+0.02512,+0.36842,+0.15012,+0.25799,+0.20649,+0.46556\}$
	$2.50$	$\{-0.06417,+0.12845,+0.04564,+0.23788,+0.11283,+0.19201,+0.27182,+0.33103\}$
	$5.00$	$\{-0.11054,+0.12032,-0.07856,+0.26166,+0.29534,+0.37376,+0.49289,+0.31190\}$
FEM	0	$\{-0.03657,+0.09307,+0.01908,+0.14870,+0.11846,+0.20004,+0.19177,+0.28882\}$
	$1.25$	$\{-0.04173,+0.15715,+0.03088,+0.35488,+0.14013,+0.24550,+0.19497,+0.45997\}$
	$2.50$	$\{-0.05352,+0.13974,+0.03518,+0.25026,+0.12340,+0.19195,+0.28454,+0.31778\}$
	$5.00$	$\{-0.12164,+0.10912,-0.06778,+0.26281,+0.30829,+0.36001,+0.49893,+0.32502\}$

and can be tested in XFOIL by the reader.

For the four optimal airfoils, we use SIMPLAS [19] with uniformly imposed velocity and initial conditions; $\mathit{v}=15\{cos\alpha ,sin\alpha \}$. Meshes are refined locally, as shown in Figure 13.

Velocity profiles for all four optimized airfoils resulting from the FEM analysis are shown in Figure 14. We can observe the aforementioned effects. For the ADF $g\left(x\right)$, Figure 15 shows the results.

7. Conclusions

From a representation point-of-view, specifically in the case of $C_L/C_D$, as well as in the shape optimization of airfoils, we achieved great success. An unconditionally stable, fully implicit time integrator was adopted. The use of Petrov–Galerkin finite elements (FE) for these applications presented no instabilities and the vortices were close to what was obtained in the finite volume analysis. The machine learning algorithm, which used 6000 epochs, was moderately successful, as not all CST versions of the airfoils in the UIUC database were reliable. Optimization produced better airfoils than those present in the UIUC database; these advancements were in line with the known details that are the rule of thumb within the recreational aircraft community.