Aeroelastic Optimization Design of the Global Stiffness for a Joined Wing Aircraft

Abstract

Due to the complexity and particularity of the joined wing layout, traditional design methods for the global stiffness of a high-aspect wing are not applicable for a joined wing. Herein, a beam-frame model and a three-dimensional wing-box model are built to solve the global stiffness aeroelastic optimization design problem for a joined wing. The goal is to minimize the weight, and the constraints are the overall aeroelastic requirements. Based on a genetic algorithm, two methods for the beam-frame model and one method for the three-dimensional model are used for comparative analysis. The results show that the optimization method for a diagonal beam section and the optimization method for an exponential/linear combination function fit are adequate for optimizing and designating the joined wing global stiffness. The distributions obtained using the two methods have good consistency and are similar to the distribution of the three-dimensional model. The stiffness distribution data and the beam section parameters can be converted from each other, which is convenient for redesigning the structure parameters using the stiffness distribution data, and is valuable for engineering applications.

1. Introduction

Aircraft design includes complex aeroelastic problems; therefore, it is necessary to use structural optimization techniques for a compromising primary design scheme, in order to solve issues related to the coupling between different disciplines, e.g., structural, aerodynamics, control, etc., and meeting aeroelastic performance requirements. A design under aeroelastic constraints is called an aeroelastic optimization. In the overall design of aeroelastic problems, the primary issue for determining the influence of structural deformation on the aerodynamic characteristics is designing a reasonable stiffness distribution for the wing structure [1], and this is also an essential basis for subsequent designs [2]. To obtain more accurate results that satisfy the aeroelastic performance requirements, it is necessary to design a wing stiffness distribution using design optimization technology, which can preferably guide the selection of the structural scheme and the arrangement of the macrostructure stiffness and mass distribution [3].

Joined-wing airplanes have been widely explored and studied in many different disciplines since Wolkovitch introduced the concept in 1976 [4]. A joined-wing airplane can be defined as an airplane that incorporates tandem wings arranged to form diamond shapes in both the plan and front views, and the fuselage can be seen as the connecting diagonal of the diamond frame. The advantages of a joined wing include light weight, high stiffness, low induced drag, good transonic area distribution, high-trimmed CLmax (maximum trimmed lift coefficient), reduced wetted area and parasite drag, direct lift control, side force control capability, excellent stability and controllability [5]. However, the connection between the front wing and the rear wing leads to dissimilarity in the structural and aerodynamic characteristics compared with a traditional layout, and the interconnected wings form a complex overconstrained system, which increases the difficulty of analysis and increases the design space of different disciplines [6,7]. This leads to integrated design changes becoming the key problem, as the authors in [8] noted that basic aeroelastic investigations must be introduced early in the design process of a joined-wing aircraft, and a stiffness distribution design is also necessary.

Aeroelasticity research has been accompanied by the development of the joined wing. Robert A. Canfield and his team performed a series of studies on the aeroelasticity of a joined wing during the first decade of the 21st century. Different structural models of joined wings were studied, including high-fidelity finite elements method (FEM)-based weight models [8], embedded antenna models [9,10], nonlinear structure models [11,12], sensor craft models [13,14], scale models [15,16,17,18], a beam model [19], and a wind tunnel model [20]. Nearly half of the research dealt with linear structural integrated design and optimization, whereas the other half dealt with the nonlinearity of the joined-wing models. Demasi et al. also performed great research on the aeroelasticity of joined wings. Their research focused mainly on the nonlinear problems of joined wings. The evaluation of aeroelastic characteristics of the joined wing was undertaken in the early design stages in [21]; both reduced-order models [22,23] and full-order models [24] were used to analyze the nonlinearity of joined wings; and post-critical phenomena were researched in [25,26]. However, the aeroelastic optimization of the models has hardly been taken into account. Preliminary analyses on joined wings showed that aerodynamic loads at the tip of the wing were sensitive to the modeling with reference to efforts [27,28]; moreover, aeroelastic analyses showed large differences in the predicted flutter speeds. Thus, it was speculated that by attributing the forces in the tip regions largely to the bending moment, the overall aeroelastic response was significantly affected by the redistribution. Flutter involving vehicle motion can be an active constraint [29].

As mentioned above, the aeroelastic analysis of a joined wing is necessary, and as stated earlier, optimization in the preliminary stage of aircraft design is also indispensable; therefore, the aeroelastic optimization of a joined wing is important. A large number of aeroelastic optimizations of joined wings have been explored in precursory studies. The rib thicknesses of a scaled joined wing was studied as a variables, and the scaled natural frequencies and scaled flutter speed of the full-scale vehicle were reproduced after optimization in [14]. An equivalent unsymmetrical beam section was used for structural optimization and was proven to adequately approximate the stiffness and deflection behavior of a real wing in [30]. Aerostructural optimization was carried out in [31,32]. A panel method and an equivalent beam finite-element model were used for aerostructural optimization with the aerodynamic design variables as the parameters of the geometric configuration, and the thicknesses of the spar wall were added as structural design variables. The application of bacterial foraging optimization to a joined wing was shown in [33] and a hybrid variant was applied to match the aeroelastic responses of a wind tunnel model with a full-scale aircraft.

Compared with traditional airplanes, structural optimization and weight estimations for the preliminary design of a joined-wing aircraft lack mature methods. In engineering, traditional wing-box design procedures and empirical formulas for traditional wing stiffness distributions have been frequently used for the preliminary design of a joined wing. However, this approach has three problems. Firstly, as it is a complex overconstrained system, similar external geometric parameters in the layout lead to completely different thicknesses of the internal structures and thus weights. Secondly, to minimize the weight by fixing the topology of the wing box, the material needs to be redistributed differently from what is usually seen in common layouts. As a result, the secondary bending moment tilts the neutral axis, and thus material is more efficiently utilized when distributed far apart to the two opposite corners. Thirdly, because of the interaction of the front wing and the rear wing, the distribution of the bending moment of the joined wing is different from traditional wings; thus, the distribution of the stiffness is different. For these reasons, optimization approaches in regard to stiffness based on traditional wings are not applicable to joined wings.

In the conceptual and preliminary design stages, there are many uncertain parameters. Stiffness design requires a large amount of relevant data, such as airfoil, load, structural layout, etc., and this information is often obtained through parallel or subsequent design processes. This lack of information will lead to a decrease in the accuracy of the stiffness design, which will require repeated modifications of the structural stiffness distribution in the subsequent design, resulting in a waste of manpower and time. Furthermore, parametric analysis of the joined wing is even more difficult due to its complex layout compared with the traditional airplane. This leads to the low design efficiency of the detailed finite-element models and nonlinear structure models. Therefore, it is necessary to introduce a simple stiffness model that can describe joined-wing structural characteristics directly and efficiently, so as to provide guidance for the subsequent structural scheme screening, macro-structural layout, and mass distribution.

It can be seen from the above that there are few direct studies on the stiffness distribution of joined-wing structures; however, it is necessary to conduct an in-depth study in the preliminary design stage, while considering the aeroelastic problem. In order to solve this problem, this paper introduces three different aeroelastic optimization methods in the preliminary stage to design the stiffness of the joined wing.

2. Materials and Methods

2.1. Beam Cross-Section Simplification

As shown in Figure 1, the out-of-plane components of a joined wing aircraft tend to bend the wings about a tilted bending axis. To resist this, the wing’s structural material must form a deep spar about this axis. This implies that the material must be concentrated near the upper leading edge and the lower trailing edge. As shown in Figure 2, to simulate the real material distribution of the beam section of a joined wing, three section shapes are simplified to facilitate calculation without losing rationality. Through a series of theoretical derivations, the best scheme for a joined wing is selected from the three section shapes. As shown in Figure 2a, section parameters $X_1$–$X_4$ are used for the design parameters in each section of the beams. The parameter $X_1$ represents the main diagonal flange size; $X_2$ represents the secondary diagonal flange size; $X_3$ represents the web thickness; $X_4$ represents the skin thickness; and $X_5$ and $X_6$ represent the width and height of the wing box, respectively, and they are fixed in each beam section. The symbol α represents the tilted angle, as shown in Figure 1 and Figure 2.

As shown in Figure 2a, the equations used to evaluate the principal inertia moments $I_{z1}$ and $I_{y1}$, the torsional stiffness coefficient $I_{t1}$, and the inertial product $I_{yz1}$ of the cross section are as follows:

$$I_{z1}=\frac{1}{2}{A}_t{X}_6^2+\frac{X_3(-{\overline{X}}^3-2X_6^3+3X_6{\overline{X}}^2)}{6}$$

(1)

$$I_{y1}=\frac{1}{2}{A}_t{X}_5^2+\frac{X_4(-{\overline{X}}^3-2X_5^3+3X_5{\overline{X}}^2)}{6}$$

(2)

$$I_{t1}=\frac{4{\mathsf{\Omega }}^2}{{\displaystyle \oint \frac{\mathrm{d}s}{\delta }}}\approx \frac{2X_3{X}_4{(X_5{X}_6)}^2}{X_4{X}_6+X_3{X}_5}$$

(3)

$$I_{yz1}=\frac{X_5{X}_6}{2}(A(X_2)-A(X_1))$$

(4)

where:

$$\begin{array}{l}A\left(X_1\right)=X_1^2,A\left(X_2\right)=X_2^2,A\left(X_3\right)=X_3{L}_w,A\left(X_4\right)=X_4{L}_s\\ A_t={\displaystyle \sum _{i=1}^{4}A\left(X_i\right)},L_s=X_5-\overline{X},L_w=X_6-\overline{X},\overline{X}=\frac{X_1+X_2}{2}\end{array}$$

(5)

In general, the flange size is larger than the size of the web thickness and skin thickness, and as mentioned earlier, to simulate the material distribution, the main diagonal flange size is larger than the secondary diagonal flange size shown in Figure 2a, so the relations of the cross section parameters are as in Equation (6):

$$\max (X_3,X_4)<X_2<X_1$$

(6)

The section area can be calculated as:

$$S=2A_{\mathrm{t}}$$

(7)

According to the rotation axis formula of the moment of inertia in Equation (8), the moment of inertia of a beam section relative to the bending axis u₁ in Figure 2a can be derived as Equation (9):

$$I_u=\frac{I_z+I_y}{2}+\frac{I_z-I_y}{2}\cos 2\alpha -I_{yz}\sin 2\alpha$$

(8)

$$I_{u1}=(1-{\alpha }^2)I_{z1}+{\alpha }^2{I}_{y1}-2\alpha I_{yz1}$$

(9)

For the case of Figure 2b:

$$I_{z2}=I_{z1}{I}_{y2}=I_{y1}{I}_{yz2}=0$$

(10)

then

$$I_{u2}-I_{u1}=2\alpha I_{yz1}<0$$

(11)

For the case of Figure 2c, the four flange areas are equal. Defining the flanges $X$ (where $X$=$X_1$=$X_2$) as variables and keeping parameters $X_3$–$X_6$ and $A_t$ unchanged, according to Equations (3)–(5):

$$I_{z3}=\frac{1}{2}{A}_t{X}_6^2+\frac{X_3(-X^3-2X_6^3+3X_6{X}^2)}{6}$$

(12)

$$I_{y3}=\frac{1}{2}{A}_t{X}_5^2+\frac{X_4(-X^3-2X_5^3+3X_5{X}^2)}{6}$$

(13)

$$I_{yz3}=0$$

(14)

$$I_{u3}=(1-{\alpha }^2)I_{z3}+{\alpha }^2{I}_{y3}-2\alpha I_{yz3}$$

(15)

By a series of derivations, as shown in Appendix A Equations (A15)–(A19), the result is:

$$\mathsf{\Delta }{I}_u=I_{u1}-I_{u3}=X_6(\alpha X_5(X_1^2-X_2^2))+X_3({\overline{X}}^2-X^2)$$

(16)

By a series of derivations, as shown in Appendix A Equations (A20)–(A25), when the Equation (17) condition is met, the result $\mathsf{\Delta }{I}_u>0$ holds.

$$N>\frac{(1-\frac{k_2^2}{4})(k+1)}{(k-1)k_2^2\alpha }$$

(17)

where:

$$N=\frac{X5}{X3},k=\frac{X1}{X2}>1,k_2=\frac{X3+\mathrm{X}4}{2X}$$

(18)

That is, if the above formula is satisfied, the moment of inertia relative to the bending axis of the case in Figure 2a is greater than that in Figure 2c when the section areas are the same. For thin-walled wing boxes, the above formula is generally satisfied, which means that when the areas of the sections are the same, i.e., when the weights of the beams are the same, the section in Figure 2a has a larger resistance to the tilted bending axis. In other words, a joined wing using the simplified beam section in Figure 2a requires less weight when the stiffness is the same; therefore, the section layout in Figure 2a is the best scheme in theory, so it is used for optimization.

2.2. Variable Model Division

Three methods are applied to solve the problem of the aeroelastic optimization of a joined wing, including two methods for a beam-frame model and one for a three-dimensional model, which are used for comparative analysis. The division of the variables of the models are described in detail below.

As shown in Figure 3, the joined wing is divided into three partitions, i.e., the front wing, the outer wing, and the rear wing. The main beam of each partition is then subdivided into several sections. More specifically, as shown in Figure 4 in different colors, there are three design sections between the fuselage and the kink point both in the front wing and the rear wing and six design sections between the kink point and the joint. The outer wing is divided into four design sections. For each section of the beam of the beam-frame model, two methods are used to address its stiffness, namely, method A and method B (The two methods will be discussed in detail in Section 2.3).

Method C is a traditional method to solve aeroelastic optimization problems with detailed scales of structure, and it is applied to the three-dimensional model, as shown in Figure 5. The wing box is divided into several sections in the same way as the beam-frame model, as shown in Figure 6. The difference is that the design variables of this method are the real structural sizes of the wing box, e.g., the thicknesses of skins, webs, and flanges.

2.3. Aeroelastic Optimization of the Global Stiffness

An aeroelastic optimization study is concerned with a standard optimization problem, which minimizes an objective function subject to constraints that limit the values of the design variables in n_d dimensional space. The problem can be written as [34]:

$$g_j(v)\le 0\left(j=1,2,\cdots ,n_{\mathrm{c}}\right)$$

(19)

$${(v_i)}_{\mathrm{lower}}\le v_i\le {(v_i)}_{\mathrm{upper}}(i=1,2,\cdots ,n_{\mathrm{d}})$$

(20)

where $g_j(\mathit{v})$ is a correlation function of the constraint index; $\mathit{v}$ is a vector set of the design variables; $v_i$ is a single design variable; ${(v_i)}_{\mathrm{lower}}$ is the lower bound of the design variables; ${(v_i)}_{\mathrm{upper}}$ is the upper bound of the design variables; n_c is the number of constraints; and n_d is the number of design variables. In aeroelastic design problems, the objective function $F\left(v\right)$ is generally mass, that is, to find the set of design variables that meet the conditions $v_i$ to minimize the total mass of the structure. The constraint conditions are generally static aeroelastic and dynamic aeroelastic constraint indexes, such as flutter speed, aileron efficiency, structural deformation and stress, which are constrained by Equation (19). The upper and lower limits of each design variable are constrained by Equation (20).

An optimization algorithm based on the basic genetic algorithm is introduced for aeroelastic optimization problems [2].

2.3.1. Method A

As shown in Figure 2a, section parameters $X_1$–$X_4$ are used for the design parameters in each section of the beams. Let the main diagonal flange size be $X_1$, the secondary diagonal flange size be $X_2$ ($X_1$>$X_2$), the web thickness be X₃, and the skin thickness be $X_4$. The width and height of the wing box are $X_5$ and $X_6$, respectively. The section parameters $X_3$ and $X_4$ of each section of the beam and the overall ratios $K_1$ and $K_2$ are the design variables as follows:

$$\begin{array}{l}{K}_1=X_1/X_2\\ K_2=2X_2/(X_3+X_4)\end{array}$$

(21)

There are 40 total design variables. Section parameters $X_5$ and $X_6$ are invariant in each section but vary with the spanwise direction of the beams. Generally, the overall ratios are larger than 1. In particular, when we make the design variable $K_1$ equal to 1, the situation then changes into a beam section, as shown in Figure 2c.

2.3.2. Method B

Referring to the stiffness distribution law of the traditional high-aspect-ratio wing, the stiffness distributions along the spanwise direction of the front wing, outer wing, and rear wing are optimized according to an exponential function as follows [2]:

$$GJ(y),EI(y)=a{e}^{by}+c$$

(22)

A set of three design variables (i.e., $a$, $b$, and $c$) can be defined for each kind of stiffness. Thirty-six design variables are required because the vertical bending, horizontal bending, torsional stiffness, and stiffness caused by section asymmetry of each partition need to be designed.

Due to the difference between a joined-wing shape and a traditional high-aspect-ratio wing, a traditional exponential stiffness distribution is not appropriate for a joined wing. As noted in previous papers, the moment of a joined wing along the spanwise direction decreases first and then increases in the middle part. To have a lighter weight, the distribution should have the same trend. Therefore, a new function needs to be found to simulate the distribution of a joined wing. As shown in Figure 6, an exponential/linear mixed function can change with its parameters and simulate the trend of the stiffness of a joined wing well, so it is used to modify the distribution as follows:

$$GJ(y),EI(y)=A{e}^{By}+C+Dy$$

(23)

A set of four design variables (i.e., $A$, $B$, $C$, and $D$) can be defined for each kind of stiffness, similarly to method B. Equation (23) is used to address the optimization problem because including a linear term Dy can simulate the transition of the stiffness curve of a joined wing well. Forty design variables are required because the vertical bending, horizontal bending, torsional stiffness, and stiffness caused by the section asymmetry of each partition need to be designed.

2.3.3. Method C

A three-dimensional model, as shown in Figure 7, is employed. The model is divided into the same three parts as the previous beam-frame model, i.e., the front wing, outer wing, and rear wing. The thickness values of the flanges, webs, and skins in each partition are used individually for the aeroelastic optimization using a genetic algorithm. The details of this method are presented in [2].

2.4. Aeroelastic Responses

The methods that deal with the aeroelastic responses involved in this paper are described as follows.

2.4.1. Equilibrium Equation for Static Aeroelastic Response Analysis

In general, the static aeroelastic equation is established under the displacement vector $f-set$, expressed as [35]:

$$[K_{ff}-\overline{q}{Q}_{ff}]u_f+M_{ff}{\ddot{u}}_f=\overline{q}{Q}_{fx}{u}_x+{\mathit{P}}_f$$

(24)

where $K_{ff}$ represents the structural stiffness matrix; $\overline{q}$ indicates the dynamic pressure; $u_f$ represents the displacement vector; $M_{ff}$ represents the structural mass matrix; $M_{ff}{\ddot{u}}_f$ represents the inertia force caused by the rigid body motions; $P_f$ represents the external load vector; $\overline{q}{Q}_{ff}{u}_f$ represents the aerodynamic increment caused by structural elastic deformation; and $\overline{q}{Q}_{fx}{u}_x$ represents the aerodynamic force caused by rudder deflection and rigid body motion of the aircraft. Superscript “..” is the quadratic derivative of time.

2.4.2. Flutter Function

There are generally three common flutter analysis methods, the $V-g$ method and $p-k$ method and the matrix eigenanalysis method [36], which are all computationally efficient and suitable for flutter analysis in optimization problems. The $p-k$ method is used in this paper. The equation is established under the modal set h [37]:

$$[{\left(\frac{V}{b}\right)}^2{p}^2{\mathit{M}}_{hh}+\frac{V}{b}p{\mathit{B}}_{hh}+{\mathit{K}}_{hh}-\frac{1}{2}\rho V^2\left({\mathit{Q}}_{hh}^R+\frac{p}{k}{\mathit{Q}}_{hh}^I\right)]{\mathit{u}}_h=0$$

(25)

where $V$ is the incoming flow velocity; $b$ is the reference half chord length; $p$ is the eigenvalue; $k$ is the reduced frequency; ${\mathit{M}}_{hh}$ is the modal mass matrix; ${\mathit{B}}_{hh}$ is the damping matrix; ${\mathit{K}}_{hh}$ is the modal stiffness matrix; ${\mathit{Q}}_{hh}$ is the aerodynamic force matrix; and the superscripts $R$ and $I$ indicate the real and imaginary parts, respectively.

3. Model Description

3.1. Beam-Frame Structure Model

A half-model of a high-aspect-ratio joined wing is used as the research object, and detailed design parameters of the wing are given in

Table 1. Aeroelastic constraint conditions.

Constraints	${\mathit{d}}_{\mathbf{t},\mathit{z}}$	${\mathit{d}}_{\mathbf{t},\mathit{x}}$	${\mathit{d}}_{\mathbf{j},\mathit{z}}$	${\mathit{d}}_{\mathbf{j},\mathit{x}}$	${\mathit{\varphi }}_{\mathit{t}}$	${\mathit{\varphi }}_{\mathit{j}}$	${\mathit{V}}_{\mathit{f}}$
Upper limit	7.5% × lt	1.5% × lt	7.5% × lj	1.5% × lj	2°	2°	-
Lower limit	-	-	-	-	−2°	−2°	90 (m/s)

where $d_{\mathrm{t},z}$ is the vertical deformation of the wing tip, $d_{\mathrm{t},x}$ is the horizontal deformation of the wing tip, $d_{\mathrm{j},z}$ is the vertical deformation of the joint, $d_{\mathrm{j},x}$ is the horizontal deformation of the joint, ${\varphi }_t$ is the twist angle of the wing tip, ${\varphi }_j$ is the twist angle of the joint, $V_f$ is the constraint flutter speed, l_t is the half wing span, and l_j is the length between the wing root and the joint.

. The structural model of the joined wing in the preliminary design stage is shown in Figure 3. It is composed of 270 nodes and 367 elements, including 93 mass point elements and 274 rod elements, which is established using Patran. The wing stiffness is simulated by an elastic axis from the wing root to the wing tip, and the mass characteristics are simulated by the lumped mass. An equivalent material was used by referring to the average values for aluminum alloys for the main beams of the wings. Its density ${\rho }_{alu}$ = 2700 kg/m³, Young’s modulus E = 72 Gpa, and Poisson’s ratio $\mu =0.33$. All structures except the beams, e.g., the fuselage, ribs, and the fin and joint plate, are rigid.

3.2. Three-Dimensional Structure Model

The wing is divided into panels, spars, ribs, and beams based on structural arrangement and dimension data. The three-dimensional structure model is shown in Figure 7. In accordance with the load-bearing characteristics of the various components, the entire wing is simulated by rod-shell elements. The upper/lower panels, wing ribs, and webs of the front/rear beams are simulated by shell elements. The spars and flanges of the front/rear beams are simulated by rod elements. The body, joint, and vertical tail are simulated by rigid beams. Although a large number of elements means that the results are more accurate, a moderate number of elements is more suitable for parameterization and optimization in the conceptual and preliminary stages, considering time-consumption. The entire wing model is composed of 8247 finite elements, which is sufficiently accurate and efficient in the current design stage.

3.3. Aerodynamic Model

The geometric shape parameters are shown in Figure 8, and the values of the parameters are shown in

Table 2. Values of the shape parameters of the joined-wing layout.

Shape Parameter	Value (m)	Shape Parameter	Value (m)	Shape Parameter	Value (°)
$l$	22.86	$b_r$	2.54	${\chi }_f$	22.0
$l_k$	4.63	$b_{rk}$	1.49	${\chi }_r$	35.0
$l_j$	15.67	$b_{rj}$	1.49	${\phi }_f$	4.0
$b_f$	3.51	$b_t$	0.75	${\phi }_r$	0.0
$b_{fk}$	2.11	$h$	1.19
$b_{fj}$	2.11

$l$ is the half span of the joined wing; $l_k$ is the length between the wing root and the kink; $l_j$ is the length between the wing root and the joint; $b_f$,$b_{fk}$, $b_{fj}$ are the chord lengths of the wing root, kink, and joint of the front wing, respectively; $b_r$, $b_{rk}$, $b_{rj}$ are the chord lengths of the wing root, kink, and joint of the rear wing, respectively; $b_t$ is the chord length of the wing tip; $h$ is the height of the joint plate; ${\chi }_f$ is the leading back swept of the front wing and outer wing; ${\chi }_r$ is the edging forward swept of the rear wing; ${\phi }_f$ is the dihedral of the front wing and outer wing; ${\phi }_r$ is the dihedral of the rear wing.

. A flat aerodynamic mesh considering curvature and the double-lattice method was applied to the aerodynamic model, as shown in Figure 9. The total area of the half-model is 72.76 m², among which the front wing is 46.51 m², whereas the rear wing is 26.25 m².

4. Results and Discussion

4.1. Aeroelastic Responses of the Optimized Models

The constraint conditions of all the methods are the same, including the static aeroelastic responses and flutter constraint. The specific constraint values are shown in Table 1.

The displacement and torsion angle of the wingtip and the joint are the most concerning static aeroelastic responses, and a comparison of those values is given in

Table 3. Aeroelastic responses of the optimized models.

Responses	Method A	Method B	Method C	Original
$d_{\mathrm{t},z}/l_{\mathrm{t}}$	7.49%	7.50%	7.50%	13.2%
$d_{\mathrm{t},x}/l_{\mathrm{t}}$	1.30%	1.15%	0.30%	2.23%
$d_{\mathrm{j},z}/l_{\mathrm{j}}$	5.64%	5.63%	6.48%	10.1%
$d_{\mathrm{j},x}/l_{\mathrm{j}}$	1.09%	1.01%	0.41%	2.05%
${\varphi }_{\mathrm{t}}$/°	0.99	0.50	0.08	1.74
${\varphi }_{\mathrm{j}}$/°	−0.05	−0.26	−0.64	−1.43
$V_F$	>90 (m/s)	>90 (m/s)	>90 (m/s)	>90 (m/s)

“Method A”, “Method B”, and “Method C” represent the aeroelastic responses of the models optimized using different optimization methods, respectively; “Original” represents the aeroelastic responses of the original model; $V_F$ is the actual flutter speed.

. The flight condition is defined at Ma = 0.2, a 1.5 g overload and 25,000 m cruise altitude, where the acoustic velocity $V_a$ = 298.4 m/s, the air density ${\rho }_{air}$ = 0.04 kg/m³. The cruise speed $V_{cr}$ = 59.7 m/s, and the dynamic pressure $q$ = 71.28 Pa. The aeroelastic deformation of the beam-frame model optimized using method A is shown in Figure 10. The deformation of the other beam-frame model optimized using method B is similar to what is shown in Figure 10, so it is not included to avoid repetition. The aeroelastic deformation of the 3-D model optimized using method C is shown in Figure 11. Figure 12 depicts one of the $V-g$ and $V-f$ diagrams obtained based on the $p-k$ method. The other flutter result diagrams are not included to avoid repetition. All the flutter speeds of the optimized models are larger than the constraint speed. The flutter constraint $V_f$ is larger than 90.0 m/s, which seems very low; however, as the designed cruise speed $V_{cr}$ is only 59.7 m/s, and $V_f$ is larger than 1.5 times $V_{cr}$, the constraint speed is reasonable.

The specific data of the aeroelastic constraints are listed in Table 3. From the table, we can see that compared with the original model, the constraints are satisfied well after optimization by all the methods, with only the vertical deformation of the wingtip reaching the constraint boundary value, and the other constraints having significant margins, so the optimization result is reasonable. It should be noted that the horizontal deformation of the wingtip is not negligible.

4.2. Results of the Optimized Parameters

The width and height of the beam sections, i.e., $X_5$ and $X_6$ respectively, are fixed and presented in Figure 13. The span station is nondimensional. The front wing and the outer wing are continuous in the span direction and the span station is from the wing root to the wing tip, whereas the span station of the rear wing is from the wing root to the joint. The descriptions of the nondimensional span station in later figures are the same.

According to Equations (A3)–(A7), the parameters of the beam sections, i.e., $X_1$–$X_6$, are positively correlated with the torsional stiffness coefficient $I_{t1}$ and inertias $I_{z1}$ and $I_{y1}$. Because the scales of the fixed parameters, i.e., $X_5$ and $X_6$, are far greater than the design parameters, i.e., $X_1$–$X_4$, the inertias are primarily determined by the fixed parameters, so when the fixed parameters change, the inertias have the same trend; thus, we mainly consider the middle parts between the kink and the joint point where the fixed parameters are invariant with the span direction.

All the optimizations were completed in 20 generations, and there were 600 individuals in each generation. In theory, as mentioned in Section 2.1, the diagonal beam section is the best-simplified beam section for the stiffness optimization of a joined wing. As shown in Figure 2a, the parameters of each beam section are designed and optimized using a genetic algorithm in method A. The beam section parameters are solved in a reverse fashion using Equations (1)–(9) in method B. On one side, the areas of the beam sections are calculated, then the volume and mass of the beam are considered the objective, and on the other side, those beam section parameters can be used as a reference for a real scale of the beam sections, which is convenient for redesign work in the primary stage of joined wing design. The parameters of method C are the scales of the skins, flanges, and webs, which are in accordance with method A.

Figure 14 and Figure 15 show the beam section parameters of the front wing and the rear wing, respectively. Based on these results, we can see that all four parameters of the three methods have the same tendency, i.e., in the inner part before the kink of both the front wing and rear wing, these parameters gradually increase with the span direction; then, in the middle part between the kink and the joint point of both the front wing and rear wing, these parameters first decrease and then increase up to the joint point; finally, in the outer part, they gradually decrease with the span direction.

Figure 16 shows two objective results; one is from the optimization using a traditional exponential function, whereas the other is from the optimization using a modified exponential/linear mixed function. We can see that the modified function results in a much lower value for the objective result, which means that when obtaining the same stiffness level, less weight is needed using the new function for optimization.

4.3. Results of the Optimized Stiffness Distribution

To test the accuracy of Methods A and B, the optimized stiffness distribution, including the vertical bending stiffness, the horizontal bending stiffness, and the torsional stiffness of both the front wing and the rear wing were calculated and the results are shown in Figure 17 and Figure 18.

We can see from the figures that the stiffness distributions obtained using the two methods have good consistency and are similar to the distribution of the three-dimensional model except for the first span station point. The reason for this error is that the first point is near the root of the wing, where the ratio of the wingspan to the width of the wing box is relatively small; therefore, the theoretical calculation of an engineering beam is not accurate. Furthermore, the rigid section assumption is adopted for the three-dimensional model, and the wings of the three-dimensional model are connected to the fuselage through rigid multipoint constraint elements. This results in a stiffness increase at the wing root. However, the results have little influence on the overall stiffness distribution and can be ignored in the preliminary design; thus, the methods are valuable for engineering applications. In practical applications, it should be noted that the design parameters of method A and method C are cross-section parameters of the wing, which are relatively intuitive and easy to adjust, whereas the design parameters of method B are parameters of mathematical expressions of stiffness, which do not directly relate to the structure, so this requires conversion between the stiffness and the actual structure.

5. Conclusions

In this study, two beam-frame model methods and one three-dimensional model method are proposed for the global stiffness design and optimization of joined wings. A beam-frame model of a joined wing, as an example for use during the preliminary design stage, is established. A three-dimensional detailed model is built as a reference and introduced for comparison. Two methods for the beam-frame model, an optimization of the diagonal beam section, and an optimization method of the exponential/linear combination function fit, defined as method A and method B, respectively, are used. The method for the three-dimensional model is a traditional optimization method based on a genetic algorithm, and is defined as method C. The static aeroelastic responses are calculated and compared under the same flight conditions at 0.2 Ma and 25 km above sea level. The main conclusions are as follows:

• The stiffness distributions of the beam-frame model obtained from method A and method B have good consistency and are similar to the distribution of the three-dimensional model of method C;
• The stiffness distribution data and the beam section parameters obtained using method B can be converted into parameter sizes of the beam section and conform with the trend of the other methods;
• Method B (the optimization method of the exponential/linear combination function fit) obtains a greater weight reduction than the traditional method (the optimization method of the exponential function fit) under the same stiffness conditions;
• Compared with method C and other detailed-model-based methods, methods A and B show high computational efficiency and are easy to implement with variable parametrical analysis of the joined wing.

In our future work, we aim to perform joint connection analysis, including the influence of fixed connections, hinged connections, and others, especially the gap of the joint between the front wing and rear wing. The gap nonlinearity of the joined wing will express many interesting problems and will have great research value.