Effect of Stiffness on Flutter of Composite Wings with High Aspect Ratio

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 thinwalled structures, the torsional stiffness and bending stiffness of the wing are obtained through calculation. We use the fluidstructure 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. +e load is calculated using a computational fluid dynamics (CFD) solver.+e 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%. +e greater the aircraft’s flight speed is, the more severe the stiffness affects the wing limit cycle oscillation (LCO) amplitude.


Introduction
e effect of stiffness needs to be factored into the design of a high aspect ratio wing (HARW). e flutter of the wing is severely restricted by the bending and torsional stiffness of the wing's cross section. rough 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 � n √ 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][4][5][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]. e 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. e 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][15][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. is 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][22][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][28][29][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. e 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. erefore, 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.

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 crosssectional area of the spar flange and the stringer to the concentrated area and equated the materials of different structures to the same material. en, 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. e 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. us, the inertia moment of the wing cross section to the inertia axis can be obtained [3].
where φ i is the reduction factor, t i is the thickness of the skin, whose number is i, i is the segment number of the simplified cross-sectional skin, and A 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 crosssectional skin of the wing. e 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. e 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. e geometric parameters of the spar on the two wings are shown in Tables 1 and 2. e wing is a full composite wing with double spars. e 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. e thicknesses of the spar at the wing root and wing tip are 3 mm and 2 mm, respectively. e 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. e composite materials used on the wing are all laminated composite materials. e 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. e 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. ere 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. e torsional stiffness increment is between 7.71 × 10 5 and 1.72 × 10 8 N * m/rad, and the bending stiffness increment is between 2.15 × 10 7 and 1.07 × 10 8 N/m.

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

Flutter Analysis Model.
We introduce the nonlinear model of [27] (Qiao) to perform the deformation analysis of the wing structure. e aerodynamic force is calculated using a CFD solver. e information exchange of the coupling surface is realized through the FLUENT external interface, which integrates self-compiled load and displacement interpolation programs. e main process of the wing flutter analysis is shown in Figure 6. e 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. e 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. e 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. e time advance in the coupling process uses a loosely coupled method.
e expression of the flutter motion equation in the modal space is as follows: where M, C, and K represent the mass matrix, damping matrix, and tangent stiffness matrix, respectively, d represents the displacement vector, and F(t) represents the load vector. e damping matrix can be expressed as follows: where α r and β r represent the damping constants corresponding to the mass and stiffness, respectively. erefore, in the modal space, the flutter motion equation can be replaced by the state equation. where where q(t) is the generalized displacement and t 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 e 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.  In this paper, the first four modes are extracted for the flutter calculation. e frequencies of the first four modes increase in sequence. e first, second, and fourth modes are bending modes, and the third is a torsion mode. e first four modes of the wing are shown in Table 4

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.

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).   Mathematical Problems in Engineering 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. e flutter analysis of two large aspect ratio composite material wings    was performed from a low speed to a subsonic speed. e 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. e 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. e 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/halfspan 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. e 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  Mathematical Problems in Engineering amplitude of the composite wing with a large aspect ratio becomes increasingly obvious.

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%. e 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) e 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. e 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. e research results

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.