# Study on a Rapid Aerodynamic Optimization Method of Flying Wing Aircraft for Conceptual Design

Authors: Chenhao Wei; Jun Huang; Lei Song
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

Searching for a design scheme satisfying the requirements in aircraft conceptual design can be a time-consuming work because of the multipeak and nonlinearity of the design space. This paper proposed a rapid aerodynamic optimization method for flying wing aircraft conceptual design. This method is aimed at reducing the induced drag at the design point by adjusting the camber and twist angle of spanwise airfoil. Firstly, the mean camber surface of the flying wing aircraft was parameterized. Secondly, the surrogate model was constructed based on the points selected by the optimization Latin square method. Thirdly, the surrogate model combined with a multi-island genetic algorithm was used for the preliminary solution of global optimum, and then, the Nonlinear Programming by Quadratic Lagrangian method combined with a vortex lattice method was used for searching the nearest exact best point from the initial best point. Finally, connecting with manual selection, the optimized flying wing layout scheme was obtained. The optimized results show that the induced drag coefficient is reduced by 10%, the pitching moment coefficient is reduced by an order of magnitude, and the lift drag ratio is increased from 26.3 to 27.3. The proposed optimization method decreases the time cost in aircraft conceptual design while achieving sufficient calculation accuracy. By using this method, the design space can be explored rapidly to search for the best design scheme satisfying the constraints.


---

## Body

## 1. Introduction

Aerodynamic shape optimization is aimed at improving the performance of aircraft including drag reduction [1], range improvement [2], and weight minimization [3], by changing its shape while satisfying specified constraints. In the early phase, the optimization mainly relies on wind tunnel experiments and flight tests. With the development of the computer technology and numerical calculation methods, computational fluid dynamic (CFD) tools (e.g., panel method and solving the Reynolds-averaged Navier-Stokes function) play a more and more important role in aerodynamic shape optimization.High-fidelity computational fluid dynamics tools usually need a long time to achieve a convergent result. Applying this type of analysis tool to the aerodynamic optimization may lead to unacceptable high time cost. The solution is to introduce an approximate model which is called a surrogate model to predict aerodynamic performance. Although the efficiency of the surrogate model will decrease with the increment of input variables and system nonlinearity [4], it still has the advantages of quickly exploring the design space and providing feasible solutions [5].The numerical optimization algorithms can be divided into two categories including derivative-based methods and derivative-free methods. Considering that the design space for the aerodynamic shape optimization is multipeak and nonlinear, a meta-heuristic algorithm belonging to a derivative-free method is more suitable. Gupta et al. evaluated the efficiency of nine meta-heuristic algorithms on solving eight mechanical design problems [6]. Meta-heuristic algorithms were widely used in multiobjective aircraft design [7, 8] and many other practical engineering problems in vehicle design [9–14]. Uncertainties exist in engineering design, and the reliability-based design optimization needs to consider these constraints. Demirci and Yıldız developed a hybrid approach to evaluate the constrain functions in the reliability-based design optimization. The robustness, simplicity, and efficiency were also analyzed [15].The aerodynamic shape of flying wing aircraft is concise. Changes of the wing surface will have a significant impact on its aerodynamic performance. Thus, numerical aerodynamic shape optimization is a fast and efficient method to improve aerodynamic performance of flying wing aircraft [16]. Recently, aerodynamic shape optimization in flying wing aircraft design mainly focused on lift-drag ratio improvement [17, 18], drag reduction in a cruise state [19–22], operating cost reduction [23], and stealth performance improvement [24]. Some optimizations were conducted based on the surrogate model [17, 18]. Both derivative-based methods and derivative-free methods have been applied in these researches [17, 19–22, 24]. Some researchers have conducted analysis on the impact of the asymmetric flight on aircraft’s minimum speed [25, 26] and the power optimization in the take-off phase [27]. This is especially problematic and risky in low-speed flight modes. In summary, little attention has been paid to the flying wing aerodynamic optimization applying the surrogate model and meta-heuristic algorithm.In this paper, a surrogate model-based aerodynamic optimization method is developed and this method is used to reduce the aerodynamic drag in a cruise state. The aerodynamic drag of aircraft in a cruise state is mainly composed of the zero-lift drag and the induced drag. Most of the zero-lift drag comes from the skin friction, and a small part comes from the pressure difference along the body caused by viscous separation. The whole zero-lift drag is approximately linearly proportional to the aircraft wetted area [28]. When the geometric relationship of the planform layout is fixed, the wetted area of the flying wing aircraft mainly changes with the maximum thickness of the airfoil, which is mainly defined by the inner structure layout. Therefore, the changing of the zero-lift drag by altering the spanwise camber and twist angle distribution is limited during the optimization. The primary purpose of this paper is to reduce the induced drag by optimizing the spanwise camber and twist angle distribution with unchanged thickness distribution.To conduct the optimization, the parametric model of the flying wing mean camber surface was established and was controlled by 16 parameters. A surrogate model was constructed based on the sample points in the design space extracted by an experimental design method. A rapid aerodynamic optimization method combining multi-island genetic algorithm (MIGA) and Nonlinear Programming by Quadratic Lagrangian method (NLPQL) was executed to reduce the induced drag under pitching moment balance constraints in flying wing aircraft conceptual design.


## 2. Materials and Methods

### 2.1. Planform Layout of the Study Object

The basic parameters of a flying wing aircraft under specific mission requirements were estimated. The preliminary general layout arrangement was further completed. According to the structural space and aerodynamic requirements, the spanwise airfoil section was manually selected. By taking the above steps, the planform of the initial design scheme was obtained and is shown in Figure1.Figure 1
Main parameters of the initial planform layout (m).


### 2.2. Vortex Lattice Method

The aerodynamic analysis is a crucial work throughout the aircraft design process. Compared with the N-S equation solving-based CFD method and the statistic-based experiential arithmetic, a vortex lattice method has a fair tradeoff between the calculation time consumption and the result accuracy. The vortex lattice method code used in this study was referred to the algorithm of tornado [29] and rewritten to achieve higher computational efficiency and accuracy, especially in the case of aircraft in sideslip. To validate the method, a flying wing aircraft [30] illustrated in Figure 2 was analyzed. Because of its concise shape, flying wing aircraft can be treated as a wing in vortex lattice method analysis.Figure 2
The planform layout and main parameters of the experimental aircraft (inch).The aerodynamic derivatives of the flying wing aircraft calculated by the vortex lattice method and NASA’s wind tunnel experiment results were compared and are shown in Figure3. It can be seen that when the angle of attack is within 10 degrees, the calculated results fit well with the experimental results. It only takes a few seconds using the vortex lattice method to calculate a state point running on a common PC. From the perspective of the computational accuracy and efficiency in an aircraft conceptual design phase, comparison results validated that the vortex lattice method in this study might be more suitable than the CFD method and wind tunnel experiment method for aerodynamic optimization.Figure 3
Comparison of experimental results with calculated results.


### 2.3. Aerodynamic Model Parameterization

The planform layout parameters of the aircraft are determined in Figure1. The aerodynamic optimization process is mainly oriented to the spanwise airfoil section. The study object needs parametric modelling before the optimization process.


#### 2.3.1. Parametric Modelling of the Airfoil Section

Based on the thin wing theory, the lift of aircraft is only related to the camber of airfoil at low speed [31]. It is further deduced that the induced drag is related to the lift distribution along the spanwise direction. Considering that the planform layout of the aircraft has been determined in this study, only aerodynamic twist and geometric twist along the spanwise direction can be optimized to reduce the induced drag. The concept of geometric twist is relatively simple. The essence of aerodynamic twist is to alter the zero-lift angle of attack of different spanwise sections by adjusting the airfoils and then influence the induced drag by changing the spanwise lift distribution.An airfoil curve was divided into two parts, including camber curve and thickness curve. This study only optimized the camber curve, and the thickness curve remains unchanged. It ensures that the optimized aerodynamic shape does not affect the preliminary general layout and structural design.In this study, the nonuniform rational B-spline curve [32] was used to parameterize the camber curve. Four parameters listed below were used to control the curve. The meaning of the parameters is shown in Figure 4. αLE is the angle between the camber line and the chord line on the leading edge. αTE is the angle between the camber line and the chord line on the trailing edge. xc is the chordwise position of the maximum camber. c is the maximum camber.Figure 4
The parameterization of the airfoil camber curve.


#### 2.3.2. Parametric Modelling of the Aircraft Mean Camber Surface

The main purpose of this section is to parameterize the aircraft main camber surface. A total of five control sections were set along the spanwise direction. The inner, middle, and outer sections are the main control airfoil sections, and each wing section was described by four parametersc, xc, αLE, and αTE. To adjust the transition surface between the main control airfoil sections, two hybrid airfoil sections were added between the main control airfoil sections. The camber parameters of each hybrid airfoil section were controlled by an independent parameter called the shape similarity ratio k (k∈0,1). When k<0.5, the hybrid airfoil is closer to the inner airfoil; otherwise, it is closer to the outer airfoil. The parametric description equations of hybrid airfoil are as follows.
(1)c=kcin+1−kcout,xc=kxcin+1−kxcout,αLE=kαLE.in+1−kαLE.out,αTE=kαTE.in+1−kαTE.out.Figure5 shows the spanwise control section distribution of the aircraft mean camber surface. The twist angle was defined as the angle that the section chord deflects around the line which passes through its leading-edge point and parallels to the Y axis. The twist angle of the aircraft was divided into two parts. The twist angle of the central section of the fuselage was defined as 0, the first twist angle was defined as the twist angle of the middle airfoil section, the second twist angle was defined as the outer airfoil section, and the twist angle performed linear transition among these sections.Figure 5
Control sections along the spanwise direction.


### 2.4. Aerodynamic Optimization Framework

#### 2.4.1. Optimization Variables and Constraints

The optimization variables in this paper were concentrated in the spanwise airfoil section since the planform layout has been determined. The main content for the general wing aerodynamic optimization is to reduce the induced drag in the cruise state by optimizing the spanwise airfoil sections and spanwise twist angle distribution.The parameters for the design point in cruise state were determined. The cruise altitude is 13000 m. The cruise speed is 180 m/s, and the lift coefficient is 0.3. The objective function is to minimize the induced drag coefficient. The optimization variables included camber parameters of the three main control airfoil sections (c, xc, αLE, and αTE), the shape similarity ratio of two hybrid airfoil sections k, and the twist angle of the middle airfoil section and outer airfoil section t. The value ranges of above optimization variables are shown in Table 1.Table 1
The range of optimization variables.
VariableMinimumMaximumc1-0.030.03c2-0.030.03c3-0.050.05αLE.1-10 deg10 degαLE.2-15 deg15 degαLE.3-20 deg20 degαTE-10 deg10 degxc0.20.8k0.20.8t-15 deg15 degFor the flying wing aircraft involved in this study, two more constraint factors must be considered in the optimization process. The first one is the longitudinal trim constraint, and the second one is stall angle of attack constraint.The pitching control of flying wing aircraft mainly depends on the lifting aileron at the trailing edge of the wing. Considering that its force arm is much smaller than that of the elevator at the trailing edge of the tail, it leads to the problem that the longitudinal trim ability is weak for flying wing aircraft. Thus, the longitudinal trim should be considered in the optimization process to ensure that the aircraft can be trimmed with control surface not deflected under the design condition.To ensure that there is a sufficient reserve lift coefficient for the aircraft to maneuver, enough margin should be kept between the angle of attack in the cruise state and the stall angle of attack. Thus, trimming angle of attack constraints should also be considered in the optimization process.Considering the longitudinal trim constraint and stall angle of attack constraint mentioned above, the constraints for the optimization are listed below:The pitching moment coefficient at the design point:−0.0005<Cm<0.0005.The angle of attack at the design point:α≤8°.


#### 2.4.2. Optimization Process Establishment

The aircraft aerodynamic optimization design space has obvious multipeak characteristics. In order to obtain the global optimum of the whole design space with computational efficiency, two core processes were executed for the aerodynamic optimization.The first step was the experimental design. The flow chart is illustrated in Figure6. The optimization Latin square method [33, 34] was used to extract 6000 points from the design space. The aerodynamic performance of these sample points was calculated. According to the results, the influence of various parameters on aerodynamic performance can be analyzed from a global perspective. The aerodynamic performance data can also support the construction of the surrogate model.Figure 6
Experimental design flow chart.The second step was the aerodynamic optimization, and it can be divided into two parts. The flow chart of the aerodynamic optimization is shown in Figure7. The RBF neural network surrogate model [35, 36] was used as an aerodynamic solver in the former part. The multi-island genetic algorithm [37–39] was used to execute global optimization. The results of the former part may be an approximate solution close to the real optimal solution and was treated as the start point in the latter part. The Nonlinear Programming by Quadratic Lagrangian method [40, 41] was used for further accurate optimization. The vortex lattice method was used as an aerodynamic solver in the latter part.Figure 7
Flow chart of the aerodynamic optimization.Because of the randomness of the multi-island genetic algorithm and the multipeak characteristics of design space, the results of each optimization round may be different. In this study, the aerodynamic optimization process circled 50 times. The optimal result among the 50 rounds was taken as the final result of the whole optimization process.


## 2.1. Planform Layout of the Study Object

The basic parameters of a flying wing aircraft under specific mission requirements were estimated. The preliminary general layout arrangement was further completed. According to the structural space and aerodynamic requirements, the spanwise airfoil section was manually selected. By taking the above steps, the planform of the initial design scheme was obtained and is shown in Figure1.Figure 1
Main parameters of the initial planform layout (m).


## 2.2. Vortex Lattice Method

The aerodynamic analysis is a crucial work throughout the aircraft design process. Compared with the N-S equation solving-based CFD method and the statistic-based experiential arithmetic, a vortex lattice method has a fair tradeoff between the calculation time consumption and the result accuracy. The vortex lattice method code used in this study was referred to the algorithm of tornado [29] and rewritten to achieve higher computational efficiency and accuracy, especially in the case of aircraft in sideslip. To validate the method, a flying wing aircraft [30] illustrated in Figure 2 was analyzed. Because of its concise shape, flying wing aircraft can be treated as a wing in vortex lattice method analysis.Figure 2
The planform layout and main parameters of the experimental aircraft (inch).The aerodynamic derivatives of the flying wing aircraft calculated by the vortex lattice method and NASA’s wind tunnel experiment results were compared and are shown in Figure3. It can be seen that when the angle of attack is within 10 degrees, the calculated results fit well with the experimental results. It only takes a few seconds using the vortex lattice method to calculate a state point running on a common PC. From the perspective of the computational accuracy and efficiency in an aircraft conceptual design phase, comparison results validated that the vortex lattice method in this study might be more suitable than the CFD method and wind tunnel experiment method for aerodynamic optimization.Figure 3
Comparison of experimental results with calculated results.


## 2.3. Aerodynamic Model Parameterization

The planform layout parameters of the aircraft are determined in Figure1. The aerodynamic optimization process is mainly oriented to the spanwise airfoil section. The study object needs parametric modelling before the optimization process.


### 2.3.1. Parametric Modelling of the Airfoil Section

Based on the thin wing theory, the lift of aircraft is only related to the camber of airfoil at low speed [31]. It is further deduced that the induced drag is related to the lift distribution along the spanwise direction. Considering that the planform layout of the aircraft has been determined in this study, only aerodynamic twist and geometric twist along the spanwise direction can be optimized to reduce the induced drag. The concept of geometric twist is relatively simple. The essence of aerodynamic twist is to alter the zero-lift angle of attack of different spanwise sections by adjusting the airfoils and then influence the induced drag by changing the spanwise lift distribution.An airfoil curve was divided into two parts, including camber curve and thickness curve. This study only optimized the camber curve, and the thickness curve remains unchanged. It ensures that the optimized aerodynamic shape does not affect the preliminary general layout and structural design.In this study, the nonuniform rational B-spline curve [32] was used to parameterize the camber curve. Four parameters listed below were used to control the curve. The meaning of the parameters is shown in Figure 4. αLE is the angle between the camber line and the chord line on the leading edge. αTE is the angle between the camber line and the chord line on the trailing edge. xc is the chordwise position of the maximum camber. c is the maximum camber.Figure 4
The parameterization of the airfoil camber curve.


### 2.3.2. Parametric Modelling of the Aircraft Mean Camber Surface

The main purpose of this section is to parameterize the aircraft main camber surface. A total of five control sections were set along the spanwise direction. The inner, middle, and outer sections are the main control airfoil sections, and each wing section was described by four parametersc, xc, αLE, and αTE. To adjust the transition surface between the main control airfoil sections, two hybrid airfoil sections were added between the main control airfoil sections. The camber parameters of each hybrid airfoil section were controlled by an independent parameter called the shape similarity ratio k (k∈0,1). When k<0.5, the hybrid airfoil is closer to the inner airfoil; otherwise, it is closer to the outer airfoil. The parametric description equations of hybrid airfoil are as follows.
(1)c=kcin+1−kcout,xc=kxcin+1−kxcout,αLE=kαLE.in+1−kαLE.out,αTE=kαTE.in+1−kαTE.out.Figure5 shows the spanwise control section distribution of the aircraft mean camber surface. The twist angle was defined as the angle that the section chord deflects around the line which passes through its leading-edge point and parallels to the Y axis. The twist angle of the aircraft was divided into two parts. The twist angle of the central section of the fuselage was defined as 0, the first twist angle was defined as the twist angle of the middle airfoil section, the second twist angle was defined as the outer airfoil section, and the twist angle performed linear transition among these sections.Figure 5
Control sections along the spanwise direction.


## 2.3.1. Parametric Modelling of the Airfoil Section

Based on the thin wing theory, the lift of aircraft is only related to the camber of airfoil at low speed [31]. It is further deduced that the induced drag is related to the lift distribution along the spanwise direction. Considering that the planform layout of the aircraft has been determined in this study, only aerodynamic twist and geometric twist along the spanwise direction can be optimized to reduce the induced drag. The concept of geometric twist is relatively simple. The essence of aerodynamic twist is to alter the zero-lift angle of attack of different spanwise sections by adjusting the airfoils and then influence the induced drag by changing the spanwise lift distribution.An airfoil curve was divided into two parts, including camber curve and thickness curve. This study only optimized the camber curve, and the thickness curve remains unchanged. It ensures that the optimized aerodynamic shape does not affect the preliminary general layout and structural design.In this study, the nonuniform rational B-spline curve [32] was used to parameterize the camber curve. Four parameters listed below were used to control the curve. The meaning of the parameters is shown in Figure 4. αLE is the angle between the camber line and the chord line on the leading edge. αTE is the angle between the camber line and the chord line on the trailing edge. xc is the chordwise position of the maximum camber. c is the maximum camber.Figure 4
The parameterization of the airfoil camber curve.


## 2.3.2. Parametric Modelling of the Aircraft Mean Camber Surface

The main purpose of this section is to parameterize the aircraft main camber surface. A total of five control sections were set along the spanwise direction. The inner, middle, and outer sections are the main control airfoil sections, and each wing section was described by four parametersc, xc, αLE, and αTE. To adjust the transition surface between the main control airfoil sections, two hybrid airfoil sections were added between the main control airfoil sections. The camber parameters of each hybrid airfoil section were controlled by an independent parameter called the shape similarity ratio k (k∈0,1). When k<0.5, the hybrid airfoil is closer to the inner airfoil; otherwise, it is closer to the outer airfoil. The parametric description equations of hybrid airfoil are as follows.
(1)c=kcin+1−kcout,xc=kxcin+1−kxcout,αLE=kαLE.in+1−kαLE.out,αTE=kαTE.in+1−kαTE.out.Figure5 shows the spanwise control section distribution of the aircraft mean camber surface. The twist angle was defined as the angle that the section chord deflects around the line which passes through its leading-edge point and parallels to the Y axis. The twist angle of the aircraft was divided into two parts. The twist angle of the central section of the fuselage was defined as 0, the first twist angle was defined as the twist angle of the middle airfoil section, the second twist angle was defined as the outer airfoil section, and the twist angle performed linear transition among these sections.Figure 5
Control sections along the spanwise direction.


## 2.4. Aerodynamic Optimization Framework

### 2.4.1. Optimization Variables and Constraints

The optimization variables in this paper were concentrated in the spanwise airfoil section since the planform layout has been determined. The main content for the general wing aerodynamic optimization is to reduce the induced drag in the cruise state by optimizing the spanwise airfoil sections and spanwise twist angle distribution.The parameters for the design point in cruise state were determined. The cruise altitude is 13000 m. The cruise speed is 180 m/s, and the lift coefficient is 0.3. The objective function is to minimize the induced drag coefficient. The optimization variables included camber parameters of the three main control airfoil sections (c, xc, αLE, and αTE), the shape similarity ratio of two hybrid airfoil sections k, and the twist angle of the middle airfoil section and outer airfoil section t. The value ranges of above optimization variables are shown in Table 1.Table 1
The range of optimization variables.
VariableMinimumMaximumc1-0.030.03c2-0.030.03c3-0.050.05αLE.1-10 deg10 degαLE.2-15 deg15 degαLE.3-20 deg20 degαTE-10 deg10 degxc0.20.8k0.20.8t-15 deg15 degFor the flying wing aircraft involved in this study, two more constraint factors must be considered in the optimization process. The first one is the longitudinal trim constraint, and the second one is stall angle of attack constraint.The pitching control of flying wing aircraft mainly depends on the lifting aileron at the trailing edge of the wing. Considering that its force arm is much smaller than that of the elevator at the trailing edge of the tail, it leads to the problem that the longitudinal trim ability is weak for flying wing aircraft. Thus, the longitudinal trim should be considered in the optimization process to ensure that the aircraft can be trimmed with control surface not deflected under the design condition.To ensure that there is a sufficient reserve lift coefficient for the aircraft to maneuver, enough margin should be kept between the angle of attack in the cruise state and the stall angle of attack. Thus, trimming angle of attack constraints should also be considered in the optimization process.Considering the longitudinal trim constraint and stall angle of attack constraint mentioned above, the constraints for the optimization are listed below:The pitching moment coefficient at the design point:−0.0005<Cm<0.0005.The angle of attack at the design point:α≤8°.


### 2.4.2. Optimization Process Establishment

The aircraft aerodynamic optimization design space has obvious multipeak characteristics. In order to obtain the global optimum of the whole design space with computational efficiency, two core processes were executed for the aerodynamic optimization.The first step was the experimental design. The flow chart is illustrated in Figure6. The optimization Latin square method [33, 34] was used to extract 6000 points from the design space. The aerodynamic performance of these sample points was calculated. According to the results, the influence of various parameters on aerodynamic performance can be analyzed from a global perspective. The aerodynamic performance data can also support the construction of the surrogate model.Figure 6
Experimental design flow chart.The second step was the aerodynamic optimization, and it can be divided into two parts. The flow chart of the aerodynamic optimization is shown in Figure7. The RBF neural network surrogate model [35, 36] was used as an aerodynamic solver in the former part. The multi-island genetic algorithm [37–39] was used to execute global optimization. The results of the former part may be an approximate solution close to the real optimal solution and was treated as the start point in the latter part. The Nonlinear Programming by Quadratic Lagrangian method [40, 41] was used for further accurate optimization. The vortex lattice method was used as an aerodynamic solver in the latter part.Figure 7
Flow chart of the aerodynamic optimization.Because of the randomness of the multi-island genetic algorithm and the multipeak characteristics of design space, the results of each optimization round may be different. In this study, the aerodynamic optimization process circled 50 times. The optimal result among the 50 rounds was taken as the final result of the whole optimization process.


## 2.4.1. Optimization Variables and Constraints

The optimization variables in this paper were concentrated in the spanwise airfoil section since the planform layout has been determined. The main content for the general wing aerodynamic optimization is to reduce the induced drag in the cruise state by optimizing the spanwise airfoil sections and spanwise twist angle distribution.The parameters for the design point in cruise state were determined. The cruise altitude is 13000 m. The cruise speed is 180 m/s, and the lift coefficient is 0.3. The objective function is to minimize the induced drag coefficient. The optimization variables included camber parameters of the three main control airfoil sections (c, xc, αLE, and αTE), the shape similarity ratio of two hybrid airfoil sections k, and the twist angle of the middle airfoil section and outer airfoil section t. The value ranges of above optimization variables are shown in Table 1.Table 1
The range of optimization variables.
VariableMinimumMaximumc1-0.030.03c2-0.030.03c3-0.050.05αLE.1-10 deg10 degαLE.2-15 deg15 degαLE.3-20 deg20 degαTE-10 deg10 degxc0.20.8k0.20.8t-15 deg15 degFor the flying wing aircraft involved in this study, two more constraint factors must be considered in the optimization process. The first one is the longitudinal trim constraint, and the second one is stall angle of attack constraint.The pitching control of flying wing aircraft mainly depends on the lifting aileron at the trailing edge of the wing. Considering that its force arm is much smaller than that of the elevator at the trailing edge of the tail, it leads to the problem that the longitudinal trim ability is weak for flying wing aircraft. Thus, the longitudinal trim should be considered in the optimization process to ensure that the aircraft can be trimmed with control surface not deflected under the design condition.To ensure that there is a sufficient reserve lift coefficient for the aircraft to maneuver, enough margin should be kept between the angle of attack in the cruise state and the stall angle of attack. Thus, trimming angle of attack constraints should also be considered in the optimization process.Considering the longitudinal trim constraint and stall angle of attack constraint mentioned above, the constraints for the optimization are listed below:The pitching moment coefficient at the design point:−0.0005<Cm<0.0005.The angle of attack at the design point:α≤8°.


## 2.4.2. Optimization Process Establishment

The aircraft aerodynamic optimization design space has obvious multipeak characteristics. In order to obtain the global optimum of the whole design space with computational efficiency, two core processes were executed for the aerodynamic optimization.The first step was the experimental design. The flow chart is illustrated in Figure6. The optimization Latin square method [33, 34] was used to extract 6000 points from the design space. The aerodynamic performance of these sample points was calculated. According to the results, the influence of various parameters on aerodynamic performance can be analyzed from a global perspective. The aerodynamic performance data can also support the construction of the surrogate model.Figure 6
Experimental design flow chart.The second step was the aerodynamic optimization, and it can be divided into two parts. The flow chart of the aerodynamic optimization is shown in Figure7. The RBF neural network surrogate model [35, 36] was used as an aerodynamic solver in the former part. The multi-island genetic algorithm [37–39] was used to execute global optimization. The results of the former part may be an approximate solution close to the real optimal solution and was treated as the start point in the latter part. The Nonlinear Programming by Quadratic Lagrangian method [40, 41] was used for further accurate optimization. The vortex lattice method was used as an aerodynamic solver in the latter part.Figure 7
Flow chart of the aerodynamic optimization.Because of the randomness of the multi-island genetic algorithm and the multipeak characteristics of design space, the results of each optimization round may be different. In this study, the aerodynamic optimization process circled 50 times. The optimal result among the 50 rounds was taken as the final result of the whole optimization process.


## 3. Results and Discussion

### 3.1. Surrogate Model Validation

Using the optimized Latin square method, 6000 sample points were extracted from the design space and the aerodynamic performance was calculated by the vortex lattice method. The histograms are drawn in Figures8–10 for the data including the angle of attack under the design condition, the pitching moment coefficient under the design condition, the induced drag coefficient under the design condition, the slope of the lift coefficient line, and the zero-lift angle of attack.Figure 8
Distribution of the angle of attack (a) and pitching moment coefficient (b).
(a)(b)Figure 9
Distribution of the zero-lift angle of attack (a) and slope of lift coefficient line (b).
(a)(b)Figure 10
Distribution of the induced drag coefficient corresponding to design point.According to the statistical histogram, except that the induced drag coefficient presents the distribution which is high in Figure8(a) and low in Figure 8(b), the other four parameters perform the distribution that the center is higher than two sides. Taking the accumulative percentage of 50%±40% as the boundary, the aerodynamic parameter range can be determined. The angle of attack under the design condition ranges from 0.4° to 8.3°, the pitching moment coefficient under the design condition ranges from -0.089 to 0.045, the induced drag coefficient under the design condition ranges from 0.0044 to 0.0096, and the slope of the lift coefficient line ranges from 0.069/deg to 0.071/deg. The zero-lift angle of attack ranges from -3.9°to 4.0°.Although the coupling relationship between aerodynamic parameters is not considered yet, the following conclusions can be approximately drawn from the above data:(1)
The optimization objectives and boundary conditions established above are reasonable and can be satisfied(2)
The self-trimming requirements of the aircraft can be realized within 8° angle of attack through reasonable arrangement of the spanwise camber distribution and twist angle distribution(3)
It is a reasonable optimization expectation that the induced drag coefficient is about 0.0040(4)
The main effect of changing the spanwise camber distribution is to change the zero-lift angle of attack, while the slope of the lift coefficient line is basically independent of the spanwise camber distributionTo improve the calculation speed of the optimization process, the RBF neural network algorithm was used to generate the surrogate model as the aerodynamic solver. In order to verify the accuracy of the surrogate model, 50 points were selected in the design space by the optimized Latin square method. The aerodynamic performance of these 50 points was calculated by the surrogate model. These results were compared with vortex lattice method calculation results. The relative errors between the results output by the surrogate model and the results calculated by the vortex lattice method are shown in Table2.Table 2
Relative error between the results output by the surrogate model and the results calculated by the vortex lattice method.
VariableAverage errorMaximum errorSquare root of errorR squaredCDi0.018490.113530.02820.98714Cm0.007220.033250.009880.99815αCL=0.30.004140.021580.006010.99934According to Table2, the RBF neural network surrogate model can obtain high fitting accuracy. With the largest error among the three output variables, the maximum relative error of CDi is within 12%. The average relative error is less than 2%. The maximum relative error of the other two variables is kept within 4%.


### 3.2. Aerodynamic Optimization Results and Analysis

50 different optimization results were obtained after 50 rounds. This phenomenon reflects the nonlinearity and multipeak characteristic of the design space. It is found that some optimization results do not satisfy the constraint conditions. It is caused by the error of the surrogate model. Although the RBF neural network has a good fitting accuracy in the design space, the number of validation points is relatively small. It can be found that there are still large errors between the surrogate model and the vortex lattice method at some design points. This error leads to the difference between the pitching moment coefficient calculated by the surrogate model and that obtained by the vortex lattice method in the Nonlinear Programming by Quadratic Lagrangian optimization process.Combining with manual selection, the optimal point of the results satisfying the constraints was taken as the final result of the whole optimization process. The thickness distribution of airfoil followed the initial airfoil parameters. The initial airfoil and the optimized airfoil of different spanwise sections are compared in Figures11–13. The initial airfoil is represented by a dotted line, and the optimized airfoil is represented by a solid line.Figure 11
Comparison of the airfoil of the BL0 (a) and BL1625 (b) section before and after optimization.
(a)(b)Figure 12
Comparison of the airfoil of the BL3250 (a) and BL5550 (b) section before and after optimization.
(a)(b)Figure 13
Comparison of the airfoil of the BL7556 section before and after optimization.Figure14 shows the comparison of aerodynamic performance before (red line) and after optimization (blue line). The optimization has no effect on the slope of the lift coefficient line, while the zero-lift angle of attack decreases slightly. In terms of the induced drag, when the angle of attack is above 2°, the induced drag coefficient is lower than that of the initial design. In terms of the lift-drag ratio, the optimized result is bigger than the initial design when the angle of attack is above 6°.Figure 14
Comparison of aerodynamic performance before and after optimization.The comparison of aerodynamic performance parameters before and after optimization is shown in Table3. It can be found that the angle of attack under cruise condition increases by 0.6° after optimization, the induced drag coefficient decreases by about 10%, and the pitching moment coefficient decreases by an order of magnitude. Since the initial design is obtained by manually adjusting the parameters, it is difficult to achieve the balance of pitching moment at the design point. Therefore, it is necessary to deflect the lift aileron for the initial design to achieve longitudinal trimming, which may bring additional trim drag.Table 3
The comparison of aerodynamic performance parameters before and after optimization.
VariableCdiCmL/DαCL=0.3Initial0.0041-0.002926.35.4Optimized0.00370.000227.34.8


## 3.1. Surrogate Model Validation

Using the optimized Latin square method, 6000 sample points were extracted from the design space and the aerodynamic performance was calculated by the vortex lattice method. The histograms are drawn in Figures8–10 for the data including the angle of attack under the design condition, the pitching moment coefficient under the design condition, the induced drag coefficient under the design condition, the slope of the lift coefficient line, and the zero-lift angle of attack.Figure 8
Distribution of the angle of attack (a) and pitching moment coefficient (b).
(a)(b)Figure 9
Distribution of the zero-lift angle of attack (a) and slope of lift coefficient line (b).
(a)(b)Figure 10
Distribution of the induced drag coefficient corresponding to design point.According to the statistical histogram, except that the induced drag coefficient presents the distribution which is high in Figure8(a) and low in Figure 8(b), the other four parameters perform the distribution that the center is higher than two sides. Taking the accumulative percentage of 50%±40% as the boundary, the aerodynamic parameter range can be determined. The angle of attack under the design condition ranges from 0.4° to 8.3°, the pitching moment coefficient under the design condition ranges from -0.089 to 0.045, the induced drag coefficient under the design condition ranges from 0.0044 to 0.0096, and the slope of the lift coefficient line ranges from 0.069/deg to 0.071/deg. The zero-lift angle of attack ranges from -3.9°to 4.0°.Although the coupling relationship between aerodynamic parameters is not considered yet, the following conclusions can be approximately drawn from the above data:(1)
The optimization objectives and boundary conditions established above are reasonable and can be satisfied(2)
The self-trimming requirements of the aircraft can be realized within 8° angle of attack through reasonable arrangement of the spanwise camber distribution and twist angle distribution(3)
It is a reasonable optimization expectation that the induced drag coefficient is about 0.0040(4)
The main effect of changing the spanwise camber distribution is to change the zero-lift angle of attack, while the slope of the lift coefficient line is basically independent of the spanwise camber distributionTo improve the calculation speed of the optimization process, the RBF neural network algorithm was used to generate the surrogate model as the aerodynamic solver. In order to verify the accuracy of the surrogate model, 50 points were selected in the design space by the optimized Latin square method. The aerodynamic performance of these 50 points was calculated by the surrogate model. These results were compared with vortex lattice method calculation results. The relative errors between the results output by the surrogate model and the results calculated by the vortex lattice method are shown in Table2.Table 2
Relative error between the results output by the surrogate model and the results calculated by the vortex lattice method.
VariableAverage errorMaximum errorSquare root of errorR squaredCDi0.018490.113530.02820.98714Cm0.007220.033250.009880.99815αCL=0.30.004140.021580.006010.99934According to Table2, the RBF neural network surrogate model can obtain high fitting accuracy. With the largest error among the three output variables, the maximum relative error of CDi is within 12%. The average relative error is less than 2%. The maximum relative error of the other two variables is kept within 4%.


## 3.2. Aerodynamic Optimization Results and Analysis

50 different optimization results were obtained after 50 rounds. This phenomenon reflects the nonlinearity and multipeak characteristic of the design space. It is found that some optimization results do not satisfy the constraint conditions. It is caused by the error of the surrogate model. Although the RBF neural network has a good fitting accuracy in the design space, the number of validation points is relatively small. It can be found that there are still large errors between the surrogate model and the vortex lattice method at some design points. This error leads to the difference between the pitching moment coefficient calculated by the surrogate model and that obtained by the vortex lattice method in the Nonlinear Programming by Quadratic Lagrangian optimization process.Combining with manual selection, the optimal point of the results satisfying the constraints was taken as the final result of the whole optimization process. The thickness distribution of airfoil followed the initial airfoil parameters. The initial airfoil and the optimized airfoil of different spanwise sections are compared in Figures11–13. The initial airfoil is represented by a dotted line, and the optimized airfoil is represented by a solid line.Figure 11
Comparison of the airfoil of the BL0 (a) and BL1625 (b) section before and after optimization.
(a)(b)Figure 12
Comparison of the airfoil of the BL3250 (a) and BL5550 (b) section before and after optimization.
(a)(b)Figure 13
Comparison of the airfoil of the BL7556 section before and after optimization.Figure14 shows the comparison of aerodynamic performance before (red line) and after optimization (blue line). The optimization has no effect on the slope of the lift coefficient line, while the zero-lift angle of attack decreases slightly. In terms of the induced drag, when the angle of attack is above 2°, the induced drag coefficient is lower than that of the initial design. In terms of the lift-drag ratio, the optimized result is bigger than the initial design when the angle of attack is above 6°.Figure 14
Comparison of aerodynamic performance before and after optimization.The comparison of aerodynamic performance parameters before and after optimization is shown in Table3. It can be found that the angle of attack under cruise condition increases by 0.6° after optimization, the induced drag coefficient decreases by about 10%, and the pitching moment coefficient decreases by an order of magnitude. Since the initial design is obtained by manually adjusting the parameters, it is difficult to achieve the balance of pitching moment at the design point. Therefore, it is necessary to deflect the lift aileron for the initial design to achieve longitudinal trimming, which may bring additional trim drag.Table 3
The comparison of aerodynamic performance parameters before and after optimization.
VariableCdiCmL/DαCL=0.3Initial0.0041-0.002926.35.4Optimized0.00370.000227.34.8


## 4. Conclusions

In this study, the mean camber surface of flying wing aircraft was parameterized. A surrogate model was constructed based on the sample points extracted from design space. A rapid aerodynamic optimization method combining the multi-island genetic algorithm (MIGA) and Nonlinear Programming by Quadratic Lagrangian method (NLPQL) was developed and investigated. The conclusions are summarized as follows:(1)
According to the performance data of the 6000 sample points, the slope of the lift coefficient line is basically independent of the camber spanwise distribution(2)
Among the verified results calculated by the surrogate model and the vortex lattice method, the maximum relative error ofCDi is within 12%. The average relative error of CDi is less than 2%. The maximum relative error of the pitching moment coefficient and the angle of attack is kept within 4%. The RBF neural network surrogate model has good fitting accuracy in a large range. However, in some local areas of the design space, the surrogate model has insufficient accuracy. Therefore, more work will need to be done to improve the surrogate model accuracy for the whole design space(3)
After the optimization of spanwise twist angle distribution and airfoil camber distribution, the optimized induced drag coefficient is reduced to 0.0037, which is 10% lower than the initial design. The pitching moment coefficient is reduced by an order of magnitude. The lift-drag ratio is increased from 26.3 to 27.3 after the optimization(4)
The optimization method established in this paper can obtain a balance between aerodynamic performance simulation result accuracy and calculation cost. By applying it in flying wing aircraft conceptual design, this method can contribute to exploring the design space quickly and searching for the feasible solution satisfying the constraints(5)
This work is focused on the drag reduction in a cruise state. The analysis of the impact of the asymmetric flight on the aircraft safe flight boundaries is important and unavoidable. Further studies need to be carried out on this issue



---

*Source: 5775355-2022-07-16.xml*