Next Article in Journal
Deep Reinforcement Learning-Augmented Spalart–Allmaras Turbulence Model: Application to a Turbulent Round Jet Flow
Previous Article in Journal
Study of the Geometry of an Oscillating Water Column Device with Five Chambers Coupled under Regular Waves through the Constructal Design Method
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Discrete and Continuous Adjoint-Based Aerostructural Wing Shape Optimization of a Business Jet

1
Parallel CFD & Optimization Unit, School of Mechanical Engineering, National Technical University of Athens, 15772 Athens, Greece
2
Dassault Aviation, 75008 Paris, France
*
Author to whom correspondence should be addressed.
Fluids 2024, 9(4), 87; https://doi.org/10.3390/fluids9040087
Submission received: 1 March 2024 / Revised: 1 April 2024 / Accepted: 3 April 2024 / Published: 5 April 2024

Abstract

:
This article presents single- and multi-disciplinary shape optimizations of a generic business jet wing at two transonic cruise flow conditions. The studies performed are based on two high-fidelity gradient-based optimization tools, assisted by the adjoint method (following both discrete and continuous approaches). Single discipline and coupled multi-disciplinary sensitivity derivatives computed from the two tools are compared and verified against finite differences. The importance of not making the frozen turbulence assumption in adjoint-based optimization is demonstrated. Then, a number of optimization runs, ranging from a pure aerodynamic with a rigid structure to an aerostructural one exploring the trade-offs between the involved disciplines, are presented and discussed. The middle-ground scenario of optimizing the wing with aerodynamic criteria and, then, performing an aerostructural trimming is also investigated.

1. Introduction

This paper is related to the design phase of modern aircraft wings, considering both aerodynamic and structural criteria. From the aerodynamics point of view, the wing design targets minimum drag and sufficient lift for the aircraft to fulfill its flight mission with minimum fuel consumption and environmental footprint. From the structure point of view, the wing should be able to withstand aerodynamic and inertial forces while being light enough to save fuel. Dealing with wings with higher aspect ratios, the required extra structural stiffness is typically achieved by adding more structural mass, which increases the aircraft’s weight. Multi-disciplinary optimization (MDO) methods help account for contradictory objectives, such as aerodynamic performance and weight. MDO is in accordance with FlightPath 2050 [1] and Destination 2050 [2], which set a long-term vision by proposing a pathway that combines new technologies.
The adjoint method, both in its discrete and continuous variants, is used herein to support the optimization studies. Adjoint methods are very popular in computational fluid dynamics (CFD)-based optimization problems since they compute the gradient of any performance metric at a cost that is independent of the number of design variables. In the discrete adjoint approach, the flow equations are first discretized and then differentiated, while in continuous, the adjoint equations are derived in the form of partial differential equations (PDEs), then discretized and numerically solved. An important aspect of the adjoint approach, strongly affecting the accuracy of the computed gradient and, thus, the optimization path itself, is the differentiation of the turbulence model in use.
In discrete adjoint, studies on the effect of the so-called frozen turbulence assumption in an aerodynamic shape optimization can be found in [3,4,5]. Regarding continuous adjoint, the group of authors from the Parallel CFD & Optimization unit of the NTUA (PCOpt/NTUA) was the first to develop the (continuous) adjoint to the Spalart-Allmaras turbulence model [6], for incompressible flows. This work was extended to compressible flows in [7] for the same turbulence model.
Concerning aerostructural optimization, some early attempts using low-fidelity tools can be found in [8,9]. During the last two decades, high-fidelity CFD and computational structural mechanics (CSM) tools have almost exclusively been used. For instance, in [10,11] Euler CFD codes were employed, while the use of CFD models involving the solution of the Reynolds-Averaged Navier-Stokes (RANS) equations into aerostructural shape optimization problems can be found in [12,13,14,15,16,17,18]. More accurate fluid flow models, such as large-eddy and direct numerical simulations, can be used too. However, a RANS equations’ solver coupled with a finite element structural analysis represents the state-of-the-art in aerostructural wing optimization [19]. In [12], the coupled aerostructural discrete adjoint for turbulent flows was applied to optimize the flight shape of two aircraft models by decreasing the drag at constant lift. Turbulence was modeled via the Spalart-Allmaras equation. Ref. [13] presented an aerostructural discrete adjoint method developments at ONERA, applied for drag and wing weight optimizations of the Airbus wing-body XRF1 configuration. It extended an aeroelastic adjoint (neglecting variations in the wing’s structural properties and sizing) towards aerostructural adjoint for flexible wings, by simultaneously optimizing the aerodynamic shape and the primary structure geometry. The fully coupled adjoint relied on a beam model, and material stresses were aggregated to reduce the high number of structural design constraints. In [15], the aerostructural wing gradient-based optimization with discrete adjoint was applied to commercial aircraft geometry while considering a powered engine, flutter, as well as overall aircraft design constraints. A (discrete) adjoint-based aerodynamic shape optimization that coupled the RANS solver with a commercial finite element solver and a thermodynamic engine cycle analysis tool was used for the gradient-based aeroelastic optimization of a full aircraft with powered engines in [14]. In [16], a framework, within SU2, for the gradient-based aerostructural optimization of wings, assisted by algorithmic differentiation, was tested on the ONERA M6 and NASA CRM wings. The results reconfirmed the importance of including aerostructural coupling in shape optimization. In [17], single- and multi-point aerostructural wing optimizations including a flow separation constraint at low-speed, high-lift conditions were performed. The flow separation constraint resulted in a substantially different wing design with better low-speed performance and only a slight decrease in cruise performance. The potential of using unconventional, tow-steered composites in constructing aircraft wings was studied in [18]. By coupling a RANS CFD solver with a finite element structural analysis code, it was shown that aircraft wings constructed by tow-steered composites can lead to 2.4% lower fuel consumption and 24% less weight compared to wings made of conventional composites. These benefits decreased when higher aspect ratio wings were concerned.
In this paper, Dassault Aviation and PCOpt/NTUA demonstrate a series of shape optimizations of a Generic Business Jet (GBJ) wing (with an aspect ratio equal to 9) using high-fidelity adjoint-based MDO tools considering two flight cruise conditions, namely, (FC1) M i n f = 0.82 , AoA = 2.5 and (FC2) M i n f = 0.80 , AoA = 2.0 both at the same flight altitude h = 41,000 ft, with M i n f being the farfield Mach number and AoA the angle of attack. In total, nine optimization runs were performed. The first eight comprise single- and multi-point wing shape optimization studies, with and without taking wing flexibility into account, incorporating only wing shape parameters. The purpose was to minimize the aircraft’s drag coefficient with constant lift in trimmed conditions. The ninth study dealt with aerostructural wing shape optimization for minimum fuel burn, involving both shape and structural parameters. The expression for fuel burn combines aerodynamic (lift-over-drag ratio) with structural metrics (weight), and constitutes a reasonable objective to investigate the trade-offs between the involved disciplines. In these studies, two different high-fidelity CFD codes and their adjoints were employed. Namely, AETHER [20] from Dassault Aviation, the adjoint of which is based on the discrete approach, and PUMA (v21.10) [21], from PCOpt/NTUA, which includes continuous adjoint. Comparisons between the results produced by the two CFD codes are also made and presented in this paper. Concerning structural analysis, the virtual performance solution (VPS) [22], CSM tool by ESI has been coupled with both CFD codes.

2. Aerodynamic Optimization Tools

The shape parameterization models, the flow analysis tools, and their (discrete/continuous) adjoints to support gradient-based aerodynamic optimization loops are briefly presented in this section. The objective function and constraints used are the lift ( C L ), drag ( C D ), and pitching moment ( C M ) coefficients, all of which are computed based on a constant reference area and length, as well as the fuel burn ( F W , Equation (13)).

2.1. Shape Parameterization and Grid Displacement

The parameterization of the GBJ wing shape was based on the Geometry ANd Inherent MEsh DEformation (GANIMEDE) tool, Figure 1 [23], which is the in-house tool of Dassault Aviation. GANIMEDE is a CAD-based geometric modeler that may handle shape parameters such as thickness, twist, and camber of wing sections.
In the applications presented in this paper, the number ( N D ) and type of design variables used in each optimization run are selected from a super-set of 11 design variables. b 1 to b 8 and b 10 , b 11 control the wing shape. The first eight control the shape of five wing sections (Figure 2); the twist angles of these sections are controlled by b 1 to b 4 , and their trailing edge camber by b 5 to b 8 . Sweep and non-dimensional thickness of the wing are controlled by b 10 and b 11 , respectively. The horizontal tail plane (HTP) angle is controlled by b 9 , enabling aircraft trimming.
The GANIMEDE tool is directly integrated within the optimization loop using the AETHER code. After generating a new surface mesh with GANIMEDE, a Laplacian-like model is used to adapt the CFD volume mesh. This uses the preconditioned conjugate gradient method, which is robust and efficient for solving sparse linear algebra problems.
The optimization loop, which uses the flow analysis and adjoint solver of PUMA, employs a properly trained surrogate model, replacing GANIMEDE; this was developed to directly provide the aircraft surface grid for each new value set of design variables. Firstly, 2 × N D (i.e., 22) surface grids corresponding to perturbations for all design variables by ± Δ b i were generated using GANIMEDE and used as training patterns. The bounds of all design variables used and, also, imposed as constraints in the optimization runs, are given in Table 1.
Let x D be the nodal coordinates of the surface grid for the datum (D) aircraft shape. Any change ± Δ b i in the i th design variable displaces the datum surface grid coordinates by Δ x b i . During the shape optimization, for any new value set of the N D design variables, the nodal coordinates x of the GBJ surface grid result from the superposition of all Δ x b i
x = x D + i = 1 N D Δ x b i
where Δ x b i is given by the quadratic polynomial
Δ x b i = a i Δ b i + c i Δ b i 2 , i = 1 , , N D
with x b i , a i , c i R 3 . Coefficients a i and c i result from the solution of three (one for each Cartesian coordinate) 2 × 2 systems of equations for each surface grid node. Illustrative examples of how (some of) the design variables may change the GBJ surface geometry are given in Figure 3. This surrogate model (SM), to be referred to as SM P , was analytically differentiated w.r.t. b i and incorporated into the PUMA adjoint software.
After generating a new surface grid during the PUMA-based optimization, a two-step radial basis function (RBF)-based grid deformation tool [24], is used to adapt the CFD volume grid to the newly defined boundary. In the first step (predictor), grid boundary nodes agglomerate to reduce the interpolation problem size; in this step, RBF with global support is used. Since the RBF interpolation applies to the boundary nodes too, the so-displaced boundaries do not perfectly match the known boundary displacements; so, the second step (corrector) corrects the position of the boundary nodes by performing local deformations using RBF kernels with local support. The whole process is accelerated by the sparse approximate inverse preconditioner and the fast multipole method [25]. Indicatively, in the studies performed using PUMA, the RBF-based grid deformation tool takes less than 10 core minutes (on a single Intel Xeon Silver 4114 CPU at 2.20 GHz) to morph ∼8.5 Mi CFD nodes given ∼200 K nodes on the aircraft surface.

2.2. The AETHER Flow Analysis and Adjoint Tool

The Dassault aviation code AETHER solves the RANS equations on unstructured grids with tetrahedral elements. It is based on a finite element method with an entropic formulation [26,27], stabilized with the streamline upwind Petrov-Galerkin (SUPG) method. Solving the RANS equations w.r.t. the entropy variables instead of the usual conservative variables has numerous advantages (such as dimensionally correct dot products, symmetric operators with positivity properties, and efficient preconditioning), providing a strong mathematical and numerical coherence. Among the many available turbulence models in AETHER, the Spalart-Allmaras model [28], is used herein. The use of an implicit scheme and the GMRES solver [29], make AETHER an efficient tool. The code has been successfully ported on many computer architectures and is fully vectorized and parallelized for shared or distributed memory machines using the MPI message passing library or OpenMP directives for many-core processors.
The adjoint of AETHER is based on the discrete approach. Thus, the objective or constraint function J, in discrete form, is augmented by the sum of the discrete residuals of the flow equations R n F , each multiplied by the corresponding vector of adjoint variables Ψ n F , n = 1 , , M F , where M F = 6 is the number of flow (i.e., RANS and Spalart–Allmaras) equations. Boldface symbols denote vectors with a size equal to the grid nodal count. By differentiating J aug w.r.t. the design parameters b i , i = 1 , , N D , we get
δ J aug δ b i = δ J δ b i + [ Ψ n F ] T δ R n F δ b i , i = 1 , , N D
where the Einstein summation applies for repeated indices. Since both J and R n F are functions of the flow variables W m F , m = 1 , , M F and the nodal coordinates x k F , k = 1 , , 3 , not directly depending on b , Equation (2) becomes
δ J aug δ b i = J W m F + [ Ψ n F ] T R n F W m F δ W m F δ b i + J x k F + [ Ψ n F ] T R n F x k F δ x k F δ b i
To avoid computing δ W m F δ b i , the discrete adjoint equations, namely [ R n F W m F ] T Ψ n F = [ J W m F ] T must be satisfied. By doing so, the second term on the r.h.s. of Equation (3) provides δ J aug δ b i . The discrete adjoint tool in AETHER [30], was programmed using the automatic differentiation tool TAPENADE [31]. This tool provides differentiated, w.r.t. the user-defined input variables, Fortran routines.

2.3. The PUMA Flow Analysis and Adjoint Tool

The PCOpt/NTUA code PUMA [21], solves the RANS equations on unstructured grids, which may consist of tetrahedra, pyramids, prisms, and/or hexahedra. It is based on the vertex-centered finite volume approach and a multi-stage Runge-Kutta scheme with residual smoothing. The inviscid fluxes are discretized using a central scheme with dissipation, which comprises a blend of second- and fourth-order solution differences. PUMA is programmed in the CUDA-C++ environment and runs on clusters of GPUs by decomposing the flow domain into overlapping subdomains; the grid decomposition process creates disjoint groups of grid edges, giving rise to subdomains with almost the same number of nodes. Computations and communications overlap, resulting in high parallel efficiency, while the MPI protocol or shared compute node memory are used for data transfer among GPUs on different or the same compute nodes, respectively. A distinguishing feature of PUMA is the use of mixed–precision arithmetic, thanks to which the GPU memory footprint and transactions are lower than with double–precision arithmetic, without jeopardizing the solver’s accuracy. Concerning turbulence closure, the Spalart–Allmaras model is used herein.
In the continuous adjoint used in the PUMA software, J is augmented by the integrals of the flow equations R n F multiplied by the adjoint variable fields Ψ n F , all of them in continuous form. The derivatives of J a u g w.r.t. b i give
δ J aug δ b i = δ J δ b i + Ω F Ψ n F δ R n F δ b i d Ω T 1 + Ω F Ψ G 2 x k 2 δ x F δ b i d Ω T 2
where n = 1 , , M F , k , = 1 , , 3 , and Ω F is the CFD domain. In Equation (4), over and above the meanflow and turbulence model equations, δ J a u g δ b i is additionally augmented by hypothetical grid deformation equations (as many as the spatial dimension of the problem) associated with additional grid–related adjoint fields ( Ψ G ). Practically, the grid is assumed to adapt to the displaced boundaries using Laplace equations (see T 2 in Equation (4)). This is just a convenient working hypothesis that allows the formulation of adjoint displacement PDEs and by no means implies that the grid should necessarily be displaced using this tool within the optimization loop; this is further discussed in Section 2.1. By applying the divergence theorem to T 1 and T 2 , a lengthy mathematical development (that is omitted here) leads to the following symbolic expressions for δ J δ b i , as well as the continuous equivalent for terms T 1 and T 2 ,
δ J δ b i = S A / C B J , m F δ W m F δ b i d S + S b D J , k δ x k F δ b i d S T 1 = Ω F C F , m F δ W m F δ b i d Ω + Ω F C F , k G δ x k F δ b i d Ω + S B F , m F δ W m F δ b i d S + S b D F , k δ x k F δ b i d S T 2 = Ω F C G , k G δ x k F δ b i d Ω + S b D G , k δ x k F δ b i d S
where S stands for the boundary of Ω (the aircraft’s wall, symmetry, and free–stream boundary). S b stands for the parameterized part of the CFD wall boundary (i.e., the aircraft’s wing and horizontal tail) and parts of the aircraft fuselage and vertical tail located close to the wing-fuselage and horizontal-vertical tail intersections. S A / C is the entire aircraft (A/C) boundary over which C L , C D , and C M are defined, m = 1 , , M F and k = 1 , , 3 . δ x k F δ b i is non-zero only over S b . Hereafter, the following notations are used: (a) C stands for expressions contributing to the field adjoint equations, (b) B contribute to the adjoint boundary conditions and (c) D contribute to the sensitivity derivatives. In these symbols, the subscript denotes the primal equation this term comes from, while the superscript denotes the adjoint PDE this term contributes to. Letter F corresponds to flow-related PDEs; G is for grid displacement; and later in this paper, S is for structure. The expressions of the multipliers B , C , and D , along with the detailed mathematical development of the continuous adjoint equations concerning only aerodynamic shape optimization, can be found in [21]. Combining Equations (4) and (5) and satisfying (a) the adjoint flow equations C F , m F = 0 , with boundary conditions B J , m F + B F , m F = 0 , to be solved for Ψ n F and (b) the adjoint grid displacement equations C F , k G + C G , k G = 0 , with Ψ l G = 0 on the farfield boundaries, to be solved for Ψ l G , the Sensitivity Derivatives ( S D s ) of J w.r.t. b i are given by the non-vanishing surface integrals, namely
δ J δ b i = S b D J , k + D F , k + D G , k δ x k F δ b i d S
Numerical stability issues of the adjoint solver, mainly due to the presence of a strong shock wave over the aircraft wing upper surface, are circumvented using the recursive projection method (RPM) [32]. The latter identifies the unstable modes and stabilizes the fixed-point iterative solver of the discretized adjoint PDEs by projecting the adjoint linear system onto the unstable subspace and solving it through Newton’s method.

2.4. Aerodynamic Cross Comparisons

A comparison of the above-mentioned flow and adjoint codes follows. This is based on the lift and drag coefficients of the GBJ as well as their sensitivity derivatives.
The polar curves computed by the two CFD codes, for M i n f = 0.80 and M i n f = 0.82 , both at the same flight altitude, are compared in Figure 4. The resulting lift-drag polars are quite close, and small differences are due to the use of different grids (16 Mi nodes tetrahedral grid for AETHER vs. 9 Mi nodes hybrid grid for PUMA) and methods (finite elements for AETHER vs. finite volumes for PUMA). All non-dimensional distances ( y + ) of the first nodes off the solid walls were below 1. A close-up view of the CFD grid used by PUMA is shown in Figure 5.
The S D s of C L , computed by the adjoint codes of AETHER and PUMA, are in agreement with those computed by finite differences (FDs, separately based on the two flow solvers and different grids), as illustrated in Figure 6. Despite major differences between AETHER and PUMA codes (different primal and adjoint equations’ discretization, use of discrete and continuous adjoints, or even use of different shape parameterization tools), their S D s are very close to each other, and so are FDs based on the two codes. Similar conclusions can be drawn for the S D s of C D and C M (not presented herein).
The frequently made frozen turbulence assumption was assessed in this case. Ignoring changes in the turbulent viscosity due to changes in the design variables results in considerable accuracy loss in the computed S D s , especially w.r.t. the wing section twist parameters ( b 1 to b 3 ). In support of the previous statement, the frozen turbulence derivative of C L w.r.t. b 3 is occasionally wrongly signed, Figure 6. This justifies our decision to include the adjoint to the Spalart–Allmaras turbulence model in all optimization runs.

3. Aerostructural Optimization Tools

Methods and tools in Section 2 are extended to account for aerostructural optimization. The structural analysis model, its coupling with the aerodynamic model, and its adjoint to support a gradient-based optimization are briefly presented below.

3.1. Structural Analysis Model

The structural analysis of the GBJ wing under aerodynamic and gravitational loading is carried out by means of the VPS software [22], based on an implicit solution scheme. The linear elasticity equations in discrete form are:
R k S = K k U S f k S f k G = 0 , with k , = 1 , , 3
where K stands for the stiffness matrix, f k S for the aerodynamic load applied on the structure nodes, and f k G for the gravitational load due to the structure weight as well as the lumped fuel load. U k S is the displacement array of the N S structure nodes, expressed as the difference between the flight ( x k S ) and the unloaded (jig) shape coordinates ( y S ), i.e., U k S = x k S y k S . For the aerostructural optimization studies performed for this paper, the wingbox layout of Figure 7 has been used. It consists of shell (for the wing’s skin) and beam (for the spar webs and ribs) elements.

3.2. Jig Shape Computation

Given a value-set for the design vector b , the parameterization tool (i.e., GANIMEDE or the surrogate model SM P ) generates an aircraft geometry x F * . To compute the corresponding wing jig shape y S * and its structural properties, the aerodynamic loads computed at this geometry and the constant wingbox layout of Figure 7 are taken into account. This computation is accomplished by solving a least-squares optimization problem that targets minimizing the difference between x F * and the CFD nodes displaced from y F * [33]. At each cycle, the structure nodes are re-positioned, and an inner optimization targeting the minimum structure weight under stiffness, plasticity, and buckling constraints is performed. The critical load cases that the inner optimization takes into account are defined according to the CS-25 certification norm [34].
The design variables b 10 and b 11 , controlling the wing’s sweep angle and thickness, strongly affect the structural model properties (and the wing’s jig shape). On the other hand, changes in b 1 to b 9 within the ranges of Table 1 do not significantly affect the structural properties, and their influence is practically negligible. Thus, the performed aerostructural optimization studies, accounting only for b 1 to b 9 , have a frozen structural model. On the other hand, in the aerostructural optimization using PUMA, which accounts for changes in b 10 and b 11 too, the structural model properties and wing’s jig shape were re-computed at the beginning of each optimization cycle. To avoid repeating the computationally expensive iterative procedure of updating the structural properties and jig shape, a polynomial surrogate model ( SM J ) was created and used during the optimization. For this, a design of experiments (DoE) based on the upper and lower bounds of b 10 and b 11 was performed. The surrogate models SM P and SM J are depicted in a schematic view in Figure 8. The aircraft’s geometry ( x F * ), produced by the SM P for the current value–set of b is used to initialize the aerostructural loop within each optimization cycle. This corresponds to a displacement from the aircraft’s jig shape to be denoted as U S * .
Apart from the reduction in the optimization wall-clock time, another benefit of using SM J is that its differentiation w.r.t. the design variables b 10 and b 11 is straightforward. In addition, using SM J , the structural design parameters are hidden from the main optimizer, which “sees” only the shape design parameters b .

3.3. Coupled Flow and Structural Analysis Tool

The deformations computed at the CSM nodes are interpolated to the boundary CFD nodes using the RBF model
U k F , i = j N S α k , j ϕ r F i , S j 2 + β k , 0 + β k , 1 x F , i + β k , 2 y F , i + β k , 3 z F , i , i = 1 , , N F
applied for any CFD node i on the wing. Here ϕ stands for Wendland’s W33 function [35], and r F i , S j = | | x F , i x S , j | | . The unknown coefficients β k , 0 , β k , 1 , β k , 2 , β k , 3 and α k , j for k = 1 , , 3 and j = 1 , , N S should satisfy the constraints
j N S α k , j = j N S α k , j x S , j = j N S α k , j y S , j = j N S α k , j z S , j = 0
and reproduce the structural deformations U k S
U k S , i = j N S α k , j ϕ r S i , S j 2 + β k , 0 + β k , 1 x S , i + β k , 2 y S , i + β k , 3 z S , i , i = 1 , , N S
Equations (9) and (10) can be written in matrix form as
A S S w k = v k
with
A S S = 0 0 0 0 1 1 1 0 0 0 0 x S , 1 x S , 2 x S , N S 0 0 0 0 y S , 1 y S , 2 y S , N S 0 0 0 0 z S , 1 z S , 2 z S , N S 1 x S , 1 y S , 1 z S , 1 ϕ S 1 , S 1 ϕ S 1 , S 2 ϕ S 1 , S N S 1 x S , 2 y S , 2 z S , 2 ϕ S 2 , S 1 ϕ S 2 , S 2 ϕ S 2 , S N S 1 x S , N S y S , N S z S , N S ϕ S N S , S 1 ϕ S N S , S 2 ϕ S N S , S N S , w k = β α k , v k = 0 U k S and ϕ S i , S j = ϕ r S i , S j 2
Equation (11) represents three systems to be solved for the three Cartesian directions to get coefficients β k and α k . To link the interpolated deformation fields over the boundary CFD nodes directly with the deformations at the CSM nodes, the expression for computing U k F is written in matrix form as U k F = A F S w k with
A F S = 1 x F , 1 y F , 1 z F , 1 ϕ F 1 , S 1 ϕ F 1 , S 2 ϕ F 1 , S N S 1 x F , 2 y F , 2 z F , 2 ϕ F 2 , S 1 ϕ F 2 , S 2 ϕ F 2 , S N S 1 x F , N F y F , N F z F , N F ϕ F N F , S 1 ϕ F N F , S 2 ϕ F N F , S N S
and, then, since w k = A S S 1 v k , U k F can be computed as
U k F = A F S A S S 1 0 1 H U k S
giving rise to the coupling matrix H , where 0 , 1 are arrays of zeros and ones with size 4 and N S , respectively. In order to reduce the size of H , deformations at the skin element centers are interpolated rather than those at all the CSM nodes.
The CFD and CSM analysis solvers are coupled via a fixed-point iteration scheme with adaptive relaxation, according to the steps shown in Figure 9. The CFD solver computes the aerodynamic loads on the aircraft’s wing ( f F ). These are mapped onto the CSM boundary nodes ( f S ) using the transpose of the coupling matrix H as f k S = H T f k F . This ensures conservation of force, moment, and virtual work [36]. The CSM solver computes the structural node displacements ( U S ), which are adaptively relaxed ( U ˜ S ) using Aitken’s relaxation formula [37], and, then, mapped onto the CFD surface grid ( U F ) as U F = H U k S . The CFD volume grid is adapted to the new wing shape, and a new aerostructural cycle begins. In Figure 9, the “CFD” box includes either the AETHER or the PUMA flow solver and the CFD grid displacement tool. The termination criteria are related to the convergence of the C L , C D , and C M values. Upon convergence, x F corresponds to the flight shape.

3.4. Coupled Adjoint Flow and Structural Solver

In the aerostructural optimization runs, the performance metrics may combine information from both the CFD and CSM disciplines. Such a metric, to be used as an objective function, is the fuel burn ( F W ) during the flight; this is computed by the Breguet-Leduc formula [18], as
F W = Z F W · e x p R · g · T S F C V · C L C D 1
where R, g, T S F C , and V stand for the aircraft’s range, the gravitational acceleration, the aircraft’s thrust specific fuel consumption, and its velocity magnitude, respectively; these quantities remain constant during optimization. The zero fuel weight ( Z F W ), provided by the SM J , depends on b 10 or b 11 .
To compute the S D s of the objective or constraint functions ( C L , C D , C M , and F W depending on the optimization run), while accounting for wing flexibility, the mathematical developments of the adjoint method of Section 2.2 and Section 2.3, i.e., of AETHER and PUMA respectively, were properly extended.
With regard to the discrete adjoint of AETHER, Equation (2) now includes an additional term containing variations in the discrete CSM residuals ( δ R k S δ b i ) multiplied by the CSM adjoint variables Ψ k S , i.e.,
δ J aug δ b i = δ J δ b i + [ Ψ n F ] T δ R n F δ b i + [ Ψ k S ] T δ R k S δ b i , n = 1 , , M F and k = 1 , , 3
Compared to the adjoint formulation (for a rigid wing) of Section 2, herein J (such as J = F W ) may directly be affected by the design variables. So,
δ J δ b i = J b i + J W m F δ W m F δ b i + J x k F δ x k F δ b i
where J b i is computed by the differentiated SM J of Section 3.1. By differentiating the expression of the CSM residuals, the last term of Equation (14) becomes
[ Ψ k S ] T δ R k S δ b i = [ Ψ k S ] T δ K k l δ b i U l S + K k l δ U l S δ b i δ f k G δ b i [ Ψ k S ] T δ f k S δ b i
Let us refer to the multiplier of the derivative of f k S w.r.t. b i as the adjoint loads computed on the CSM surface nodes, f k S , a d j = Ψ k S . Introducing the derivatives of J, R n F , and R k S w.r.t. b i in Equation (14) and taking into account that the aerodynamic loads f k F (mapped from the CSM onto the CFD surface nodes using H T ) are functions of the flow variables W n F and the flight shape geometry x k F , one gets
δ J aug δ b i = J b i + [ Ψ k S ] T δ K k l δ b i U l S + K k l δ U l S δ b i δ H T δ b i f k F δ f k G δ b i + J W m F + [ Ψ n F ] T R n F W m F [ f k S , a d j ] T H T f k F W m F δ W m F δ b i + J x l F + [ Ψ n F ] T R n F x l F [ f k S , a d j ] T H T f k F x l F δ x l F δ b i
According to Equation (15), f k S , a d j is first mapped onto the CFD surface nodes (using matrix H ) before being multiplied with the derivatives of f k F ; mapping is performed using the expression f k F , a d j = H f k S , a d j . The flight shape geometry x k F is linked with the jig shape y k F using x k F = y k F + H U k S . By differentiating the expression of x k F , the last term of Equation (15) takes the form
J x l F + [ Ψ n F ] T R n F x l F [ f k S , a d j ] T H T f k F x l F δ x l F δ b i = J x l F + [ Ψ n F ] T R n F x l F [ f k S , a d j ] T H T f k F x l F δ y k F δ b i + δ H δ b i U k S + H δ U k S δ b i
Let us refer to the multiplier of the derivative of U k S w.r.t. b i as the adjoint displacement vector; this is computed over the CFD surface nodes and mapped onto the CSM ones using H T as
U k F , a d j = J x k F T + R n F x k F T Ψ n F f m F x k F T f m F , a d j and U k S , a d j = H T U k F , a d j
So, Equation (15) becomes
δ J aug δ b i = J b i + [ Ψ k S ] T δ K k l δ b i U l S δ H T δ b i f k F δ f k G δ b i + [ U k a d j , F ] T δ y k F δ b i + δ H δ b i U k S + J W m F + [ Ψ n F ] T R n F W m F [ f k F , a d j ] T f k F W m F δ W m F δ b i + [ Ψ k S ] T K k l + [ U l S , a d j ] T δ U l S δ b i
To avoid the computation of δ W m F δ b i and δ U l S δ b i , the adjoint CFD and CSM equations
R n F W m F T Ψ n F f k F W m F T f k F , a d j = J W m F T
K k l Ψ l S + U k S , a d j = 0
should be satisfied. In contrast to the primal aerostructural loop, Figure 9, the adjoint displacements ( U k F , a d j ) are transferred from the adjoint CFD tool to the adjoint CSM one, and the adjoint loads ( f k S , a d j ) the other way around; this is done using H T or H , respectively. The remaining terms give the S D s , namely
δ J δ b i = J b i + [ Ψ k S ] T δ K k l δ b i U l S δ H T δ b i f k F δ f k G δ b i + [ U k F , a d j ] T δ y k F δ b i + δ H δ b i U k S
where J b i , δ K k l δ b i , δ f k G δ b i , δ H δ b i , and δ y k F δ b i are computed by the differentiated structural and parameterization tools.
The development of the continuous adjoint variant of PUMA is similar. Equation (4) is further expanded by accommodating a discrete part containing the structural adjoint variables and the derivatives of the residuals of the discrete CSM equations w.r.t. b i , leading to a hybrid continuous-discrete FSI adjoint. The development starts by
δ J aug δ b i = δ J δ b i + Ω F Ψ n F δ R n F δ b i d Ω + Ω F Ψ m G 2 x k F 2 δ x m F δ b i d Ω + [ Ψ k S ] T δ R k S δ b i
Introducing the expression of the structural model residual after a lengthy mathematical development, Equation (19) takes the form
δ J a u g δ b i = S A / C B J , m F δ W m F δ b i d S + S b D J , k δ x k F δ b i d S + Ω F C F , m F δ W m F δ b i d Ω + Ω F C F , k G δ x k F δ b i d Ω + S B F , m F δ W m F δ b i d S + S b D F , k δ x k F δ b i d S + Ω F C G , k G δ x k F δ b i d Ω + S b D G , k δ x k F δ b i d S + [ Ψ k S ] T K k l δ U l S δ b i + [ Ψ k S ] T δ K k l δ b i U l S δ f k G δ b i [ Ψ k S ] T H T δ f k F δ b i
where multipliers B , C , and D are the same as in Section 2.3. It should be noted here that the CFD solid wall surface splits not only into S A / C and S b but also S F S I , denoting the fluid–structure interface (i.e., the wing).
Following the terminology used in discrete adjoint, the adjoint loads computed at the CSM surface nodes, namely f k S , a d j = Ψ k S , are mapped onto the CFD surface grid nodes over S F S I W as f k F , a d j = H f k S , a d j . Collecting, then, the multipliers of δ x k F δ b i in the surface integrals over S F S I , one gets the adjoint displacement field
U k F , a d j = D J , k + D F , k + D G , k + f m F , a d j f m F x k
where the last term comes from the differentiation of the aerodynamic forces. U F , a d j is mapped onto the CSM surface nodes as U k S , a d j = H T U k F , a d j . Satisfying (a) the adjoint CFD equations C F , m F = 0 , subject to B J , m F + B F , m F + f k F , a d j f k F W m = 0 , with f k F , a d j being non-zero only over S F S I , (b) the adjoint grid displacement equation C F , k G + C G , k G = 0 , subject to Ψ l G = 0 over the farfield boundaries and a zero–Neumann condition over S A / C , and (c) the adjoint CSM equations K k l Ψ l S + U k S , a d j = 0 , the derivatives of J w.r.t. b i are computed as
δ J δ b i = J b i + [ Ψ k S ] T δ K k l δ b i U l S δ H T δ b i f k F δ f k G δ b i + [ U k F , a d j ] T δ H δ b i U k S + S b S F S I D J , k + D F , k + D G , k δ x k F * δ b i d S + S F S I x k F , a d j δ y k F δ b i d S
Figure 10 illustrates the coupled adjoint aerostructural loop. The adjoint CFD tool computes the so-called adjoint displacements ( U F , a d j ). These are mapped onto the CSM boundary nodes and, then, scaled as in the primal problem. Then, the adjoint CSM tool computes the adjoint load ( f S , a d j ) field, which is mapped onto the CFD boundary nodes, and a new adjoint aerostructural cycle starts. In the adjoint aerostructural loop, the convergence criterion accounts for the norm of the quantities exchanged between the disciplines.
A comparison of the adjoint-based S D s of the aircraft C L with FDs at FC1 is presented in Figure 11, which corroborates the above-mentioned hybrid FSI adjoint.

4. Applications

Initially, a rigid wing was assumed, and shape optimization runs were performed at both flight conditions (FC1 and FC2) by using different sub-sets of design variables. Then, aerostructural optimizations follow. Optimization runs are abbreviated as OptRxx (xx being the case number). All of them start from a trimmed GBJ configuration, in which the HTP rotation ( b 9 ) is adapted to yield zero C M ; this will be referred to as the baseline configuration (superscript B) at the corresponding flight condition and for rigid or flexible wing. In all runs, by either code, the SLSQP [38], which is a sequential quadratic programming-based algorithm for non-linearly constrained gradient-based optimization problems, was used.

4.1. Aerodynamic Shape Optimization with Rigid Wing Structure

The first runs focused on the aerodynamic optimization of the GBJ by assuming a rigid wing structure and considering only the first 9 design variables ( N D = 9 ). Optimization runs at FC1 (OptR1) and FC2 (OptR2) were performed, the first one exclusively using PUMA and the second using both codes. In both cases, the target was to minimize the GBJ drag coefficient and maintain the lift and pitching moment values within a certain margin w.r.t. the baseline configuration, i.e.,
min . C D s . t . | C L C L B | < 10 4 | C M | < 10 4
Figure 12 presents the convergence history of the two runs using the PUMA code and its adjoint. After 29 (OptR1) or 24 (OptR2) optimization cycles, a reduction in C D by ∼4% and ∼2%, respectively, was achieved while satisfying the constraints. Since the optimizations started with trimmed configurations (baseline), the C M constraint was satisfied from the first cycle. Each optimization cycle involves one flow and three adjoint (for the three coefficients of Equation (22)) problem solutions. Given that the cost to solve each adjoint system of equations is similar to that of the flow equations, the cost per optimization cycle is considered to be 4 equivalent flow solutions (EFS). The aforementioned reductions in C D can be interpreted (Equation (13)) as ∼4.8% and ∼2.6% reduction in fuel burn, respectively.
Figure 13 and Figure 14 compare the pressure coefficient fields (left) and the fields of local contributions to C D (right) on the GBJ surface of the baseline and the optimized configurations of OptR1&2. A strong shock wave across the wing suction side can be observed in both flight conditions, for both the baseline and optimized configurations. It can also be seen that the reduction in the objective function value is (mainly) attributed to drag reduction in the area close to the wing-fuselage junction over the suction side of the wing.
Figure 15 compares the convergence histories of AETHER and PUMA in OptR2. Both codes reduced C D by ∼2% while satisfying the constraints. Despite some small differences in the optimal design variable vectors, Table 2, as computed by the two methods, the trend is the same, suggesting a negative twist for section 0 (that close to the wing-fuselage intersection) and a positive for the rest sections. A slightly negative twist of the horizontal tail is also suggested by both CFD codes for trimming the aircraft. The aforementioned small differences between the optimal design vectors computed by AETHER and PUMA are attributed to the different geometric modelers the two optimization loops use (AETHER uses GANIMEDE, replaced by a surrogate model in the loop of PUMA).
In order to assess the role of each type of design variable in the achieved improvements, two extra optimization runs (OptR3 and OptR4) were carried out at FC2, both using AETHER. In OptR3, the design variables related to the wing twist angles ( b 1 to b 4 ) remained fixed and equal to their baseline values, whereas b 5 to b 9 were allowed to change within the bounds of Table 1. In OptR4, over and above to what was done in OptR3, the AoA was allowed to change within the range [ 1.9 , 2.1 ] ; changing the AoA is almost equivalent to a uniform span-wise change of the wing twist. The achieved reductions in C D are 0.27 % in OptR3 ( N D = 5 ) and 0.30 % in OptR4 ( N D = 6 ), which are both much lower compared to the reduction achieved in OptR2. In addition, OptR4 suggests that AoA = 1.95 , which is very close to the constant AoA = 2 value of OptR2, without though a noticeable gain in C D . It is interesting to note that, despite the narrow range of AoA, the optimal value reached neither the upper nor the lower bound. This justifies the selection of such a narrow band. OptR1 to OptR4 reveal that, for this aircraft at FC2, modifying the span-wise wing twist distribution (especially the twist of the wing-fuselage intersection) is the main mechanism to reduce drag, while respecting the constraints; in other words, twist distribution is much more important than changes in the trailing edge camber.

4.2. Aerostructural Shape Optimization

The next shape optimization runs take into account the wing structure flexibility. Initially, the wing structural properties (and the corresponding finite element model) are kept constant. Then, an aerostructural optimization including structural parameters was carried out and presented next.

4.2.1. Single-Point Optimization with Fixed Structural Model

OptR1 (at FC1) and OptR2 (at FC2) are herein revisited (OptR5&6) by considering the wing structure flexibility. The structural model properties and wing’s jig shape were computed once (following the iterative procedure of Section 3.1) and kept constant during the optimization runs with the first N D = 9 design variables. The objective and constraint functions, as well as the design variable ranges, are those presented in Section 4.1, Equation (22) and Table 1.
Figure 16 presents the convergence histories of OptR5&6 using PUMA. A reduction in C D by ∼4.0% and ∼1.8% was achieved after 25 and 22 cycles, respectively, while also satisfying the imposed constraints. These reductions, using expression (13), lead to ∼3.9% and ∼2.3% reduction in F W . Recall that similar reductions, after more or less the same number of cycles, obtained from OptR1&2, see Section 4.1. However, the post-hoc aerostructural analysis of the outcomes of OptR1&2 results to solutions violating the constraints. For instance, the re-evaluated solutions have | C M | 3.5 · 10 3 , which is way higher than the used threshold of 10 4 .
The flight shapes of the optimal solutions obtained from OptR5&6 result from the combined wing deformation due to (a) the design variable changes and (b) flexibility. The former stands for the deviations between x m F * , 0 and x m F * , where x m F * , 0 is the shape obtained from the SM J using the baseline design variable set and x m F * is the corresponding shape for the optimized configuration. The latter stands for the deviations between x m F * and x m F , with x m F being the flight shape for the optimized configuration. Figure 17 shows the wing deformation fields for OptR5&6. Deformations due to the design variables changes and flexibility are of the same order, but their maximum values appear in different areas over the wing. As expected, deformations due to the wing flexibility systematically get higher values towards the tip. On the other hand, changes in the design variables mostly affect the wing shape close to the fuselage, and through this, the C D values, as in OptR1&2, too.
OptR6 was also computed using AETHER. Figure 18 compares the convergence histories of OptR2 and OptR6 using AETHER. Both runs, i.e., considering rigid and flexible wing, yielded a reduction by ∼2% in C D . As with PUMA, the aerodynamically optimized configuration (that of OptR2 with AETHER) was re-evaluated using the aerostructural analysis model and failed to satisfy the constraints. To obtain a feasible solution, an aerostructural trimming of the OptR2 optimal solution, in the form of an aerostructural optimization by varying b 9 and AoA, targeting zero C M by also constraining C L , was carried out. The trimming process was completed after 3 cycles and led to a feasible solution with a ∼1.7% reduction in C D w.r.t. the baseline configuration of OptR6 (blue filled triangle of Figure 18). For these flow conditions, the gains from the aerostructural optimization and the aerodynamic one with trimming are almost the same; nevertheless, the aerostructural optimization yielded a slightly better solution (an additional gain of 0.3 % in the C D value). This came with a higher computational cost, as it required the additional solution of the CSM problem and its adjoint.

4.2.2. Two-Point Optimization with Fixed Structural Model

Wing flexibility effects were also investigated by considering a two-point optimization scenario. This used N D = 10 design variables in total, namely the 4 twist angles, the 4 camber parameters, and the 2 HTP rotation angles (one for each point). The objective was to minimize the weighted sum of C D values at FC1 and FC2, namely, J = w 1 C D 1 + w 2 C D 2 , with ( w 1 , w 2 ) = ( 0.75 , 0.25 ) , under constraints on C L and C M (imposed at each flow condition). Two runs (OptR7&8) dealing with the aforementioned target and constraints were performed using AETHER. Wing flexibility was ignored in OptR7, while OptR8 took wing flexibility into account by also keeping the wing’s structural model fixed during the optimization. As in the single-point studies, OptR7&8 resulted in similar reductions in the objective function after, more or less, the same number of optimization cycles. Trimming the optimized solution obtained from OptR7 was, again, necessary in order to obtain a feasible solution. However, in this case, trimming resulted in a solution with a reduction in the objective function value very close to that obtained from OptR8. Thus, for the above-mentioned selection of weights linking the two flight conditions, there is almost no gain from performing an aerostructural optimization instead of an aerodynamic one followed by aerostructural trimming.

4.2.3. Optimization with Varying Structural Model

PUMA was also used for an aerostructural shape optimization (OptR9) with a varying structural model, which was updated during the optimization using the SM J . This way, the structural design parameters were hidden from the main optimizer, which handled the N D = 11 shape design variables of Table 1. OptR9 target was to minimize F W (Equation (13)) at FC2, with constraints C L C L B 10 4 , | C M | < 10 4 .
The evolution of the design variables is presented in Figure 19, while the optimized and baseline wing shapes are in Figure 20. Figure 21 shows the computed pressure coefficient fields on the baseline and the optimized GBJ surface. It can be seen that the shock wave on the suction side became stronger, increasing C D . At the same time, on the pressure side, the pressure close to the trailing edge increased. The latter, combined with the pressure reduction on the suction side, resulted in a higher C L .
Figure 22 presents the convergence history of the aerostructural optimization. After ∼30 cycles, the fuel burn was reduced by more than 10 % while satisfying the C L and C M constraints. However, it should be noted that this came at the expense of a drag increase and a significantly higher lift (leading to an overall increase of 8 % in C L C D ). Aside from the fact that performing the mission at a higher drag is incompatible with fuel reduction, the underlying assumptions upon which the Breguet-Leduc formula is based, i.e., quasi-level flight and quasi-steady flight to yield a realistic approximation for the fuel burn, are not met to an important extent after a single optimization run. A more complex iterative process on top of the existing loop would be required to converge to a feasible solution, but this lies beyond the scope of this study.

5. Discussion and Conclusions

Nine optimization runs (OptR1-9) related to the aerodynamic and aerostructural shape optimization of a generic business jet at two flight conditions were presented. Studies were based on two different codes (a finite-element and a finite-volume one) with different adjoint formulations (continuous and discrete adjoint, respectively). Another difference between the two CFD tools, as used within the optimization loops, is that the first code (AETHER) relied upon a built-in CAD-based parameterization tool (GANIMEDE), whereas the second code (PUMA) implemented a surrogate model built using a set of shapes generated by GANIMEDE for distinct variable combinations. For the structural analysis and optimization, the VPS tool was coupled with both CFD codes. The paper presented the mathematical development of the (aerodynamic and aerostructural) adjoint method for both codes. A first interesting outcome is the absolutely satisfactory matching of the gradients computed by the two adjoint codes, as well as finite differences, despite the aforesaid differences. It was also demonstrated that the frozen turbulence assumption, in this particular application, computes wrong, sometimes even wrongly signed, sensitivity derivatives.
Three interesting findings from these studies are summarized below:
(a) The idea of replacing the (computationally expensive) aerostructural shape optimization (OptR7&8) with a pure aerodynamic (OptR1&2), followed by a post-hoc aerostructural evaluation of the optimized shapes (less expensive overall), was tried at both flight conditions. Though, in either optimization, the gain in C D was practically the same, the less expensive approach failed to satisfy the imposed constraints on C L and C M when wing flexibility was post-hoc accounted for. To meet these two constraints, a final trimming was necessary; the latter stands for an aerostructural optimization with fewer design variables (the aircraft’s HTP rotation and AoA, only). This very last step led to a feasible solution but showed less gain in C D than the expensive approach. This finding was reconfirmed using either PUMA or AETHER and a fixed structural model.
(b) A similar treatment was considered in a case in which the objective function is the weighted sum of the GBJ performances at the two flight conditions. In this case, the optimized solution obtained by the less expensive tool (the solution of OptR7), post-hoc evaluated by the expensive tool, also violated the constraints. The final trimming led to a feasible solution that performed as well as the one obtained by the expensive approach (OptR8).
(c) The last optimization run (OptR9), based on an aerostructural tool with varying structural models (and its adjoint), was quite expensive but opened the way for a high-fidelity trade-off that could be used in a broader MDO context during the design process.
These findings highlight the importance of taking wing flexibility into account in the GBJ design, either by running an aerostructural optimization loop or, at least, by aerostructurally trimming the outcome of a pure aerodynamic optimization. Depending on the case, the latter may reduce the overall cost even though, in some cases, sub-optimal solutions may result.

Author Contributions

Conceptualization, G.R. and K.G.; methodology, K.T., X.T., G.R., S.K. and L.M.; software, K.T., X.T., V.A., G.R., S.K., L.M. and S.J.; writing—original draft preparation, X.T., V.A. and K.G.; writing—review and editing, K.T., X.T., V.A., K.G., G.R., S.J., L.M. and S.K.; supervision, G.R. and K.G. All authors have read and agreed to the published version of the manuscript.

Funding

This work was part of the European Union’s Horizon 2020 research and innovation program under grant agreement No 769025 (MADELEINE).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The datasets presented in this article are not available due to commercial restrictions.

Acknowledgments

The authors would like to thank Gabriel Fougeron (ESI) and Olivier Amoignon for their valuable contributions in the development of the wing structure finite element model and the RBF-based non-matching interface tool.

Conflicts of Interest

Authors Gilbert Rogé, Sarah Julisson, Ludovic Martin and Steven Kleinveld were employed by the company Dassault Aviation. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
CFDComputational Fluid Dynamics
CSMComputational Structural Mechanics
FDsFinite Differences
GBJGeneric Business Jet
HTPHorizontal Tail Plane
MDOMulti-Disciplinary Optimization
RANSReynolds-Averaged Navier-Stokes
RBFRadial Basis Function
SMSurrogate Model
SDsSensitivity Derivatives

References

  1. European Commission. Flightpath 2050—Europe’s Vision for Aviation; European Commission: Brusssels, Belgium, 2011. [Google Scholar]
  2. van der Sman, E.; Peerlings, B.; Kos, J.; Lieshout, R.; Boonekamp, T. Destination 2050—A route to net zero European aviation. NLR-CR-2020-510, 2020. Available online: https://www.destination2050.eu/wp-content/uploads/2021/03/Destination2050_Report.pdf (accessed on 2 February 2024).
  3. Dwight, R.; Brezillon, J. Effect of approximations of the discrete adjoint on gradient-based optimization. AIAA J. 2006, 44, 3022–3031. [Google Scholar] [CrossRef]
  4. Kim, C.; Kim, C.; Rho, O. Feasibility study of constant eddy-viscosity assumption in gradient-based design optimization. J. Aircr. 2003, 40, 1168–1176. [Google Scholar] [CrossRef]
  5. Martin, L.; Forestier, N.; Colo, L.; Billard, F.; Chalot, F.; Johan, Z.; Mallet, M. Extension of Linearized CFD Methods for Complex Aerodynamic Flows and Application to Unsteady Load Evaluations. In Proceedings of the International Forum on Aerolasticity and Structural Dynamics, IFASD 2022, Madrid, Spain, 13–17 June 2022. [Google Scholar]
  6. Zymaris, A.; Papadimitriou, D.; Giannakoglou, K.; Othmer, C. Continuous Adjoint Approach to the Spalart-Allmaras Turbulence Model for Incompressible Flows. Comput. Fluids 2009, 38, 1528–1538. [Google Scholar] [CrossRef]
  7. Bueno-Orovio, A.; Castro, C.; Palacios, F.; Zuazua, E. Continuous adjoint approach for the Spalart-Allmaras model in aerodynamic optimization. AIAA J. 2012, 50, 631–646. [Google Scholar] [CrossRef]
  8. Haftka, R. Optimization of flexible wing structures subject to strength and induced drag constraints. AIAA J. 1977, 15, 1101–1106. [Google Scholar] [CrossRef]
  9. Sobieszczanski-Sobieski, J. Structural Shape Optimization in Multidisciplinary System Synthesis; NASA Technical Memorandum 100538; Springer: Dordrecht, The Netherlands, 1988. [Google Scholar]
  10. Martins, J.R.R.A.; Alonso, J.; Reuther, J. High-Fidelity Aerostructural Design Optimization of a Supersonic Business Jet. J. Aircr. 2004, 41, 523–530. [Google Scholar] [CrossRef]
  11. Kenway, G.K.W.; Kennedy, G.J.; Martins, J.R.R.A. Scalable Parallel Approach for High-Fidelity Steady-State Aeroelastic Analysis and Adjoint Derivative Computations. J. Comput. Phys. 2014, 52, 935–951. [Google Scholar] [CrossRef]
  12. Abu-Zurayk, M.; Brezillon, J. Shape Optimization Using the Aero-structural Coupled Adjoint Approach for Viscous Flows. In Proceedings of the EUROGEN 2011, Capua, Italy, 14–16 September 2011. [Google Scholar]
  13. Ghazlane, I.; Carrier, G.; Dumont, A.; Desideri, J.A. Aerostructural Adjoint Method for Flexible Wing Optimization. In Proceedings of the 53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Honolulu, HI, USA, 23–26 April 2012. [Google Scholar]
  14. Merle, A.; Ilic, C.; Abu-Zurayk, M.; Häßy, J.; Becker, R.; Schulze, M.; Klimmek, T. High-Fidelity Adjoint-based Aircraft Shape Optimization with Aeroelastic Trimming and Engine Coupling. In Proceedings of the EUROGEN 2019, Guimarães, Portugal, 12–14 September 2019. [Google Scholar]
  15. Abu-Zurayk, M.; Merle, A.; Ilic, C.; Keye, S.; Goertz, S.; Schulze, M.; Klimmek, T.; Kaiser, C.; Quero, D.; Häßy, J.; et al. Sensitivity-based Multifidelity Multidisciplinary Optimization of a Powered Aircraft Subject to a Comprehensive Set of Loads. In Proceedings of the AIAA AVIATION 2020 Forum, Virtual Event, 15–19 June 2020. [Google Scholar]
  16. Bombardieri, R.; Cavallaro, R.; Sanchez, R.; Gauger, N. Aerostructural wing shape optimization assisted by algorithmic differentiation. Struct. Multidiscip. Optim. 2021, 64, 739–760. [Google Scholar] [CrossRef]
  17. Bons, N.; Martins, J.R.R.A. Aerostructural Design Exploration of a Wing in Transonic Flow. Aerospace 2020, 7, 118. [Google Scholar] [CrossRef]
  18. Brooks, T.R.; Martins, J.R.R.A.; Kennedy, G.J. High-fidelity aerostructural optimization of tow-steered composite wings. J. Fluids Struct. 2019, 88, 122–147. [Google Scholar] [CrossRef]
  19. Coder, J.G.; Pulliam, T.H.; Hue, D.; Kenway, G.K.W.; Sclafani, A.J. Contributions to the 6th AIAA CFD Drag Prediction Workshop Using Structured Grid Methods. In AIAA SciTech Forum; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2017. [Google Scholar]
  20. Chalot, F.; Mallet, M.; Ravachol, M. A comprehensive finite element Navier-Stokes solver for low- and high-speed aircraft design. In 32nd Aerospace Sciences Meeting and Exhibit, AIAA 94-0814; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 1994. [Google Scholar]
  21. Trompoukis, X.; Tsiakas, K.; Asouti, V.; Kontou, M.; Giannakoglou, K. Optimization of an Internally Cooled Turbine Blade—Mathematical Development and Application. Int. J. Turbomach. Propuls. Power 2021, 6, 20. [Google Scholar] [CrossRef]
  22. ESI-Group. Available online: https://www.esi-group.com/products/virtual-performance-solution/ (accessed on 1 February 2024).
  23. Kleinveld, S.; Rogé, G.; Daumas, L.; Dinh, Q. Differentiated parametric CAD used within the context of automatic aerodynamic design optimization. In Proceedings of the 12th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Victoria, BC, Canada, 10–12 September 2008. [Google Scholar]
  24. Gagliardi, F.; Giannakoglou, K. A Two-Step Radial Basis Function-Based CFD Mesh Displacement Tool. Adv. Eng. Softw. 2019, 128, 86–97. [Google Scholar] [CrossRef]
  25. Fong, W.; Darve, E. The black-box fast multipole method. J. Comput. Phys. 2009, 228, 8712–8725. [Google Scholar] [CrossRef]
  26. Hughes, T.; Franca, L.; Mallet, M. A New Finite Element Formulation for Computational Fluid Dynamics: I. Symmetric Forms of the Compressible Euler and Navier-Stokes Equations and the Second Law of Thermodynamics. Comput. Methods Appl. Mech. Eng. 1986, 54, 223–234. [Google Scholar] [CrossRef]
  27. Chalot, F. Industrial Aerodynamics; John Wiley & Sons: Hoboken, NJ, USA, 2004. [Google Scholar]
  28. Spalart, P.; Allmaras, S. A one-equation turbulence model for aerodynamic flows. Rech. Aérospatiale 1992, 1, 5–21. [Google Scholar]
  29. Saad, Y.; Schultz, M.H. GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SISC 1986, 7, 856–869. [Google Scholar] [CrossRef]
  30. Martin, L.; Rogé, G. Calcul de la sensibilité d’ordre deux d’une observation aérodynamique. ESAIM Procs 2009, 27, 138–155. [Google Scholar] [CrossRef]
  31. Hascoët, L. TAPENADE: A Tool for Automatic Differentiation of programs. In Proceedings of the 4th European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS), Jyvaskyla, Finland, 24–28 July 2004. [Google Scholar]
  32. Shroff, G.; Herbert, B. Stabilization of unstable procedures: The Recursive Projection Method. SIAM J. Numer. Anal. 1993, 40, 1099–1120. [Google Scholar] [CrossRef]
  33. Brooks, T.R.; Kenway, G.J.; Martins, J.R.R.A. Benchmark Aerostructural Models for the Study of Transonic Aircraft Wings. AIAA J. 2018, 56, 2840–2855. [Google Scholar] [CrossRef]
  34. European Aviation Safety Agency. Available online: https://www.easa.europa.eu/en/certification-specifications/cs-25-large-aeroplanes (accessed on 13 January 2024).
  35. Jakobsson, S.; Amoignon, O. Mesh deformation using Radial Basis Functions for gradient-based aerodynamic shape optimization. Comput. Fluids 2007, 36, 1119–1136. [Google Scholar] [CrossRef]
  36. Farhat, C.; Lesoinne, M.; LeTallec, P. Load and motion transfer algorithms for fluid/structure interaction problems with non-matching discrete interfaces: Momentum and energy conservation, optimal discretization and application to aeroelasticity. Comput. Methods Appl. Mech. Eng. 1998, 157, 95–114. [Google Scholar] [CrossRef]
  37. Degroote, J.; Haelterman, R.; Annerel, S.; Bruggeman, P.; Vierendeels, J. Performance of partitioned procedures in fluid-structure interaction. Comput. Struct. 2010, 88, 446–457. [Google Scholar] [CrossRef]
  38. Kraft, D. Algorithm 733: TOMP–Fortran modules for optimal control calculations. ACM Trans. Math. Softw. 1994, 20, 262–281. [Google Scholar] [CrossRef]
Figure 1. Twist angle (left) and trailing edge (TE) camber (right) effects. In both images, the tip section represents a cut at the winglet junction. Red indicates the baseline wing and green the deformed one.
Figure 1. Twist angle (left) and trailing edge (TE) camber (right) effects. In both images, the tip section represents a cut at the winglet junction. Red indicates the baseline wing and green the deformed one.
Fluids 09 00087 g001
Figure 2. Sections of the wing and HTP directly controlled by the design variables. b 1 , b 2 , and b 3 control the twist, and b 5 , b 6 , and b 7 control the trailing edge camber of sections 0, 1, and 2, respectively. The twist and trailing edge camber of sections 3 and 4 are controlled by the same design variables, namely b 4 and b 8 . b 9 controls the rotation of both (0 and 1) HTP sections.
Figure 2. Sections of the wing and HTP directly controlled by the design variables. b 1 , b 2 , and b 3 control the twist, and b 5 , b 6 , and b 7 control the trailing edge camber of sections 0, 1, and 2, respectively. The twist and trailing edge camber of sections 3 and 4 are controlled by the same design variables, namely b 4 and b 8 . b 9 controls the rotation of both (0 and 1) HTP sections.
Fluids 09 00087 g002
Figure 3. Indicative changes in the GBJ shape resulting from changing (a) b 1 , (b) b 8 and (c) b 9 to their upper values. Blue contours correspond to the datum values of the design variables.
Figure 3. Indicative changes in the GBJ shape resulting from changing (a) b 1 , (b) b 8 and (c) b 9 to their upper values. Blue contours correspond to the datum values of the design variables.
Fluids 09 00087 g003
Figure 4. Comparison of the polar curves computed by AETHER and PUMA, at M i n f = 0.8 and 0.82 .
Figure 4. Comparison of the polar curves computed by AETHER and PUMA, at M i n f = 0.8 and 0.82 .
Fluids 09 00087 g004
Figure 5. Cut across the 3D grid generated around the GBJ baseline geometry and a close-up view of the nacelle–pylon–fuselage area, as used by PUMA. Hexahedra and prisms with high aspect ratios are used close to the solid walls to capture boundary layer physics. Tetrahedra fill the rest of the domain, while pyramids are used to connect elements of different types.
Figure 5. Cut across the 3D grid generated around the GBJ baseline geometry and a close-up view of the nacelle–pylon–fuselage area, as used by PUMA. Hexahedra and prisms with high aspect ratios are used close to the solid walls to capture boundary layer physics. Tetrahedra fill the rest of the domain, while pyramids are used to connect elements of different types.
Fluids 09 00087 g005
Figure 6. S D s of C L w.r.t. the five design variables ( b 1 to b 4 and b 9 ) controlling the wing twist angle distribution and the HTP rotation (left) and those ( b 5 to b 8 ) controlling the wing’s trailing edge camber distribution (right). The S D s computed by the adjoint of PUMA and AETHER are compared with FDs at FC1.
Figure 6. S D s of C L w.r.t. the five design variables ( b 1 to b 4 and b 9 ) controlling the wing twist angle distribution and the HTP rotation (left) and those ( b 5 to b 8 ) controlling the wing’s trailing edge camber distribution (right). The S D s computed by the adjoint of PUMA and AETHER are compared with FDs at FC1.
Fluids 09 00087 g006
Figure 7. The GBJ wing model made of 240 skin panels, 25 ribs, and 48 spars.
Figure 7. The GBJ wing model made of 240 skin panels, 25 ribs, and 48 spars.
Fluids 09 00087 g007
Figure 8. Schematic representation of the inputs and outputs of the two surrogate models ( SM P and SM J ). Purple circles are inputs, while red ones are outputs. Z F W stands for the zero fuel weight of the GBJ. SM J also builds matrix H , which links the displacements of the structural model nodes ( U k S ) and those on the CFD surface grid nodes ( U k F ), as U k F = H U k S [35].
Figure 8. Schematic representation of the inputs and outputs of the two surrogate models ( SM P and SM J ). Purple circles are inputs, while red ones are outputs. Z F W stands for the zero fuel weight of the GBJ. SM J also builds matrix H , which links the displacements of the structural model nodes ( U k S ) and those on the CFD surface grid nodes ( U k F ), as U k F = H U k S [35].
Fluids 09 00087 g008
Figure 9. Schematic representation of the CFD-CSM workflow. Purple circles are inputs whereas white circles indicate intermediate quantities. f G , K and H are computed by SM J .
Figure 9. Schematic representation of the CFD-CSM workflow. Purple circles are inputs whereas white circles indicate intermediate quantities. f G , K and H are computed by SM J .
Fluids 09 00087 g009
Figure 10. Schematic representation of the adjoint CFD-CSM workflow.
Figure 10. Schematic representation of the adjoint CFD-CSM workflow.
Fluids 09 00087 g010
Figure 11. Aerodynamic sensitivities (computed by PUMA at FC1) of C L w.r.t. b 1 to b 4 and b 9 controlling the wing twist angle distribution and the HTP rotation (left) and w.r.t. b 5 to b 8 controlling the wing’s trailing edge camber distribution (right).
Figure 11. Aerodynamic sensitivities (computed by PUMA at FC1) of C L w.r.t. b 1 to b 4 and b 9 controlling the wing twist angle distribution and the HTP rotation (left) and w.r.t. b 5 to b 8 controlling the wing’s trailing edge camber distribution (right).
Fluids 09 00087 g011
Figure 12. Aerodynamic optimization with rigid wing structure using PUMA, for OptR1 (left) and OptR2 (right): Convergence history of the objective and constraint functions.
Figure 12. Aerodynamic optimization with rigid wing structure using PUMA, for OptR1 (left) and OptR2 (right): Convergence history of the objective and constraint functions.
Fluids 09 00087 g012
Figure 13. Aerodynamic optimization with rigid wing structure (OptR1) using PUMA: Pressure coefficient (left) and drag coefficient integrand (right) fields, computed for the baseline and the optimized configurations.
Figure 13. Aerodynamic optimization with rigid wing structure (OptR1) using PUMA: Pressure coefficient (left) and drag coefficient integrand (right) fields, computed for the baseline and the optimized configurations.
Fluids 09 00087 g013
Figure 14. Aerodynamic optimization with rigid wing structure (OptR2) using PUMA: Pressure coefficient (left) and drag coefficient integrand (right) fields, computed for the baseline and the optimized configurations.
Figure 14. Aerodynamic optimization with rigid wing structure (OptR2) using PUMA: Pressure coefficient (left) and drag coefficient integrand (right) fields, computed for the baseline and the optimized configurations.
Fluids 09 00087 g014
Figure 15. Aerodynamic optimization with rigid wing structure (OptR2): Comparison of the convergence histories of C D (left) and C L (constraint, right) during the optimization, using AETHER (red) and PUMA (blue).
Figure 15. Aerodynamic optimization with rigid wing structure (OptR2): Comparison of the convergence histories of C D (left) and C L (constraint, right) during the optimization, using AETHER (red) and PUMA (blue).
Fluids 09 00087 g015
Figure 16. Aerostructural optimization with fixed structural model (OptR5: left and OptR6: right) using PUMA: Convergence history of the objective and constraint functions.
Figure 16. Aerostructural optimization with fixed structural model (OptR5: left and OptR6: right) using PUMA: Convergence history of the objective and constraint functions.
Fluids 09 00087 g016
Figure 17. Aerostructural optimization with fixed structural model (OptR5: top and OptR6: bottom), using PUMA: Deformation fields plotted on the optimized shape; the wing deformation due to design variable changes (from baseline to optimized x F * , left) is superimposed to the wing deformation due to flexibility (from the optimized x F * to the optimized flight shape x F , center), forming the total deformation field (from baseline x F * to optimized x F , right).
Figure 17. Aerostructural optimization with fixed structural model (OptR5: top and OptR6: bottom), using PUMA: Deformation fields plotted on the optimized shape; the wing deformation due to design variable changes (from baseline to optimized x F * , left) is superimposed to the wing deformation due to flexibility (from the optimized x F * to the optimized flight shape x F , center), forming the total deformation field (from baseline x F * to optimized x F , right).
Fluids 09 00087 g017
Figure 18. Aerostructural optimization with fixed structural model (OptR6), using AETHER: Comparison of the convergence histories of OptR2 (purple) and OptR6 (green) wing structure. The blue square corresponds to the aerostructurally trimmed configuration (OptR6B) starting from the aerodynamically optimized one.
Figure 18. Aerostructural optimization with fixed structural model (OptR6), using AETHER: Comparison of the convergence histories of OptR2 (purple) and OptR6 (green) wing structure. The blue square corresponds to the aerostructurally trimmed configuration (OptR6B) starting from the aerodynamically optimized one.
Fluids 09 00087 g018
Figure 19. Aerostructural optimization with varying structural model (OptR9), using PUMA: Evolution of the design variables ( b 1 to b 11 ) during the optimization.
Figure 19. Aerostructural optimization with varying structural model (OptR9), using PUMA: Evolution of the design variables ( b 1 to b 11 ) during the optimization.
Fluids 09 00087 g019
Figure 20. Aerostructural optimization with varying structural model (OptR9), using PUMA: Comparison of the wing shape of the baseline GBJ and the outcome of OptR9.
Figure 20. Aerostructural optimization with varying structural model (OptR9), using PUMA: Comparison of the wing shape of the baseline GBJ and the outcome of OptR9.
Fluids 09 00087 g020
Figure 21. Aerostructural optimization with varying structural model (OptR9), using PUMA: Pressure coefficient fields computed on the suction (left) and pressure (right) side for the baseline and the aerostructurally optimized aircraft.
Figure 21. Aerostructural optimization with varying structural model (OptR9), using PUMA: Pressure coefficient fields computed on the suction (left) and pressure (right) side for the baseline and the aerostructurally optimized aircraft.
Fluids 09 00087 g021
Figure 22. Aerostructural optimization with varying structural model (OptR9), using PUMA: Convergence history of the fuel weight and the ZFW.
Figure 22. Aerostructural optimization with varying structural model (OptR9), using PUMA: Convergence history of the fuel weight and the ZFW.
Fluids 09 00087 g022
Table 1. Design variable bounds used for generating the surrogate parameterization tool; these are also imposed as constraints during the optimization runs. b i + , b i are the upper and lower bounds, respectively, and b i D are the tabulated datum values.
Table 1. Design variable bounds used for generating the surrogate parameterization tool; these are also imposed as constraints during the optimization runs. b i + , b i are the upper and lower bounds, respectively, and b i D are the tabulated datum values.
b i D b i b i +
b 1 to b 4 0 2.0 2.0
b 5 to b 8 0.0 0.02 0.02
b 9 0 2.0 2.0
b 10 0 2.0 2.0
b 11 1.0 0.9 1.1
Table 2. Aerodynamic optimization with rigid wing structure (OptR2): Optimal values of the design variables computed by the AETHER- and PUMA-based runs.
Table 2. Aerodynamic optimization with rigid wing structure (OptR2): Optimal values of the design variables computed by the AETHER- and PUMA-based runs.
Var ID b 1 b 2 b 3 b 4 b 5 b 6 b 7 b 8 b 9
AETHER−1.5791.2900.4090.454−0.002080.004220.00506−0.0176−0.0374
PUMA−2.0001.1290.3180.178−0.006450.003750.00446−0.0171−0.1752
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.

Share and Cite

MDPI and ACS Style

Tsiakas, K.; Trompoukis, X.; Asouti, V.; Giannakoglou, K.; Rogé, G.; Julisson, S.; Martin, L.; Kleinveld, S. Discrete and Continuous Adjoint-Based Aerostructural Wing Shape Optimization of a Business Jet. Fluids 2024, 9, 87. https://doi.org/10.3390/fluids9040087

AMA Style

Tsiakas K, Trompoukis X, Asouti V, Giannakoglou K, Rogé G, Julisson S, Martin L, Kleinveld S. Discrete and Continuous Adjoint-Based Aerostructural Wing Shape Optimization of a Business Jet. Fluids. 2024; 9(4):87. https://doi.org/10.3390/fluids9040087

Chicago/Turabian Style

Tsiakas, Konstantinos, Xenofon Trompoukis, Varvara Asouti, Kyriakos Giannakoglou, Gilbert Rogé, Sarah Julisson, Ludovic Martin, and Steven Kleinveld. 2024. "Discrete and Continuous Adjoint-Based Aerostructural Wing Shape Optimization of a Business Jet" Fluids 9, no. 4: 87. https://doi.org/10.3390/fluids9040087

APA Style

Tsiakas, K., Trompoukis, X., Asouti, V., Giannakoglou, K., Rogé, G., Julisson, S., Martin, L., & Kleinveld, S. (2024). Discrete and Continuous Adjoint-Based Aerostructural Wing Shape Optimization of a Business Jet. Fluids, 9(4), 87. https://doi.org/10.3390/fluids9040087

Article Metrics

Back to TopTop