Next Article in Journal
Numerical Study of Laminar Flow and Vortex-Induced Vibration on Cylinder Subjects to Free and Forced Oscillation at Low Reynolds Numbers
Next Article in Special Issue
Bridging Large Eddy Simulation and Reduced-Order Modeling of Convection-Dominated Flows through Spatial Filtering: Review and Perspectives
Previous Article in Journal
Predicting Wall Pressure in Shock Wave/Boundary Layer Interactions with Convolutional Neural Networks
Previous Article in Special Issue
Numerical Simulations of Scalar Transport on Rough Surfaces
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Gradient-Based Aero-Stealth Optimization of a Simplified Aircraft

by
Charles Thoulon
1,2,*,
Gilbert Roge
1 and
Olivier Pironneau
2
1
Dassault-Aviation, 92210 St. Cloud, France
2
Laboratoire Jacques-Louis Lions, Sorbonne Université, 75005 Paris, France
*
Author to whom correspondence should be addressed.
Fluids 2024, 9(8), 174; https://doi.org/10.3390/fluids9080174
Submission received: 20 May 2024 / Revised: 30 June 2024 / Accepted: 25 July 2024 / Published: 30 July 2024
(This article belongs to the Special Issue Recent Advances in Fluid Mechanics: Feature Papers, 2024)

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   GHz to 12   GHz , corresponding to wavelengths between 0.025   m and 0.3   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) σ defined as follows:
σ ( θ , ϕ ) = lim R R 2 | E P | 2 | E i | 2 ,
where E P is the scattered electric fields at a point P ( R , θ , ϕ ) in polar coordinates (defined in Figure 2) and E i is the incident electric field at the target, whose norm is considered to be 1 V / m in practice. Usually, the RCS is expressed in decibels:
σ dB = 10 log 10 ( σ ) .
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 E P at a point P can be expressed as follows:
E P = j k Z 0 4 π e j k R R S u s × ( u s × J ) exp j k R · ( u i u s ) ,
with j = 1 , k the wave number, Z 0 the impedance of free space, S the surface of the target, u i and u s the unit direction vectors for the incident and scattered fields, respectively, R = O P , and R = | R | , and where J is the surface current density vector. In the PO approximation,
J = 2 Z 0 n × u i × E i in the illuminated zone , 0 in the shadow zone ,
where n is the local normal unit vector to the surface.
Note that by (2),
E P = O 1 R as R .
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 x e , the coordinates of the vertices of the surface grid used for the electromagnetism problem.
If
F e ( ν ) = F ^ e ( σ ( x e ( ν ) ) )
is a function depending explicitly only on the RCS σ , then the directional derivative of F e with respect to a set of geometric variables ν (e.g., sweep angle, chord length) is
d F e = F ^ e ( σ ) · σ ( x e ) · x e ( ν ) · d ν .
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 d F e / d σ and d σ / d 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 d ν 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 F e · d ν = d ν · x e ( ν ) · σ ( x e ) · F ^ e ( σ ) . When F e is scalar-valued and the dimension of ν is large, and assuming that the storage capacity is available, it is theoretically more efficient to compute the Jacobian matrix F e of F e row by row and then compute d F e = F e d ν 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 d x e / d ν , 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 d F e / d ν 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, ε ), (k, ω ), (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
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 d F a d ξ and d F a d ν with respect to a set of aerodynamic variables ξ (e.g., the angle of attack) and geometric variables ν (e.g., the wing sweep angle), respectively. The variables ν and ξ 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 x a is subjected to a deformation owing to the CAD modeler, the deformation is propagated to the 3D grid of coordinates X a with the help of an elliptic operator described by the system:
L ( X a , x a ) = 0 .
Let us write F a and R as functions of ξ , X a , and W. Note that ξ is an independent variable but X a depends on ν and W depends on both ξ and ν .

2.4.1. Gradient with Respect to Aerodynamic Variables

First, the derivative of F a and R with respect to ξ can be written as
d F a = F a ξ W Ψ T R ξ W d ξ .
where Ψ is the solution to the aerodynamic adjoint system:
R W ξ T Ψ = F a W ξ T .

2.4.2. Gradient with Respect to Geometric Variables

The derivative of R and F a with respect to ν is
d F a = F a X a W Ψ T R X a W X a ν d ν ,
where Ψ is the same as defined in Equation (8).
In this case, it is not enough, though, as X a ν is expensive to compute. In order to find a more convenient expression, we must write an adjoint problem of the system, which finally yields
d F a = Φ T L x a X a x a ν d ν ,
where Φ is the solution of the mesh adjoint equation system:
L X a x a T Φ = F a X a x a T R X a x a T Ψ .

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 λ , then the RCS σ converges to the section area of the sphere as the wavelength decreases:
lim λ 0 σ = π 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   M Hz and 20   G 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 ( 2 π f R 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 ( 2 π f R 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 d x d ν 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 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 m , the frequency should be higher than 25 MHz in order to meet the criterion of Section 3:
50 < 2 π f L r e f 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 GHz and spanning 100 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 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:
F = w F a + ( 1 w ) F e ,
with 0 w 1 .

5.1. Aerodynamic Criterion

For the aerodynamic optimization, we minimize the drag coefficient ( C D ):
F a = C D C D r e f .
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 ( θ 0 , θ 1 ) × ( ϕ 0 , ϕ 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:
F e = σ ˜ σ ˜ r e f ,
where σ ˜ is a weighed average of a function of σ dB :
σ ˜ = θ 0 θ 1 ϕ 0 ϕ 1 Γ d ( θ ) δ s ( σ dB ) sin ( ϕ ) d θ d ϕ .
The function δ 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:
δ s ( x ) = x s 1 + e τ ( x s ) ,
where τ is a parameter that allows us to control the sharpness of the peak. The threshold s = 5 dB was chosen.
The function Γ 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:
Γ d ( θ ) = 1 1 2 i = 0 4 6 i 10 exp θ d 10 i τ 2 .
Graphs of the function Γ d are shown in Figure 11 for various values of τ . For this study, we chose, somewhat arbitrarily, a peak direction at a bearing of d = 30 ° , and τ = 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 ( C L ).
Furthermore, we want to preserve the reference surface ( S ref ) of the wing. S 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 ν and ν + ϵ d ^ (where 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 its derivative, and Δ ϵ f and δ ϵ f its derivative approximated with forward and centered FDs, respectively, with a step size ϵ . It is well known that
Δ ϵ f f = o ( ϵ ) as ϵ 0 ,
δ ϵ f f = o ϵ 2 as ϵ 0 .
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. C 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:
| δ F a / δ ν F a / ν | F a / ν .
When the gradient is too small, the relative error may be misleading, so the absolute error is used instead:
| δ F a / δ ν F a / ν | .
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 LE 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 ± 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 C 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 [ 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.

Funding

This research received no external funding.

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
ADAutomatic Differentiation
AETHERIn-house CFD solver
CADComputer-Assisted Design
CFDComputational Fluid Dynamics
COPOLAIn-house differentiated PO and PTD solver
GANIMEDEIn-house CAD modeler
GFAGeneric Fighter Aircraft (baseline geometry)
POPhysical Optics
PTDPhysical Theory of Diffraction
RCSRadar Cross-Section
TAPENADEAutomatic Differentiation tool developed by INRIA
DSMDesign Structure Matrix

References

  1. Knott, E.F.; Schaeffer, J.F.; Tulley, M.T. Radar Cross Section; SciTech Publishing: Raleigh, NC, USA, 2004. [Google Scholar]
  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+. In Proceedings of the AIAA AVIATION 2022 Forum, Chicago, IL, USA, 27 June–1 July 2022; p. 3322. [Google Scholar]
  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. In Proceedings of the 45th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 8 January 2007–11 January 2007; p. 36. [Google Scholar]
  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. In Proceedings of the 15th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Atlanta, GA, USA, 16–20 June 2014; p. 2304. [Google Scholar]
  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, 274–284. [Google Scholar]
  6. Parr, J.; Holden, C.M.; Forrester, A.I.; Keane, A.J. Review of efficient surrogate infill sampling criteria with constraint handling. In Proceedings of the 2nd International Conference on Engineering Optimization, Lisbon, Portugal, 6–10 September 2011. [Google Scholar]
  7. Lyu, Z.; Xu, Z.; Martins, J. Benchmarking optimization algorithms for wing aerodynamic design optimization. In Proceedings of the 8th International Conference on Computational Fluid Dynamics, Chengdu, China, 11–13 November 2014; Citeseer: Princeton, NJ, USA, 2014; Volume 11, p. 585. [Google Scholar]
  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, 1–43. Available online: https://gitlab.inria.fr/tapenade/tapenade (accessed on 29 July 2024). [CrossRef]
  10. Lyu, Z.; Martins, J.R. Aerodynamic design optimization studies of a blended-wing-body aircraft. J. Aircr. 2014, 51, 1604–1617. [Google Scholar] [CrossRef]
  11. Lyu, Z.; Kenway, G.K.; Martins, J.R. Aerodynamic shape optimization investigations of the common research model wing benchmark. AIAA J. 2015, 53, 968–985. [Google Scholar] [CrossRef]
  12. Bons, N.P.; He, X.; Mader, C.A.; Martins, J.R. Multimodality in aerodynamic wing design optimization. AIAA J. 2019, 57, 1004–1018. [Google Scholar] [CrossRef]
  13. Garnier, E.; Langlet, S.; Klotz, P.; Cadillon, J.; Simon, J.; Castelli, J.C.; Levadoux, D. Surrogate-based conception for aerodynamic/stealth compromise. In Proceedings of the AVT-324 Specialists’ Meeting on Multidisciplinary Design Approaches and Performance Assessment of Future Combat Aircraft, virtuel, Italy, 28–30 September 2020. [Google Scholar]
  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. [Google Scholar]
  15. Liu, Z.; Song, W.; Han, Z.; Wang, Y. Multifidelity Aerodynamic/Stealth Design Optimization Method for Flying Wing Aircraft. In Proceedings of the 32th Congress of the International Council of the Aeronautical Sciences, ICAS, Shanghai, China, 6–10 September 2021. [Google Scholar]
  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. [Google Scholar] [CrossRef]
  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. [Google Scholar] [CrossRef]
  18. Pan, Y.; Huang, J.; Li, F.; Yan, C. Integrated design optimization of aerodynamic and stealthy performance for flying wing aircraft. In Proceedings of the International MultiConference of Engineers and Computer Scientists, Hong Kong, 15–17 March 2017; Volume 2, pp. 15–17. [Google Scholar]
  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, 156–169. [Google Scholar] [CrossRef]
  20. 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; p. 5952. [Google Scholar]
  21. Chalot, F. Industrial Aerodynamics. In Encyclopedia of Computational Mechanics; Stein, E., de Borst, R., Hughes, T.J.R., Eds.; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2004; Volume 3: Fluids, Chapter 12. [Google Scholar]
  22. Barclay, A. SQP Methods for Large-Scale Optimization; University of California: San Diego, CA, USA, 1999. [Google Scholar]
  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. [Google Scholar] [CrossRef]
  25. Michaeli, A. Elimination of infinities in equivalent edge currents, Part I: Fringe current components. IEEE Trans. Antennas Propag. 1986, 34, 912–918. [Google Scholar] [CrossRef]
  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. [Google Scholar] [CrossRef]
  27. Alfonsi, G. Reynolds-averaged Navier–Stokes equations for turbulence modeling. Appl. Mech. Rev. 2009, 62, 040802. [Google Scholar] [CrossRef]
  28. Spalart, P.R.; Allmaras, S.R. One-equation turbulence model for aerodynamic flows. La Rech. Aerosp. 1994, 1, 5–21. [Google Scholar]
  29. Martin, L.; Roge, G. Calcul de la sensibilité d’ordre deux d’une observation aérodynamique. In Proceedings of the ESAIM: Proceedings; EDP Sciences: Les Ulis, France, 2009; Volume 27, pp. 138–155. [Google Scholar]
  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. [Google Scholar]
  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. In Proceedings of the 32nd Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 10–13 January 1994; p. 814. [Google Scholar]
  33. Jenn, D. Radar and Laser cross Section Engineering; American Institute of Aeronautics and Astronautics, Inc.: Reno, NV, USA, 2005. [Google Scholar]
  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. [Google Scholar]
Figure 1. DSM representation of the optimization process.
Figure 1. DSM representation of the optimization process.
Fluids 09 00174 g001
Figure 2. Definition of the bearing ( θ ) and elevation ( ϕ ) angles.
Figure 2. Definition of the bearing ( θ ) and elevation ( ϕ ) angles.
Fluids 09 00174 g002
Figure 3. The RCS of the sphere of radius R = 1   m . (a) The RCS function of the bearing angle. f = 5   G Hz . (b) Dimensionless RCS function of the radius-to-wavelength ratio.
Figure 3. The RCS of the sphere of radius R = 1   m . (a) The RCS function of the bearing angle. f = 5   G Hz . (b) Dimensionless RCS function of the radius-to-wavelength ratio.
Fluids 09 00174 g003
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 4. Comparison of the RCS function of the bearing angle predicted by the COPOLA and BEM codes on a GFA geometry (see Figure 5).
Fluids 09 00174 g004
Figure 5. The Generic Fighter Aircraft. (a) 3D-view. (b) Bottom view with the leading and trailing edge sweep angles.
Figure 5. The Generic Fighter Aircraft. (a) 3D-view. (b) Bottom view with the leading and trailing edge sweep angles.
Fluids 09 00174 g005
Figure 6. Distribution of local Mach number around the fuselage of the GFA in the symmetry plane.
Figure 6. Distribution of local Mach number around the fuselage of the GFA in the symmetry plane.
Fluids 09 00174 g006
Figure 7. Distribution of local Mach number around the cross-section of the wing of the GFA at Y / b = 0.45 .
Figure 7. Distribution of local Mach number around the cross-section of the wing of the GFA at Y / b = 0.45 .
Fluids 09 00174 g007
Figure 8. Distribution of local Mach number around the cross-section of the wing of the GFA at Y / b = 0.8 .
Figure 8. Distribution of local Mach number around the cross-section of the wing of the GFA at Y / b = 0.8 .
Fluids 09 00174 g008
Figure 9. Profile of RCS with respect to the bearing angle on the GFA (elevation: ϕ = 20 ° ). (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 9. Profile of RCS with respect to the bearing angle on the GFA (elevation: ϕ = 20 ° ). (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.
Fluids 09 00174 g009
Figure 10. The wing parameterized with the control sections (red lines) and control points (red dots).
Figure 10. The wing parameterized with the control sections (red lines) and control points (red dots).
Fluids 09 00174 g010
Figure 11. The function Γ d ( θ ) with d = 30 ° for various values of τ .
Figure 11. The function Γ d ( θ ) with d = 30 ° for various values of τ .
Fluids 09 00174 g011
Figure 12. Definition of S ref .
Figure 12. Definition of S ref .
Fluids 09 00174 g012
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 ( ϵ ).
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 ( ϵ ).
Fluids 09 00174 g013
Figure 14. Comparison between the value of the gradient of the aerodynamic drag with respect to the geometric variables ν 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 14. Comparison between the value of the gradient of the aerodynamic drag with respect to the geometric variables ν 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.
Fluids 09 00174 g014
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 10 4 . (a) Aerodynamic optimization history. (b) Stealth optimization history. (c) w = 0.8 shape optimization history. (d) w = 0.75 shape optimization history.
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 10 4 . (a) Aerodynamic optimization history. (b) Stealth optimization history. (c) w = 0.8 shape optimization history. (d) w = 0.75 shape optimization history.
Fluids 09 00174 g015
Figure 16. Optimal shape from a purely aerodynamic standpoint.
Figure 16. Optimal shape from a purely aerodynamic standpoint.
Fluids 09 00174 g016
Figure 17. Comparison of the distributions of local Mach number around the cross-section of the wing profile at Y / b = 0.45 between the baseline (top) and the optimized aerodynamic shape (bottom).
Figure 17. Comparison of the distributions of local Mach number around the cross-section of the wing profile at Y / b = 0.45 between the baseline (top) and the optimized aerodynamic shape (bottom).
Fluids 09 00174 g017
Figure 18. Comparison of the distributions of local Mach number around the cross-section of the wing profile at Y / b = 0.8 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 Y / b = 0.8 between the baseline (top) and the optimized aerodynamic shape (bottom).
Fluids 09 00174 g018
Figure 19. Detailed analysis of the distribution of aerodynamic coefficients along their respective axes.
Figure 19. Detailed analysis of the distribution of aerodynamic coefficients along their respective axes.
Fluids 09 00174 g019
Figure 20. RCS of the optimized aerodynamic shape compared to that of the baseline.
Figure 20. RCS of the optimized aerodynamic shape compared to that of the baseline.
Fluids 09 00174 g020
Figure 21. Optimal shape from a purely low-observability standpoint.
Figure 21. Optimal shape from a purely low-observability standpoint.
Fluids 09 00174 g021
Figure 22. RCS of the stealth-optimal shape compared to that of the baseline.
Figure 22. RCS of the stealth-optimal shape compared to that of the baseline.
Fluids 09 00174 g022
Figure 23. Comparison of the distributions of the local Mach number around the cross-section of the wing profile at Y / b = 0.45 between the baseline (top) and the stealth-optimal shape (bottom).
Figure 23. Comparison of the distributions of the local Mach number around the cross-section of the wing profile at Y / b = 0.45 between the baseline (top) and the stealth-optimal shape (bottom).
Fluids 09 00174 g023
Figure 24. Comparison of the distributions of the local Mach number around the cross-section of the wing profile at Y / b = 0.8 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 Y / b = 0.8 between the baseline (top) and the stealth-optimal shape (bottom).
Fluids 09 00174 g024
Figure 25. Pareto front.
Figure 25. Pareto front.
Fluids 09 00174 g025
Figure 26. Shape corresponding to the w = 0.8 point.
Figure 26. Shape corresponding to the w = 0.8 point.
Fluids 09 00174 g026
Figure 27. Shape corresponding to the w = 0.75 point.
Figure 27. Shape corresponding to the w = 0.75 point.
Fluids 09 00174 g027
Figure 28. Comparison of the distributions of local Mach number around the cross-section of the wing profile at Y / b = 0.45 between the baseline (top) and the w = 0.8 shape (bottom).
Figure 28. Comparison of the distributions of local Mach number around the cross-section of the wing profile at Y / b = 0.45 between the baseline (top) and the w = 0.8 shape (bottom).
Fluids 09 00174 g028
Figure 29. Comparison of the distributions of local Mach number around the cross-section of the wing profile at Y / b = 0.8 between the baseline (top) and the w = 0.8 shape (bottom).
Figure 29. Comparison of the distributions of local Mach number around the cross-section of the wing profile at Y / b = 0.8 between the baseline (top) and the w = 0.8 shape (bottom).
Fluids 09 00174 g029
Figure 30. Comparison of the distributions of the local Mach number around the cross-section of the wing profile at Y / b = 0.45 between the baseline (top) and the w = 0.75 shape (bottom).
Figure 30. Comparison of the distributions of the local Mach number around the cross-section of the wing profile at Y / b = 0.45 between the baseline (top) and the w = 0.75 shape (bottom).
Fluids 09 00174 g030
Figure 31. Comparison of the distributions of local Mach number around the cross-section of the wing profile at Y / b = 0.8 between the baseline (top) and the w = 0.75 shape (bottom).
Figure 31. Comparison of the distributions of local Mach number around the cross-section of the wing profile at Y / b = 0.8 between the baseline (top) and the w = 0.75 shape (bottom).
Fluids 09 00174 g031
Figure 32. RCS of the w = 0.8 shape compared to that of the baseline.
Figure 32. RCS of the w = 0.8 shape compared to that of the baseline.
Fluids 09 00174 g032
Figure 33. RCS of the w = 0.75 shape compared to that of the baseline.
Figure 33. RCS of the w = 0.75 shape compared to that of the baseline.
Fluids 09 00174 g033
Table 1. Summary of the values of the design variables for a few chosen optimal shapes.
Table 1. Summary of the values of the design variables for a few chosen optimal shapes.
Design VariableUnitSectionIndexStep SizeMinMaxBaselineStealth-Optimal ShapeAero-Optimal Shape
Sweep angle ° 210.052858483058
420.05−10100010
Chord lengthmm23 10 3 59917863748970597863
44 10 3 26434845440512114845
65 10 3 9431238117911231238
Spanmm2–66 10 3 48646675623066745737
Twist angle ° 270.05020−0.5−0.3
380.05−220−0.8−0.6
590.05−220−0.6−1.8
6100.05−220−0.8−1.9
Leading edge camber%311 10 3 −110−1−0.69
chord512 10 3 −11000.64
length613 10 3 −11010.63
AoA ° 140.011322.32.7
Table 2. Summary of the cost and constraint values for various shapes in the Pareto front.
Table 2. Summary of the cost and constraint values for various shapes in the Pareto front.
w C D C L F RCS σ ^ S ref N Cycle
1.0346.97 112.5 · 10 3 186.57 · 10 3 8.4 · 10 3 68.641
0.995352.16 112.5 · 10 3 144.41 · 10 3 6.9 · 10 3 68.663
0.99354.34 112.5 · 10 3 104.80 · 10 3 6.0 · 10 3 68.662
0.95365.58 112.5 · 10 3 94.15 · 10 3 6.3 · 10 3 68.640
0.9367.49 112.5 · 10 3 84.99 · 10 3 7.8 · 10 3 68.659
0.8372.88 112.5 · 10 3 78.02 · 10 3 25.3 · 10 3 68.647
0.75380.85 112.5 · 10 3 40.75 · 10 3 81.6 · 10 3 68.669
0.52381.72 112.5 · 10 3 40.31 · 10 3 84.4 · 10 3 68.665
0.51426.99 112.5 · 10 3 2.24 · 10 3 20.9 · 10 3 68.691
0.5427.22 112.5 · 10 3 2.23 · 10 3 21.0 · 10 3 68.698
0.0432.54 112.5 · 10 3 1.74 · 10 3 23.3 · 10 3 68.696
baseline387.46 112.6 · 10 3 82.27 · 10 3 7.4 · 10 3 68.6
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

Thoulon, C.; Roge, G.; Pironneau, O. Gradient-Based Aero-Stealth Optimization of a Simplified Aircraft. Fluids 2024, 9, 174. https://doi.org/10.3390/fluids9080174

AMA Style

Thoulon C, Roge G, Pironneau O. Gradient-Based Aero-Stealth Optimization of a Simplified Aircraft. Fluids. 2024; 9(8):174. https://doi.org/10.3390/fluids9080174

Chicago/Turabian Style

Thoulon, Charles, Gilbert Roge, and Olivier Pironneau. 2024. "Gradient-Based Aero-Stealth Optimization of a Simplified Aircraft" Fluids 9, no. 8: 174. https://doi.org/10.3390/fluids9080174

APA Style

Thoulon, C., Roge, G., & Pironneau, O. (2024). Gradient-Based Aero-Stealth Optimization of a Simplified Aircraft. Fluids, 9(8), 174. https://doi.org/10.3390/fluids9080174

Article Metrics

Back to TopTop