Effect of Aeroelastic Tailoring Design on Wing Mode

,


Introduction
Aeroelastic problems have plagued aircraft engineers since the birth of the first aircraft. Aeroelasticity of aircraft refers to the interaction of aerodynamics, elastic forces, and inertia forces of the structure. With the rapid development of material science, the research and development of composite materials make the aeroelasticity problem of aircraft have a new development direction. The successful test flight of the X-29 demonstrator also marks the successful application of composite design technology. Compared with traditional metal materials, composite materials have many advantages such as high specific strength, high specific stiffness, good corrosion resistance, and directional designability. Therefore, aircraft designers can obtain a stronger and lighter structure by changing the layering design of composite materials according to the actual situation. Becky Aircraft firm made the world's first all-composite aircraft. It uses a heat-resistant carbon fiber layer sandwiched with an epoxide. A protective layer of graphite and epoxide surrounds a honeycomb material. The composite structure is half as light as the alloys of aluminum, steel, and titanium commonly used today and has nearly the same strength and heat resistance. At the same time, the aeroelastic deformation problem can be greatly improved by the design of the aircraft. According to statistics, the composite structure of large aircraft in the civil aviation market accounts for about 50% of the body structure [1]. In military aircraft, the amount of composite materials in the world's advanced military aircraft accounts for about 20% to 50% of the total structural weight of the aircraft. The main components of composite materials include fairing, flat tail, vertical tail, wings, and middle and front fuselage. In 1967, the U.S. Navy's F-14A heavy carrier-borne fighter aircraft made up less than 1% of its composite structure, while the U.S. Air Force's B-2 strategic bomber now makes up 50% of its composite structure. In the rapidly growing UAV industry, the U.S. Air Force MQ-1 Predator uses composite materials for 92% of the total weight of the structure, and the Israeli-developed Shadow multipurpose UAV uses composite materials for 95% of the structure [2]. To sum up, composite structure design has become an increasingly important part of modern aircraft design.
At present, a lot of basic research has been done by designers engaged in aerospace engineering and laminate design field in composite structure design. Zhou et al., Wan et al., and Liang et al. from Beihang University [3][4][5] took layering angle, layering ratio, layering unbalanced coefficient, and layering order as research targets to describe the aeroelastic effects of the above structural variables on high-aspect-ratio wings. Belbachir et al. [6] used the four-variable optimal plate theory and proposed an improved theory of laminated plates to study the bending problem of antisymmetric cross-laminated plates under nonlinear thermal and force loads. They reduced the original variables and derived the analytical solution of simply supported plates and verified it. Li et al. [7] discussed the relationship between the distribution law of static aeroelastic deformation and the unbalanced coefficient of high-aspectratio wing with layering unbalanced coefficient as a single variable. Xu et al. and Mostefa and Yousef [8,9] studied the deformed wing by numerical simulation and experiment, taking into account the influence of layering angle, layer thickness and position, and other variables according to the finite element analysis method. Zhu et al. and Xu et al. [10,11] focused on the aeroelastic tailoring and structural weight reduction of composite wings, which effectively controlled the deformation of wings and reduced the structural weight.
In terms of structural layering optimization, Lei et al. [12] used commercial ANSYS software to optimize the layering mode of a high-aspect-ratio wing. Basri et al. [13] studied the layering combination structure of different angles and further determined the optimal layering combination of the nodule wing leading edge. Bai et al. [14] proposed a three-stage optimization algorithm. With strength and deformation as constraints, the wing structure mass and flutter speed were significantly changed by the three-stage optimization algorithm. Liu and Xiang [15] considered nonlinear aerodynamic factors and took flutter velocity as the objective function to study the relationship between flutter velocity and layering angle of composite wings with different sections and improved the rear wing flutter velocity by 22 [16][17][18][19] performed a nonlinear dynamics study of variable stiffness composites and optimized them for aeroelastic problems.
Structural vibration is an important engineering issue that must be properly addressed when designing and developing aircraft or wind turbine equipment. In terms of modal studies of wings and fan blades, Liu et al. from Northwestern Polytechnical University and Huang and Yao from Nanjing University of Aeronautics and Astronautics [20,21] have explored the thermal modal effects of wings. In terms of research methods, the former mainly uses computational fluid dynamics and numerical heat transfer theory to carry out theoretical analysis and simulation, and the latter uses finite element method and experimental means to correct and observe the wing thermal modal effect. Zhong and Xu [22] described a modal method for calculating the mesh deformation of the flow around the wing based on the dynamic grid foundation and verified the method by using AGARD 445.6 wing. The numerical results obtained were in good agreement with the experimental results, and the calculation time was greatly reduced. Qiu et al. [23] used PATRAN/NASTRAN to analyze a certain composite wing model and obtained the natural frequency of the structure through the vibration mode experiment of the trailing edge of the wing, providing reliable data for the structural design of the wing. The team of An et al. and Zhang et al. from North China Electric Power University [24,25] took the wind turbine blade as the research object and analyzed the influence rule of layering parameters on the modal of the fan blade by simulation and experiment methods.
Among the many representative studies listed above, there are three main categories of research directions. The first is the structural design research represented by layering sequence, layering thickness, layering proportion, and many other parameters. The main purpose of this study is to determine the nfluence of the distribution of parameters on the aerodynamic or structural performance of wing. The second type is multidisciplinary comprehensive optimization of aircraft performance. Based on sensitivity algorithm, genetic algorithm, or hybrid algorithm, it can improve flutter speed or reduce structure weight by optimizing layering variables. The third is to study the dynamic response of the change of layering parameters to the structure. However, less attention has been paid to the unbalanced coefficients of laminates in the above studies. The static and dynamic problems caused by the change of the unbalanced coefficient of the laminate are still a subject to be studied.
In this work, a certain type of high-aspect-ratio wing is taken as the research object, and the influence of layingangle, laying-unbalanced coefficient, and laying-reference direction of composite material on the modal characteristics of high-aspect-ratio wing is investigated.

Computational Theory
2.1. Numerical Calculation Process. The calculation methods of aircraft static aeroelastic coupling include full coupling, strong coupling, and loose coupling. Among them, there is no time lag in the process of full coupling solution, and the solution accuracy is very high. However, this method has a large number of equations in the calculation, and the calculation is complicated and tedious, which is difficult to implement in practical applications. The strongly coupled calculation, on the other hand, is an interleaved coupling algorithm that introduces a prediction-multiple iterations within the same coupling time step. This algorithm requires repeated iterations in each physical time step, which brings a very high computational cost although the solution accuracy is relatively high. In contrast, the loosely coupled solution approach first solves the fluid control equations and then passes the solution of the fluid equations to the structural equations by interpolation after obtaining the aerodynamic loads and subsequently obtains the various physical parameters of the structural equations. Finally, the physical parameters of the structure are reinterpolated into the flow field and iteratively calculated. This way of solving makes the equations smaller in size and facilitates the saving of computational resources while satisfying the computational accuracy. Therefore, this paper solves the hydrostatic 2 International Journal of Aerospace Engineering aeroelastic problem using a loose coupling approach. The coupled calculation process is shown in Figure 1.
In order to accurately describe the vibration and inherent frequency of a structure under load, a nonlinear structural analysis must first be performed to obtain the correct stiffness matrix and then extract its eigenvalues to obtain the exact inherent frequency and vibration. For solving iterative methods, the Newton-Raphson method requires updating the stiffness matrix and inverse matrix solution for each iteration, which takes more time for each iteration, but converges relatively fast and requires fewer iterations. In contrast, the modified Newton method can iterate several times and then regenerate the stiffness matrix and inverse matrix, which takes less time for each iteration, but convergence is relatively slow and requires more iterations, and sometimes it is even difficult to converge. Therefore, the traditional Newton-Raphson method is chosen for this paper ( Figure 2).
In the numerical solution, the initial values of the structural nodal displacements and the cell stiffness matrix in the local coordinate system are first found using linear analysis of the grid cells, and then, the local equations are integrated into the structural stiffness matrix in the overall coordinate system. The whole process is iterated until the structure converges.

Aerodynamic Solution.
The Navier-Stoke equation in the Cartesian coordinate system can be expressed as [27] where Ω is an arbitrary fixed control body in the flow field where p is the mesh center density; Re is the Reynolds number; u, v, and w are the velocities in X, Y, and Z directions; ρ is the pressure; E is the total energy; τ is the stress tensor; ϕ is heat flux; and H is the enthalpy.

Laminate Structure Solution.
According to the classical laminate theory, for the whole laminate, the internal force and internal moment acting on the plate surface can be expressed as [28] As shown in Figure 3, the following formula can be derived:    Among them, where A ij is the tensile stiffness matrix, B ij is the coupling stiffness matrix, and D ij is the bending stiffness matrix. We use ω for out of plane; K x = ∂ 2 ω/∂x 2 , K y = ∂ 2 ω/∂y 2 , and K xy = ∂ 2 ω/∂x∂y; the displacement is called the bending curvature of the laminated plate.

Geometric Nonlinearity and Modal Solution.
The effect of structural geometric nonlinearity on the aeroelastic characteristics of a large aspect ratio wing is primarily reflected in the change of structural stiffness characteristics with load conditions and the effect of wing nonlinear deformation on the aerodynamic distribution of the wing. The structural strain in the geometric nonlinear large deformation of the wing is still small at this point, and the stress-strain of the structural material still satisfies the linear constitutive relationship, but the structure's strain-displacement is nonlinear. The tangent stiffness matrix must redefine the structural stiffness characteristics as follows [29]: where K inc is the principal tangential stiffness matrix, K u is the large displacement stiffness matrix, and K a is the initial load matrix. In structural calculation, its balance equation can be expressed as [24] M€ u where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, u is the displacement vector, _ u quarter is the velocity vector, € u is an acceleration vector, and F is the force vector.
The analysis of structural dynamic characteristics considering geometric nonlinearity shows that when the wing  Figure 4: HIRENASD wing model parameters.  5 International Journal of Aerospace Engineering vibrates slightly near the equilibrium position of large deformation, the vibration can be considered to be linear. The vibration equation is In the formula, x is the displacement array of the wing structure deviating from the equilibrium position of large deformation, M is the mass matrix of wing structure, and K is the tangent stiffness matrix of the wing structure at the equilibrium position of large deformation.     International Journal of Aerospace Engineering The wing microvibration is assumed to be simple harmonic vibration.
Substituting into equation (7) gives that where ω 2 is the eigenvalue and X is the eigenvector, corresponding to the quasilinear modal frequency and modal shape of the wing at the large deformation equilibrium position.

Method Verification
The HIRENASD wing (shown in Figure 4) was used as an example to verify the above calculation method. The longitudinal deformation of the leading edge of the wing at point A for different angles of attack was calculated and compared with the wind tunnel data. The specific incoming flow conditions are shown in Table 1.
In order to capture the air additive effects on the airfoil more accurately, the slit meshes at the leading and trailing edges of the airfoil and the airfoil junction are processed by the proximity method, the rest of the airfoil and airframe are processed by the curvature method, and the external flow field body mesh is divided using unstructured tetrahedral mesh cells. Finally, a total of five sets of meshing schemes are generated to verify the mesh independence of the flow field computational domain. The calculated and experimental values of the lift coefficient C L and drag coefficient C D of the wing are compared during the validation process. The details are shown in Table 2.
By comparing with the lift and drag coefficients of wind tunnel experiments, it can be seen that the lift and drag coefficients of the wing-body assembly will no longer change significantly when the minimum number of grid cells reaches 1.63 million or below. Therefore, in the subsequent numerical simulation, a discrete grid distribution scheme with 290,000 nodal cells and 1.63 million grid cells is chosen in this paper.
For the flow field calculation, the turbulence model uses the Spalart-Allmaras equation, and in the spatial discretization term, the viscous flux vector is selected in the central format and the convective flux vector is selected in the Roe-FDS format. The LU-SGS implicit time discretization method is used to advance the solution in the time discretization term [30,31]. The root plane of the wing-body assembly is set as the symmetric boundary condition, both the wing and the fuselage are set as the object boundary condition, and the rest are set as the pressure far-field boundary condition. In the structural calculation, the wing structure is made of 18-nickel martensitic aging steel, and the finite element modeling of the wing structure uses solid cells.
The tetrahedral mesh method is used for meshing. The number of structural nodes is 90,000, and the number of structural meshes is 51,000. The wing root is fixedly constrained, and the upper and lower surfaces of the wing are defined as fluid-solid coupling surfaces.
The final meshing and calculation results are shown in Figure 5. The results show that the numerical simulation results are in good agreement with the experimental results under the small angle of attack variation [32,33]. In summary, the method has high simulation accuracy within a certain angle of attack calculation range, so the method can be used to investigate the regularity of the static aeroelasticity of the wing with large span ratio.
Then, modal analysis of the HIRENASD wing model is carried out. Table 3 shows the comparison of the first six modes of the HIRENASD validated airfoil with the experimental values. The deviation is less than 5%, and the deviation between simulation value and experimental value is very small. Therefore, this method can be used for modal analysis of high-aspect-ratio wing structures.

Conceptual Design
The pneumatic-structural model takes a certain type of high span ratio UAV wing as the research object, as shown in Figure 6. On the premise of not changing the overall stiffness and affecting the calculation accuracy, other structures of the wing are simplified accordingly. Only the main load-bearing structures such as skin, wing beam and wing rib are retained, so as to achieve the purpose of improving the calculation speed. The specific geometric parameters of the wing model are as follows: the half-spread length is 16 m, the spreading chord ratio λ is 20.65, the root tip ratio η is 3.4, the leading edge swept back angle χ Θ is 6°, and the airfoil type is NACA63212 laminar flow airfoil type. The internal structure of the wing is a normal double wing beam with multiple   Table 4.
The initial reference direction of the ply is taken as the 0°r eference direction in the tangential direction of the wing root plane. The overall wing is paved with variable thickness, and the thickness decreases from the wing root to the wing tip. In the ply area, the starting point is the wing root, and 40% of the half span along the wing span is defined as ply area I; then, 20% of the wing span is defined as ply area II, and the rest of the wing to the wing tip is defined as ply area III. The reference direction is taken as the direction normal to the wing root plane. The specific division of the layup area is shown in Figure 7. At the same time, the thickness of single-layer board is defined as 0.2 mm.
In order to explore the influence of layering angle on wing modal under preload, layering angle selection is carried out in the first stage. Define the following: layer area I layer mode for½α°2/0°2/α°/0°/α°/0°2/α°2 s , a total of 22 layers, and thickness of 4.4 mm; layer area II layer mode for½α°/0°2/α°/0°/α°/0°2/α° s , a total of 18 layers, and thickness    In each layup, the same layup angle is adopted for α°in the three layup areas. Therefore, six layup schemes are produced by positive layup angle. Subsequently, in order to observe and compare whether there is a difference between positive and negative angles, reverse layering is carried out based on the above angles. That is, the layering angle is -15°, -30°, -45°, -60°, -75°, and -90°.
In the first stage, the appropriate laying-angle α°is determined, and in the second stage, the influence of the change of the unbalanced coefficient on the structural modal of the wing is explored. The unbalanced coefficient ranges from 0.2 to 0.8. Taking the unbalanced coefficient of 0.2 as an example, the sequence of single-angle layering is shown in Table 5.
Subsequently, another layup angle β°and a layup angle α°are added to form a hybrid layup scheme, where the layup angle β°is the layup angle that has a secondary effect on the mode in the previous phase of the angle screening exercise. The proportion of β°and 0°in the control layering scheme is 50% of the total layering thickness, respectively. The modal changes of wing structure after mixed layer were observed. The two layering angles of α°and β°were inverted to compare the modal changes caused by the change of unbalanced coefficient in mixed layering. Taking the unbalanced coefficient of 0.2 and layering angle α°as an example, the order of mixed layering is shown in Table 6.
The final is the third stage of the layering scheme design. From the second stage, the layup scheme with the minimum modal frequency is determined, and based on the layup sequence, the influence of layup reference direction on the modal of the structure is studied. The reference direction is defined as a gradual change, with a change of 15°. It gradually changes from 0°until the off-axis angle of the reference  9 International Journal of Aerospace Engineering direction of layering is 180°. In the stage design of the above three layering schemes, the layering angle, layering unbalanced coefficient, and layering reference direction are considered comprehensively. Based on this, the structural modal analysis of high-aspect-ratio wing under preload is carried out.
Taking the chord direction of the wing as the X-axis and the spreading direction as the Y-axis, the angle of change of the layup reference direction in the XY plane can be expressed as A second reference direction needs to be introduced in the spatial coordinate system to describe the spatial coordinates, and for convenience, the coordinate system of both directions is used in the form of a matrix at the same time. At this point, the angles in the first, second, third, and fourth quadrants can all be expressed as Based on the above theory, all the required variations in the layup reference direction angles can be deduced. The specific variation principle is shown in Figure 8. The final direction of each reference direction change vector is shown in Table 7.
After verification of grid and iteration step independence, the obtained finite element model of high-aspectratio wing is shown in Figure 9.

Effect of Layering Angle on Modal Frequency.
Vibration mode refers to the inherent vibration form of elastic body or elastic system. For the modal study of high-aspect-ratio wing structures, since higher-order modes require more energy injection, and lower-order vibrations are more likely  Figure 11: Different angle laminate structure diagram. 10 International Journal of Aerospace Engineering   11 International Journal of Aerospace Engineering to cause structural problems, lower-order modes are generally more concerned. Through the above calculation, the distribution of the first six modes of the high-aspect-ratio composite wing structure is obtained, as shown in Figure 10. It is mainly manifested as the first-order vertical bending mode, the second-order vertical bending mode, the first-order horizontal torsion mode, the third-order vertical bending mode, the fourth-order vertical bending mode, and the second-order horizontal torsion mode. This means that the main characteristics of the vibration characteristics of the wing of the structural layout type are in the form of bending and torsion coupling.
When the layup angle of the wing is changed, the fundamental vibration pattern remains the same. Table 8 shows the modal frequencies of the wing structure at each layup angle. It can be found that the modal frequencies of the first four orders of the airfoil structure show an increasing trend with the increase of the layup angle, and the larger the layup angle, the larger the incremental change of the modal frequencies. The increments of the first-order modal frequencies are 0.011 and 0.066 from 15°to 45°, but 0.455 and 0.490 when the layup angle is changed from 60°to 90°. The rest of the orders are similar. On the other hand, the frequency results for the positive and negative layup angles are essentially the same and do not appear to be very different. This implies that in the composite layup design of the wing structure, when the layup ratio is constant, the positive and negative angle layups show irrelevance, both of which can make the bendingtorsion coupling type wing stiffness consistent and can obtain high modal frequencies.
The reason for this increase in frequency can be broadly attributed to the fact that the stiffness of the structure in a certain direction changes when the angle of the ply is changed. The modulus of elasticity of the laminate structure decreases with increasing angle, while the shear modulus increases with increasing angle. In addition, as the angle of the ply increases, the directional range of stiffness changes and is no longer a single range. This results in a gradual increase in the bending and torsional resistance of the wing structure, while the resistance to bending deformation decreases to a certain extent. On the other hand, the modulus of elasticity in the main direction of the positive and negative laminates is the same during the layup process, but the direction is reversed. Therefore, in the wing ply design, the structural frequency does not change after the positive and negative angle plies are laid. The structural parameters of the laminate for each layup angle are shown in Figure 11.

Effect of Layering Unbalanced Coefficient on Modal
Frequency. The effect of different unbalanced coefficients on the modal frequency of the wing structure is shown in Figure 12. When there are only two angular layups, the results are shown in Figures 12(a)-12(f). It is found that the unbalanced coefficient does not affect the modal frequencies of the wing structure significantly. Although the value of the frequency increases slightly with the increase of the unbalanced coefficient, the amount of this change is extremely subtle. The layup scheme consisting of 0°and 90°layup angles results in a greater modal frequency of the wing structure compared to the combination of only 0°and 45°layups, which means that the combination of 0°and 90°layup schemes has a greater structural stiffness compared to the other layup schemes. The reason for this phenomenon can be attributed to the fact that when the layup angles are only 0°and 45°, the axial stiffness carrying capacity of the layup structure is smaller, while the in-plane shear stiffness carrying capacity is relatively larger. Such a pavement structure mainly bears the torsional influence caused by external load, which is called shear skin. With 90°pavement angle, the 90°axial stiffness of the skin structure is strengthened and the tensile load carrying capacity is stronger, but the shear strength is reduced. This kind of skin is also called hard skin. Another noteworthy phenomenon is that the fifthand sixth-order modal frequencies show an anomalous

12
International Journal of Aerospace Engineering trend, no longer is the frequency of 0°and 45°layup combinations less than that of 0°and 90°layup combinations. Rather, the opposite is true: the frequency of the 0°and 45°ply combinations is greater than that of the 0°and 90°ply combinations. This may be influenced by the higher order deformation, which makes the torsional properties of the wing enhanced.
On the other hand, the comparison of the chordwise deformation of the three regions of the wing is shown in Figure 13. The change of the unbalanced coefficient has a large effect on the chordal deformation of the wing. When the unbalanced coefficient is 0.2, the chordal deformation is linearly distributed, while when the unbalanced coefficient increases gradually, the chordal deformation is nonlinearly distributed, which is different from the usual conclusion that the chordal deformation of the wing is rigid.
The frequencies for the two-angle pavement and threeangle pavement cases were then compared, as shown in Figure 14. As can be seen from the figure (observe the column), the first six orders of modal frequencies for the threeangle hybrid paving scheme are better than those for the      two-angle hybrid paving scheme, especially the first-order frequencies and the second-order frequencies are particularly enhanced; i.e., the addition of the 90°paving angle enhances the vibration characteristics of the lower order. Secondly, observing the same line diagram, it can be seen that the modal frequency of the mixed layup scheme consisting of three layup angles is greater than that of the scheme consisting of two layup angles. Compared with the fourth-order frequency of the laminate with an unbalanced factor of 0. 8 It is also found that the modal frequencies of the pavement design scheme consisting of 0°and 45°as the dominant and 90°as the supplementary are smaller than those of the pavement design scheme consisting of 0°and 90°as the dominant and 45°as the supplementary. In terms of firstand second-order modal frequencies, the hybrid ply scheme with [0°/90°/45°] angles is more likely to have higher frequencies at unbalanced coefficients of 0.3 to unbalanced coefficients of 0.5. This is because the laminate evolves gradually from a low-equilibrium laminate to a quasiequilibrium laminate as the unbalanced factor increases. This evolution leads to a change in the mechanical and coupling properties of the laminate, thus changing the vibration characteristics at the macroscopic level.
Taking the unbalanced factor of 0.5 as an example, Figure 15 gives a microscopic comparison of the four layup schemes in terms of laminate structure. It can be seen that when the [0°/45°] ply scheme is used, the axial stiffness is smaller than that of the [0°/90°] ply scheme, but the relative stiffness bearing range is wider; when the [0°/45°/90°] ply scheme is used, the shear bearing capacity is much larger than that of the [0°/90°/45°] ply scheme, but at the same time, the axial stiffness also decreases sharply. In the [0°/45°/90°] layup scheme, the laminates have mixed layup angles and gradually evolve into a quasi-isotropic plate.

Effect of Reference Direction of Layering on Modal
Frequency. After completing the exploration of the unbalanced coefficient in the second stage, the third stage of the scheme focuses on the analysis of the effect of the layup reference direction on the modal frequencies of the wing structure with an unbalanced coefficient of 0.5 (i.e., symmetrical layup) as an example. Figure 16 gives the frequencies of the first six orders of the modalities of the wing structure under the change of the layup reference direction when the unbalanced coefficient is 0.5. According to the information in the figure, with the change of the layup reference direction, the lowest frequency generally exists in the 0°reference direction except for the sixth-order modal frequency, and one to two peaks and valleys appear in each order modal frequency, but there is no obvious regularity in the location of the peaks and valleys to be found. The specific frequency peaks and valleys should vary with the changing force characteristics of the wing structure. For the first-order mode,

16
International Journal of Aerospace Engineering the maximum frequency appears in the 165°reference direction; for the second-order mode, the maximum frequency appears in the 120°reference direction; for the third-order mode, the maximum frequency appears in the 165°reference direction; for the fourth-order mode, the maximum frequency appears near the 120°reference direction; for the fifth-order mode, the maximum frequency appears near the 120°reference direction; for the sixth-order mode, the maximum frequency appears near the 75°reference direction. The above results imply that the second quadrant-based layup reference direction angle may be helpful in the selection of the layup reference direction for the enhancement of the modal frequency and structural stiffness. The above results also illustrate that the different angles of reference direction selection can, to some extent, make the structural frequencies change significantly.

Conclusions
Based on CFD/CSD coupling method, by designing the composite structure of the high-aspect-ratio wing, the influence of layering angle, layering unbalanced coefficient, and layering reference direction on the modal frequency of wing structure is discussed. The following conclusions are drawn:  (1) The first sixth-order modal frequency of the wing structure increases with the increase of the ply angle, and the larger the ply angle, the larger the incremental change of its frequency; the frequency results of the positive ply angle and the negative ply angle are basically the same and do not appear to be very different. This is because the modulus of elasticity in the main direction is the same for positive and negative laminates during the paving process, but the direction is reversed. Therefore, in the design of the wing ply, it will appear that the structural frequency does not change basically after the positive and negative angle plies (2) The modal frequency of the wing structure is not very sensitive to the change of the unbalanced coefficient. The modal frequencies obtained in the mixed ply scheme with three angles are significantly greater than those in the ply scheme with only two angles. In this example, the laminate evolves gradually from a low-equilibrium laminate to a quasiequilibrium laminate as the unbalanced coefficient increases, so it will lead to a higher frequency at the firstand second-order modal frequencies for the hybrid paving scheme with [0°/90°/45°] angles at the unbalanced coefficient of 0.3 to the unbalanced coefficient of 0.5 (3) The lowest frequencies are generally present in the 0°reference direction as the layup reference direction changes. Several extremes appear in each order of modal frequencies, but there is no obvious regularity to find where the extremes appear; the angle of the reference direction dominated by the second quadrant is beneficial to enhance the modal frequencies. In summary, after determining the layup angle, layup order, and layup ratio, proper adjustment of the layup reference direction according to the actual design situation will be beneficial to change the vibration characteristics of the airfoil structure

Data Availability
Some of the data and models used to support the results of this study can be obtained from the authors upon reasonable request. All simulation data were used only during the study and do not involve any commercial profit situation, which includes the Excel worksheet used for documentation. The experimental data in this paper include calculated stress strain, deformation and failure data, design solutions for the layups, and the wing model parameters used. Information on the computational grid and initial computational conditions is also available through the authors.

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