Gradient-Based Aero-Stealth Optimization of a Simplified Aircraft

Abstract

Modern fighter aircraft increasingly need to conjugate aerodynamic performance and low observability. In this paper, we showcase a methodology for a gradient-based bidisciplinary aero-stealth optimization. The shape of the aircraft is parameterized with the help of a CAD modeler, and we optimize it with the SLSQP algorithm. The drag, computed with the help of a RANS method, is used as the aerodynamic criterion. For the stealth criterion, a function is derived from the radar cross-section in a given cone of directions and weighed with a function whose goal is to cancel the electromagnetic intensity in a given direction. Stealth is achieved passively by scattering back the electromagnetic energy away from the radar antenna, and no energy is absorbed by the aircraft, which is considered as a perfect conductor. A Pareto front is identified by varying the weights of the aerodynamic and stealth criteria. The Pareto front allows for an easy identification of the CAD model corresponding to a chosen aero-stealth trade-off.

1. Introduction

With the omnipresence of radar technology, from early detection to target acquisition by a hostile air defense system to active radar homing for missiles, modern fighter aircraft design needs to take stealth into consideration. The designer needs to maximize aerodynamic performances while keeping a low Radar Cross-Section (RCS).

The frequency bands commonly used by modern radars are the L, S, C, and X bands [1]. Their working frequencies range from $1 \mathrm{GHz}$ to $12 \mathrm{GHz}$, corresponding to wavelengths between $0.025 \mathrm{m}$ and $0.3 \mathrm{m}$, which are much smaller than typical aircraft dimensions. Therefore, the radar waves are considered high-frequency. This means that the scattering characteristics of the aircraft are heavily influenced by its geometric details. On the other hand, the aerodynamic qualities of an aircraft can be very sensitive to even small changes to its geometry. There is then a need for methodologies able to efficiently design shapes that balance aerodynamic and stealth performances.

Some studies such as [2] deal with the contradictory effects of shape on aerodynamic efficiency and stealth, but few of them use integrated aero-stealth optimization methods. Among them, Lee et al. [3] used a robust evolutionary algorithm for the wing planform and sections design. In order to mitigate the large computation costs due to the number of objective function evaluations incurred by the evolutionary algorithm, they used a potential flow and physical optics solver. Later, Wu et al. [4] used the adaptive metamodel method in order to optimize a high aspect ratio wing for the aerodynamics, stealth, and structure, while Jiang et al. [5] used the full factorial design method in combination with radial basis functions in order to optimize the aero-stealth of the rotor airfoil of a helicopter.

All those studies used the gradient-free optimization method, which has the advantage of not necessitating the costly computation of derivatives, which can be hard to derive or implement. However, they need a high number of objective function evaluations, whose computation may be prohibitive when using high-fidelity solvers. In an industrial design process, gradient-free algorithms are typically used for preliminary design with low-fidelity models, while gradient-based methods are more adapted to large-scale aircraft shape design with mid- to high-fidelity models. Furthermore, surrogate-based algorithms are known to struggle with constraints handling, satisfying them with a poor accuracy and struggling to meet the stop criteria [6].

1.1. Gradient-Based Method

It is well known that gradient-based methods usually necessitate significantly less function evaluations than their gradient-free counterpart. For example, Lyu et al. [7] benchmarked several optimization algorithms, both gradient-free and gradient-based, on three different problems with increasing complexity. They showed that the number of evaluations required by gradient-free methods tend to grow as the square or the cube of the number of optimization parameters, while the gradient-based methods keep a linear behavior. It is clear form their work that as soon as the design space has a higher dimension than ten, gradient-free algorithms require tens or hundreds of thousands of function evaluation while gradient-based algorithms using an analytic method such as the adjoint formulation to compute the gradient require less than a hundred.

The adjoint formulation of the problem allows for an efficient computation of the gradient of the objective and constraints with respect to the optimization parameters. For the applications of this study, the cost of the gradient derivation is roughly linear with respect to the number of design variables and the size of the mesh. Various publications have already covered the application of the adjoint method to gradient-based optimization. For instance, Martin applied the adjoint method to robust airfoil optimization [8]. Martin used the Automatic Differentiation (AD) tool TAPENADE [9] to develop his discrete adjoint solver which drastically reduced the implementation time. Lyu et al. applied their own adjoint solver to the optimization of a blended-wing-body aircraft [10] and to the Common Research Wing [11].

However, it is well known that gradient-based algorithms have a lower probability of converging towards a global minimum as opposed to a local minimum than global, gradient-free algorithms and are very sensitive to the starting point in the case where the response surface in the design space is multimodal. Furthermore, there is increasing evidence that the aerodynamic shape optimization is multimodal [12], thus making gradient-based algorithms ineffective at exploring unusual yet potentially revolutionary designs. In an industrial framework, this difficulty is solved by a back-and-force between the early design phase where global algorithms are used in conjunction with lower-fidelity methods in order to select a limited number of interesting designs, which are then investigated further during the preliminary design phase with a higher fidelity method and local, gradient-based optimization algorithms. In this paper, we study a use-case representative of the preliminary design phase.

1.2. Aero-Stealth Optimization

As for aero-stealth optimization, gradient-free algorithm have proven very popular in the past few years. For instance, Garnier et al. [13] used an algorithm based on the coupling of the NSGA-II genetic algorithm and the multi-objective efficient global optimization approach [14] to optimize an UCAV planform for aero-stealth. Liu et al. [15] used hierarchical kriging models and other surrogate models for the aero-stealth optimization of a flying-wing aircraft. Zhou and Huang [16] used a two-particle search algorithm to optimize the aero-stealth performances of a vertical tail. Finally, Taj et al. [17] used a machine learning algorithm called a Gaussian process regression to perform the design exploration of aerodynamics and RCS for a fighter aircraft and a genetic algorithm to minimize drag and RCS. Because they used the gradient-free algorithm, these works are restricted to using surrogate models or low-fidelity methods that are not suited for the level of precision required at the stage of preliminary design in an industrial process.

Few studies, however, have dealt with the gradient-based aero-stealth optimization. Pan et al. [18] performed a CAD-based, gradient-based aero-stealth optimization of a flying-wing aircraft, solving the Prandtl–Glauert equation to compute the pressure coefficients of the aircraft and a Physical Optics (PO) method to predict RCS. Li et al. [19] used a RANS and a PO solver as well as a free-form deformation algorithm in order to perform a CAD-free, gradient-based aero-stealth optimization of a flying-wing aircraft. They performed various optimizations using the aerodynamic and stealth criteria alternatively as the objective and a constraint and showed that the better-performing optimization result is obtained when using an objective function built as a trade-off of the aerodynamic and stealth criteria.

1.3. Contribution

The purpose of this study is to showcase a gradient-based optimization framework in an industrial context to efficiently explore aero-stealth design trade-offs using our in-house CAD modeler GANIMEDE [20], a RANS solver [21], a physical optics solver, and the Sequential Least Square Quadratic Programming (SLSQP) algorithm [22]. We used a realistic criterion as a PO objective that is representative of modern stealth shape design consideration. We demonstrate our methodology on a simplified aircraft shape with 14 design variables. Because we use a RANS solver, which is a medium-fidelity solver whose computation time is about 14 min in the case we consider, a design space with more than ten design variable is too expensive for a gradient-free algorithm to explore. We identify a Pareto front to showcase various optimal shapes corresponding to various aero-stealth trade-offs.

The use of CAD-based optimization offers a lot of advantages in an industrial framework. First, it greatly reduces the number of design variables and constraints needed to parameterize and manipulate a realistic aircraft shape and allows us to use variables that are independent of one another. It also allows us to easily specify and compute geometric constraints such as the wing surface reference (see Section 4.4) or the wing absolute thickness. Furthermore, it is possible to define design variable that make sense from an aeronautic designer’s point of view, such as the leading edge sweep angle or the trailing edge camber. Lastly, the CAD geometry is an important piece of the aircraft design process, being a tool of communication between the different disciplines at each step of the design process, from the early design step to the manufacturing step.

This paper is organized as follows. In Section 2, we present the gradient-based aero-stealth coupled optimization process. Figure 1 shows the process as a Design Structure Matrix (DSM). We also describe some numerical methods used throughout this study. Next, the results of the code validation are shown in Section 3 and, in Section 4 and Section 5, we introduce the problem formulation, baseline geometry of the flying wing, meshes, and shape design variables. In Section 6, the accuracy of the sensitivities is tested. After that, in Section 7, a Pareto front is showcased and particularly interesting geometries are highlighted. Finally, in Section 8, we discuss the conclusions based on the results.

2. Methodology

2.1. The Physical Optics Approximation

In an aircraft, stealthiness can mean a number of different things. Here, we understand it as the capacity of the aircraft to scatter a limited amount of energy back into the direction of the emitter. We consider a monostatic radar, which represents the large majority of the systems in use. This means that the emitting antenna and the receiving antenna are collocated. A popular measure of the scattered electromagnetic energy is the radar cross-section (RCS) $\sigma$ defined as follows:

(1)$\sigma (\theta ,\varphi )=\underset{R\to \infty }{lim}{R}^2{\displaystyle \frac{|{\mathbf{E}}_P{|}^2}{|{\mathbf{E}}_i{|}^2}},$

where ${\mathbf{E}}_P$ is the scattered electric fields at a point $P(R,\theta ,\varphi )$ in polar coordinates (defined in Figure 2) and ${\mathbf{E}}_i$ is the incident electric field at the target, whose norm is considered to be 1 $\mathrm{V}/\mathrm{m}$ in practice. Usually, the RCS is expressed in decibels:

${\sigma }_{\mathrm{dB}}=10log10\left(\sigma \right).$

We restrict this study to radar bands L (1–5 GHz, used for moderate-distance surveillance) through X (around 10 GHz, used for missile guidance systems). These frequency bands produce on a target the size of an aircraft of predominantly high-frequency scattering [1]. Thus, an asymptotically exact method such as a Boundary Element Method (BEM) or Fast Multipole Method (FMM) would be too expensive for an optimization process. PO is one of the most common high-frequency approximate methods for RCS prediction. It is much cheaper than BEM, and it gives a better prediction than a purely optical approximation [23].

For such frequencies, the size of the aircraft is large with respect to the wavelength. Another reasonable assumption is that the distance between the aircraft and the radar antenna is also large with respect to the wavelength, so that the incident electric field can be considered a plane wave, and with respect to the size of the aircraft, so the far-field Green’s function approximation to the Stratton–Chu [24] integral equation is valid. Together with a tangent plane approximation, these can be simplified, and the scattered electric field ${\mathbf{E}}_P$ at a point P can be expressed as follows:

(2)${\mathbf{E}}_P=-j{\displaystyle \frac{k{Z}_0}{4\pi }}{\displaystyle \frac{e^{-jkR}}{R}}\underset{S}{\int }{\mathbf{u}}_s\times ({\mathbf{u}}_s\times \mathrm{J})exp\left(jk\mathbf{R}·({\mathbf{u}}_i-{\mathbf{u}}_s)\right),$

with $j=\sqrt{-1}$, k the wave number, $Z_0$ the impedance of free space, S the surface of the target, ${\mathbf{u}}_i$ and ${\mathbf{u}}_s$ the unit direction vectors for the incident and scattered fields, respectively, $\mathbf{R}=\overrightarrow{OP}$, and $R=\left|\mathbf{R}\right|$, and where J is the surface current density vector. In the PO approximation,

(3)$\mathrm{J}=\left\{\begin{array}{cc}{\displaystyle \frac{2}{Z_0}}\mathbf{n}\times {\mathbf{u}}_i\times {\mathbf{E}}_i& \hfill \mathrm{in}\mathrm{the}\mathrm{illuminated}\mathrm{zone},\\ \mathbf{0}& \hfill \mathrm{in}\mathrm{the}\mathrm{shadow}\mathrm{zone},\end{array}\right.$

where $\mathbf{n}$ is the local normal unit vector to the surface.

Note that by (2),

${\mathbf{E}}_P=\mathcal{O}\left({\displaystyle \frac{1}{R}}\right)\mathrm{as}R\to \infty .$

Therefore, the RCS does not depend on the range.

Although this approximate method does not account for diffraction, and in particular wedge diffraction, extensions exist to correct the induced current along wedges. One such extension is called Partial Theory of Diffraction (PTD) [25] and has been added to PO to give a better prediction of the RCS of an angular shape, such as the trailing edge of a wing.

The solver including PO and PTD along with their differentiation, described on the next section, is dubbed COPOLA.

2.2. Automatic Differentiation of Physical Optics

If the characteristics of the radar wave (wavelength, polarization, and direction) are given, the RCS depends only on ${\mathbf{x}}_e$, the coordinates of the vertices of the surface grid used for the electromagnetism problem.

If

$F_e\left(\nu \right)={\widehat{F}}_e\left(\sigma \left({\mathbf{x}}_e\left(\nu \right)\right)\right)$

is a function depending explicitly only on the RCS $\sigma$, then the directional derivative of $F_e$ with respect to a set of geometric variables $\nu$ (e.g., sweep angle, chord length) is

(4)$\mathrm{d}{F}_e={\widehat{F}}_e^{\prime }\left(\sigma \right)·{\sigma }^{\prime }\left({\mathbf{x}}_e\right)·{\mathbf{x}}_e^{\prime }\left(\nu \right)·\mathrm{d}\nu .$

We used the Automatic Differentiation (AD) tool TAPENADE version 3.16 [9] in order to generate a differentiated program from the original code, which is able to compute the terms $\mathrm{d}{F}_e/\mathrm{d}\sigma$ and $\mathrm{d}\sigma /\mathrm{d}{\mathbf{x}}_e$. With this tool, it is possible to save a lot of development and debugging time. TAPENADE offers two modes: tangent (or direct) and reverse. The tangent mode computes the directional derivative of $F_e$ in the direction $\mathrm{d}\nu$ using the chain rule, from left to right. This mode is immediate to use as it only needs to create a copy of the input code, interleaved with the derivative instructions. The reverse mode is an efficient way of computing the gradient of $\nabla F_e·\mathrm{d}\nu =\mathrm{d}{\nu }^{\ast }·{\mathbf{x}}_e^{\prime }{\left(\nu \right)}^{\ast }·{\sigma }^{\prime }{\left({\mathbf{x}}_e\right)}^{\ast }·\widehat{F}{e}^{\prime }{\left(\sigma \right)}^{\ast }$. When $F_e$ is scalar-valued and the dimension of $\nu$ is large, and assuming that the storage capacity is available, it is theoretically more efficient to compute the Jacobian matrix $\nabla F_e$ of $F_e$ row by row and then compute $\mathrm{d}{F}_e=\nabla F_e\mathrm{d}\nu$ than to compute it by the chain rule with the tangent mode.

Here, however, all the steps of the chain rule computation are cheap to compute since they do not require us to solve a system. Therefore, the actual gain in efficiency from the reverse mode is expected to be small, and since it is much more difficult to implement than the tangent mode, we did not deem it worth the trouble.

The term $\mathrm{d}{\mathbf{x}}_e/\mathrm{d}\nu$, which is the derivative of the vertex coordinates of the surface grid with respect to a set of geometric variables, is computed by the CAD modeler (see Section 4.1) and is given as an input to the differentiated PO program.

This allows us to compute the cost function $F_e$ as well as its derivative $\mathrm{d}{F}_e/\mathrm{d}\nu$ in a single call of the PO code.

2.3. The Aerodynamic Solver

Our in-house Navier–Stokes AETHER code [21,26] solves the 2D, axisymmetric, and 3D compressible Navier–Stokes equations. It uses a finite element approach based on a symmetric form of the equations written in terms of entropy variables. The advantages of this change of variables are numerous: in addition to the strong mathematical and numerical coherence they provide such as dimensionally correct dot product, symmetric operators with positive properties and efficient preconditioning, entropy variables yield further improvements over the usual conservation variables, particularly in the context of chemically reacting flows.

The code can handle the unstructured mixture of numerous types of elements (triangles and quadrilaterals in 2D; tetrahedral, bricks, and prisms in 3D).

Several one- and two-equation Reynolds-averaged turbulence models are available: Spalart–Allmaras, (k, $\epsilon$), (k, $\omega$), (k, kl), DRSM… Additional equations devoted to the evaluation of the laminar–turbulence transition can be activated. These models are integrated down to the wall. Extensions to DES and LES are also available.

The linearization of Navier–Stokes equations and turbulence equations requires the differentiation of all the routines involved in the aerodynamic computation. These differentiated routines were obtained using TAPENADE, an Automatic Differentiation tool. This tool provides the differentiated Fortran routines with respect to input variables prescribed by the user.

We used a compressible Reynolds-Averaged Navier–Stokes [27] method with the Spalart–Allmaras [28] turbulence model in order to compute the aerodynamic criterion and constraints.

2.4. Adjoint Approach for CFD

Let W be a vector state solution of the non-linear RANS equation system which we note as

(5)$\mathcal{R}=0.$

Finally, let $F_a$ be an observation function of the field W. For the process of optimization, we need to compute the derivatives ${\displaystyle \frac{\mathrm{d}{F}_a}{\mathrm{d}\xi }}$ and ${\displaystyle \frac{\mathrm{d}{F}_a}{\mathrm{d}\nu }}$ with respect to a set of aerodynamic variables $\xi$ (e.g., the angle of attack) and geometric variables $\nu$ (e.g., the wing sweep angle), respectively. The variables $\nu$ and $\xi$ are independent of one another.

In order to mitigate the cost of computing the gradient of the aerodynamic criterion with respect to many variables, we use a double adjoint formulation [29].

When a surface grid of coordinates ${\mathbf{x}}_a$ is subjected to a deformation owing to the CAD modeler, the deformation is propagated to the 3D grid of coordinates ${\mathbf{X}}_a$ with the help of an elliptic operator described by the system:

(6)$\mathcal{L}({\mathbf{X}}_a,{\mathbf{x}}_a)=0.$

Let us write $F_a$ and $\mathcal{R}$ as functions of $\xi$, ${\mathbf{X}}_a$, and W. Note that $\xi$ is an independent variable but ${\mathbf{X}}_a$ depends on $\nu$ and W depends on both $\xi$ and $\nu$.

2.4.1. Gradient with Respect to Aerodynamic Variables

First, the derivative of $F_a$ and $\mathcal{R}$ with respect to $\xi$ can be written as

(7)$\mathrm{d}{F}_a=\left({\left.{\displaystyle \frac{\partial F_a}{\partial \xi }}\right|}_W-{\mathsf{\Psi }}^{\mathrm{T}}{\left.{\displaystyle \frac{\partial \mathcal{R}}{\partial \xi }}\right|}_W\right)\mathrm{d}\xi .$

where $\mathsf{\Psi }$ is the solution to the aerodynamic adjoint system:

(8)${\left({\left.{\displaystyle \frac{\partial \mathcal{R}}{\partial W}}\right|}_{\xi }\right)}^{\mathrm{T}}\mathsf{\Psi }={\left({\left.{\displaystyle \frac{\partial F_a}{\partial W}}\right|}_{\xi }\right)}^{\mathrm{T}}.$

2.4.2. Gradient with Respect to Geometric Variables

The derivative of $\mathcal{R}$ and $F_a$ with respect to $\nu$ is

(9)$\mathrm{d}{F}_a=\left({\left.{\displaystyle \frac{\partial F_a}{\partial {\mathbf{X}}_a}}\right|}_W-{\mathsf{\Psi }}^{\mathrm{T}}{\left.{\displaystyle \frac{\partial \mathcal{R}}{\partial {\mathbf{X}}_a}}\right|}_W\right){\displaystyle \frac{\partial {\mathbf{X}}_a}{\partial \nu }}\mathrm{d}\nu ,$

where $\mathsf{\Psi }$ is the same as defined in Equation (8).

In this case, it is not enough, though, as ${\displaystyle \frac{\partial {\mathbf{X}}_a}{\partial \nu }}$ is expensive to compute. In order to find a more convenient expression, we must write an adjoint problem of the system, which finally yields

(10)$\mathrm{d}{F}_a=-{\Phi }^{\mathrm{T}}{\left.{\displaystyle \frac{\partial \mathcal{L}}{\partial {\mathbf{x}}_a}}\right|}_{{\mathbf{X}}_a}{\displaystyle \frac{\partial {\mathbf{x}}_a}{\partial \nu }}\mathrm{d}\nu ,$

where $\Phi$ is the solution of the mesh adjoint equation system:

(11)${\left({\left.{\displaystyle \frac{\partial \mathcal{L}}{\partial {\mathbf{X}}_a}}\right|}_{{\mathbf{x}}_a}\right)}^{\mathrm{T}}\Phi ={\left({\left.{\displaystyle \frac{\partial F_a}{\partial {\mathbf{X}}_a}}\right|}_{{\mathbf{x}}_a}\right)}^T-{\left({\left.{\displaystyle \frac{\partial \mathcal{R}}{\partial {\mathbf{X}}_a}}\right|}_{{\mathbf{x}}_a}\right)}^T\mathsf{\Psi }.$

3. Validation

3.1. Validation of the COPOLA Code

We used a sphere geometry as a test case for the validation of COPOLA. According to Mie’s theory, if a sphere of radius R is illuminated by a wave of wavelength $\lambda$, then the RCS $\sigma$ converges to the section area of the sphere as the wavelength decreases:

$\underset{\lambda \to 0}{lim}\sigma =\pi R^2.$

Moreover, due to the symmetric properties of the sphere, its RCS does not depend on the bearing angle.

Since the PO method is a high-frequency approximation, it should predict an RCS close to the analytical value. Figure 3 shows the PO prediction of the RCS of a sphere of radius 1 m with respect to the bearing angle (Figure 3a) and frequency (Figure 3b). The plot in Figure 3b is computed with 2000 frequencies between $10 \mathrm{M}\mathrm{Hz}$ and $20 \mathrm{G}\mathrm{Hz}$.

In Figure 3a, the RCS is independent of the bearing angle, as it should be, and the error of the analytical value is around 0.15%. Figure 3b shows that when the frequency is too low (${\displaystyle \frac{2\pi fR}{c}}<50$, where c is the speed of light), the predicted RCS is very far from the analytical value, and the PO approximation is not trustworthy. When the frequency is too high (${\displaystyle \frac{2\pi fR}{c}}>100$), the prediction starts to deteriorate again due to the mesh being too coarse for the higher frequencies.

Figure 4 shows the comparison between an RCS simulation and a BEM computation on the GFA geometry (described in Section 4). The BEM computation was performed owing to the SPECTRE solver [30,31].

PO is able to predict accurately both the height and the width of the two main peaks. Away from of these peaks, the RCS is very low (around $-20$ dB). The precision of PO is not as good in these zones because low-frequency effects such as creeping waves are predominant, but since the RCS is so much lower than that of the peaks, it is negligible for the purpose of this study.

We present the validation of the sensitivities in Section 6.

3.2. Validation of the Aerodynamic Tool

AETHER is the solver used by Dassault-Aviation in an industrial context. It has been validated extensively both with numerical benchmark cases and against experimental data. For example, see [32] or [21].

4. Baseline Description

The geometry we use is a simplified shape of a Generic Fighter Aircraft (GFA, see Figure 5). It features the fuselage without air intake and a pair of wings. The wings have a trapezoidal shape without leading edge extension. It does not have a tail fin, horizontal stabilizer, or canard. The front of the fuselage has a ridge on both sides. We did not use a realistic geometry for the rear of the fuselage, so the physics of the wake is not accurately represented. Since we only performed supersonic simulations, this does not disturb the flow upstream. This also allows us to be lenient with meshing the downstream region, thus saving computation time.

4.1. The CAD Modeler

In order to model and modify the geometry, we used the CAD tool GANIMEDE [20]. In this study, we used some functionalities of GANIMEDE that have not yet been differentiated, which prevented us from using the exact differentiation presented in [20]. The surface is represented by a set of NURBS grids and parameterized with control points. GANIMEDE uses both global and local variables as input. Local variables are the osculating parameters (i.e., position, tangent, and curvature) at a control point. Global variables allow us to modify usual aeronautical design variables such as sweep and twist angles, or the chord of a wing. In this study, we used only global variables.

If a surface mesh is associated with the baseline geometry, GANIMEDE projects the mesh onto the deformed geometry and computes the derivatives of the mesh node coordinates with respect to the design variables ${\displaystyle \frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\nu }}$ using centered finite difference. Because of the smoothing and other refinement operations performed by the CAD modeler on the NURBS geometry, too small a step size can lead to a bad gradient, so we have to find an optimal step size for each design variable. The optimal step sizes are determined on an individual basis by analyzing the field of the second derivative of the displacement of the geometry with respect to each design variable. This field should preserve the spatial symmetries or anti-symmetries of the original design variable and should not be too distorted. The final step sizes that we determined for each design variable are presented in Table 1, which summarizes the chosen design variables.

4.2. Aerodynamic Analysis of the Baseline

We performed a CFD simulation of the GFA geometry. The upstream Mach number is 1.8, with an angle of attack of $2°$. The reference surface used for the computation of the aerodynamic coefficients is $68{\mathrm{m}}^2$. The overall drag coefficient is a $387.5$ drag count and the lift coefficient is $0.113$. These values are presented in Table 2.

The field of the local Mach number in the symmetry plane and in the planes parallel to the symmetry plane at $Y/b=0.45$ and $Y/b=0.8$ are displayed in Figure 6, Figure 7, and Figure 8, respectively.

In the symmetry plane (Figure 6), we can clearly see the weak, attached shock on the nose of the GFA. A large detached area is present at the rear of the aircraft, which is expected since no attention was directed toward this area.

In Figure 7 and Figure 8, directly upstream of the wing leading-edge, we can see a weak bow shock.

4.3. Stealth Analysis of the Baseline

The span of the baseline shape being $6\mathrm{m}$, the frequency should be higher than $25\mathrm{MHz}$ in order to meet the criterion of Section 3:

$50<{\displaystyle \frac{2\pi f{L}_{\mathrm{r}ef}}{c}}<100$

The COPOLA simulation was performed on a range of bearing angles from 0 to $70°$ with a step of $0.011°$ and on a range of elevation angles from −20 to $-10°$ with the same step in order to simulate an aircraft reflecting hostile waves from a radar on the ground far in front of the plane. For each line of sight, the simulation is run for 11 frequencies around a main frequency at $1\mathrm{GHz}$ and spanning $100\mathrm{MHz}$. The step for the bearing and elevation angles has been chosen according to the Shannon criterion.

In Figure 9a,b, a peak is visible around $48°$ of bearing. This corresponds to the direction the leading edge is facing. This peak is the specular reflection of the incident wave onto the leading edge. A second peak appears at around $-15°$. This peak has a much lower intensity compared to the specular reflection peak due to the non-continuity of the normal at the trailing edge wedges. Physically, this means that the trailing edge wedges diffract the incident wave. However, the peak appears to have the same intensity no matter the polarization, which, given the orientation of the trailing edge (see Figure 5) should not be the case: the trailing edges being horizontal, the diffraction should be negligible when observed with a horizontal polarization [33]. We conclude that PO alone is not enough to compute the trailing edge diffraction, and we need to add the PTD extension to our simulation to adequately predict the RCS. When taking diffraction into account, the second peak appears much brighter with the vertical polarization (Figure 9c), only a few dB less intense than the specular reflection peak, and it is almost absent with the horizontal polarization (Figure 9d).

Over the rest of the observation window, the RCS is very low, around −30 to −20 dB, which is practically negligible but increases as the bearing angle comes close to $70°$. This is likely due to the influence of the fuselage, but as it is far away from the axis of the plane, it is not relevant to this study. One might note, however, that the value around $70°$ is not the same with or without the PTD extension. This is likely due to a diffraction by the wing tip of wedges on the fuselage that are overestimated by PO alone.

Since the vertical polarization is strictly richer than the horizontal polarization in our case, it is the one we used for the RCS simulation in the rest of this study.

4.4. Optimization Parameters

The wing of the GFA is parameterized with the help of 6 control sections, with 15 control points along each of them (see Figure 10). All the sections are and stay parallel to the symmetry plane. Section number 1 is at the symmetry plane and is not modified. Section 2 is a few centimeters away from the root of the wing. At Section 4, the continuity of the tangent is not ensured, allowing for a sharp angle on the wing. The fuselage is not modified.

The geometry is parameterized with the help of 13 global variables:

• The leading edge sweep angle in Section 2 and Section 4;

• The chord of Section 2, Section 4 and Section 6;

• The distance between Section 2 and Section 6;

• The twist angle of Section 2, Section 3, Section 5 and Section 6;

• The leading edge camber of Section 3, Section 5 and Section 6.

The chords of Section 3 and Section 5 are interpolated between those of Section 2, Section 4, and Section 6, respectively. The twist angle for a given section is a rotation along an axis perpendicular to the section going through a point at 25% chord and 50% thickness. Both the leading edge camber and the twist angle in Section 4 are interpolated between Section 3 and Section 5.

Finally, the angle of attack is the 14th optimization parameter. It is private to the aerodynamic process.

In Table 1, we summarize the design variable as well as various pieces of information such as their index, their bounds, the step size used for the computation of the gradient of the displacement, their baseline values, and their value for the monodiscipline optimized shapes.

5. Cost and Constraints

In order to explore multiple trade-offs between aerodynamic and stealth optima, we define the cost function as a weighed sum of an aerodynamic criterion and a stealth criterion:

(12)$F=w{F}_a+(1-w)F_e,$

with $0\le w\le 1$.

5.1. Aerodynamic Criterion

For the aerodynamic optimization, we minimize the drag coefficient (${\mathrm{C}}_{\mathrm{D}}$):

(13)$F_a={\displaystyle \frac{C_D}{{C_D}^{ref}}}.$

This is somewhat simplistic and will lead to unrealistic shapes. Practically, a criterion based on the pitching moment coefficient should be added to the drag criteria, as well as a criterion based on the rolling stability or maneuverability.

However, since the goal of this study is to showcase the methodology rather than find a realistic shape, we settled for this simplified criterion.

5.2. Stealth Criterion

Typically, a stealth aircraft shape is designed to scatter radar wave away from any hostile radar antenna, which is usually located roughly in front of the aircraft and slightly below. For example, the upward facing planes of the F117 scatter the radar wave into the sky, while more modern aircraft, like the B2, are built to concentrate the reflected radar wave only as a very narrow signal in a few specific directions (“peak directions”). The stealth criterion is built to reflect and favor this behavior: we want to minimize the RCS in every direction in a given sector $({\theta }_0,{\theta }_1)\times ({\varphi }_0,{\varphi }_1)$, except for a limited interval around a peak direction where the RCS can be high. We built an electromagnetic criterion in order to reflect this design methodology as follows.

In order to manage dimensionless quantities, the electromagnetic criterion is defined as follows:

(14)$F_e={\displaystyle \frac{\tilde{\sigma }}{{\tilde{\sigma }}^{ref}}},$

where $\tilde{\sigma }$ is a weighed average of a function of ${\sigma }_{\mathrm{dB}}$:

(15)$\tilde{\sigma }={\int }_{{\theta }_0}^{{\theta }_1}{\int }_{{\varphi }_0}^{{\varphi }_1}{\mathsf{\Gamma }}_d\left(\theta \right){\delta }_s\left({\sigma }_{\mathrm{dB}}\right)sin\left(\varphi \right)\mathrm{d}\theta \mathrm{d}\varphi .$

The function ${\delta }_s$ is a differentiable threshold function built to cut the low-intensity variations that we note in Figure 9, which have a negative impact on the computation of the gradient of the cost function:

(16)${\delta }_s\left(x\right)={\displaystyle \frac{x-s}{1+e^{-\tau (x-s)}}},$

where $\tau$ is a parameter that allows us to control the sharpness of the peak. The threshold $s=-5\mathrm{dB}$ was chosen.

The function ${\mathsf{\Gamma }}_d$ is a weighting function that we chose in order to favor the peak directions. It is close to 1 everywhere, except in the neighborhood of the direction d where it decreases quickly:

(17)${\mathsf{\Gamma }}_d\left(\theta \right)=1-{\displaystyle \frac{1}{2}}\sum _{i=0}^4{\displaystyle \frac{6-i}{10}}exp\left(-{\left({\displaystyle \frac{\theta -d}{{10}^i\tau }}\right)}^2\right).$

Graphs of the function ${\mathsf{\Gamma }}_d$ are shown in Figure 11 for various values of $\tau$. For this study, we chose, somewhat arbitrarily, a peak direction at a bearing of $d=30°$, and $\tau =0.1$.

5.3. Constraints

In order to ensure that the reduced drag coefficient does not come at the cost of a lower lift, we constrain the optimization to work at a constant lift coefficient (${\mathrm{C}}_{\mathrm{L}}$).

Furthermore, we want to preserve the reference surface (${\mathrm{S}}_{\mathrm{ref}}$) of the wing. ${\mathrm{S}}_{\mathrm{ref}}$ is defined in Figure 12.

6. Sensitivity Validation

6.1. RCS with Respect to Geometric Variables

We validated the RCS gradient with respect to the 13 shape design variables by comparing its finite differences (FDs) for various step sizes ranging from 1 to ${10}^{-12}$ in the unit of each individual variable. For the sake of simplicity, each FD gradient is computed with the same step size for each design variable. In practice, of course, we should fine tune a different step size for each individual variable in order to get a more precise gradient.

If the gradient is properly computed, we expect that, if the step size is small enough, the error decreases with the step size, up until a point where the numerical error on the function is of the same order of magnitude as the difference between the values of the RCS function at $\nu$ and $\nu +ϵ\widehat{d}$ (where $\widehat{d}$ is an arbitrary direction unit vector). After this critical point, the error should increase roughly linearly as the step size continues decreasing.

Let f denote a real-valued function that is at least three times differentiable, $f^{\prime }$ its derivative, and ${\mathsf{\Delta }}_{ϵ}f$ and ${\delta }_{ϵ}f$ its derivative approximated with forward and centered FDs, respectively, with a step size $ϵ$. It is well known that

(18)$\begin{array}{cc}\hfill {\mathsf{\Delta }}_{ϵ}f-f^{\prime }=o\left(ϵ\right)& \mathrm{as}ϵ\to 0,\hfill \end{array}$

(19)$\begin{array}{cc}\hfill {\delta }_{ϵ}f-f^{\prime }=o\left({ϵ}^2\right)& \mathrm{as}ϵ\to 0.\hfill \end{array}$

Hence, the logarithm of the error between the exact gradient and the FD should be linear with respect to the step size if the step size is small enough but bigger than the critical step size, and its slope should be 1 for the forward FD and 2 for the centered FD. We consider that the minimal error should be at most ${10}^{-3}$ for the forward FD and ${10}^{-6}$ for the centered FD in order to deem the gradient acceptable.

Figure 13 shows the norm of the difference (the solid line) between the gradient computed by COPOLA and the FD gradients function of the step size for both forward (in red) and centered (in green) differences, as well as the slope of the logarithm of the error (the dashed line).

The forward FD behaves very well: it is linear with a slope of 1 for a wide range of step sizes, and the minimal error is a few ${10}^{-4}$. The centered difference is also acceptable with a minimal error of around ${10}^{-6}$ and a slope of 2 right before the critical step size.

6.2. ${\mathrm{C}}_{\mathrm{D}}$ with Respect to Geometric and Aerodynamic Variables

Figure 14 shows the error between the adjoint gradient and the finite differences, both centered and forward, for the optimal step size determined for the CAD gradient. When the value of the adjoint gradient is high enough (more than ${10}^{-3}$), the relative error is used:

${\displaystyle \frac{|\delta F_a/\delta \nu -\partial F_a/\partial \nu |}{\partial F_a/\partial \nu }}.$

When the gradient is too small, the relative error may be misleading, so the absolute error is used instead:

$|\delta F_a/\delta \nu -\partial F_a/\partial \nu |.$

Most of the points show an accurate gradient, with an error of at most ${10}^{-3}$. The only exceptions are the gradients with respect to the leading edge camber of the 3rd and 5th sections, but the objective function is less sensitive to these variables. In these cases, the forward and central finite differences are not in good agreement with each other, which indicates that the step size imposed by the CAD gradient is not small enough to compute an accurate aerodynamic gradient. In future developments, the new features of the CAD modeler will be differentiated analytically, and the inaccuracy due to the CAD finite differences will not plague us anymore.

A more accurate validation of the aerodynamic gradient would be achievable with a methodology similar to that used in Section 6.1 in order to free ourselves from the step size determined by the CAD gradient, but this was not feasible for the CFD code.

7. Optimization Results

In this section, we present the results of the various optimizations we performed. We analyze the qualities of both the optimized aerodynamic and stealth shapes with respect to the baseline and each other. We then draw a Pareto front to illustrate the various trade-offs between stealth and aerodynamics. Furthermore, we present two shapes corresponding to two trade-offs in particular.

Figure 15 shows the optimization history of those four optimizations, with the evolution of the objective function and the normalized constraints with the iterations. The aerodynamic optimization (Figure 15a) needs 18 iterations to converge, meaning 18 gradient evaluations and 22 objective function evaluations, which is standard for this type of optimization. Because the objective varies a lot with a small shape variation, the stealth optimization (Figure 15b) needs more than double that number, with 40 iterations or 40 gradient evaluations and 56 objective function evaluations. This is still an acceptable performance for a gradient-based process. The trade-offs shape fall in between, as could be expected. The $w=0.8$ shape (Figure 15c) took 16 gradient evaluations and 31 objective function evaluations (16 iterations) to converge while the $w=0.75$ shape (Figure 15d) needed 21 gradient evaluations and 48 objective function evaluations (21 iterations).

In each case, the constraints are satisfied to a precision that is an order of magnitude better than the constraint that we had set (${10}^{-4}$).

7.1. Aerodynamic Optimization

First, we performed an optimization with a purely aerodynamic objective function ($w=1$). The planform of the optimal shape compared to that of the baseline is presented in Figure 16. The corresponding values of the design variables can be found in Table 1. Most notably, the leading edge sweep angles reach their respective upper-bounds, leading to an unrealistic shape. As discussed in Section 5.1, this was expected because of a simplistic criterion, relying only on the drag coefficient rather than dealing with any lateral stability criterion. Since the simulation is supersonic, the main contributor to the drag is the wave drag due to the shock waves developing on the nose of the aircraft and the wing leading edge. Increasing the leading edge sweep angle decreases the shock strength and the wave drag.

Figure 17 and Figure 18 show a comparison of the fields of local Mach number around the cross-section of the wing of the baseline (top) and the aerodynamic optimal shape (bottom) at $Y/b=0.45$ and $0.8$, respectively. The decreasing strength of the bow shockwave can be clearly seen at both positions.

The reduction in wave drag owing to the weakening of the bow shock leads to a reduction of 40 drag counts, a little more than 10% of the baseline. The lift coefficient stays the same, at 0.0113, as per the constraint.

Figure 19a shows the integration of the total drag coefficients along the X axis for different shapes. The fuselage is unmodified by the optimization process, so no difference can be seen for $X<X_{\mathrm{LE}}^{\mathrm{root}}$, but the drag of the optimized aerodynamic shape (red line) starts increasing later and more gradually than for the other shapes, which translates to the gain in wave drag because of the higher leading sweep angles. Figure 19b shows the distribution of the wing lift coefficient with respect to the normalized Y coordinate for different shapes. The lift is better distributed along the wing.

On the other hand, the multiple values of leading and trailing edge sweep angles have an adverse effect on the stealth criterion. This is because this leads to multiple RCS peaks as well as a wider spread of the RCS due to wedge diffraction, where the stealth criterion privileges a few, very sharp peaks.

In Figure 20, three peaks are visible. The very wide peak around $5°$ corresponds to the diffraction on the first section of the trailing edge wedge, and the peak around $36°$ corresponds to the second section of the trailing edge. The more intense peak around $58°$ corresponds to the specular reflection on the first section of the leading edge. The peak corresponding to the specular reflection on the second section of the leading would be in the $63°$ direction and is not visible in this figure. This yields a value of the stealth criterion more than twice that of the baseline.

Both the drag coefficient and the stealth criterion value can be read in Table 2 along with other values for various shapes.

7.2. Stealth Optimization

We then performed an optimization with a purely stealth cost function ($w=0$). The planform of the optimal shape compared to that of the baseline is presented in Figure 21. The corresponding values of the design variables can be found in Table 1. The leading and trailing edges are aligned with a sweep angle of $\pm 30°$, which is the result expected from the choice of the stealth criterion. This allows for the RCS to be concentrated in a sharp peak around a bearing angle of $30°$, as seen in Figure 22. As a result, the stealth criterion value for the optimized shape is almost 50 times lower than that of the baseline (see Table 2).

However, the decrease in leading edge sweep angle is very detrimental to the aerodynamic drag and particularly to the wave drag: the drag count is 55 points greater than that of the baseline (see Table 2). This can be seen clearly in Figure 23 and Figure 24, where a strong shock appears just upstream of the leading edge, which was absent from the other shapes. This is confirmed in Figure 19a, which shows a much steeper and earlier increase in ${\mathrm{C}}_{\mathrm{D}}$ than for the other shapes just downstream of the leading edge position. The pressure drag of the stealth-optimal shape alone is greater than the total drag on most of the baseline wing and on the totality of the aerodynamic optimal shape. The contribution of the pressure drag to the total drag is also much more important for the stealth-optimal shape than the other shapes.

In Figure 19b, we can see that the lift is more concentrated next to the root for the stealth-optimal shape than for the other shapes. On a realistic aircraft, such a behavior would be detrimental to stalling resistance at great AoA maneuvers, for example.

7.3. Pareto Optimality

In the two previous sections, we show that the requirements for a good aerodynamic shape and a good stealth shape are mainly conflicting. Therefore, the need for compromise arises. We explore different combinations and draw a Pareto front in Figure 25. We present a few selected shapes. The values of the criteria and constraints for these shapes are summarized in Table 2.

Note that the Pareto front features three main attraction points, around $w=0.9$, $w=0.75$, and $w\in [0,0.5]$. These attractors are probably due to the steepness of the weighting used in the RCS cost function (see Figure 11). The shapes around $w=0.5$ and below are essentially the same as the purely optimized stealth shape (Figure 21).

The shape corresponding to $w=0.75$ and $w=0.8$ is showcased in Figure 26 and Figure 27, respectively.

The drag difference between these shapes is mild (8 drag count) and is largely due to the sweep angle as can be seen in the difference in pressure drag in Figure 19a. This can be explained by the variation in intensity of the bow shock at the leading of the wings, because of the different sweep angles (see Figure 28, Figure 29, Figure 30 and Figure 31). The difference in RCS cost function is more notable. Although the RCS cost function is lower for the $w=0.75$ shape than the $w=0.8$, its average RCS is higher. This is because, on the $w=0.8$, the trailing edge exhibits two different sweep angles, thus creating a secondary diffraction peak at $20°$ (see Figure 32), while there is only one trailing edge sweep angle on the $w=0.75$ shape and thus only one diffraction peak at $30°$ (see Figure 33). These two shapes are also better than the baseline both in terms of aerodynamic drag (7 and 15 drag counts lower than the baseline respectively) and RCS function, though the $w=0.8$ shape is barely better than the baseline (see Table 2).

8. Conclusions

In this paper, we present a constrained, gradient-based aero-stealth optimization framework in order to determine an optimal shape for an aircraft. For the sake of keeping computation time low, we showcased our method on a simplified shape with fourteen parameters, but it is able to deal with a more complex shape parameterized with hundreds of design variables.

We computed the aerodynamic gradient through an adjoint formulation and the RCS gradient with the help of the AD tool TAPENADE. We validated these gradients by comparing them to finite differences. The results showed an acceptable difference.

After developing and validating the PO solver, we optimized a simplified aircraft geometry, with a weighted sum of the aerodynamic drag and a criterion based on the RCS distribution as the cost function. We varied the weights in order to identify a Pareto front. We extracted four of those optimal shapes and analyzed them.

The algorithm optimized the RCS mainly by tweaking the leading edge sweep angle and the chord length, which had an impact on the trailing edge sweep angle. This is because the stealth criterion depended heavily on the angular distribution of the RCS rather than merely its intensity.

When the weight is 1, the cost function is purely the aerodynamic drag. The optimal shape features a 40 drag count decrease but at the cost of spreading the RCS across the observation window, which is detrimental to the stealth criterion. On the other hand, when the weight is 0, the cost function is purely the stealth criterion. The algorithm is able to improve the RCS of the shape across most of the observation window by concentrating the scattered wave in a single direction, but this has a dramatic impact on the wave drag of the aircraft.

With weights between 0.5 and 0.9, the optimizer was able to find trade-offs that were strictly better than the baseline shape.

The framework allowed for efficient optimization, each one converging in less than 100 objective function and gradient evaluations and taking less than 10,000 CPU h to complete, which is consistent with the standards of industrial aerodynamic design. Furthermore, the use of a CAD facilitates interdisciplinary communication and is easy to set up by aerodynamic engineers, even when they are not shape parameterization experts. This makes it very useful in an industrial context.

Going forward, a more advanced method to identify the Pareto front, such as the one proposed by Desideri [34], would be beneficial.

Author Contributions

Conceptualization, C.T., G.R. and O.P.; Methodology, C.T. and G.R.; Software, C.T.; Validation, C.T.; Writing—original draft preparation, C.T.; Writing—review, G.R. and O.P.; Supervision, O.P. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

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

Acknowledgments

Our thanks to Quentin Carayol, Hervé Steve, and Nicolas Réau for their help in modifying COPOLA; to Steven Kleinveld for his help with GANIMEDE; and to Laurent Hascoët from INRIA for his help with TAPENADE. Finally, we are grateful to Sarah Julisson for maintaining and adapting the optimization system of Dassault-Aviation to our needs.

Conflicts of Interest

Authors Charles Thoulon and Gilbert Roge were employed by the company Dassault-Aviation. The remaining author declares that the research was conducted in the absence of any commercial or financial relationship that could be construed as a potential conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript

AD	Automatic Differentiation
AETHER	In-house CFD solver
CAD	Computer-Assisted Design
CFD	Computational Fluid Dynamics
COPOLA	In-house differentiated PO and PTD solver
GANIMEDE	In-house CAD modeler
GFA	Generic Fighter Aircraft (baseline geometry)
PO	Physical Optics
PTD	Physical Theory of Diffraction
RCS	Radar Cross-Section
TAPENADE	Automatic Differentiation tool developed by INRIA
DSM	Design Structure Matrix

Footnotes

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

Figures and Tables

Figure 1. DSM representation of the optimization process.

Figure 2. Definition of the bearing ([Forumla omitted. See PDF.]) and elevation ([Forumla omitted. See PDF.]) angles.

Figure 3. The RCS of the sphere of radius [Forumla omitted. See PDF.]. (a) The RCS function of the bearing angle. [Forumla omitted. See PDF.]. (b) Dimensionless RCS function of the radius-to-wavelength ratio.

Figure 4. Comparison of the RCS function of the bearing angle predicted by the COPOLA and BEM codes on a GFA geometry (see Figure 5).

Figure 5. The Generic Fighter Aircraft. (a) 3D-view. (b) Bottom view with the leading and trailing edge sweep angles.

Figure 6. Distribution of local Mach number around the fuselage of the GFA in the symmetry plane.

Figure 7. Distribution of local Mach number around the cross-section of the wing of the GFA at [Forumla omitted. See PDF.].

Figure 8. Distribution of local Mach number around the cross-section of the wing of the GFA at [Forumla omitted. See PDF.].

Figure 9. Profile of RCS with respect to the bearing angle on the GFA (elevation: [Forumla omitted. See PDF.]). (a) Purely PO simulation with horizontal polarization. (b) Purely PO simulation with vertical polarization. (c) PO simulation with PTD extension with horizontal polarization. (d) PO simulation with PTD extension with vertical polarization.

Figure 10. The wing parameterized with the control sections (red lines) and control points (red dots).

Figure 11. The function [Forumla omitted. See PDF.] with [Forumla omitted. See PDF.] for various values of [Forumla omitted. See PDF.].

Figure 12. Definition of [Forumla omitted. See PDF.].

Figure 13. Error between the gradient of the stealth cost function and the same computed by the FD (both forward and central) function of the FD step size ([Forumla omitted. See PDF.]).

Figure 14. Comparison between the value of the gradient of the aerodynamic drag with respect to the geometric variables [Forumla omitted. See PDF.] as computed by the adjoint formulation and using the forward (Fwd) and centered (Ctd) FDs for the optimal step size determined for the CAD gradient. The variable indices and step sizes are summarized in Table 1.

Figure 15. Optimization history. The cost is read on the left axis. The normalized constraints are read on the right axis. The tolerance for the constraints is [Forumla omitted. See PDF.]. (a) Aerodynamic optimization history. (b) Stealth optimization history. (c) [Forumla omitted. See PDF.] shape optimization history. (d) [Forumla omitted. See PDF.] shape optimization history.

Figure 16. Optimal shape from a purely aerodynamic standpoint.

Figure 17. Comparison of the distributions of local Mach number around the cross-section of the wing profile at [Forumla omitted. See PDF.] between the baseline (top) and the optimized aerodynamic shape (bottom).

Figure 18. Comparison of the distributions of local Mach number around the cross-section of the wing profile at [Forumla omitted. See PDF.] between the baseline (top) and the optimized aerodynamic shape (bottom).

Figure 19. Detailed analysis of the distribution of aerodynamic coefficients along their respective axes.

Figure 20. RCS of the optimized aerodynamic shape compared to that of the baseline.

Figure 21. Optimal shape from a purely low-observability standpoint.

Figure 22. RCS of the stealth-optimal shape compared to that of the baseline.

Figure 23. Comparison of the distributions of the local Mach number around the cross-section of the wing profile at [Forumla omitted. See PDF.] between the baseline (top) and the stealth-optimal shape (bottom).

Figure 24. Comparison of the distributions of the local Mach number around the cross-section of the wing profile at [Forumla omitted. See PDF.] between the baseline (top) and the stealth-optimal shape (bottom).

Figure 25. Pareto front.

Figure 26. Shape corresponding to the [Forumla omitted. See PDF.] point.

Figure 27. Shape corresponding to the [Forumla omitted. See PDF.] point.

Figure 28. Comparison of the distributions of local Mach number around the cross-section of the wing profile at [Forumla omitted. See PDF.] between the baseline (top) and the [Forumla omitted. See PDF.] shape (bottom).

Figure 29. Comparison of the distributions of local Mach number around the cross-section of the wing profile at [Forumla omitted. See PDF.] between the baseline (top) and the [Forumla omitted. See PDF.] shape (bottom).

Figure 30. Comparison of the distributions of the local Mach number around the cross-section of the wing profile at [Forumla omitted. See PDF.] between the baseline (top) and the [Forumla omitted. See PDF.] shape (bottom).

Figure 31. Comparison of the distributions of local Mach number around the cross-section of the wing profile at [Forumla omitted. See PDF.] between the baseline (top) and the [Forumla omitted. See PDF.] shape (bottom).

Figure 32. RCS of the [Forumla omitted. See PDF.] shape compared to that of the baseline.

Figure 33. RCS of the [Forumla omitted. See PDF.] shape compared to that of the baseline.

Table 1

Summary of the values of the design variables for a few chosen optimal shapes.

Design Variable	Unit	Section	Index	Step Size	Min	Max	Baseline	Stealth-Optimal Shape	Aero-Optimal Shape
Sweep angle	$°$	2	1	0.05	28	58	48	30	58
		4	2	0.05	−10	10	0	0	10
Chord length	mm	2	3	${10}^{-3}$	5991	7863	7489	7059	7863
		4	4	${10}^{-3}$	2643	4845	4405	1211	4845
		6	5	${10}^{-3}$	943	1238	1179	1123	1238
Span	mm	2–6	6	${10}^{-3}$	4864	6675	6230	6674	5737
Twist angle	$°$	2	7	0.05	0	2	0	−0.5	−0.3
		3	8	0.05	−2	2	0	−0.8	−0.6
		5	9	0.05	−2	2	0	−0.6	−1.8
		6	10	0.05	−2	2	0	−0.8	−1.9
Leading edge camber	%	3	11	${10}^{-3}$	−1	1	0	−1	−0.69
	chord	5	12	${10}^{-3}$	−1	1	0	0	0.64
	length	6	13	${10}^{-3}$	−1	1	0	1	0.63
AoA	$°$	–	14	0.01	1	3	2	2.3	2.7

Table 2

Summary of the cost and constraint values for various shapes in the Pareto front.

w	${\mathbf{C}}_{\mathbf{D}}$	${\mathbf{C}}_{\mathbf{L}}$	${\mathit{F}}_{\mathbf{RCS}}$	$\widehat{\mathit{\sigma }}$	${\mathbf{S}}_{\mathbf{ref}}$	N Cycle
1.0	346.97	$112.5·{10}^{-3}$	$186.57·{10}^{-3}$	$8.4·{10}^{-3}$	68.6	41
0.995	352.16	$112.5·{10}^{-3}$	$144.41·{10}^{-3}$	$6.9·{10}^{-3}$	68.6	63
0.99	354.34	$112.5·{10}^{-3}$	$104.80·{10}^{-3}$	$6.0·{10}^{-3}$	68.6	62
0.95	365.58	$112.5·{10}^{-3}$	$94.15·{10}^{-3}$	$6.3·{10}^{-3}$	68.6	40
0.9	367.49	$112.5·{10}^{-3}$	$84.99·{10}^{-3}$	$7.8·{10}^{-3}$	68.6	59
0.8	372.88	$112.5·{10}^{-3}$	$78.02·{10}^{-3}$	$25.3·{10}^{-3}$	68.6	47
0.75	380.85	$112.5·{10}^{-3}$	$40.75·{10}^{-3}$	$81.6·{10}^{-3}$	68.6	69
0.52	381.72	$112.5·{10}^{-3}$	$40.31·{10}^{-3}$	$84.4·{10}^{-3}$	68.6	65
0.51	426.99	$112.5·{10}^{-3}$	$2.24·{10}^{-3}$	$20.9·{10}^{-3}$	68.6	91
0.5	427.22	$112.5·{10}^{-3}$	$2.23·{10}^{-3}$	$21.0·{10}^{-3}$	68.6	98
0.0	432.54	$112.5·{10}^{-3}$	$1.74·{10}^{-3}$	$23.3·{10}^{-3}$	68.6	96
baseline	387.46	$112.6·{10}^{-3}$	$82.27·{10}^{-3}$	$7.4·{10}^{-3}$	68.6	–

References

1. Knott, E.F.; Schaeffer, J.F.; Tulley, M.T. Radar Cross Section; SciTech Publishing: Raleigh, NC, USA, 2004.

2. Aliaga, C.; Kopp, M.; Salui, K.B.; Mahapatra, D.; Sardesai, A. Aerodynamic and stealth studies of canard-wing Configurations at Transonic Speeds using Ansys Fluent & Ansys HFSS SBR+. Proceedings of the AIAA AVIATION 2022 Forum; Chicago, IL, USA, 27 June–1 July 2022; 3322.

3. Lee, D.; Gonzalez, L.F.; Srinivas, K.; Auld, D.; Periaux, J. Multi-objective/multidisciplinary design optimisation of blended wing body UAV via advanced evolutionary algorithms. Proceedings of the 45th AIAA Aerospace Sciences Meeting and Exhibit; Reno, NV, USA, 8 January 2007–11 January 2007; 36.

4. Wu, D.; Long, T.; Li, Y.; Jiang, M.; Huang, B. Aero-structure-stealth coupled optimization for high aspect ratio wing using adaptive metamodeling method. Proceedings of the 15th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference; Atlanta, GA, USA, 16–20 June 2014; 2304.

5. Jiang, X.; Zhao, Q.; Zhao, G.; Meng, C. Numerical Optimization on Aerodynamic/Stealth Characteristics of Airfoil Based on CFD/CEM Coupling Method. Trans. Nanjing Univ. Aeronaut. Astronaut.; 2016; 33, pp. 274-284.

6. Parr, J.; Holden, C.M.; Forrester, A.I.; Keane, A.J. Review of efficient surrogate infill sampling criteria with constraint handling. Proceedings of the 2nd International Conference on Engineering Optimization; Lisbon, Portugal, 6–10 September 2011.

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

8. Martin, L. Conception Aérodynamique Robuste. Ph.D. Thesis; Université de Toulouse, Université Toulouse III-Paul Sabatier: Toulouse, France, 2010; Available online: https://core.ac.uk/download/pdf/12095692.pdf (accessed on 29 July 2024).

9. Hascoët, L.; Pascual, V. The Tapenade Automatic Differentiation tool: Principles, Model, and Specification. ACM Trans. Math. Softw.; 2013; 39, pp. 1-43. Available online: https://gitlab.inria.fr/tapenade/tapenade (accessed on 29 July 2024). [DOI: https://dx.doi.org/10.1145/2450153.2450158]

10. Lyu, Z.; Martins, J.R. Aerodynamic design optimization studies of a blended-wing-body aircraft. J. Aircr.; 2014; 51, pp. 1604-1617. [DOI: https://dx.doi.org/10.2514/1.C032491]

11. Lyu, Z.; Kenway, G.K.; Martins, J.R. Aerodynamic shape optimization investigations of the common research model wing benchmark. AIAA J.; 2015; 53, pp. 968-985. [DOI: https://dx.doi.org/10.2514/1.J053318]

12. Bons, N.P.; He, X.; Mader, C.A.; Martins, J.R. Multimodality in aerodynamic wing design optimization. AIAA J.; 2019; 57, pp. 1004-1018. [DOI: https://dx.doi.org/10.2514/1.J057294]

13. Garnier, E.; Langlet, S.; Klotz, P.; Cadillon, J.; Simon, J.; Castelli, J.C.; Levadoux, D. Surrogate-based conception for aerodynamic/stealth compromise. Proceedings of the AVT-324 Specialists’ Meeting on Multidisciplinary Design Approaches and Performance Assessment of Future Combat Aircraft; virtuel, Italy, 28–30 September 2020.

14. Guerra, J. Optimisation Multi-Objectif Sous Incertitudes de Phénomènes de Thermique Transitoire. Ph.D. Thesis; Institut Superieur de L’aeronautique et de L’espace (ISAE): Toulouse, France, 2016.

15. Liu, Z.; Song, W.; Han, Z.; Wang, Y. Multifidelity Aerodynamic/Stealth Design Optimization Method for Flying Wing Aircraft. Proceedings of the 32th Congress of the International Council of the Aeronautical Sciences, ICAS; Shanghai, China, 6–10 September 2021.

16. Zhou, Z.; Huang, J. Quantitative weight and two-particle search algorithm to optimize aero-stealth performance of a backward inclined vertical tail. Aerospace; 2023; 10, 345. [DOI: https://dx.doi.org/10.3390/aerospace10040345]

17. Taj, Z.U.D.; Bilal, A.; Awais, M.; Salamat, S.; Abbas, M.; Maqsood, A. Design exploration and optimization of aerodynamics and radar cross section for a fighter aircraft. Aerosp. Sci. Technol.; 2023; 133, 108114. [DOI: https://dx.doi.org/10.1016/j.ast.2023.108114]

18. Pan, Y.; Huang, J.; Li, F.; Yan, C. Integrated design optimization of aerodynamic and stealthy performance for flying wing aircraft. Proceedings of the International MultiConference of Engineers and Computer Scientists; Hong Kong, 15–17 March 2017; Volume 2, pp. 15-17.

19. Li, M.; Bai, J.; Li, L.; Meng, X.; Liu, Q.; Chen, B. A gradient-based aero-stealth optimization design method for flying wing aircraft. Aerosp. Sci. Technol.; 2019; 92, pp. 156-169. [DOI: https://dx.doi.org/10.1016/j.ast.2019.05.067]

20. Kleinveld, S.; Rogé, G.; Daumas, L.; Dinh, Q. Differentiated parametric CAD used within the context of automatic aerodynamic design optimization. Proceedings of the 12th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference; Victoria, BC, Canada, 10–12 September 2008; 5952.

21. Chalot, F. Industrial Aerodynamics. Encyclopedia of Computational Mechanics; Stein, E.; de Borst, R.; Hughes, T.J.R. John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2004; Volume 3: Fluids, Chapter 12.

22. Barclay, A. SQP Methods for Large-Scale Optimization; University of California: San Diego, CA, USA, 1999.

23. Laval, D. Application de Méthodes Asymptotiques à la Simulation de la Diffraction électromagnétique par un Corps Régulier. Ph.D. Thesis; Université Sciences et Technologies-Bordeaux I: Bordeaux, France, 2006; Available online: https://theses.hal.science/tel-00197165/document (accessed on 29 July 2024).

24. Stratton, J.A.; Chu, L. Diffraction theory of electromagnetic waves. Phys. Rev.; 1939; 56, 99. [DOI: https://dx.doi.org/10.1103/PhysRev.56.99]

25. Michaeli, A. Elimination of infinities in equivalent edge currents, Part I: Fringe current components. IEEE Trans. Antennas Propag.; 1986; 34, pp. 912-918. [DOI: https://dx.doi.org/10.1109/TAP.1986.1143913]

26. Chalot, F.; Johan, Z.; Mallet, M.; Billard, F.; Martin, L.; Barré, S. Extension of methods based on the stabilized finite element formulation for the solution of the Navier–Stokes equations and application to aerodynamic design. Comput. Methods Appl. Mech. Eng.; 2023; 417, 116425. [DOI: https://dx.doi.org/10.1016/j.cma.2023.116425]

27. Alfonsi, G. Reynolds-averaged Navier–Stokes equations for turbulence modeling. Appl. Mech. Rev.; 2009; 62, 040802. [DOI: https://dx.doi.org/10.1115/1.3124648]

28. Spalart, P.R.; Allmaras, S.R. One-equation turbulence model for aerodynamic flows. La Rech. Aerosp.; 1994; 1, pp. 5-21.

29. Martin, L.; Roge, G. Calcul de la sensibilité d’ordre deux d’une observation aérodynamique. Proceedings of the ESAIM: Proceedings; EDP Sciences: Les Ulis, France, 2009; Volume 27, pp. 138-155.

30. Poggio, A.J.; Miller, E.K. Integral Equation Solutions of Three-Dimensional Scattering Problems; MB Assoc.: Winnipeg, MB, Canada, 1970; Chapter 4 pp. 195-264.

31. Carayol, Q. Développement et Analyse D’une Méthode Multipôle Multiniveau Pour L’électromagnétisme. Ph.D. Thesis; Université Paris 6: Paris, France, 2002; Available online: https://theses.fr/2002PA066485 (accessed on 19 May 2024).

32. Chalot, F.; Mallet, M.; Ravachol, M. A comprehensive finite element Navier-Stokes solver for low-and high-speed aircraft design. Proceedings of the 32nd Aerospace Sciences Meeting and Exhibit; Reno, NV, USA, 10–13 January 1994; 814.

33. Jenn, D. Radar and Laser cross Section Engineering; American Institute of Aeronautics and Astronautics, Inc.: Reno, NV, USA, 2005.

34. Désidéri, J.A. MGDA Variants for Multi-Objective Optimization; Institut National de Recherche en Informatique et en Automatique: Sophia-Antipolis, France, 2012.