#### Abstract

High aspect ratio wing (HARW) structures will deform greatly under aerodynamic loads, and changes in the stiffness will have a great impact on the flutter characteristics of such wings. Based on this, this paper presents an effective method to determine the effect of the stiffness on the flutter characteristics of HARWs. Based on the calculation theory of the mechanical profile of thin-walled structures, the torsional stiffness and bending stiffness of the wing are obtained through calculation. We use the fluid-structure coupling method to analyze the flutter characteristics of the wing, and we use our research results based on the corotational (CR) method to perform structural calculations. The load is calculated using a computational fluid dynamics (CFD) solver. The results show that, compared with the original wing, when the bending stiffness and torsional stiffness of the wing along the spanwise direction increase by 8.28% and 5.22%, respectively, the amplitude of the flutter decreases by approximately 30%. Increasing the stiffness in the range of 0.4 to 0.6 Mach has a greater impact on the flutter critical velocity, which increases by 12.03%. The greater the aircraft’s flight speed is, the more severe the stiffness affects the wing limit cycle oscillation (LCO) amplitude.

#### 1. Introduction

The effect of stiffness needs to be factored into the design of a high aspect ratio wing (HARW). The flutter of the wing is severely restricted by the bending and torsional stiffness of the wing’s cross section. Through a large number of analyses and experimental results on a two-dimensional airfoil, it has been shown that if the torsional and bending stiffnesses of the wing are simultaneously increased by *n* times, the flutter critical velocity of the wing will be accordingly increased by times [1]. Researchers both domestically and abroad have done little research on the effect of stiffness on the aeroelasticity of the wing and have not conducted any systematic research. G. Ortiz-Torres et al. [2] carried out systematic research on the safety of the vertical take-off and landing UAV systems.

In current stiffness computations, many researchers first simplify the model using the related theory of thin-walled structural mechanics and then program the simplified model to calculate the wing cross-sectional stiffness center [3–6], which is efficient and useful in engineering designs. Others have studied the effect of the position of the elastic axis on aerodynamics [7]. Based on the optimization design, a blade cross-sectional internal layout for helicopters was designed to maximize the static twist actuation while satisfying the locations of the center of the elastic axis [8, 9]. The effect of the elastic axis was investigated, with the elastic axis close to the airfoil midpoint, and the pitch amplitude increases slightly, while the plunge amplitude fluctuates with an extremum [10].

In terms of flutter calculations, many use nonlinear computational structure dynamics (CSD) solvers to perform the structural calculations, and Dang [11] used the loose coupling computational fluid dynamics (CFD)/CSD method to analyze the static aeroelastic characteristics of wings with an HARW. The research results show that the effect of geometric nonlinearity on such wings cannot be ignored. Similarly, taking the HARW as the research object and using the loose coupling CFD/CSD method, Wang et al. [12] established an efficient model for the analysis of aeroelastic characteristics. Demasi et al. [13] established an unsteady aerodynamic model in the frequency domain and combined the updated Lagrange (UL) method to conduct aeroelastic coupling research on composite thin-walled structures. Some researchers have taken the delta wing as the research object, using the large-deformation structure model and the linear fluid model to conduct limit cycle oscillation (LCO) phenomenon research [14–16]. Based on the structural dynamics equations and Navier-Stokes (N-S) equations, Yang et al. [17] developed a fluid-structure coupling solver whose purpose is to study the effect of aerodynamics between the structure and the aerodynamic model. Fairuz et al. [18] used a mesh-based parallel code coupling interface (MPCCI) as the coupling platform and the CFD and CSD solvers that came with the platform to study the effect of the interaction between the fluid and structure on the flap wing deformation. Based on the same time step and combined with the nonlinear finite element method, some researchers have developed an aeroelastic coupling calculation program [19]. Some researchers have introduced a coupling method to predict the response of a three-dimensional wing. This coupling method includes modal methods, structural modal equations, and N-S equations [20]. To conduct an aeroelastic analysis and consider the computational efficiency at the structural level, a 3D shell element was established by introducing the CR method [21–23].

Although the above research results have achieved their expected purpose, when considering the application of these research results in actual engineering, it is urgent to further consider the issue of computational efficiency. To account for computational efficiency, the CR method has been developed. Currently, to conduct aeroelastic analyses of HARW, many researchers have developed nonlinear spar [24, 25], rod [26], and shell [27–30] elements using the CR method. Kirsch et al. [31] also developed a calculation program specifically for the aeroelastic analysis of flexible aircraft. Barnes et al. [32] conducted a study on the effect of stiffness on the laminar separation flutter based on an NACA0012 airfoil and considered the effect of the Reynolds number. The results show that the change in stiffness at any Reynolds number will have a greater impact on the laminar separation chatter, resulting in even more nonlinear aeroelastic responses. After studying the effective stiffness calculation method and flutter calculation method separately, we also considered the serious effect of stiffness on the flutter of HARWs. Therefore, in this paper, an effective and practical stiffness calculation method and an efficient nonlinear flutter calculation method are used to study the effect of stiffness on the flutter characteristics of an HARW. An effective method to determine the effect of stiffness on the flutter of HARWs is then proposed.

#### 2. Stiffness Calculation

Wing stiffness centers were calculated for different cross sections selected from root to tip. However, to use the stiffness calculation formula to calculate the wing cross-sectional rigidity, the cross-sectional area and material needed to be equivalent, as shown in [3].

In the stiffness calculation process, to facilitate programming, we used equivalent methods to process the area and material of the cross section of the wing according to the force characteristics of the structure. We equated the cross-sectional area of the spar flange and the stringer to the concentrated area and equated the materials of different structures to the same material. Then, the moment of inertia, the static distance, and the rigidity coordinates of the equivalent cross section were sequentially calculated. Finally, the stiffness calculation was performed according to the bending and torsional stiffness calculation formulas. The simplified cross section of the wing is shown in Figure 1.

According to coordinate transfer theory, stiffness center computation is conducted in the inertia axis. We define an inertial coordinate system for the simplified wing cross section, in which the wing chord length direction is defined as the *y*-axis, and the wing thickness direction is defined as the *z*-axis. Thus, the inertia moment of the wing cross section to the inertia axis can be obtained [3].where is the reduction factor, is the thickness of the skin, whose number is , is the segment number of the simplified cross-sectional skin, and is the sum area of all component cross sections. To calculate the cross-sectional area of the skin, *s* is used here to represent the length of the cross-sectional skin of the wing.

The centroid axis coordinate system is rotated around the origin; when it rotates to a certain angle, making the product of inertia equal to 0, the centroid axis coordinate system at this time is called the inertial axis coordinate system. Equation (1) is the moment of the reduced cross sections in the inertia axis; then, according to the related theory of thin-walled structural mechanics, the stiffness center of any section of the wing can be calculated according to the following equation:where represents twice the cross-sectional area and is the distance between the centroid and tangent of each component.

According to the geometric parameters and the number of spars, the wing is divided into three parts, as shown in Figure 2. The cross-sectional shape of the spar flange is a rectangle, and we use W and H to denote the width and height of the rectangle, respectively, as shown in Figure 3. The geometric parameters of the spar on the two wings are shown in Tables 1 and 2. The wing is a full composite wing with double spars. The unmanned aerial vehicle (UAV) has an aspect ratio of 12 and a wingspan greater than 17 meters, and the root chord length is more than twice the wing tip chord length. The thicknesses of the spar at the wing root and wing tip are 3 mm and 2 mm, respectively. The material and thickness of the skin between two adjacent ribs are different. From the wing root to the wing tip, the thickness of the skin and Poisson’s ratio of the material gradually decrease. The composite materials used on the wing are all laminated composite materials. The number of layers and fiber laying angles of the laminated composite material between two adjacent ribs are different.

From the wing root to the wing tip, the geometric parameters of the spar are unevenly distributed, so it is divided into three parts. The material parameters of the spar are shown in Table 3. When the geometric parameters of the three-part spars were increased by 16.7%, 25%, and 25%, the bending and torsional stiffness distributions of the two wings along the span were obtained as shown in Figures 4 and 5. There are 12 wing ribs from wing root to wing tip. We choose the cross-sectional position used to calculate the stiffness based on the cross-sectional position of the 12 ribs. When the geometric parameters of the spars increase simultaneously according to the ratio, although the wing bending and torsional stiffness distributions in the span direction are not uniform, there are relatively large changes in the aileron section. The torsional stiffness increment is between and N∗m/rad, and the bending stiffness increment is between and N/m.

#### 3. Flutter Analysis

The flutter calculation is performed for the wing with an enlarged stiffness and compared with the original wing flutter results. The degree of the change in stiffness’s effect on the flutter characteristics of the wing is studied.

##### 3.1. Flutter Analysis Model

We introduce the nonlinear model of [27] (Qiao) to perform the deformation analysis of the wing structure. The aerodynamic force is calculated using a CFD solver. The information exchange of the coupling surface is realized through the FLUENT external interface, which integrates self-compiled load and displacement interpolation programs. The main process of the wing flutter analysis is shown in Figure 6.

The flutter calculation is performed based on the CFD/CSD coupling method. In the calculation process, we use the FLUENT interface to integrate the User-Defined Function (UDF) program to realize the transfer of node information between the aerodynamic model and the structural model. The load interpolation program is then used to transfer the nodal loads on the aerodynamic model to the nodes of the structural model, and the displacement interpolation program is used to transfer the nodal displacements on the structural model to the nodes on the coupling surface of the aerodynamic model. We use the boundary layer grid function in GEMBIT software to address the boundary layer problems. We define the distance between the wall and the second layer as 0.00005 meters, and the gradient change from the wall to the outer boundary of the flow field is defined as 1.2. The CFD software FLUENT is used for the aerodynamic analysis. A coupled-implicit Spalart-Allmaras (S-A) turbulence model is used to perform CFD, and the N-S equations are used for the aerodynamic model. The time advance in the coupling process uses a loosely coupled method.

The expression of the flutter motion equation in the modal space is as follows:where , , and represent the mass matrix, damping matrix, and tangent stiffness matrix, respectively, represents the displacement vector, and represents the load vector. The damping matrix can be expressed as follows:where and represent the damping constants corresponding to the mass and stiffness, respectively.

Therefore, in the modal space, the flutter motion equation can be replaced by the state equation.wherewhere is the generalized displacement and is the time.

With this equation, we can use the fourth-order Runge-Kutta method to perform the calculations more conveniently; the calculation results are more accurate and can be implemented more conveniently by using programs. Finally, a generalized displacement calculation formula is formed.where

The finite element model and the flow field calculation model of the wing are shown in Figures 7 and 8, respectively, and the coupling surface grid model is shown in Figure 9. The entire computation zone containing 291,991 nodes and 1,724,960 elements is made of unstructured and tetrahedral elements for the aerodynamic analysis. The dimensions of the domain of the HARW are 100 m × 60 m × 60 m. There are 28679 nodes and 38968 elements within the structural model.

##### 3.2. Flutter Analysis of the Original Wing

Based on the modal method, the flutter calculation is performed under different flight speed conditions, and the calculation results are reflected by the history curve of generalized displacement with time. By comparing the results at different flight speeds, the flutter characteristics of HARW are studied.

In this paper, the first four modes are extracted for the flutter calculation. The frequencies of the first four modes increase in sequence. The first, second, and fourth modes are bending modes, and the third is a torsion mode. The first four modes of the wing are shown in Table 4. Using the halving method, through multiple calculations, it was determined that the time step used in this flutter calculation was 0.04 s [27]. Finally, the flutter calculations were conducted at flighting speeds of 0.4 Ma and 0.8 Ma. The response histories of the generalized displacements with time are shown in Figures 10 and 11.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

It can be seen in Figures 10 and 11 that, compared with the flutter response of the first-order mode, the effect of other third-order modes on flutter can be ignored, and the first-order mode is the most important mode for the flutter analysis. With an incoming Mach number of 0.4 Ma, the generalized displacement response of the wing over time reaches a state of constant amplitude oscillation when the flight speed is 215.3 m/s. Therefore, the flight speed of 215.3 m/s is the flutter critical speed of the wing at 0.4 Ma. Correspondingly, the flutter critical velocity of the wing at 0.8 Ma is 144.8 m/s. When the flight speed is less than the flutter critical speed, the generalized displacement amplitude of the wing gradually decays over time. In contrast, when the flight speed is greater than the flutter critical speed, the generalized displacement amplitude of the wing diverges with time.

##### 3.3. Flutter Analysis of Wings with an Enlarged Stiffness

To study the effect of stiffness on flutter, we increased the stiffness of the original wing. We still use the flutter analysis model in Section 3.1 to calculate the flutter of the wings with an increased stiffness. The flutter results at Mach numbers of 0.5 Ma, 0.8 Ma, and 0.9 Ma are shown in Figures 12–14, respectively.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

When the incoming flow Mach number is 0.5 Ma, the generalized displacement response of the wing over time reaches a state of constant amplitude oscillation when the flight speed is 228.9 m/s. Therefore, the flight speed of 228.9 m/s is the flutter critical speed of the wing at 0.5 Ma. Correspondingly, when the Mach numbers are 0.8 Ma and 0.9 Ma, the flutter critical speeds of the wing are 155.7 m/s and 145.8 m/s, respectively. Similarly, when the flight speed is less than the flutter critical speed, the generalized displacement amplitude of the wing gradually decays over time. In contrast, when the flight speed is greater than the flutter critical speed, the generalized displacement amplitude of the wing diverges with time.

#### 4. Effect of Stiffness on the Flutter Characteristics

Compared with the original wing, when the wing bending stiffness and torsional stiffness of each section along the span increase by an average of 8.28% and 5.22%, the flutter amplitude value decreases by approximately 30%. Under the condition that the Mach number is 0.8 Ma, the greater the stiffness is, the greater the flutter critical speed is, and the increase rate reaches 7.53% (the flutter critical speed is increased by 10.9 m/s).

To finally determine the change trend of the flutter critical speed of the large aspect ratio composite wing from the low speed to the subsonic speed range, we performed a comparative analysis between the initial stiffness wing and the enlarged stiffness wing flutter critical speed. The flutter analysis of two large aspect ratio composite material wings was performed from a low speed to a subsonic speed. The curve of the change of the flutter critical speed of these two wings with the flight speed is shown in Figure 15.

It can be seen from Figure 15 that, in the low speed to subsonic speed range, the flutter critical speed of the two wings changes slowly with the flight speed, but increasing the stiffness in the range of 0.4 to 0.6 Mach has a greater impact on the flutter critical velocity. The flutter critical speed increased by 25.9 m/s, and the amplitude increased by 12.03%. As the flight speed continued to increase, the magnitude of the reduction in the flutter critical speed of the two wings was increased compared with that at low speed, and the minimum flutter critical speed was reached at a flight speed of 0.9 Ma. In general, increasing the stiffness of the wing increases the flutter critical speed. The flutter critical speed of the two stiffness wings is closest in the range of 0.8 Ma to 0.9 Ma, and the flutter critical speed is increased by 7.53%.

Based on the established flutter analysis model considering the geometric nonlinearity of the structure, the numerical calculation of the limit cycle oscillation (LCO) characteristics of the two wings under different flight speed conditions was performed, and the variation curve of the LCO amplitude of the two wings with the flight speed is shown in Figure 16.

Figure 16 shows that the overall increase in the wing stiffness from the wing root to the wing tip reduces the LCO amplitude. At low flight speeds, the LCO amplitude/half-span lengths of the two wings are relatively close; the values differ by 0.015. As the flight speed continued to increase, the difference between the LCO amplitude/half-span lengths of the two wings gradually increased, and the two wings reached a maximum difference of 0.048 at 0.9 Ma. The increase in the flight speed increases the vertical displacement of the wing, which is in line with reality. As the flight speed increases, the effect of the increasing stiffness on the LCO amplitude of the composite wing with a large aspect ratio becomes increasingly obvious.

#### 5. Conclusions

Under different flight speeds, the flutter characteristics of the original wing and the enlarged stiffness wing were considered, particularly the geometrical nonlinear characteristics of the structure. By comparing the calculation results of the flutter critical speed, the degree of the effect of the stiffness on the flutter characteristics of a high aspect ratio composite wing was determined. Some useful conclusions are as follows:(1)In terms of the overall impact, the stiffness has a greater effect on the wing flutter critical speed. Increasing the wing stiffness overall increases the wing flutter critical speed.(2)Compared with the original wing, when the wing bending stiffness and torsional stiffness of each section along the span increase by averages of 8.28% and 5.22%, the flutter amplitude value decreases by approximately 30%. The flutter critical speed increased by 25.9 m/s, and the amplitude increased by 12.03%. Reasonably increasing the wing stiffness can effectively reduce the wing flutter, for example, increasing the size of the main bearing structure near the control surface.(3)The flutter critical speed is more affected by the lower flight speed, and the degree of the effect decreases slightly as the flight speed continues to increase. The flutter critical speed of the two stiffness wings is closest in the range of 0.8 Ma to 0.9 Ma, and the flutter critical speed is increased by 7.53%. With the continuous increase in the flight speed, the effect of stiffness on the flutter of the wing is reduced, but the effect on the LCO amplitude is greater.(4)In the structural design process, attention should be given to the above effects to achieve the purpose of reducing the structural quality. The research results provide effective key technical support for the comprehensive design of aeroelasticity and stiffness of composite wings with large aspect ratios.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was funded by the Scientific Research Foundation for Doctors, the Natural Science Basic Research Plan in Shaanxi Province of China (no. 2019JQ-912), and the Natural Science Foundation of Xi’an Aeronautical University (no. 2018KY1226).