1. Introduction
With the omnipresence of radar technology, from early detection to target acquisition by 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 low Radar Cross-Section (RCS).
The frequency bands commonly used by the modern radars are L, S, C and X bands [
1]. Their working frequencies range from 1GHz to 12GHz, corresponding to wavelength between 0.025 m and 0.3 m, much smaller than a 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 wing planform and sections design. In order to mitigate the large computation cost 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 aerodynamic, 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 gradient-free optimization method, which have 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 o high-fidelity models.
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. [
6] 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.
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 already cover the application of the adjoint method to gradient-based optimization. For instance, Martin applied the adjoint method to robust airfoil optimization [
7]. Martin used the Automatic Differentiation (AD) tool TAPENADE [
8] 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 [
9] and to the Common Research Wing [
10].
1.2. Aero-Stealth Optimization
As for the aero-stealth optimization, Garnier et al. [
11] used an algorithm based on the coupling of the NSGA-II genetic algorithm and the multi-objective efficient global optimization approached [
12] to optimize an UCAV planform for aero-stealth, Liu et al. [
13] used hierarchical kriging models and surrogate models for aero-stealth optimization of a flying-wing aircraft and Zhou and Huang [
14] used a two-particle search algorithm to optimize the aero-stealth performances of a vertical tail.
Few studies however have dealt with the gradient-based aero-stealth optimization. Pan et al. [
15] performed a CAD-based, gradient-based aero-stealth optimization of a flying-wing aircraft, solving Prandtl-Glauert equation to compute the pressure coefficients of the aircraft and a PO method to predict RCS. Li et al. [
16] 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 aero and stealth criteria alternatively as the objective and a constraint and showed that the most performant optimization result is obtained when using an objective function built as a trade-off of teh aero and stealth criteria.
1.3. Contribution
The purpose of this study is to develop a methodology for gradient-based aero-stealth design using our in-house CAD modeler GANIMEDE [
17], a RANS solver [
18], a physical optics solver and the Sequential Least Square Quadratic Programming (SLSQP) algorithm [
19] in an industrial framework. We identify a Pareto front to showcase various optimal shapes corresponding to various aero-stealth trade-offs.
This paper is organized as follows. In
Figure 1, the flow chart of the gradient-based aero-stealth coupled optimization is demonstrated, and some numerical methods used in this work are described. Next, the results of the code validation are shown in
Section 3 and the problem formulation, baseline geometry of the flying wing, meshes and shape design variables are introduced in
Section 4 and
Section 5. 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 will 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 :
where
is the scattered electric fields at a point
in polar coordinates and
is the incident electric field at the target, whose norm is considered to be 1
n practice. Usually, the RCS is expressed in decibels :
We will 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 will produce on a target the size of an aircraft a 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. Physical Optics (PO) is one of the most common high frequency approximate method for RCS prediction. It is much cheaper than BEM, and it gives a better prediction than a purely optical approximation [
20].
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 [
21] integral equation is valid. Together with a tangent plane approximation, those can be simplified and the scattered electric field
at a point
P can be expressed as :
with
,
k the wave number,
the impedance of free space,
S the surface of the target,
and
the unit direction vectors for the incident and scattered fields respectively,
and
, and where
J is the surface current density vector. In the PO approximation,
where
is the local normal unit vector to the surface.
Note that by (
2) that
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) [
22] 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 the coordinates of the vertices of the surface grid used for the electromagnetism problem.
If
is a function depending explicitly only on the RCS
, then the directional derivative of
with respect to a set of geometric variables
(e.g. sweep angle, chord length...) is :
We used the Automatic Differentiation (AD) tool TAPENADE [
8] in order to generate a differentiated program from the original code which is able to compute the terms
and
. With this tool, it is possible to save a lot of development and debugging time. TAPENADE offers two modes : tangent (or direct) and reverse.
Tangent mode compute the directional derivative of in the direction 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. Reverse mode is an efficient way of computing the gradient of . When 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 of row by row and then compute than to compute it by the chain rule with the tangent mode.
Here, however, all the step of the chain rule computation are cheap to compute since they do not require 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
, that is the derivative of the vertex coordinates of the surface grid with respect to a set of geometric variable, is computed by the CAD modeler (see §
4.1) and is given as an input to the differentiated PO program.
This allows us to compute the cost function as well as its derivative in a single call of the PO code.
2.3. The aerodynamic solver
Our in-house Navier-Stokes AETHER code [
18,
23] 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 positivity properties and efficient preconditioning, entropy variables yield further improvements over the usual conservation variables, in particular 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 have been obtained using an Automatic Differentiation tool, TAPENADE. This tool provides the differentiated Fortran routines with respect to input variables prescribed by the user.
We used a compressible Reynolds-Averaged Navier Stokes [
24] method with the Spalart-Allmaras [
25] 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
Finally, let
be an observation function of the field
W. For the process of optimization, we need to compute the derivatives
and
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 an adjoint approach [
26].
When a surface grid of coordinates
is subjected to a deformation thanks to the CAD modeler, the deformation is propagated to the 3D grid of coordinates
with the help of an elliptic operator described by the system :
Let us write and as functions of , and W. Note that is an independent variable but depends on and W depends on both and .
2.4.1. Gradient with Respect to Aerodynamic Variables
First, the derivative of
and
with respect to
can be written as
and
The term
in (
6) is not easy to compute in itself, but it can be expressed differently thanks to (
7) :
Consequently, (
6) can then be written as :
Where
is the solution to the aerodynamic adjoint system :
2.4.2. Gradient with Respect to Geometric Variables
The derivative of
and
with respect to
are :
and
As in
Section 2.4.1, we can express the term
of equation (
11) thanks to (
12) as :
and then :
where
is the same as defined in equation (
10).
The work is not yet done, though, as the term
is not to compute. In order to find a more convenient expression, we must write the derivative of
with respect to
:
and then once again we can express the term
as a function of more convenient derivatives :
which finally yields :
where
is solution of the mesh adjoint equation system :
2.4.3. Summary
If
is an observation depending on
W, then the derivative of
with respect to a set of geometric variables
and respectively to a set of aerodynamic variables
can be written as :
where
is the coordinates of the 2D mesh for the aerodynamic problem,
is the coordinates of the 3D aerodynamic mesh and
is the mesh adjoint, solution of :
with
the aerodynamic adjoint, solution of :
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 :
Moreover, due to the symmetric properties of the sphere, its RCS does not depend on the bearing angle.
Since 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 has been computed with 2000 frequencies between
and
.
In
Figure 3a, the RCS is independent on the bearing angle, as it should be, and the error to the analytical value is around 0.15%.
Figure 3b shows that, when the frequency is too low (
, with
c 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 (
), 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 on
Section 4). The BEM computation was performed thanks to the SPECTRE solver [
27,
28].
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 ). 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.
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 side. We did not use a realistic geometry for the rear of the fuselage, so the physics of the wake will not be accurately represented. Since we only performed supersonic simulations, this will 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 use the CAD tool GANIMEDE [
17]. 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 [
17]. The surface is represented by a set of NURBS grids and parametrized 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 variable 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 to the baseline geometry, GANIMEDE will project the mesh onto the deformed geometry, and compute the derivatives of the mesh node coordinates with respect to the design variables 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.
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 overall drag coefficient is
drag count and the lift coefficient is
. These values are gathered in
Table 1.
The field of local Mach number in the symmetry plane and in the planes parallel to the symmetry plane at
and
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 towards 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
, the frequency should be higher than
in order to meet the criterion of
Section 3.1:
The COPOLA simulation is done on a range of bearing angles from 0 to with a step of and on a range of elevation angles from -20 to 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 } and spanning }. The step for the bearing and elevation angles has been chosen according to the Shannon criterion.
In
Figure 9b and
Figure 9a, a peak is visible around
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 [
29]. 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 . 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 is not the same with or without the PTD extension. This is likely due to a diffraction by the wing tip of by 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 will be using for RCS simulation in the rest of this study.
4.4. Optimization Parameters
The wing of the GFA is parametrized with the help of 6 control sections, with 15 control points along each of them (see
Figure 10). All the section are and will stay parallel to the symmetry plane. Section number 1 is at the symmetry plane and won’t be modified.
Section 2 is a few centimeters away from the root of the wing. At section 4, continuity of the tangent is not ensured, allowing for a sharp angle on the wing. The fuselage is not modified.
The geometry is parametrized with the help of 13 global variables :
the leading edge sweep angle the section 2 and 4,
the chord of section 2, 4 and 6,
the distance between section 2 and 6,
the twist angle of section 2, 3, 5 and 6,
the leading edge camber of section 3, 5 and 6.
The chords of section 3 and 5 are interpolated between that of section 2 and 4 and 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 at section 4 are interpolated between sections 3 and 5.
Finally, the angle of attack is the 14 optimization parameter. It is a private to the aerodynamic process.
Figure 11.
Definition of Sref
Figure 11.
Definition of Sref
5. Cost and Constraints
In order to explore multiple tradeoffs between aerodynamic and stealth optima, we define the cost function as a weighed sum of an aerodynamic criterion and a stealth criterion :
with
.
5.1. Aerodynamic Criterion
For the aerodynamic optimization, we minimize the drag coefficient (
) :
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 will settle 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, 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 , except for a limited interval around a peak direction where the RCS can be high.
In order to manage dimensionless quantities, the electromagnetic criterion is defined as follows :
where
is a weighed average of a function of
:
The function
is a differentiable threshold function built to cut the low intensity variations that we noted in
Figure 9, that have a negative impact on computation of the gradient of the cost function :
with
a parameter that allows to control the sharpness of the peak. The threshold
was chosen.
The function
is a weighting function chosen to favor the peak directions. It is close to 1 everywhere, except in the neighborhood of the direction
d where it decreases quickly :
Graphs of the function
are shown in
Figure 12 for various values of
. For this study, we have chosen, somewhat arbitrarily, a peak direction at a bearing of
, and
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 ().
Furthermore, we want to preserve the reference surface (
) of the wing.
is defined in
Figure 11.
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 it finite differences (FD) for various step sizes ranging from 1 to , in the unit of each individual variables. 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 (where 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, at least three times differentiable,
its derivative and
and
its derivative approximated with respectively forward and centered FD with a step size
. It is well known that
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 forward FD and 2 for centered FD. We consider that the minimal error should be at most for forward FD and for 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 DF gradients function of the step size, for both forward (in red) and centered (in green) difference, 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 . The centered difference is also acceptable with a minimal error of around and a slope of 2 right before the critical step size.
6.2. 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
), the relative error is used :
When the gradient is too small, the relative error maybe misleading, so the absolute error is used instead :
Most of the points show a pretty accurate gradient, with an error of at most . The only exceptions are the gradient with respect to the leading edge camber of the 3 and 5 sections. 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.
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
7.1. Aerodynamic Optimization
First, we performed an optimization with a purely aerodynamic objective function (
.). The planform of the optimal shape compared to that of the baseline is presented of
Figure 15. Most notably, the leading edge sweep angles have reached 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 will decrease the shock strength and the wave drag.
Figure 16 and
Figure 17 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
and
respectively. The decreasing strength of the bow shockwave can be clearly seen at both position.
Figure 23 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
, but the drag of the aero optimal shape (red line) starts increasing later and more gradually than for the other shapes, which translates the gain in wave drag because of the higher leading sweep angles.
Figure 24 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 has 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 18, 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 on this figure.
7.2. Stealth Optimization
We then performed an optimization with a purely stealth cost function (
). The planform of the optimal shape compared to that of the baseline is presented of
Figure 19. The leading and trailing edge are aligned with a sweep angle of
, which is the result expected from the choice of the stealth criterion. This allows the RCS to be concentrated in a sharp peak around a bearing angle of
, as seen in
Figure 20.
However, the decrease in leading edge sweep angle is very detrimental to the aerodynamic drag and particularly to the wave drag, as can be seen in
Figure 21 and
Figure 22, where a strong shock has appeared just upstream of the leading edge, which was absent from the other shapes. This is confirmed in
Figure 23 which shows a much steeper and earlier increase in
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 24, we can see than 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.
Figure 23.
Integration of CD along the X axis.
Figure 23.
Integration of CD along the X axis.
Figure 24.
Distribution of the wing CL.
Figure 24.
Distribution of the wing CL.
Figure 25.
Comparison of the aerodynamic polars of the wings of the different shapes.
Figure 25.
Comparison of the aerodynamic polars of the wings of the different shapes.
Table 2.
Gradient-Based Method
Table 2.
Gradient-Based Method
7.3. Pareto Optimality
In the two previous section, we have shown that the requirements for a good aerodynamic shape and a good stealth shape are mainly in competition. Therefore, the need for compromise arise. We have explored different combination and drawn a Pareto front (
Figure 34). We will present a few selected shapes. The values of the criteria and constraints for these shapes are summarized in
Table 3.
Note that the Pareto front features three main attraction points, around
,
and
. Those attractors are probably due to the steepness of the weighting used in the RCS cost function (see
Figure 12). The shapes around
and below are essentially the same as the purely RCS optimal shape (
Figure 19). The shape corresponding to
and
are showcased in
Figure 27 and
Figure 26 respectively. The drag difference between those shapes is pretty mild (8 drag count) and is largely due to the sweep angle as can be seen on the difference in pressure drag in
Figure 23. The difference in RCS cost function is more notable. Though the RCS cost function is lower for the
shape than the
, its average RCS is higher. This is due to the secondary reflection peak at
aligning with the primary peak at
. Those 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
shape is barely better than the baseline.
Figure 26.
Shape corresponding to the w = 0.8 point.
Figure 26.
Shape corresponding to the w = 0.8 point.
Figure 27.
Shape corresponding to the w = 0.75 point.
Figure 27.
Shape corresponding to the w = 0.75 point.
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).
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).
Figure 30.
omparison 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.75 shape (bottom).
Figure 30.
omparison 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.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).
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 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.
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.
8. Conclusion
In this study, we presented a constrained, gradient-based aero-stealth optimization method 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 fifteen parameters, but it is able to deal with a more complex shape parametrized 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 those gradients by comparing them to finite differences. The results show an acceptable difference.
After developing and validating the OP solver, we optimized a simplified aircraft geometry, with a weighted sum of the aerodynamic drag and a criterion based on the RCS repartition 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.
Each optimization takes less than 100 cycles to reach an optimal solution, a cycle being the evaluation of either the cost function and the constraints or their gradients. This is efficient, and it allows optimizing more complex shapes or use higher fidelity methods while keeping a reasonable computation time.
Going forward, a more advanced method to identify the Pareto front, such as the one proposed by Desideri [
30] would be beneficial.
Author Contributions
Conceptualization, Charles Thoulon, Gilbert Roge and Olivier Pironneau ; Methodology, Charles Thoulon and Gilbert Roge ; Software, Charles Thoulon ; Validation, Charles Thoulon ; Writing — original draft preparation, Charles Thoulon ; Writing — review, Gilbert Roge and Olivier Pironneau ; Supervision, Olivier Pironneau. All authors have read and agreed to the published version of the manuscript.
Funding
This research received no external funding
Acknowledgments:
Our thanks to Quentin Carayol, Hervé Steve and Nicolas Réau for their help to modify 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 declare that the research was conducted in the absence of any commercial or financial relationship that could be construed as a potential conflict of interest
References
- Knott, E.F.; Schaeffer, J.F.; Tulley, M.T. Radar cross section; SciTech Publishing, 2004. [Google Scholar]
- 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+. AIAA AVIATION 2022 Forum, 2022, p. 3322. [CrossRef]
- 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. 45th AIAA Aerospace Sciences Meeting and Exhibit, 2007, p. 36. [CrossRef]
- Wu, D.; Long, T.; Li, Y.; Jiang, M.; Huang, B. Aero-structure-stealth coupled optimization for high aspect ratio wing using adaptive metamodeling method. 15th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, 2014, p. 2304. [CrossRef]
- Jiang, X.; Zhao, Q.; Zhao, G.; Meng, C. Numerical Optimization on Aerodynamic/Stealth Characteristics of Airfoil Based on CFD/CEM Coupling Method. Transactions of Nanjing University of Aeronautics & Astronautics 2016, 33, 274–284. [Google Scholar]
- 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, Sichuan, China. Citeseer, 2014, Vol. 11, p. 585.
- Martin, L. Conception aérodynamique robuste. PhD thesis, Université de Toulouse, Université Toulouse III-Paul Sabatier, 2010. [Google Scholar]
- Hascoët, L.; Pascual, V. The Tapenade Automatic Differentiation tool: Principles, Model, and Specification. ACM Transactions On Mathematical Software 2013, 39. [Google Scholar] [CrossRef]
- Lyu, Z.; Martins, J.R. Aerodynamic design optimization studies of a blended-wing-body aircraft. Journal of Aircraft 2014, 51, 1604–1617. [Google Scholar] [CrossRef]
- Lyu, Z.; Kenway, G.K.; Martins, J.R. Aerodynamic shape optimization investigations of the common research model wing benchmark. AIAA journal 2015, 53, 968–985. [Google Scholar] [CrossRef]
- Garnier, E.; Langlet, S.; Klotz, P.; Cadillon, J.; Simon, J.; Castelli, J.C.; Levadoux, D. Surrogate-based conception for aerodynamic/stealth compromise. AVT-324 Specialists’ Meeting on Multidisciplinary Design Approaches and Performance Assessment of Future Combat Aircraft, 2020.
- Guerra, J. Optimisation multi-objectif sous incertitudes de phénomènes de thermique transitoire. PhD thesis, INSTITUT SUPERIEUR DE L’AERONAUTIQUE ET DE L’ESPACE (ISAE), 2016. [Google Scholar]
- Liu, Z.; Song, W.; Han, Z.; Wang, Y. Multifidelity Aerodynamic/Stealth Design Optimization Method for Flying Wing Aircraft. ICAS 2021. [Google Scholar]
- 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]
- 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 2017, 2, 15–17. [Google Scholar]
- Li, M.; Bai, J.; Li, L.; Meng, X.; Liu, Q.; Chen, B. A gradient-based aero-stealth optimization design method for flying wing aircraft. Aerospace Science and Technology 2019, 92. [Google Scholar] [CrossRef]
- Kleinveld, S.; Rogé, G.; Daumas, L.; Dinh, Q. Differentiated parametric CAD used within the context of automatic aerodynamic design optimization. 12th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, 2008, p. 5952. [CrossRef]
- Chalot, F. Industrial Aerodynamics. In Encyclopedia of Computational Mechanics; Fluids, chapter 12; Stein, E., de Borst, R., Hughes, T.J.R., Eds.; John Wiley & Sons, Ltd., 2004; Vol. 3. [Google Scholar]
- Barclay, A. SQP methods for large-scale optimization; University of California: San Diego, 1999. [Google Scholar]
- Laval, D. Application de méthodes asymptotiques à la simulation de la diffraction électromagnétique par un corps régulier. PhD thesis, Université Sciences et Technologies-Bordeaux I, 2006. [Google Scholar]
- Stratton, J.A.; Chu, L. Diffraction theory of electromagnetic waves. Physical Review 1939, 56, 99. [Google Scholar] [CrossRef]
- Michaeli, A. Elimination of infinities in equivalent edge currents, Part I: Fringe current components. IEEE Transactions on Antennas and Propagation 1986, 34, 912–918. [Google Scholar] [CrossRef]
- 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. Computer Methods in Applied Mechanics and Engineering 2023, 417, 116425. [Google Scholar] [CrossRef]
- Alfonsi, G. Reynolds-averaged Navier–Stokes equations for turbulence modeling 2009.
- Spalart, P.R.; Allmaras, S.R. One-equation turbulence model for aerodynamic flows. La Recherche Aerospatiale 1994, 1, 5–21. [Google Scholar]
- Martin, L.; Roge, G. Calcul de la sensibilité d’ordre deux d’une observation aérodynamique. ESAIM: Proceedings. EDP Sciences 2009, 27, 138–155. [Google Scholar] [CrossRef]
- Poggio, A.J.; Miller, E.K. Integral equation solutions of three-dimensional scattering problems; MB Assoc., 1970. [Google Scholar]
- Carayol, Q. Développement et analyse d’une méthode multipôle multiniveau pour l’électromagnétisme. PhD thesis, Paris 6, 2002. [Google Scholar]
- Jenn, D. Radar and laser cross section engineering; American Institute of Aeronautics and Astronautics, Inc., 2005. [Google Scholar]
- Désidéri, J.A. MGDA Variants for Multi-Objective Optimization. 2012.
Figure 1.
Flowchart of the automatic optimization process chain
Figure 1.
Flowchart of the automatic optimization process chain
Figure 2.
Definition of the bearing
(θ) and elevation (ϕ) angles.
Figure 2.
Definition of the bearing
(θ) and elevation (ϕ) angles.
Figure 3.
RCS of the sphere of radius R = 1m.
Figure 3.
RCS of the sphere of radius R = 1m.
Figure 4.
Comparison of the RCS function of the bearing angle predicted by COPOLA and BEM code on a GFA geometry (see
Figure 5).
Figure 4.
Comparison of the RCS function of the bearing angle predicted by COPOLA and BEM code on a GFA geometry (see
Figure 5).
Figure 5.
The Generic Fighter Aircraft.
Figure 5.
The Generic Fighter Aircraft.
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.
Figure 7.
istribution of local Mach number around the cross-section of the wing of the GFA, at Y/b = 0.45.
Figure 7.
istribution of local Mach number around the cross-section of the wing of the GFA, at Y/b = 0.45.
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.
Figure 9.
Profile of RCS with respect to the bearing angle on the GFA (elevation : ϕ = −20°).
Figure 9.
Profile of RCS with respect to the bearing angle on the GFA (elevation : ϕ = −20°).
Figure 10.
The wing parametrized with the control sections and control points.
Figure 10.
The wing parametrized with the control sections and control points.
Figure 12.
The function Γd(θ) with d = 30° for various values of τ.
Figure 12.
The function Γd(θ) with d = 30° for various values of τ.
Figure 13.
Error between the gradient of the stealth cost function and the same computed by 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 FD (both forward and central) function of the FD step size (ϵ).
Figure 14.
Error between the gradient of the aerodynamic drag and the same computed by FD (both forward and centered) for the optimal step size determined for the CAD gradient.
Figure 14.
Error between the gradient of the aerodynamic drag and the same computed by FD (both forward and centered) for the optimal step size determined for the CAD gradient.
Figure 15.
Optimal shape from a purely aerodynamic standpoint.
Figure 15.
Optimal shape from a purely aerodynamic standpoint.
Figure 16.
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 aero-optimal shape (bottom).
Figure 16.
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 aero-optimal shape (bottom).
Figure 17.
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 aero-optimal shape (bottom).
Figure 17.
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 aero-optimal shape (bottom).
Figure 18.
RCS of the aero-optimal shape compared to that of the baseline.
Figure 18.
RCS of the aero-optimal shape compared to that of the baseline.
Figure 19.
Optimal shape from a purely low-observability standpoint.
Figure 19.
Optimal shape from a purely low-observability standpoint.
Figure 20.
RCS of the stealth-optimal shape compared to that of the baseline.
Figure 20.
RCS of the stealth-optimal shape compared to that of the baseline.
Figure 21.
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 stealth-optimal shape (bottom).
Figure 21.
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 stealth-optimal shape (bottom).
Figure 22.
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 stealth-optimal shape (bottom).
Figure 22.
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 stealth-optimal shape (bottom).
Table 1.
Summary of the cost and contraints values for various shapes in the Pareto front.
Table 1.
Summary of the cost and contraints values for various shapes in the Pareto front.
Table 3.
Summary of the values of the design variables for a few chosen optimal shapes.
Table 3.
Summary of the values of the design variables for a few chosen optimal shapes.
|
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. |
© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).