# Advanced Aerostructural Optimization Techniques for Aircraft Design

Authors: Yingtao Zuo; Pingjian Chen; Lin Fu; Zhenghong Gao; Gang Chen
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Traditional coupled aerostructural design optimization (ASDO) of aircraft based on high-fidelity models is computationally expensive and inefficient. To improve the efficiency, the key is to predict aerostructural performance of the aircraft efficiently. The cruise shape of the aircraft is parameterized and optimized in this paper, and a methodology named reverse iteration of structural model (RISM) is adopted to get the aerostructural performance of cruise shape efficiently. A new mathematical explanation of RISM is presented in this paper. The efficiency of RISM can be improved by four times compared with traditional static aeroelastic analysis. General purpose computing on graphical processing units (GPGPU) is adopted to accelerate the RISM further, and GPU-accelerated RISM is constructed. The efficiency of GPU-accelerated RISM can be raised by about 239 times compared with that of the loosely coupled aeroelastic analysis. Test shows that the fidelity of GPU-accelerated RISM is high enough for optimization. Optimization framework based on Kriging model is constructed. The efficiency of the proposed optimization system can be improved greatly with the aid of GPU-accelerated RISM. An unmanned aerial vehicle (UAV) is optimized using this framework and the range is improved by 4.67% after optimization, which shows effectiveness and efficiency of this framework.


---

## Body

## 1. Introduction

Coupled ASDO has reached a high academic level and now it is widely utilized in enterprises and universities. Many researchers construct complex multifidelity optimization framework to satisfy the various requirements of different phases of aircraft design [1, 2]. Various coupled ASDO frameworks are developed in the past decades [3–5]. Typical applications include the detailed wing optimization, winglet design [6], supersonic aircraft design [7, 8], and so forth.Low-fidelity models are widely used in coupled ASDO [9–11]. Optimization based on low-fidelity models can reflect the complex interaction of aerodynamics and structure. It is widely used in conceptual design to reach the optimal solution space rapidly. Optimization based on high-fidelity models is necessary in detail design phase to enhance the optimization quality further. High-fidelity models include the compressible Euler/NS equations and structural finite element method. The huge expense of optimization based on high-fidelity models encumbers its engineering application. Therefore, one of the main tasks is to increase the optimization efficiency. Numerous efforts have been made, and these efforts include the following:(1)
Tightly coupled solving strategy is adopted to solve the static aeroelastic problem, which can reduce the computational expense [12].
(2)
Many optimization frameworks [13, 14] based on various kinds of surrogate model are developed to reduce the computational expense.
(3)
Gradient-based optimization [15–18] is utilized to improve the optimization efficiency.
(4)
GPGPU is adopted to speed up the optimization.The most popular aerostructural prediction method is loosely coupled aeroelastic analysis for its simplicity. It solves the aerodynamic problem until it converges and then the structural problem and then the aerodynamic problem until it converges again. This procedure continues until the deformation of wing converges. It needs five iterations of repeated full-cycle aerodynamic/structural analysis at least, which is very time-consuming and inefficient. On the contrary, the tightly coupled analysis method is very efficient [17]. The flow-solid system consists of structure and flow subproblems, and they are successively solved at each time step. With time marching, the flow and deformation converge and the aerostructural performance is achieved. Tightly coupled strategy is relatively complex for implementation and is rarely found in commercial software packages. So it is seldom used in engineering.Gradient-based optimization is very efficient. The key is to get the gradient of objective with efficiency. Usually, this is accomplished by adjoint method, in which the expense is independent on the number of design variables. Adjoint method is widely used in aerodynamic design optimization. However, the derivation of adjoint equations in aerostructural design optimization is very complex and the convergence of the adjoint equations needs to be improved greatly [18] before it is widely used.Driven by the insatiable market demand of calculation, GPGPU has experienced a fast development in recent years. Many applications have been done in computational fluid dynamics (CFD) [19–24], and twenty times of acceleration have been achieved at least compared with conventional computation based on CPU [25]. Generally speaking, the computational expense of ASDO is much larger than that of the conventional aerodynamic optimization. Thus, optimization based on GPU provides a very bright prospect for ASDO, and ASDO based on GPU architecture is an innovative research area. Some work has been done. Typical work includes Patterson’s conceptual design studies of electric aircraft [24], with a higher order vortex lattice solver run on GPU to predict the aerodynamic performance. No researches have been done in coupled CFD/CSD design optimization up to now.This paper aims to construct an efficient multidiscipline design optimization framework based on CFD/CSD coupling. Two measures are taken to ensure the efficiency. Firstly, we propose to directly optimize the cruise shape instead of the jig shape. A novel methodology named RISM is adopted to evaluate the stress and strain of the aircraft [26], and a new mathematical explanation of RISM is presented. The advantage of RISM is that it avoids the repeated and time-consuming aerodynamic and structural iteration, and the efficiency can be improved greatly. Secondly, a GPU-accelerated CFD solver is developed, which can greatly raise the efficiency. Based on this novel methodology and GPU-accelerated CFD solver, a high-efficiency aerostructural optimization framework is constructed.The following part of this paper is divided into 4 sections. Section2 provides an introduction of the tools and methods used in the paper, and the novel aerostructural performance prediction methodology is described in detail. Section 3 provides optimization framework used in this work. Section 4 validates the proposed methodologies and an aerostructural optimization design case is done. Finally, Section 5 gives the conclusions.


## 2. Methodology

### 2.1. Flow Solver

In this context, a multiblock viscous flow solver named LMNS3D [25] for steady and unsteady turbulent flows based on GPU parallel methodology under the finite volume frame is employed. The nondimensional Navier-Stokes equations are solved on body-fitted structured mesh. All the variables are stored and operated in the cell center.Several approximate Riemann solvers, for example, ROE, Harten-Lax-van Leer with contact discontinuities (HLLC), and Simple Low-dissipation AUSM (SLAU), are available for inviscid flux calculation [27]. SLAU is mainly exploited due to its parameter-free property and all-speed simulation capability. In order to obtain high order, third-order Monotone Upwind Schemes for Scalar Conservation Laws (MUSCL) and fifth-order Weighted Essentially Nonoscillatory (WENO) scheme are utilized to reconstruct primitive variables on the cell interface. Second-order and fourth-order central scheme can be chosen to evaluate the viscous flux in a conservative way. Considering the turbulent model, one-equation Spalart-Allmaras (SA) and two-equation Menter’s k
-
w SST models are adopted for steady RANS simulations. These eddy turbulent models are discretized and updated in a loosely coupled way from the mean governing equation. Furthermore, scale-resolving methods, for example, Detached Eddy Simulation (DES) and Scale-Adaptive Simulation (SAS), are constructed for high-fidelity turbulent flow simulations.Data-Parallel Lower–Upper Relaxation (DP-LUR) [28] method is chosen as the implicit time marching method as it eliminates the data dependency in the LU-SGS method by replacing forward and backward sweeps as point-wise iterations. This key property enables highly efficient GPU cards based parallel simulations. Second-order temporal accuracy is enforced by imposing dual time-stepping strategy. The heterogeneous multiple CPU + GPU coprocessing system is established by implementing NVIDIA’s CUDA as well as Message Passing Interface (MPI) programming models. The data transfer time is hidden by a carefully designed concurrent coping and execution algorithm. On the other hand, multithreaded OpenMp parallel technique is also developed for shared-memory simulations. Through extensive numerical experiments, the robustness, accuracy, and efficiency of current in-house solver are well validated based on consumer-market oriented GPU cards and professional GPU cards. Double precision arithmetic is always kept through the entire residual computations with the help of the latest GPU hardware and careful design of CFD codes. The aerodynamic performance calculated by the GPU-accelerated solver is the same as the result calculated by the corresponding CPU solver.


### 2.2. The Proposed Aerostructural Performance Prediction Methodology

Jig shape is the aircraft shape that does not undergo aerodynamic load, gravity, and the propulsive force. It deforms into cruise shape at cruise condition by applied force. Traditional coupled aerostructural design optimization deals with the jig shape of the aircraft. The jig shape is parameterized and static aeroelastic analysis is adopted to get the corresponding aerostructural performance, which is used to evaluate the objective function in multidisciplinary design optimization. We also get the corresponding cruise shape after static aeroelastic analysis. Usually, the loosely coupled CFD/CSD simulation is adopted in static aeroelastic analysis [29]. The procedure of loosely coupled aeroelastic analysis is presented in Figure 1.Figure 1
Procedure of loosely coupled aeroelastic analysis.Cruise shape can be corrected into jig shape. The procedure to get the jig shape from cruise shape is called jig shape correction and it is widely used in aircraft design. It is notable that the jig shape correction procedure provides not only the jig shape, but also the distribution of the stress and strain. In light of it, we put forward to optimize the cruise shape of aircraft directly. For a given cruise shape, the aerodynamic performance of the aircraft at cruise condition can be achieved directly by aerodynamic analysis. The repeated aerodynamic/structural iterations to the jig shape as done in traditional coupled aerostructural design are avoided. Now, what we need to do is to get the stress and strain of the aircraft efficiently.We can apply the jig shape correction to get the structural performance of the aircraft. Aly’s methodology to get initial jig shape is considered firstly [30]. The procedure of Aly’s method is listed as follows:(1)
Get the aerodynamic load of cruise shape.
(2)
Applying the above aerodynamic load and gravity in the reverse direction to the cruise shape, conduct structural analysis to get the displacement of the aircraft.
(3)
Add the displacement to the cruise shape to get the deformed aircraft.If finite element method is adopted, the structural nodes coordinateX
j
i
g of jig shape can be achieved by solving the following: (1)
X
d
=
K
-
1
F
,
(2)
X
j
i
g
=
X
c
+
X
d
,where [
K
] represents the stiffness matrix of the cruise shape, X
d means the unknown vector of the structural nodes displacement, and F means the forces including the aerodynamic force and force of gravity acting on the structural nodes. X
c is the structural nodes coordinates of cruise. The structural performance of the aircraft can also be achieved by solving (1).If we apply the aerodynamic load of cruise shape and force of gravity to the jig shape, the deflected jig shape has some difference with the cruise shape as pointed by Aly [30]. The distribution of stress and strain also has some difference with the former structural performance. Alternatively, if static aeroelastic analysis is done to the above jig shape, the deflected jig shape does not coincide with the cruise shape too. This is mainly caused by the difference of the stiffness matrices of these two configurations. Inaccurate jig shape leads to inaccurate structural performance of the aircraft. To get the correct structural performance of the cruise shape, improvements must be done.The general improved jig shape correction or the conventional jig shape correction (CJSC) is precise. It gives the correct jig shape as well as the accurate structural performance. The main improvement of CJSC compared with Aly’s method is that static aeroelastic analysis is performed to the jig shape, and the jig shape is adjusted according to the difference between the deflected jig shape and cruise shape. Usually, CFD and CSD method is adopted in CJSC. The procedure can be summarized in detail as follows [31]:(1)
Call the CFD solver to feature the aerodynamic characteristics of the cruise shape and we get the aerodynamic load of the aircraft at cruise condition.
(2)
Get the initial jig shapeX
j
i
g by Aly’s methodology.
(3)
Aeroelastic analysis as mentioned above is conducted to get the deflected jig shapeX
d
i
s.
(4)
Compare the configurations of coordinatesX
d
i
s and X
c. Their difference is achieved by Δ
X
=
X
c
-
X
d
i
s.
(5)
If the 2-norm ofΔ
X is small enough, finish. Otherwise, the jig shape is updated by X
j
=
X
j
+
ω
Δ
X, where ω is a factor between 0 and 1. Go to (
3
).The procedure can be represented by Figure2.Figure 2
Conventional jig shape correction method.The aeroelastic analysis is needed only once to get the aerostructural performance of an aircraft in conventional aerostructural optimization. However, the above method involves several times of static aeroelastic analysis to get the aerostructural performance of an aircraft. Obviously, it is inefficient and not suitable for use in aerostructural optimizations.It is noticed in the above procedure that as the deflected jig shape converges to the cruise shape, the aerodynamic load acting on the jig shape converges to the aerodynamic load of the cruise shape. In the last iteration, the aerodynamic load acting on the jig shape equals the aerodynamic load of the cruise shape. This prompts a way of using the loads of the cruise shape to get the deflected jig shape directly without iteratively calling the CFD solver. Here, a novel aerostructural performance analysis methodology named RISM is adopted as follows [26]:(1)
Call the GPU-accelerated CFD solver to obtain the aerodynamic characteristics of the cruise shape and we get the aerodynamic load of the aircraft.
(2)
Get the initial jig shapeX
j
i
g by Aly’s methodology.
(3)
Apply the aerodynamic forces of the cruise shape and the force of gravity in the right direction to the jig shape to get the deflected jig shapeX
d
i
s
′.
(4)
Compare the displacement of every structural node ofX
c and X
d
i
s
′. Calculate the coordinate difference vector of the corresponding structural nodes and we get Δ
X, where Δ
X
=
X
c
-
X
d
i
s
′. Add ω
Δ
X to the structural nodes of the jig shape and we get the updated jig shape.
(5)
Go to step(
3
) unless the 2-norm of Δ
X is small enough.About 15 iterations are needed in this procedure throughout which the aerodynamic load is invariable. This procedure can be presented by Figure3.Figure 3
Aerostructural analysis of aircraft.As can be seen in Figure3, the CFD solver is called only once to get the aerodynamic performance of the cruise shape, and CSD solver is called iteratively to get the structural performance. Generally, structural models composed of shell and beam elements are used to represent a stiffened aircraft wing in aerostructural optimization. The expense to conduct one time of structural analysis is negligible compared with that of aerodynamic analysis based on CFD. The CFD solver is called at least five times in general loosely coupled static aeroelastic analysis. Therefore, the efficiency to get the aerostructural performance by RISM can be improved by at least four times compared with the loosely coupled aeroelastic analysis.The coordinate difference vector of finite element mesh of the jig shape and cruise shape is denoted byU. Now that the jig shape will deform into the cruise shape under the aerodynamic load of cruise shape and force of gravity, the strain of the jig shape should be U. The coordinate vector of the finite element mesh of the jig shape is X
c
-
U. The deformation of the jig shape under the aerodynamic load of cruise shape and force of gravity can be modeled as follows:(3)
K
X
c
,
U
U
-
F
=
0
,where K
(
X
c
,
U
) is the structural stiffness matrix of the jig shape. It depends on the finite element mesh of the cruise shape X
c, the displacement U between the jig shape and cruise shape, and the property of material. F is the resultant force of the aerodynamic load of the cruise shape and force of gravity.Equations (3) are nonlinear algebraic equations. The overrelaxation iteration method is adopted to solve them. Firstly, an appropriate initial value U
0 is provided by Aly’s method. Then, it can be solved by the following iterations:(4)
U
k
+
1
=
1
-
ω
U
k
+
ω
K
X
c
,
U
k
-
1
F
k
=
0,1
,
…
,where ω is the relaxation factor between 0 and 1.There lies an interesting contrast between the loosely coupled static aeroelastic analysis and RISM. The finite element mesh is fixed in loosely coupled static aeroelastic analysis, and the CFD grid as well as aerodynamic loads of FEM updates during iterations. But in RISM, the CFD grid and aerodynamic loads are fixed, and the finite element mesh updates during iterations.To accelerate the above procedures, the former GPU-accelerated CFD solver can be adopted and the GPU-accelerated RISM is formed.


### 2.3. Geometry Parameterization

The free form deformation (FFD) method is used to parameterize the cruise shape in this paper. The FFD was proposed by Sederberg and Parry [32] in 1986 and then used in graphics processing, and now it is widely used by CAD and cartoon. Later, it was introduced into aircraft design [33, 34]. FFD constructs a R
3
→
R
′
3 mapping function X
=
f
(
X
′
) from physical space to parameter space, where X
′ is logic coordinate of parameter space and X is coordinate of physical space. The deformation of aircraft in physical space is controlled by the movement of the control points of the parameters s, t, and u by the following formula:(5)
X
s
,
t
,
u
+
Δ
X
s
,
t
,
u
=
∑
i
=
1
l
∑
j
-
1
m
∑
k
=
1
n
B
l
-
1
i
-
1
s
B
m
-
1
j
-
1
t
B
n
-
1
k
-
1
u
P
i
,
j
,
k
+
Δ
P
i
,
j
,
k
,where P
i
,
j
,
k and Δ
P
i
,
j
,
k are the matrixes representing the original coordinates and the displacements of the node point (i
,
j
,
k) of the control box, respectively. The number of control points of the hexahedral control box is l
×
m
×
n. (s
,
t
,
u) is the local curvilinear coordinates mapped into the control box, and it is also called lattice coordinates. B
l
-
1
i
-
1
(
s
) is the (
i
-
1
)th Bernstein polynomial of degree l
-
1 defined as follows:(6)
B
l
-
1
i
-
1
s
=
l
-
1
!
i
-
1
!
l
-
i
!
s
i
-
1
1
-
s
l
-
i
.In the matrix form, equation can be written as follows:(7)
Δ
X
=
B
s
,
t
,
u
·
Δ
P
,
δ
x
1
δ
x
2
δ
x
3
=
B
1,1
,
1
⋯
B
l
,
m
,
n
Δ
P
1,1
,
1,1
Δ
P
1,1
,
1,2
Δ
P
1,1
,
1,3
⋮
⋮
⋮
Δ
P
l
,
m
,
n
,
1
Δ
P
l
,
m
,
n
,
2
Δ
P
l
,
m
,
n
,
3
,where(8)
B
i
,
j
,
k
=
B
l
-
1
i
-
1
s
B
m
-
1
j
-
1
t
B
n
-
1
k
-
1
u
.The procedures of FFD are listed as follows:(1)
Construct control lattice around the objective to be parameterized.
(2)
Calculate the logic coordinate.
(3)
Move the control point of lattice.
(4)
Calculate the deformation of the aircraft with (5).


## 2.1. Flow Solver

In this context, a multiblock viscous flow solver named LMNS3D [25] for steady and unsteady turbulent flows based on GPU parallel methodology under the finite volume frame is employed. The nondimensional Navier-Stokes equations are solved on body-fitted structured mesh. All the variables are stored and operated in the cell center.Several approximate Riemann solvers, for example, ROE, Harten-Lax-van Leer with contact discontinuities (HLLC), and Simple Low-dissipation AUSM (SLAU), are available for inviscid flux calculation [27]. SLAU is mainly exploited due to its parameter-free property and all-speed simulation capability. In order to obtain high order, third-order Monotone Upwind Schemes for Scalar Conservation Laws (MUSCL) and fifth-order Weighted Essentially Nonoscillatory (WENO) scheme are utilized to reconstruct primitive variables on the cell interface. Second-order and fourth-order central scheme can be chosen to evaluate the viscous flux in a conservative way. Considering the turbulent model, one-equation Spalart-Allmaras (SA) and two-equation Menter’s k
-
w SST models are adopted for steady RANS simulations. These eddy turbulent models are discretized and updated in a loosely coupled way from the mean governing equation. Furthermore, scale-resolving methods, for example, Detached Eddy Simulation (DES) and Scale-Adaptive Simulation (SAS), are constructed for high-fidelity turbulent flow simulations.Data-Parallel Lower–Upper Relaxation (DP-LUR) [28] method is chosen as the implicit time marching method as it eliminates the data dependency in the LU-SGS method by replacing forward and backward sweeps as point-wise iterations. This key property enables highly efficient GPU cards based parallel simulations. Second-order temporal accuracy is enforced by imposing dual time-stepping strategy. The heterogeneous multiple CPU + GPU coprocessing system is established by implementing NVIDIA’s CUDA as well as Message Passing Interface (MPI) programming models. The data transfer time is hidden by a carefully designed concurrent coping and execution algorithm. On the other hand, multithreaded OpenMp parallel technique is also developed for shared-memory simulations. Through extensive numerical experiments, the robustness, accuracy, and efficiency of current in-house solver are well validated based on consumer-market oriented GPU cards and professional GPU cards. Double precision arithmetic is always kept through the entire residual computations with the help of the latest GPU hardware and careful design of CFD codes. The aerodynamic performance calculated by the GPU-accelerated solver is the same as the result calculated by the corresponding CPU solver.


## 2.2. The Proposed Aerostructural Performance Prediction Methodology

Jig shape is the aircraft shape that does not undergo aerodynamic load, gravity, and the propulsive force. It deforms into cruise shape at cruise condition by applied force. Traditional coupled aerostructural design optimization deals with the jig shape of the aircraft. The jig shape is parameterized and static aeroelastic analysis is adopted to get the corresponding aerostructural performance, which is used to evaluate the objective function in multidisciplinary design optimization. We also get the corresponding cruise shape after static aeroelastic analysis. Usually, the loosely coupled CFD/CSD simulation is adopted in static aeroelastic analysis [29]. The procedure of loosely coupled aeroelastic analysis is presented in Figure 1.Figure 1
Procedure of loosely coupled aeroelastic analysis.Cruise shape can be corrected into jig shape. The procedure to get the jig shape from cruise shape is called jig shape correction and it is widely used in aircraft design. It is notable that the jig shape correction procedure provides not only the jig shape, but also the distribution of the stress and strain. In light of it, we put forward to optimize the cruise shape of aircraft directly. For a given cruise shape, the aerodynamic performance of the aircraft at cruise condition can be achieved directly by aerodynamic analysis. The repeated aerodynamic/structural iterations to the jig shape as done in traditional coupled aerostructural design are avoided. Now, what we need to do is to get the stress and strain of the aircraft efficiently.We can apply the jig shape correction to get the structural performance of the aircraft. Aly’s methodology to get initial jig shape is considered firstly [30]. The procedure of Aly’s method is listed as follows:(1)
Get the aerodynamic load of cruise shape.
(2)
Applying the above aerodynamic load and gravity in the reverse direction to the cruise shape, conduct structural analysis to get the displacement of the aircraft.
(3)
Add the displacement to the cruise shape to get the deformed aircraft.If finite element method is adopted, the structural nodes coordinateX
j
i
g of jig shape can be achieved by solving the following: (1)
X
d
=
K
-
1
F
,
(2)
X
j
i
g
=
X
c
+
X
d
,where [
K
] represents the stiffness matrix of the cruise shape, X
d means the unknown vector of the structural nodes displacement, and F means the forces including the aerodynamic force and force of gravity acting on the structural nodes. X
c is the structural nodes coordinates of cruise. The structural performance of the aircraft can also be achieved by solving (1).If we apply the aerodynamic load of cruise shape and force of gravity to the jig shape, the deflected jig shape has some difference with the cruise shape as pointed by Aly [30]. The distribution of stress and strain also has some difference with the former structural performance. Alternatively, if static aeroelastic analysis is done to the above jig shape, the deflected jig shape does not coincide with the cruise shape too. This is mainly caused by the difference of the stiffness matrices of these two configurations. Inaccurate jig shape leads to inaccurate structural performance of the aircraft. To get the correct structural performance of the cruise shape, improvements must be done.The general improved jig shape correction or the conventional jig shape correction (CJSC) is precise. It gives the correct jig shape as well as the accurate structural performance. The main improvement of CJSC compared with Aly’s method is that static aeroelastic analysis is performed to the jig shape, and the jig shape is adjusted according to the difference between the deflected jig shape and cruise shape. Usually, CFD and CSD method is adopted in CJSC. The procedure can be summarized in detail as follows [31]:(1)
Call the CFD solver to feature the aerodynamic characteristics of the cruise shape and we get the aerodynamic load of the aircraft at cruise condition.
(2)
Get the initial jig shapeX
j
i
g by Aly’s methodology.
(3)
Aeroelastic analysis as mentioned above is conducted to get the deflected jig shapeX
d
i
s.
(4)
Compare the configurations of coordinatesX
d
i
s and X
c. Their difference is achieved by Δ
X
=
X
c
-
X
d
i
s.
(5)
If the 2-norm ofΔ
X is small enough, finish. Otherwise, the jig shape is updated by X
j
=
X
j
+
ω
Δ
X, where ω is a factor between 0 and 1. Go to (
3
).The procedure can be represented by Figure2.Figure 2
Conventional jig shape correction method.The aeroelastic analysis is needed only once to get the aerostructural performance of an aircraft in conventional aerostructural optimization. However, the above method involves several times of static aeroelastic analysis to get the aerostructural performance of an aircraft. Obviously, it is inefficient and not suitable for use in aerostructural optimizations.It is noticed in the above procedure that as the deflected jig shape converges to the cruise shape, the aerodynamic load acting on the jig shape converges to the aerodynamic load of the cruise shape. In the last iteration, the aerodynamic load acting on the jig shape equals the aerodynamic load of the cruise shape. This prompts a way of using the loads of the cruise shape to get the deflected jig shape directly without iteratively calling the CFD solver. Here, a novel aerostructural performance analysis methodology named RISM is adopted as follows [26]:(1)
Call the GPU-accelerated CFD solver to obtain the aerodynamic characteristics of the cruise shape and we get the aerodynamic load of the aircraft.
(2)
Get the initial jig shapeX
j
i
g by Aly’s methodology.
(3)
Apply the aerodynamic forces of the cruise shape and the force of gravity in the right direction to the jig shape to get the deflected jig shapeX
d
i
s
′.
(4)
Compare the displacement of every structural node ofX
c and X
d
i
s
′. Calculate the coordinate difference vector of the corresponding structural nodes and we get Δ
X, where Δ
X
=
X
c
-
X
d
i
s
′. Add ω
Δ
X to the structural nodes of the jig shape and we get the updated jig shape.
(5)
Go to step(
3
) unless the 2-norm of Δ
X is small enough.About 15 iterations are needed in this procedure throughout which the aerodynamic load is invariable. This procedure can be presented by Figure3.Figure 3
Aerostructural analysis of aircraft.As can be seen in Figure3, the CFD solver is called only once to get the aerodynamic performance of the cruise shape, and CSD solver is called iteratively to get the structural performance. Generally, structural models composed of shell and beam elements are used to represent a stiffened aircraft wing in aerostructural optimization. The expense to conduct one time of structural analysis is negligible compared with that of aerodynamic analysis based on CFD. The CFD solver is called at least five times in general loosely coupled static aeroelastic analysis. Therefore, the efficiency to get the aerostructural performance by RISM can be improved by at least four times compared with the loosely coupled aeroelastic analysis.The coordinate difference vector of finite element mesh of the jig shape and cruise shape is denoted byU. Now that the jig shape will deform into the cruise shape under the aerodynamic load of cruise shape and force of gravity, the strain of the jig shape should be U. The coordinate vector of the finite element mesh of the jig shape is X
c
-
U. The deformation of the jig shape under the aerodynamic load of cruise shape and force of gravity can be modeled as follows:(3)
K
X
c
,
U
U
-
F
=
0
,where K
(
X
c
,
U
) is the structural stiffness matrix of the jig shape. It depends on the finite element mesh of the cruise shape X
c, the displacement U between the jig shape and cruise shape, and the property of material. F is the resultant force of the aerodynamic load of the cruise shape and force of gravity.Equations (3) are nonlinear algebraic equations. The overrelaxation iteration method is adopted to solve them. Firstly, an appropriate initial value U
0 is provided by Aly’s method. Then, it can be solved by the following iterations:(4)
U
k
+
1
=
1
-
ω
U
k
+
ω
K
X
c
,
U
k
-
1
F
k
=
0,1
,
…
,where ω is the relaxation factor between 0 and 1.There lies an interesting contrast between the loosely coupled static aeroelastic analysis and RISM. The finite element mesh is fixed in loosely coupled static aeroelastic analysis, and the CFD grid as well as aerodynamic loads of FEM updates during iterations. But in RISM, the CFD grid and aerodynamic loads are fixed, and the finite element mesh updates during iterations.To accelerate the above procedures, the former GPU-accelerated CFD solver can be adopted and the GPU-accelerated RISM is formed.


## 2.3. Geometry Parameterization

The free form deformation (FFD) method is used to parameterize the cruise shape in this paper. The FFD was proposed by Sederberg and Parry [32] in 1986 and then used in graphics processing, and now it is widely used by CAD and cartoon. Later, it was introduced into aircraft design [33, 34]. FFD constructs a R
3
→
R
′
3 mapping function X
=
f
(
X
′
) from physical space to parameter space, where X
′ is logic coordinate of parameter space and X is coordinate of physical space. The deformation of aircraft in physical space is controlled by the movement of the control points of the parameters s, t, and u by the following formula:(5)
X
s
,
t
,
u
+
Δ
X
s
,
t
,
u
=
∑
i
=
1
l
∑
j
-
1
m
∑
k
=
1
n
B
l
-
1
i
-
1
s
B
m
-
1
j
-
1
t
B
n
-
1
k
-
1
u
P
i
,
j
,
k
+
Δ
P
i
,
j
,
k
,where P
i
,
j
,
k and Δ
P
i
,
j
,
k are the matrixes representing the original coordinates and the displacements of the node point (i
,
j
,
k) of the control box, respectively. The number of control points of the hexahedral control box is l
×
m
×
n. (s
,
t
,
u) is the local curvilinear coordinates mapped into the control box, and it is also called lattice coordinates. B
l
-
1
i
-
1
(
s
) is the (
i
-
1
)th Bernstein polynomial of degree l
-
1 defined as follows:(6)
B
l
-
1
i
-
1
s
=
l
-
1
!
i
-
1
!
l
-
i
!
s
i
-
1
1
-
s
l
-
i
.In the matrix form, equation can be written as follows:(7)
Δ
X
=
B
s
,
t
,
u
·
Δ
P
,
δ
x
1
δ
x
2
δ
x
3
=
B
1,1
,
1
⋯
B
l
,
m
,
n
Δ
P
1,1
,
1,1
Δ
P
1,1
,
1,2
Δ
P
1,1
,
1,3
⋮
⋮
⋮
Δ
P
l
,
m
,
n
,
1
Δ
P
l
,
m
,
n
,
2
Δ
P
l
,
m
,
n
,
3
,where(8)
B
i
,
j
,
k
=
B
l
-
1
i
-
1
s
B
m
-
1
j
-
1
t
B
n
-
1
k
-
1
u
.The procedures of FFD are listed as follows:(1)
Construct control lattice around the objective to be parameterized.
(2)
Calculate the logic coordinate.
(3)
Move the control point of lattice.
(4)
Calculate the deformation of the aircraft with (5).


## 3. Multidisciplinary Design Optimization Framework

Direct optimization based on GPU-accelerated RISM is still time-consuming despite its high efficiency. Thus, surrogate model is adopted in this paper, which is computationally cheap. The widely used optimization framework based on the surrogate model and the genetic algorithm is adopted. In this framework, most high-fidelity simulations are replaced by recursively updating surrogate models. Latin hypercube is selected as the sampling method since it is the most popular version of stratified sampling ones. The design space for each factor is uniformly divided by this technique. The optimization is organized as follows:(1)
Initial samples are generated by Latin hypercube sampling method. The responses of these sample points are evaluated by high-fidelity models such as GPU-accelerated RISM.
(2)
Construct Kriging surrogate model based on the sample points and corresponding responses.
(3)
Search the model to get the optimum by genetic algorithm. Validate the optimum by high-fidelity models.
(4)
If the variation of objective function is small enough, stop; otherwise, the new results are added to the sample dataset and go to step(
2
).The procedure can be represented by Figure4.Figure 4
Optimization based on Kriging model.


## 4. Testing and Analysis

### 4.1. Validation of the CFD Solver

The numerical accuracy of in-house solver LMNS3D is validated very carefully in [25], and the computational efficiency is also tested on consumer-market oriented GPU cards. Now, it is validated for efficiency in a new workstation provided by Beijing Rongtian Huihai Technology Company Limited. The workstation is a multi-GPU cluster. It contains three NVIDIA Tesla K10 GPUs and two six-core, Intel Xeon Processor E5-2620 CPUs. The memory is 40 GB. For benchmarking, the NS solver was firstly run on CPU, and then it was run on GPU.The same DLR F6-WBNP configuration and grid used in [25] is tested at first. The Mach number is 0.75, the angle of attack is 1.0 degree, and the Reynolds number is 4.0 × 106. The grid is shown in Figure 5, and the grid size is 8.2 million. WENO 5th reconstruction scheme coupled with SLAU flux scheme for inviscid terms and 4th-order central differencing scheme for viscous terms are adopted. Menter’s k
-
w SST model is adopted. The DP-LUR method is used for time advancement.Figure 5
Surface grid of DLR F6-WBNP.Table1 shows the elapsed time and speedup factor for 500 time steppings. As can be seen, the program is very effective. The CFD results are in line with that of [25]. The speedup factor is different because the test is performed on a different platform. The platform in [25] was a personal computer with NVIDIA GeForce 560Ti graphic card. It is a workstation with NVIDIA Tesla K10 GPU in this paper.Table 1
Elapsed wall-time and speedup factor for 500 time steppings.
Single core (CPU)
Single K10
Triple K10
Time (minutes)
501.15
18.39
8.57
Speedup factor
1.0
27.24
58.48A high altitude long endurance UAV is tested, too. The flight speed is 0.6 Mach. The Reynolds number is 1.5 × 106. The size of the CFD grid is 3.1 million. The classical SA model is adopted to describe turbulence. WENO 5th reconstruction scheme coupled with SLAU flux scheme for inviscid terms and 4th-order central differencing scheme for viscous terms are adopted. Figure 6 displays part of the CFD grid of the aircraft. Table 2 shows the elapsed time and speedup factor for 500 time steppings. As can be seen from this table, 27.8× speedup is achieved by single GPU platform compared to single core (CPU).Table 2
Elapsed wall-time and speedup factor for 500 time steppings.
Single core (CPU)
Single K10
Triple K10
Time (minutes)
200.2
7.20
3.35
Speedup factor
1.0
27.80
59.73Figure 6
CFD grid of UAV.


### 4.2. Validation of RISM Methodology

The prediction of aerostructural performance of the aircraft is the basis of the multidisciplinary design optimization. Therefore, the effectiveness of RISM should be carefully validated.The function of GPU cards is to speed up the CFD calculation. It does not affect the final aerodynamic force, and the results given by GPU-accelerated RISM are the same as that of RISM alone. So we do not compare the results given by GPU-accelerated RISM with that of RISM alone in the following of this paper.The above UAV is used to validate the precision and efficiency of the proposed methodology. The semispan of the UAV is 8.0 m, and the aspect ratio of this UAV is 17.6. The flight altitude of this UAV is 20 Km.In [26], the surface grid of CFD grid did not coincide with that of the FEM Mesh. Therefore, fluid-structure interaction (FSI) was adopted and error was introduced. This makes it difficult to judge the accuracy of RISM. To avoid this problem, the CFD surface grid of the UAV is directly extracted to construct the surface grid of FEM in this paper as shown in Figure 7. The main load-carrying components of the wing box are considered, including skins, ribs, wing spars, and stringers. The front and rear spar are defined at 15% and 65% along the chord, respectively. The spars and ribs are assumed to be “T-beams.” The aluminium alloy adopted has the elasticity modulus of 70 GPa and Poisson’s ratio of 0.33.Figure 7
Wing structural model.The aerostructural performance prediction of this UAV is done by the GPU-accelerated RISM and the jig shape is achieved at the same time. The result is presented in Table3. Two measures are taken to check the result. Firstly, the CJSC is accomplished, and we get the corresponding jig shape and the deflected jig shape (Config1). Secondly, static aeroelastic analysis (SAA) of the jig shape achieved by GPU-accelerated RISM is conducted, and the deflected jig shape (Config2) is achieved. Table 3 presents all the structural performance and wing mass of these right-hand side jig shapes achieved by all these methods. Table 4 shows the twist angle of the wing tip section. Figure 8 shows the twist angle of jig shape predicted by different methods along the spanwise. Table 5 presents the aerodynamic performance of different configurations, and the pressure contour and sectional pressure distribution are compared in Figure 9.Table 3
Structural performance of the UAV predicted by different methodologies.
Prediction method
Maximum displacement, m
Mass of the wing, Kg
Maximum stress, MPa
CJSC
0.286
81.5623
93.24
Aly’s methodology
0.291
81.6183
86.56
GPU-accelerated RISM
0.291
81.6074
95.64
SAA
0.2871
81.6074
95.10Table 4
Wing tip twist angle of jig shape predicted by different methods.
Prediction method
CJSC
Aly’s methodology
GPU-accelerated RISM
Twist (degree)
−0.3865
−0.4050
−0.3813Table 5
Aerodynamic characteristics of different configurations.
Configuration
C
L
C
D
C
m
Cruise shape
0.550
0.028735
−0.0804
Config1
0.550
0.028740
−0.08035
Config2
0.5495
0.028727
−0.08033Figure 8
Wing twist angle predicted by different method.Figure 9
Pressure contour and sectional pressure distribution of the cruise shape and the deflected jig shape.As can be seen from Table3, the maximum displacement difference between the jig shapes achieved by GPU-accelerated RISM and CJSC is 0.005 m. The maximum stress predicted by RISM is 2.5% larger than that of the CJSC. The maximum stress predicted by Aly’s method is 7.2% smaller than that of the CJSC. Table 4 shows that the twist angles’ difference of jig shapes predicted by CJSC and RISM is 0.0052°, which is negligible, while the difference of jig shape predicted by CJSC and Aly’s method reaches 0.0185°. It can be seen from Table 5 that the aerodynamic performance of the aircraft predicted by different methods is almost the same.With the aid of GPU architecture, the efficiency can be improved by about 110 times (with single K10 GPU) or 239 times (with triple K10 GPU) compared with loosely coupled static aeroelastic analysis running on a single CPU core. This result is impressive. The objective function used in aerostructural optimization consists of the mass of the wing and the drag of the cruise shape. The maximum stress and displacement are usually used as constraints. As can be seen from the above, the mass and drag of the aircraft predicted by RISM are very precise, indicating that the objective function given by RISM will be very precise. The error of maximum stress is relatively large, but it is not as critical for optimization as the objective function. All these results show that the proposed methodology can be adopted in aerostructural optimization.Figure10 shows convergence history of RISM with different ω, where the longitudinal axis represents the maximum distance between the corresponding structural nodes of deflected jig shape and the cruise shape. Figure 11 shows the cruise shape and jig shape achieved by RISM.Figure 10
Convergence history of RISM with differentω.Figure 11
Comparison of the jig shape and the cruise shape.


### 4.3. UAV Optimization

FFD approach is used to parameterize the above UAV. All the parameters such as the weight of the UAV and reference area refer to that of the right-hand side of the UAV in the following, for simplicity. The control framework is shown in Figure12. The fuselage is fixed during optimization. The wing geometry variables include chord of root, chord of tip, semispan, and twist of the wing. Twist is defined at four sections across the wing. A graphical representation of the geometric design variables is shown in Figure 12. The number of aerodynamic design variables is 7. In this work, the topology of the wing structure remains unchanged, which means that the number of spars and ribs and their planform-view locations are all fixed. The FEM is divided into three segments along the spanwise direction as shown in Figure 7. The thickness of skins, area of spar cap of the front spar, and rear spar of each segment are selected as design variables. Therefore, the number of structural variables is 9, and the total number of design variables is 16.Figure 12
FFD control framework.The range is chosen to be the objective function, which is given by Breguet equation. This equation considers the tradeoff between drag and structural weight very well. It can be written as follows:(9)
R
=
V
C
L
D
l
n
W
1
W
2
,where R is the range, V is the cruise velocity, C is the specific fuel consumption of the powerplant, L
/
D is the lift-to-drag ratio, and W
1 and W
2 are the initial and final weights of the UAV during cruise, respectively. The initial weight of the aircraft consists of the structural weight, fuel weight, and so on. It is fixed to be 12357.8 N in this optimization. The final weight is simply the structural weight of wing plus a given constant weight. We try to lighten the weight of the wing to increase the fuel (W
1-W
2). The actual lift-to-drag ratio varies over the different cruise stages due to aerostructural effects produced by changes in aircraft weight and the inertial fuel load distribution. In order to simplify the analysis, we maintain a constant lift-to-drag ratio and C. The cruise velocity and the flight altitude are fixed, so we only need to maximize L
/
D
ln
⁡
(
W
1
/
W
2
).Constraints must be enforced in optimization. Here, the lift in cruise condition is fixed to be 11132.8 N. The reference area to calculate the lift coefficient and the drag coefficient is fixed to be 7.27 m2 for convenience, which is the projected wing area of the initial right-hand side wing. Thus, the lift coefficient is fixed to be 0.55 in optimization, which can be achieved by adjusting the angle of attack periodically during aerodynamic calculation. The maximum von Mises stress of the UAV must be below an allowable value at maneuver and cruise condition. The load at cruise condition is multiplied by a load factor equal to 2.5, to consider the peak loads encountered during various flight maneuvers or caused by turbulent air. The operating condition is listed in Table 6.Table 6
Operating conditions.
State
Mach
Altitude, (Km)
Lift coefficient
Load factor
Cruise
0.60
20
0.55
1.0
Maneuver
0.60
10
0.326
2.5Now, we can summarize the optimization problem as follows:(10)
M
a
x
i
m
i
z
e
:
F
=
L
D
ln
⁡
W
1
W
2
S
u
b
j
e
c
t
t
o
:
C
L
=
0.55
6.3
m
2
<
S
w
i
n
g
<
6.7
m
2
δ
max
⁡
1
<
95
M
P
a
δ
max
⁡
2
<
225
M
P
a
,where δ
m
a
x
1 represents the maximum von Mises stress of the wing at cruise condition. δ
m
a
x
2 represents the maximum von Mises stress of the wing at maneuver condition. S
w
i
n
g is the projected area of the right-hand side wing. The UAV is optimized by the above optimization framework, and the number of initial samples is 280. Three K10 GPU cards are adopted in this optimization. The aerostructural performance at cruise condition is predicted by GPU-accelerated RISM, and the aerostructural performance at maneuver condition is predicted by GPU-accelerated loosely coupled static aeroelastic analysis.Totally, 295 times of GPU-accelerated RISM and static aeroelastic analysis are called, respectively, during this optimization. A comparison of the geometric design variables, the projected wing area, and aerostructural performance of the UAV before and after optimization is listed in Table7. As can be seen, the drag is reduced by 2.1 counts and the mass of the wing is reduced by 4.35 Kg. The range of the UAV is increased by 4.67%.Table 7
Optimization results.
Initial
Optimized
Twist 1 (degree)
0
−0.653
Twist 2 (degree)
0
−0.7634
Twist 3 (degree)
0
−1.57
Twist 4 (degree)
0
−1.59
Chord root (m)
1.210
1.290
Chord tip (m)
0.556
0.533
Semispan (m)
8.0
7.602
S
wing (m2)
6.576
6.43
C
D
0.0287
0.02849
δ
max
⁡
1 (MPa)
93.24
76.42
δ
max
⁡
2 (MPa)
225.0
183.0
Mass of wing (Kg)
81.6
77.25
F
4.235
4.433To check the optimization results, similar validation is conducted as in Section4.2. Table 8 lists the results of RISM, CJSC, and static aeroelastic analysis of the jig shape achieved by RISM. As can be seen, the difference is negligible. The correctness of RISM is validated again.Table 8
Optimization results validation.
CJSC
GPU-accelerated RISM
SAA
C
L
0.550
0.550
0.4982
C
D
0.02849
0.02849
0.02846
C
m
−0.0642
−0.0642
−0.0641
Maximum stress (MPa)
76.25
77.76
77.43
Maximum displacement (m)
0.2178
0.2198
0.2162The pressure contour of the wing before and after optimization is given in Figure13. Figures 14 and 15 show the comparisons of the von Mises stress and planform before and after optimization separately. The convergence history is shown in Figure 16. A comparison of the deflections of the initial and optimized design of the aircraft at cruise condition and a 2.5 g maneuver condition is given in Figure 17. All these results validate the proposed optimization framework.Figure 13
Pressure contour comparison before and after optimization.Figure 14
Von Mises stress before and after optimization.Figure 15
Planform comparison before and after optimization.Figure 16
Convergence history of the optimization.Figure 17
Front view showing comparison of the jig shape and configuration at maneuver condition before and after optimization.


## 4.1. Validation of the CFD Solver

The numerical accuracy of in-house solver LMNS3D is validated very carefully in [25], and the computational efficiency is also tested on consumer-market oriented GPU cards. Now, it is validated for efficiency in a new workstation provided by Beijing Rongtian Huihai Technology Company Limited. The workstation is a multi-GPU cluster. It contains three NVIDIA Tesla K10 GPUs and two six-core, Intel Xeon Processor E5-2620 CPUs. The memory is 40 GB. For benchmarking, the NS solver was firstly run on CPU, and then it was run on GPU.The same DLR F6-WBNP configuration and grid used in [25] is tested at first. The Mach number is 0.75, the angle of attack is 1.0 degree, and the Reynolds number is 4.0 × 106. The grid is shown in Figure 5, and the grid size is 8.2 million. WENO 5th reconstruction scheme coupled with SLAU flux scheme for inviscid terms and 4th-order central differencing scheme for viscous terms are adopted. Menter’s k
-
w SST model is adopted. The DP-LUR method is used for time advancement.Figure 5
Surface grid of DLR F6-WBNP.Table1 shows the elapsed time and speedup factor for 500 time steppings. As can be seen, the program is very effective. The CFD results are in line with that of [25]. The speedup factor is different because the test is performed on a different platform. The platform in [25] was a personal computer with NVIDIA GeForce 560Ti graphic card. It is a workstation with NVIDIA Tesla K10 GPU in this paper.Table 1
Elapsed wall-time and speedup factor for 500 time steppings.
Single core (CPU)
Single K10
Triple K10
Time (minutes)
501.15
18.39
8.57
Speedup factor
1.0
27.24
58.48A high altitude long endurance UAV is tested, too. The flight speed is 0.6 Mach. The Reynolds number is 1.5 × 106. The size of the CFD grid is 3.1 million. The classical SA model is adopted to describe turbulence. WENO 5th reconstruction scheme coupled with SLAU flux scheme for inviscid terms and 4th-order central differencing scheme for viscous terms are adopted. Figure 6 displays part of the CFD grid of the aircraft. Table 2 shows the elapsed time and speedup factor for 500 time steppings. As can be seen from this table, 27.8× speedup is achieved by single GPU platform compared to single core (CPU).Table 2
Elapsed wall-time and speedup factor for 500 time steppings.
Single core (CPU)
Single K10
Triple K10
Time (minutes)
200.2
7.20
3.35
Speedup factor
1.0
27.80
59.73Figure 6
CFD grid of UAV.


## 4.2. Validation of RISM Methodology

The prediction of aerostructural performance of the aircraft is the basis of the multidisciplinary design optimization. Therefore, the effectiveness of RISM should be carefully validated.The function of GPU cards is to speed up the CFD calculation. It does not affect the final aerodynamic force, and the results given by GPU-accelerated RISM are the same as that of RISM alone. So we do not compare the results given by GPU-accelerated RISM with that of RISM alone in the following of this paper.The above UAV is used to validate the precision and efficiency of the proposed methodology. The semispan of the UAV is 8.0 m, and the aspect ratio of this UAV is 17.6. The flight altitude of this UAV is 20 Km.In [26], the surface grid of CFD grid did not coincide with that of the FEM Mesh. Therefore, fluid-structure interaction (FSI) was adopted and error was introduced. This makes it difficult to judge the accuracy of RISM. To avoid this problem, the CFD surface grid of the UAV is directly extracted to construct the surface grid of FEM in this paper as shown in Figure 7. The main load-carrying components of the wing box are considered, including skins, ribs, wing spars, and stringers. The front and rear spar are defined at 15% and 65% along the chord, respectively. The spars and ribs are assumed to be “T-beams.” The aluminium alloy adopted has the elasticity modulus of 70 GPa and Poisson’s ratio of 0.33.Figure 7
Wing structural model.The aerostructural performance prediction of this UAV is done by the GPU-accelerated RISM and the jig shape is achieved at the same time. The result is presented in Table3. Two measures are taken to check the result. Firstly, the CJSC is accomplished, and we get the corresponding jig shape and the deflected jig shape (Config1). Secondly, static aeroelastic analysis (SAA) of the jig shape achieved by GPU-accelerated RISM is conducted, and the deflected jig shape (Config2) is achieved. Table 3 presents all the structural performance and wing mass of these right-hand side jig shapes achieved by all these methods. Table 4 shows the twist angle of the wing tip section. Figure 8 shows the twist angle of jig shape predicted by different methods along the spanwise. Table 5 presents the aerodynamic performance of different configurations, and the pressure contour and sectional pressure distribution are compared in Figure 9.Table 3
Structural performance of the UAV predicted by different methodologies.
Prediction method
Maximum displacement, m
Mass of the wing, Kg
Maximum stress, MPa
CJSC
0.286
81.5623
93.24
Aly’s methodology
0.291
81.6183
86.56
GPU-accelerated RISM
0.291
81.6074
95.64
SAA
0.2871
81.6074
95.10Table 4
Wing tip twist angle of jig shape predicted by different methods.
Prediction method
CJSC
Aly’s methodology
GPU-accelerated RISM
Twist (degree)
−0.3865
−0.4050
−0.3813Table 5
Aerodynamic characteristics of different configurations.
Configuration
C
L
C
D
C
m
Cruise shape
0.550
0.028735
−0.0804
Config1
0.550
0.028740
−0.08035
Config2
0.5495
0.028727
−0.08033Figure 8
Wing twist angle predicted by different method.Figure 9
Pressure contour and sectional pressure distribution of the cruise shape and the deflected jig shape.As can be seen from Table3, the maximum displacement difference between the jig shapes achieved by GPU-accelerated RISM and CJSC is 0.005 m. The maximum stress predicted by RISM is 2.5% larger than that of the CJSC. The maximum stress predicted by Aly’s method is 7.2% smaller than that of the CJSC. Table 4 shows that the twist angles’ difference of jig shapes predicted by CJSC and RISM is 0.0052°, which is negligible, while the difference of jig shape predicted by CJSC and Aly’s method reaches 0.0185°. It can be seen from Table 5 that the aerodynamic performance of the aircraft predicted by different methods is almost the same.With the aid of GPU architecture, the efficiency can be improved by about 110 times (with single K10 GPU) or 239 times (with triple K10 GPU) compared with loosely coupled static aeroelastic analysis running on a single CPU core. This result is impressive. The objective function used in aerostructural optimization consists of the mass of the wing and the drag of the cruise shape. The maximum stress and displacement are usually used as constraints. As can be seen from the above, the mass and drag of the aircraft predicted by RISM are very precise, indicating that the objective function given by RISM will be very precise. The error of maximum stress is relatively large, but it is not as critical for optimization as the objective function. All these results show that the proposed methodology can be adopted in aerostructural optimization.Figure10 shows convergence history of RISM with different ω, where the longitudinal axis represents the maximum distance between the corresponding structural nodes of deflected jig shape and the cruise shape. Figure 11 shows the cruise shape and jig shape achieved by RISM.Figure 10
Convergence history of RISM with differentω.Figure 11
Comparison of the jig shape and the cruise shape.


## 4.3. UAV Optimization

FFD approach is used to parameterize the above UAV. All the parameters such as the weight of the UAV and reference area refer to that of the right-hand side of the UAV in the following, for simplicity. The control framework is shown in Figure12. The fuselage is fixed during optimization. The wing geometry variables include chord of root, chord of tip, semispan, and twist of the wing. Twist is defined at four sections across the wing. A graphical representation of the geometric design variables is shown in Figure 12. The number of aerodynamic design variables is 7. In this work, the topology of the wing structure remains unchanged, which means that the number of spars and ribs and their planform-view locations are all fixed. The FEM is divided into three segments along the spanwise direction as shown in Figure 7. The thickness of skins, area of spar cap of the front spar, and rear spar of each segment are selected as design variables. Therefore, the number of structural variables is 9, and the total number of design variables is 16.Figure 12
FFD control framework.The range is chosen to be the objective function, which is given by Breguet equation. This equation considers the tradeoff between drag and structural weight very well. It can be written as follows:(9)
R
=
V
C
L
D
l
n
W
1
W
2
,where R is the range, V is the cruise velocity, C is the specific fuel consumption of the powerplant, L
/
D is the lift-to-drag ratio, and W
1 and W
2 are the initial and final weights of the UAV during cruise, respectively. The initial weight of the aircraft consists of the structural weight, fuel weight, and so on. It is fixed to be 12357.8 N in this optimization. The final weight is simply the structural weight of wing plus a given constant weight. We try to lighten the weight of the wing to increase the fuel (W
1-W
2). The actual lift-to-drag ratio varies over the different cruise stages due to aerostructural effects produced by changes in aircraft weight and the inertial fuel load distribution. In order to simplify the analysis, we maintain a constant lift-to-drag ratio and C. The cruise velocity and the flight altitude are fixed, so we only need to maximize L
/
D
ln
⁡
(
W
1
/
W
2
).Constraints must be enforced in optimization. Here, the lift in cruise condition is fixed to be 11132.8 N. The reference area to calculate the lift coefficient and the drag coefficient is fixed to be 7.27 m2 for convenience, which is the projected wing area of the initial right-hand side wing. Thus, the lift coefficient is fixed to be 0.55 in optimization, which can be achieved by adjusting the angle of attack periodically during aerodynamic calculation. The maximum von Mises stress of the UAV must be below an allowable value at maneuver and cruise condition. The load at cruise condition is multiplied by a load factor equal to 2.5, to consider the peak loads encountered during various flight maneuvers or caused by turbulent air. The operating condition is listed in Table 6.Table 6
Operating conditions.
State
Mach
Altitude, (Km)
Lift coefficient
Load factor
Cruise
0.60
20
0.55
1.0
Maneuver
0.60
10
0.326
2.5Now, we can summarize the optimization problem as follows:(10)
M
a
x
i
m
i
z
e
:
F
=
L
D
ln
⁡
W
1
W
2
S
u
b
j
e
c
t
t
o
:
C
L
=
0.55
6.3
m
2
<
S
w
i
n
g
<
6.7
m
2
δ
max
⁡
1
<
95
M
P
a
δ
max
⁡
2
<
225
M
P
a
,where δ
m
a
x
1 represents the maximum von Mises stress of the wing at cruise condition. δ
m
a
x
2 represents the maximum von Mises stress of the wing at maneuver condition. S
w
i
n
g is the projected area of the right-hand side wing. The UAV is optimized by the above optimization framework, and the number of initial samples is 280. Three K10 GPU cards are adopted in this optimization. The aerostructural performance at cruise condition is predicted by GPU-accelerated RISM, and the aerostructural performance at maneuver condition is predicted by GPU-accelerated loosely coupled static aeroelastic analysis.Totally, 295 times of GPU-accelerated RISM and static aeroelastic analysis are called, respectively, during this optimization. A comparison of the geometric design variables, the projected wing area, and aerostructural performance of the UAV before and after optimization is listed in Table7. As can be seen, the drag is reduced by 2.1 counts and the mass of the wing is reduced by 4.35 Kg. The range of the UAV is increased by 4.67%.Table 7
Optimization results.
Initial
Optimized
Twist 1 (degree)
0
−0.653
Twist 2 (degree)
0
−0.7634
Twist 3 (degree)
0
−1.57
Twist 4 (degree)
0
−1.59
Chord root (m)
1.210
1.290
Chord tip (m)
0.556
0.533
Semispan (m)
8.0
7.602
S
wing (m2)
6.576
6.43
C
D
0.0287
0.02849
δ
max
⁡
1 (MPa)
93.24
76.42
δ
max
⁡
2 (MPa)
225.0
183.0
Mass of wing (Kg)
81.6
77.25
F
4.235
4.433To check the optimization results, similar validation is conducted as in Section4.2. Table 8 lists the results of RISM, CJSC, and static aeroelastic analysis of the jig shape achieved by RISM. As can be seen, the difference is negligible. The correctness of RISM is validated again.Table 8
Optimization results validation.
CJSC
GPU-accelerated RISM
SAA
C
L
0.550
0.550
0.4982
C
D
0.02849
0.02849
0.02846
C
m
−0.0642
−0.0642
−0.0641
Maximum stress (MPa)
76.25
77.76
77.43
Maximum displacement (m)
0.2178
0.2198
0.2162The pressure contour of the wing before and after optimization is given in Figure13. Figures 14 and 15 show the comparisons of the von Mises stress and planform before and after optimization separately. The convergence history is shown in Figure 16. A comparison of the deflections of the initial and optimized design of the aircraft at cruise condition and a 2.5 g maneuver condition is given in Figure 17. All these results validate the proposed optimization framework.Figure 13
Pressure contour comparison before and after optimization.Figure 14
Von Mises stress before and after optimization.Figure 15
Planform comparison before and after optimization.Figure 16
Convergence history of the optimization.Figure 17
Front view showing comparison of the jig shape and configuration at maneuver condition before and after optimization.


## 5. Conclusions

An efficient multidisciplinary design optimization framework is proposed in this paper. Message Passing Interface (MPI) and CUDA were used to accelerate the flow simulation in this framework. GPU-accelerated RISM was proposed to predict the aerodynamic and structural performance of the aircraft. The efficiency to predict the aerostructural performance of the aircraft at cruise condition is raised 4 times by RISM alone compared with conventional loosely coupled static aeroelastic analysis. It can be further improved by about 110 times (with single K10 GPU) or 239 times (with triple K10 GPU) with the aid of GPU and MPI. Tests show that RISM has almost the same fidelity as that of the conventional method. Optimization of an UAV proved the effectiveness of the optimization framework.



---

*Source: 753042-2015-10-13.xml*