Nonintrusive Aerodynamic Shape Optimisation with a POD-DEIM Based Trust Region Method

Abstract

This work presents a strategy to build reduced-order models suitable for aerodynamic shape optimisation, resulting in a multifidelity optimisation framework. A reduced-order model (ROM) based on a discrete empirical interpolation (DEIM) method is employed in lieu of computational fluid dynamics (CFD) solvers for fast, nonlinear, aerodynamic modelling. The DEIM builds a set of interpolation points that allows it to reconstruct the flow fields from sets of basis obtained by proper orthogonal decomposition of a matrix of snapshots. The aerodynamic reduced-order model is completed by introducing a nonlinear mapping function between surface deformation and the DEIM interpolation points. The optimisation problem is managed by a trust region algorithm linking the multiple-fidelity solvers, with each subproblem solved using a gradient-based algorithm. The design space is initially restricted; as the optimisation trajectory evolves, new samples enrich the ROM. The proposed methodology is evaluated using a series of transonic viscous test cases based on wing configurations. Results show that for cases with a moderate number of design variables, the approach proposed is competitive with state-of-the-art gradient-based methods; in addition, the use of trust region methodology mitigates the likelihood of the optimiser converging to, shallower, local minima.

1. Introduction

The pressure to reduce design cycles, increase performance or further explore design spaces keeps on demanding more efficient analysis tools. Typical engineering problems, such as aerodynamic shape optimisation, require sophisticated analysis tools to solve parametric partial differential equations (PDEs). This is a challenging process due to the high computational cost associated with having to interrogate a large and complex model multiple times. The advent of adjoint methods to compute gradients at a cost independent of the number of design variables opened up the perspective of affordable design whilst retaining a high-fidelity model [1,2,3,4,5,6]. The solution of the CFD adjoint equations requires similar computational resources to those used to solve the original flow equations. Furthermore, adjoint methods require the modification of the original software code, which may or may not be available to users. Therefore, the use of high-fidelity models in design optimisation is feasible, but not necessarily trivial or affordable.

Several strategies are available to accelerate the time to solution required for a given analysis, in particular for problems such as optimisation or uncertainty quantification, where it is necessary to interrogate a large and complex model multiple times. One common approach is to build surrogate models and devise a multifidelity strategy to compute the quantity of interest efficiently. Surrogates are typically based on interpolation methods such as response surfaces or kriging [7,8,9,10] or neural networks [11,12]; it is also possible to resort to lower-fidelity models that remain representative of some relevant physical phenomenon to act as an effective surrogate [13,14]. The obvious penalty of such strategies is on how to guarantee that the method converges to a local optimum common with the high-fidelity model; in addition, users face challenges in associating models and respective parameterisations that use very different representations of the same problem.

A slightly different alternative is to build a ROM by examining the governing equations and performing some type of model reduction that is suitable for the physics of the problem. Here, the common approach is to assume that the parametric behaviour of the high-fidelity model can be approximated by a small number of modes or basis, typically obtained by methods such as proper orthogonal decomposition (POD) [15], balance POD [16,17] or proper generalised decomposition [18], among others—extensive reviews on the wider subject of model reduction can be found in [19,20,21,22,23,24]. The exploitation of ROMs for aerodynamic shape optimisation was investigated by LeGresley and Alonso [25] for inverse aerofoil design, where a POD ROM was used to approximate the gradient of the objective function by finite differencing. Investigators have also been successful in deriving efficient ROMs for shape optimisation problems described by elliptical PDEs, such as those found in marine or biological applications [26,27,28]. The problem of aerofoil shape optimisation subject to compressible flows was revisited by Zahr and Farhat; the authors assumed a monotonic nonlinear trust region optimisation method that updates the reduced basis as the optimisation progresses when the residual fails to reach a required threshold [29]. Carlberg et al. exploit a least-squares Petrov–Galerkin projection [30] to reduce the state equations and respective sensitivities; hence, each snapshot involved concatenating samples of the fluid-state variables and sensitivities with respect to the design variables. Yao et al. employed a similar strategy, but for gradient-based optimisation for aerofoils in compressible and transonic flows, reducing the time to solution by 30–50% with respect to an adjoint-based optimisation [31].

One particular challenge with the aforementioned ROM strategies is that they are intrusive, i.e., they require the manipulation of the original CFD solver, which may not be available or desirable. In reference [32], the authors investigate the ability of the DEIM [33] of reconstructing unsteady flowfields for the prediction of dynamic aeroelastic instabilities. Inspired by these findings, and to overcome the difficulties of intrusive methods whilst obtaining high computational efficiency, this paper proposes a nonintrusive strategy suitable for aerodynamic shape optimisation; this is achieved by a new nonlinear mapping between the surface DEIM points and flow quantities of interest, allowing the DEIM ROM to reconstruct surface flow fields for new shapes. The paper will also show how a trust region model management (TRMM) ROM is able to drive the aerodynamic shape to new local minima at a reduced cost with respect to adjoint-based optimisation, while treating the high-fidelity model as a black box.

The remainder of the paper provides details of the DEIM formulation and the TRMM algorithm. This is followed by the description and analysis of two transonic, viscous wing test cases. The paper is completed with a conclusion section.

2. Methodology

2.1. The General Optimisation Problem

A constrained aerodynamic shape optimisation can be represented by [31]

$$\begin{array}{cc}\hfill \underset{\mu \in \mathcal{D}}{\mathrm{minimise}}& \mathcal{F}\left(\mathbf{w}\right(\mu ),\mu )\hfill \\ \hfill \mathrm{subject}\mathrm{to}& \mathbf{R}\left(\mathbf{w}\right(\mu ),\mu )=0,\hfill \\ \hfill & {\mathcal{C}}_{\mathcal{E}}(\mathbf{w}\left(\mu \right),\mu )=0,\hfill \\ \hfill & {\mathcal{C}}_{\mathcal{I}}(\mathbf{w}\left(\mu \right),\mu )\le 0.\hfill \end{array}$$

(1)

In Equation (1), the objective function $\mathcal{F}$ depends on the flow variables $\mathbf{w}=[\rho ,\rho \mathbf{U},\rho E]$ and parameter set $\mu$. Typically, $\mathcal{F}$ denotes an aerodynamic coefficient such as $C_L$ or $C_D$; therefore, the nonintrusive ROM aims to reconstruct the surface pressure and wall shear stress fields to compute the required aerodynamic coefficients. This leads to a data-driven ROM, based on samples of the aforementioned quantities from the CFD full-order model (FOM), which are reduced by POD and a DEIM algorithm; this allows the flow field of every new shape produced by the optimiser to be reconstructed. The parameter set $\mu$ represents the design variables manipulated by the optimiser to control the aerodynamic shape. The different components required are outlined are described in the following sections.

2.2. Nonintrusive Reduced-Order Modelling

The DEIM aims to approximate $\mathbf{u}$ as a nonlinear function of $\mu$, where $\mathbf{u}=\mathbf{u}(\mathbf{w},\mu )$, by projecting it on the subspace spanned by the basis ${\mathbf{T}}_M$:

$$\mathbf{u}\left(\mu \right)\approx {\mathbf{T}}_M\mathbf{c}\left(\mu \right),$$

(2)

where $\mathbf{c}\left(\mu \right)$ represents a set of coefficients and the matrix ${\mathbf{T}}_M$ is obtained by POD of flow snapshots of the quantities represented by $\mathbf{u}$, which can be p, ${\tau }_w$, etc. The task now is how to select the M distinct rows or interpolation indices that allow a unique solution of $\mathbf{c}\left(\mu \right)$ to be obtained. This can be achieved by following the algorithm in [33]; the interpolation indices are coded into the matrix ℘, which can be obtained using Algorithm 1, and the result is the DEIM approximation:

$$\mathbf{u}\left(\mu \right)\approx {\mathbf{T}}_M{\left({\wp }^T{\mathbf{T}}_M\right)}^{-1}{\wp }^T\widehat{\mathbf{u}}\left(\mu \right).$$

(3)

Algorithm 1 DEIM method for interpolation index selection [33]

Input: Subspace ${\mathbf{T}}_M\in {\mathbb{R}}^{N\times n_r}$  Output: Interpolation indices $\wp ={[{\wp }_1,\dots ,{\wp }_{n_r}]}^T$ 1: $[i_{max},{\wp }_1]=[argmax{\left({\mathbf{T}}_{\mathbf{1}}\right)}^a,max\left(\left|{\mathbf{T}}_{\mathbf{1}}\right|\right)]$ 2: for $i=2$ to $n_r$ do 3:    $\mathbf{X}=\mathbf{T}(:,i)$ ; 4:    solve $\left({\wp }^T\mathbf{T}\right)\mathbf{c}\left(\mu \right)={\wp }^T\mathbf{X}$ for $\mathbf{c}\left(\mu \right)$ ; 5:    $\mathbf{r}=\mathbf{X}-\mathbf{T}\mathbf{c}\left(\mu \right)$ ; 6:    $[i_{max},{\wp }_i]=[argmax\left(\mathbf{r}\right),max\left(\left|\mathbf{r}\right|\right)]$ ; 7:    $\wp ={[\wp ,{\wp }_i]}^T$ ; 8: end for

For cases such as those presented here, the objective function is obtained through integrating the fluid forces over surfaces; hence, the size of the ROM can be further reduced by letting $\widehat{\mathbf{u}}\left(\mu \right)$ represent only the values of the surface pressure and shear stress corresponding to the nonzero indices in ℘, i.e., the M elements required to recover the full-order vector $\mathbf{u}\left(\mu \right)$.

Thus, for viscous problems, the pressure and each component of the wall shear stress are individually reconstructed using Equation (3) and a respective set of bases, ${\mathbf{T}}_M$.

2.3. Proper Orthogonal Decomposition

The basis ${\mathbf{T}}_M$ can be built by taking the POD modes of the flow field. This requires obtaining a set of snapshots and assembling the matrix $\mathbf{A}$, where each column represents the flow solution vector for one particular set of parameter values $\mu$:

$$\mathbf{A}=\left[{\mathbf{u}}_1,{\mathbf{u}}_2,\dots ,{\mathbf{u}}_{n_s}\right],$$

(4)

where $n_s$ is the number of samples. Since $n_s$ is much smaller than the dimension of $\mathbf{u}$, the POD modes can be obtained by performing singular value decomposition (SVD) to ${\mathbf{A}}^T\mathbf{A}$, where ${\mathbf{A}}^T\mathbf{A}\in {\mathbb{R}}^{n_s\times n_s}$:

$${\mathbf{A}}^T\mathbf{A}=\mathbf{T}\mathbf{S}{\mathbf{V}}^T,$$

(5)

where $\mathbf{S}$ is a diagonal matrix containing the eigenvalues and $\mathbf{T},\mathbf{V}$ are the left and right singular vectors, respectively. The POD bases, ${\mathbf{T}}_M\in {\mathbb{R}}^{n_s\times N}$, are obtained with proper truncation according to the energy rank given by the amplitude of the first M eigenvalues:

$${\mathbf{T}}_M=\mathbf{A}{\mathbf{V}}_M{\mathbf{S}}_M^{-1/2}$$

(6)

2.4. Interpolation Point Estimation

To employ the DEIM in lieu of a high-fidelity model, it is necessary to obtain $\widehat{\mathbf{u}}\left(\mu \right)$ for any shape produced by the optimiser. In this work, a nonlinear mapping between the design parameters $\mu$ and the flow values at the interpolation indices needs to be built. Note that $\mu$ and $\widehat{\mathbf{u}}$ are vectors of distinct sizes, which results in a multi−input/multi−output problem. Any suitable interpolation method can be used here; following from reference [32], an exact radial basis function (RBF) neural network using the MatLab function newrbe [34] was built using the respective input/output sample pairs $({\mu }_i,{\widehat{\mathbf{u}}}_i),i=1,\dots ,n_s$. This interpolation method ensures that the function values interpolated are exact at the sample points.

2.5. Trust Region Model Management

Obtaining a set of global bases is usually not achievable, or would require and excessive number of snapshots, rendering the ROM approach unattractive. An alternative is to start with a small number of samples and augment the snapshot matrix as the optimisation progresses, as demonstrated in [29,31]. Trust region methods are provable to be convergent and able to manage the exchanges between different fidelity models. At each major iteration k, an optimisation subproblem is defined on the trust region centred at ${\mu }_k$ and radius ${\Delta }_k$:

$$\begin{array}{cc}\hfill \underset{s\in {\mathfrak{B}}_k}{\mathrm{minimise}}& \widehat{\mathcal{F}}(\mathbf{u}\left({\mu }_k\right),\mu )\hfill \\ \hfill \mathrm{subject}\mathrm{to}& {\widehat{\mathcal{C}}}_{\mathcal{E}}(\mathbf{u}\left({\mu }_k\right),{\mu }_k+s)=0,\hfill \\ \hfill & {\widehat{\mathcal{C}}}_{\mathcal{I}}(\mathbf{u}\left({\mu }_k\right),{\mu }_k+s)\le 0,\hfill \\ \hfill & lb\le ({\mu }_k+s)\le ub,\hfill \\ \hfill & {\left|\right|s\left|\right|}_{\infty }\le {\Delta }_k,\hfill \end{array}$$

(7)

where the “$\widehat{}$” symbol indicates quantities computed using Equation (3); s is the optimisation step size; and ${\mathfrak{B}}_k=\{\mu \in {\mathbb{R}}^{n_p}:\left|\right|\mu -{\mu }_k\left|\right|\le {\Delta }_k\}$. The optimisation subproblem was solved using MatLab’s fmincon with the gradients obtained by finite differences using the DEIM ROM.

The effectiveness of the trust region step is evaluated by the ratio of the actual improvement over the improvement predicted by the ROM, $\rho$:

$${\rho }_k={\displaystyle \frac{\Delta \mathcal{F}}{\Delta \widehat{\mathcal{F}}}},$$

(8)

where $\Delta \mathcal{F}=\mathcal{F}\left({\mu }_k\right)-\mathcal{F}({\mu }_k+s)$ and $\Delta \widehat{\mathcal{F}}=\widehat{\mathcal{F}}\left({\mu }_k\right)-\widehat{\mathcal{F}}({\mu }_k+s)$. The step is accepted, except if $\Delta \mathcal{F}$ is nonpositive. In addition, for values of $\rho <{\delta }_1$ or $\rho >{\delta }_2$, the trust region radius is maintained if $\Delta \mathcal{F}>0$, and decreased if $\Delta \mathcal{F}\le 0$; if ${\delta }_1<\rho <{\delta }_2$, the trust region size is increased. In this work, the trust region is set up with ${\delta }_1=0.80$ and ${\delta }_2=1.20$; the trust region size, ${\Delta }_{k+1}$, is then reduced by a factor of $0.75$, maintained or increased by a factor $1.25$. Numerical experiments show these values to perform adequately over a range of test cases. The sequence of steps is outlined in Figure 1.

3. Results

The framework described was demonstrated using two wing test cases and different drag minimisation problems. The first test case involved the Onera M6 wing, test 2308, as described in [35]. This is a widely used test case by the CFD community [36]; hence, it acts as the first benchmark for the method proposed in this work. The second test case employed the NASA Common Research Model (CRM) wing [37]. The FOMs were built using the $S{U}^2$ software v7.07 [38]; in all cases, the compressible RANS equations were solved employing the Spallart–Almaras turbulence model; the full details of the flow equation descriptions and implementation, including the turbulence model, are given in reference [38].

3.1. Onera M6 Test Case

A baseline mesh for the optimisation study was selected following the grid convergence results shown in Figure 2. This study was performed at $M_{\infty }=0.84$, $Reynoldsn.=11.7\times {10}^6$ and $\alpha =3.06°$. The meshes included a fine mesh containing about 2.2 million elements, a medium mesh with approximately 315 thousand cells and a coarse mesh made up with approximately 140 thousand elements. The medium mesh, shown in Figure 3a, was selected for the baseline optimisation case as it gives adequate results without demanding excessive computational resources (mesh obtained from: https://github.com/su2code/TestCases/tree/master/optimization_rans/steady_oneram—accessed on 3 June 2022).

The wing is parameterised using 60 control points equally distributed over the top and lower surfaces of the wing. The surface mesh shape is controlled by the position of the control points and deformed by an RBF morphing method [39]. The surface mesh and associated control points are illustrated in Figure 3b.

To build the ROM, four quantities are required from the FOM: pressure, and three Cartesian components of the shear stress projected on the surface mesh. Each quantity is used to build a set of bases and interpolation point matrix, ℘; therefore, each function evaluation using the ROM requires reconstructing a matrix with $(N_{Surf}\times 4)$ elements. The reconstruction of the surface values at the sample points was exact; as an example, considering the test case conditions described above, Figure 4a–c shows two samples and resultant DEIM points (the DEIM points locations are highlighted by the black markers), the $C_p$ ROM values are plotted against the FOM results for each respective sample in Figure 4b–d, illustrating the exact reconstruction of the surface field.

The result of applying the ROM built from these two samples to a new wing shape is compared in Figure 5. As expected, the ROM predictions degraded. The level of degradation was mitigated by controlling the trust region radius. This shows that in principle, an optimisation problem can be initiated with any number of samples, including a single point; however, numerical experiments show that this is inefficient, as due to the low accuracy during the initial steps, 5−10 initial samples tend to produce faster convergences.

To facilitate benchmarking, the optimisation problem set up targets the wing at $M_{\infty }=0.84$, $Reynoldsn.=11.7\times {10}^6$ and $\alpha =3.06°$, and aims to minimise $C_D$ whilst maintaining the $C_L$ above or on the value corresponding to the baseline wing:

$$\begin{array}{cc}\hfill \underset{\mu \in \mathcal{D}}{\mathrm{minimise}}& C_D\hfill \\ \hfill \mathrm{subject}\mathrm{to}& C_L>0.25,\hfill \end{array}$$

(9)

where $\mu$ defines the control point positions that manipulate the wing shape.

The ROM at the start of the optimisation was built with ten snapshots. The problem was run on 4 Intel Core 5 CPUs, with both FOM and fmincon gradient calculations running in parallel. Note that the gradients were obtained by finite differences using the ROM; this operation scales linearly with the number of CPUs available. The optimisation convergence is shown in Figure 6, where both constraint and objective functions are normalised with respect to the baseline values for $C_D$ and $C_L$, which are 0.018 and 0.25, respectively. The right vertical axis shows the normalised value of the trust region radius, $\Delta$, with respect to its initial value; the convergence is plotted against the total number of FOM calls, i.e., the number of trust region updates plus the FOM calls used to build the initial ROM. The cost of the ROM calls is about $10\%$ of a single FOM evaluation, per trust region iteration; this cost scales with the number of design variables, but is also dependent on the convergence of each optimisation subproblem.

The solution terminates after 150 iterations and results in a drag reduction of about $36\%$. Note that the trust region radius increases significantly as new snapshots enrich the POD basis, meaning that a larger portion of the design space is explored by each optimisation subproblem. In the end, the large trust region radius resulted in increasingly thinner wing sections, producing inviable geometries, which led to the optimiser stopping after 150 iterations. The same optimisation problem was solved using a CAD-based adjoint method in reference [40], and resulted in $14\%$ performance improvement.

The surface flow field for the optimum obtained is compared against the original wing in Figure 7, showing how the shock on the upper surface is nearly eliminated everywhere.

The pressure coefficient distribution at two span stations is compared in Figure 8. Results suggest that the optimiser shifts the section loading further aft, resulting in a more uniform pressure distribution along the chord, akin to supercritical sections.

Following these results, the problem was modified by requiring the wing thickness to remain above $95\%$ of the original. The optimisation convergence history is shown in Figure 9. The result was an improvement of approximately $14\%$ in $C_D$, but unlike the previous case, the trust region radius declines rapidly and remains relatively small throughout the optimisation trajectory. In this case, as the wing remains relatively thick, there is more potential for shape perturbations to generate strong shocks that the ROM fails to capture adequately, leading to larger discrepancies between the ROM and FOM predictions. The upper surface pressure coefficient of the final shape is compared in Figure 10, showing a more limited and uniform pressure drop over most of the upper surface. It is clear that a mild shock remains on the wing upper surface. The surface pressure coefficient distributions of the inboard and outboard sections of the wing are given in Figure 11, showing a reduction on the leading edge suction peak, which mitigates the downstream shock, pushing it further downstream as the tip is reached.

3.2. CRM Wing Test Case

Recall that for the DEIM to reconstruct the surface variables, it requires the solution at the DEIM points, which are interpolated based on their respective positions. Hence, the next test case assesses the methodology using a wing twist parameterisation. The test case is based on the AIAA Aerodynamic Design Optimization Discussion Group case 4.1 (data available online at: https://drive.google.com/file/d/107HHAUioOdiyFG_e63_S6R3TwlwsVM47/view, retrieved on 27 February 2023), which makes use of the wing geometry from the Common Research Model wing body configuration proposed for the AIAA Fifth Drag Prediction Workshop; full details of the model are given by Vassberg [37]. The result is a wing geometry representative of modern wide-body commercial airliner, with a span of approximately 59 m and a mean aerodynamic chord of 7 m. The flow problem set up corresponds to a flight condition with $M_{\infty }=0.85$ at an angle of attack $\alpha =2.2°$, with a corresponding Reynolds number based on the mean aerodynamic chord of $5\times {10}^6$. This test case is investigated in depth by Lyu et al. [41] using an adjoint solver; the authors provide and assess a series of meshes, which are also used in this study (meshes obtained from https://data.mendeley.com/datasets/7jnyjdbvf9/1, accessed on 10 June 2020).

Table 1. SU2 grid convergence study and comparison with respect to reference [41].

Mesh Level	Mesh Size	${\mathit{C}}_{\mathit{D}}$—[41]	${\mathit{C}}_{\mathit{D}}$—SU2
L0	28,835,840	0.01997	0.02007
L1	3,604,480	0.02017	0.02034
L2	450,560	0.02111	0.02169

shows a grid convergence study comparison between the SU2 results obtained as part of this work and those from Lyu et al., showing very good agreement. Lyu et al. show that the coarse mesh captures the design space trends adequately; hence, only that mesh is used in this optimisation study.

The selected mesh contains approximately 450 thousand elements, of which about 11 thousand are on the surface. The control points are organised in seven sets, equally spread along the span of the wing, resulting in seven design variables to control the local wing twist at each station. The normalised surface model, including the sections used as design variables is shown in Figure 12.

The optimisation problem definition is identical to the previous test case:

$$\begin{array}{cc}\hfill \underset{\mu \in \mathcal{D}}{\mathrm{minimise}}& C_D\hfill \\ \hfill \mathrm{subject}\mathrm{to}& C_L>0.50,\hfill \end{array}$$

(10)

where $\mu$ represents the different wing twist angles at each span location and $C_L$ corresponds to the CRM design lift coefficient. The ROM was initialised with 10 random samples, generated using Latin hypercube sampling from a uniform distribution, with bounds $[-0.5°,0.5°]$.

The optimisation convergence history is given in Figure 13; results show a decrease in drag of just over $8\%$. Similarly to the last case, the trust region radius decreased rapidly, however, with new samples, the accuracy of the ROM increased, leading to further gains, and eventually, the improvements became less than a set threshold, and the optimiser stopped after 51 iterations.

An overview of the surface pressure coefficient is shown in Figure 14; results indicate a slight attenuation of the shock pattern on the upper surface, particularly in the inboard and wing tip regions. This is consistent with the the local effects of the optimisation shown in Figure 15.

Of course, wave drag is not the only contributor to the wing’s aerodynamic performance; Figure 16 shows a comparison between the span loadings of the original and optimised wings against an elliptical load distribution. Results indicate that the baseline wing should reduce the inboard lift and increase the lift produced outboard, which is what the optimised wing loading does. However, the new loading outboard exceeds the ideal case and the lift inboard mostly undershoots the elliptical loading, suggesting that any improvements to the induced drag are marginal at best. However, the overall effect is a drag reduction of about $8\%$.

4. Conclusions

A new surrogate-based optimisation strategy is proposed using a POD-DEIM ROM to replace full-order CFD simulations. The ROM uses the POD-DEIM to reconstruct surface flowfields of the quantities of interest; the reconstructed variable values at the DEIM points required for the field reconstructions were interpolated and obtained using an RBF neural network based on the position of the DEIM points on the shapes proposed by the optimiser. The resultant framework is able to approximate flowfields for transonic, viscous problems with a relative small number of POD bases, i.e., requiring a limited number of snapshots. The optimisation framework was completed with a TRMM algorithm that manages the optimisation and triggers new samples that are used to enrich the ROM as the optimisation progresses. The framework was tested against two transonic wing cases. The results obtained demonstrate the ROM-based surrogate model to be a practical, nonintrusive method for problems with a moderate number of design variables, without compromising the critical physics of the problem or the performance of the optimiser.