# Effects of Unbalanced Lamination Parameters on the Static Aeroelasticity of a High Aspect Ratio Wing

Authors: JinYang Li; JunLi Wang; ZhiGui Ren; WeiFeng Wei
Journal: International Journal of Aerospace Engineering
 (2021)

---

## Abstract

In this paper, in order to understand the influence of the unbalanced coefficient of composite laminates on the static aeroelasticity of high aspect wings, a series of numerical simulation calculations were carried out, and this work wants to provide some reference for the structural design of aircraft. Considering the influence of geometric nonlinearity, the unidirectional fluid-solid coupling calculation method based on loose coupling is used to control the change of unbalanced coefficient of laminates on the basis of layering angle, layering thickness, and layering region, so as to observe the changes caused to the wings. The relationship between the unbalanced coefficient and the constant thickness layup and the variable thickness layup with 0° and ±45° layup angles was studied, respectively. Then, the layup angle of 90° was added to study the influence of the unbalanced coefficient on the static aeroelasticity of the wing structure with the change of the layup angle and the different choice of layup region. The results show that the deformation is the smallest when the unbalanced coefficient is 0.5, and the deformation trend is evenly distributed along both sides when the unbalanced coefficient is 0.5. When the unbalanced coefficient is changed, adding the 90° layup angle can significantly reduce the overall deformation of the wing and show different sensitivity characteristics to different layup areas. The increase of the unbalanced coefficient makes the chordal displacement gradually change from linear distribution to nonlinear distribution along the spread direction, and the displacement will gradually decrease.


---

## Body

## 1. Introduction

Advanced composite material (ACM) was born in the late 1960s. In the early 1970s, the U.S. Air Force Materials Research Institute applied composite materials to the structure of fighter aircraft in order to improve the performance of military fighter aircraft. Since 1980, the US Navy has applied a large number of advanced composite materials to the latest fighter aircraft, starting with the F/A-18 carrier-based multirole fighter aircraft. For example, the composite material on the USAF’s F-22 “Raptor” fighter aircraft occupies 25% of the weight of the structure, the composite material on the B-2 “Ghost” bomber aircraft occupies 50% of the weight of the structure, and the composite material consumption of the “Predator” UAV even reaches 92% of the total weight of the structure [1, 2]. Similarly, in the civil aviation market, B787 and A380 of the world’s top two airlines, led by Boeing and Airbus, use composite materials for 50% and 52%, respectively [3, 4]. Composite materials have gradually replaced metal and alloy to some extent and become the “new favorite” in aviation design and manufacturing. It can be said that the application of composite material science and technology has promoted the rapid development of modern aircraft structural design. Compared with the traditional metal materials, the composite material has the characteristics of high specific strength, high specific stiffness, strong fatigue corrosion resistance, and low thermal expansion coefficient, so it is easier to design the body structure with light weight, large bearing force, and small deformation. In modern aircraft design with high aspect ratio, aeroelastic design of composite wing is of great significance due to the characteristics of large wing flexibility and small stiffness [5]. In order to optimize the flight performance and aeroelastic characteristics of aircraft with high aspect ratio, aeroelastic tailoring technology is often needed, and the comprehensive design of layering number, angle, and sequence is usually carried out.With the increasing maturity of composite material technology and aerodynamic tailoring technology, coupled with its extensive application in aircraft design engineering practice, more and more researchers are devoted to the research of composite material aircraft layup and structure. In the aspect of aeroelastic tailoring design, the main work at present is to study the structural parameters such as laminate order, laminate thickness, and laminate proportion of composite materials and explore the influence of parameter allocation on aeroelastic characteristics of the wings. Wan et al. [6] studied the aeroelastic effects of unevenly laid wings with high aspect ratio by unevenly laid wings. Liang et al. [7] carried out aeroelastic optimization design for composite wing panels of large aircraft and obtained reasonable distribution of lamination proportion. Zhou et al. [8] studied the influence of composite lamination parameters on aeroelasticity optimization of the wing with high aspect ratio, and the results showed that the skin-unbalanced lamination mainly affected the aeroelasticity of the wing. Xu et al. [9] changed the symmetry and equilibrium of the laminated plate by adding the laminated layer on top of the symmetrical laminated layer and obtained the ideal position of the added layer and laying angle. Xingyin et al. and Zhang et al. [10, 11] focused on the coupling analysis of bending and torsion deformation of asymmetric/unbalanced laminates by using experimental methods. Wang et al. [12] explored and optimized the layup angle so as to reduce the deformation of the wing with high aspect ratio and effectively improve the aeroelastic characteristics of the wing. Basri et al. [13] studied layup combination structures with different angles and further determined the optimal layup combination for the leading edge of the wing with nodules. Koohi et al. and Xie et al. [14, 15] carried out aeroelastic analysis on the composite wing beam model with high aspect ratio on the basis of considering geometric large deformation and achieved a good calculation accuracy. Kirsch et al. [16] developed the GEBTAero code to evaluate the aeroelastic tailoring effect and carry out numerical calculation. Bourchak and Dobah [17] designed a modified wing with good performance by simulation and experiment on the deformed wing from the perspective of layup, nonequilibrium, layup number, and layup position. Dillinger et al. [18] proposed a modified model for static aeroelastic aerodynamic load of composite wing and successfully applied it to stiffness optimization of forward swept wing. Bramsiepe et al. and Yu et al. [19, 20] carried out numerical simulation of maneuver load and gust load of composite flying wing UAV and carried out parameter optimization design. Shrivastava et al. [21] adopted the combination of classical genetic algorithm and CAE to carry out weight optimization design for a transport aircraft wing and achieved the expected calculation results. DLR-AE [22] at the Aero-Elastic Research Institute of the German Aerospace Center calculated the aeroelastic response of the composite wing under multiple working conditions using internal development tools and compared it with that of the aluminum alloy structure wing.Some other researches mainly focus on the mathematical model construction of laminates; the purpose of which is to describe the mechanical properties of laminates more accurately and provide some theoretical basis for engineering design. Shukla and Singh [23] used the radial basis function to analyze the flexure of the laminates and derived the governing differential equation using the energy principle. Then, the accuracy of the method was verified by using MATLAB software to write a program. Kwon [24] studied the failure of composite structures under hydrodynamic loads by means of fluid-solid coupling numerical simulation and then discussed the multiscale modeling method of composite laminates. Guillén-Rujano et al. [25] derived a new calculation formula for composite laminates by using the first-order shear deformation theory and the minimum potential energy theorem. This formula can be used to calculate the bending and torsion problems of composite laminates well, and the calculation accuracy is in good agreement with the experimental data. Thakur et al. [26] proposed a finite element model of C0 laminated plates, which used von Karman and green-Lagrangian strain-displacement relation to simulate geometric nonlinearity and then derived the nonlinear motion governing equation by using Hamiltonian principle. Finally, the Newton-Raphson method was used to solve the problem. The calculated results are reliable and efficient, which can solve the free vibration and forced vibration of composite laminates. Yang and Yang [27] analyzed the dynamics of laminated plates under nonlinear elastic constraints, established different binding models, and obtained the transverse displacement of laminated plates through transient analysis. Using the analysis principle of micromechanics, Nallim et al. [28] proposed a variational formula for the analysis of the free vibration of the elastic constraint asymmetric composite laminates and gave the influence of different mechanical and geometric parameters on the dynamic properties of different laminates. Hafizah et al. [29] studied the free vibration problem of laminated plates with antisymmetric angles of varying thickness. By introducing the theory of high-order shear deformation plates, the shear correction factor was successfully eliminated and the accuracy of transverse shear stress was improved. Liu et al.’s [30] laminated plate is studied in the nonlinear response of a subsonic condition, using the Hamilton’s principle, the composite laminated plate is established in subsonic aerodynamic force motion control of nonlinear partial differential equation, and the frequency response curve under different parameters is obtained; thus, the exciting force and laminated plate structure provide reference for the coupling problem.In summary, the above studies are mainly divided into two aspects. On the one hand, the performance design of structural parameters such as laying angle, laying order, and asymmetry/unbalanced of high aspect ratio composite wing was carried out, or the sensitivity algorithm and genetic algorithm were used to optimize the aeroelasticity of composite wing. On the other hand, a more accurate and efficient calculation method is obtained by studying the mathematical model of laminates and deriving their constitutive equations. However, in the current design of composite materials for aircraft, few studies have paid attention to the changes in the mechanical properties of laminates caused by the changes in unbalanced coefficients, and few studies have involved the effects of the differences in unbalanced coefficients on the static aeroelasticity of aircraft with high aspect ratio. The theoretical system of composite wing structure design still has some shortcomings. Therefore, based on the calculation method of fluid-structure interaction, this paper focuses on the analysis of the static aeroelasticity of high aspect ratio wing under different unbalanced coefficients (layer disequalization coefficient, i.e., the ratio of the number of layers on the upper/lower 45° wing surface to the number of layers on the upper/lower 45° wing surface and the number of layers on the upper/lower 45° wing surface and -45° wing surface), the layer thickness, angle, and the influence of each layer segment on the static aeroelasticity.In the first part, the theoretical solution method of static aeroelasticity is introduced. In the second part, the establishment of the finite element model of the wing with high aspect ratio and the design of the layup thickness, angle, and wing segment are described in detail. In the third part, the effects of various variables on the aerostatic elasticity of the wings with high aspect ratio are analyzed under different unbalanced coefficients.


## 2. Computational Theory

### 2.1. Calculation Method

Compared with the vortex lattice method in the common static aeroelastic analysis, the unidirectional fluid-solid coupling method based on loose coupling has the advantages of fast calculation speed and relatively high accuracy. In the numerical simulation, the Navier-Stokes equation is used as the governing equation to solve the aerodynamic load by the finite volume method. If the aerodynamic load converges, the interpolation method is used to insert it into the structural grid nodes. Finally, the finite element method and Newton-Raphson iterative method are used to solve the structural equation under aerodynamic force, and the deformation and stress-strain values of the wing structure with high aspect ratio are obtained. The specific implementation process is shown in Figure1.Figure 1
Calculation process.


### 2.2. Aerodynamic Calculation

Using the Euler description method, the differential form of three-dimensional conservation N-S equation in Cartesian coordinate system is as follows [30–33]:
(1)∂W∂t+∂f∂x+∂g∂y+∂q∂z=∂R∂x+∂S∂y+∂T∂z,where W is the conserved variable vector, i.e.,
(2)W=ρρuρvρwρE.f, g, and q are the convective flux vectors, i.e.,
(3)f=ρuρu2+pρuvρuwρuE+up,g=ρvρvuρv2+pρvwρvE+vp,q=ρwρwuρwvρw2+pρwE+wp.R, S, and T are viscous flux vectors, i.e.,
(4)R=0τxxτyxτzxφx,S=0τxyτyyτzyφy,T=0τxzτyzτzzφz.Among them,(5)φx=uτxx+vτxy+wτxz+k∂T∂x,φy=uτyx+vτyy+wτyz+k∂T∂yφz=uτzx+vτzy+wτzz+k∂T∂z.,For Newtonian fluid, the formula for calculating the viscous stressτij is as follows:
(6)τxx=2μux−23μux+vy+wz,τxy=τyx=μuy+vx,τyy=2μvy−23μux+vy+wz,τxz=τzx=μuz+wx,τzz=2μwz−23μux+vy+wz,τyz=τzy=μvz+wy.In order to make the N-S equation closed, the state equation is introduced as follows:(7)p=γ−1ρE−12u2+v2+w2,where p is the fluid pressure, ρ is the fluid density, E is the total fluid energy per unit mass, u, v, and w are the three velocity components of fluid velocity V in the Cartesian coordinate system, τij is the component of the viscous stress term τ¯¯, γ is the specific heat ratio, T is the fluid temperature, k is the fluid conductivity coefficient, and μ is the viscosity coefficient of the fluid.The viscosity coefficient of air is generally calculated according to Sutherland formula, and its expression is as follows:(8)μμ0=T0+CT+CTT01.5.In the equation,T0 and μ0 are the temperature and viscosity coefficient of the standard atmosphere at sea level, respectively. T0=288.15K and μ0=1.7894×10−5N·s/m2. The constant C=110.4k.The relationship between the heat conduction coefficient and viscosity coefficient of a complete gas is as follows:(9)k=cpμPr,where cp is the isobaric specific heat capacity, Pr is Prandtl number, and Pr of laminar air is 0.72.


### 2.3. Structural Calculation

In structural calculation, its equilibrium equation can be expressed as follows [33–37]:
(10)Mu¨+Cu̇+Ku=F.When the geometrically nonlinear problem of large deformation is involved, the stiffness matrixK can be expressed as follows [38]:
(11)K=Kinc+Ku−Ka,where M is the mass matrix,C is the damping matrix,K is the stiffness matrix,uis the displacement vector,u is the velocity vector,u¨is the acceleration vector, Ft is the force vector,Kinc is the main tangential stiffness matrix,Ku is the large displacement stiffness matrix, and Ka is the initial load matrix.


### 2.4. Coupled Data Transfer

In order to improve the computational accuracy of geometrical nonlinearity of structures, the interpolation accuracy of fluid-structure interface will play an important role in the calculation. The conserved interpolation method is used to satisfy the data exchange of the coupling surface. In the case of upper-solid coupling of coupling interface, it can be expressed as follows [39]:
(12)uf=Hus,δuf·fa=δusT·fs,whereuf is the pneumatic point displacement, us is the structural point displacement, H is the interpolation matrix, fa is the aerodynamic load, fs is the equivalent load of structure, and δuf and δus are the virtual displacements of aerodynamic points and structural points, respectively.IfH is independent of displacement, then
(13)δuf=H·δus.Then, the relationship between equivalent load of structure and aerodynamic load is obtained, i.e.,(14)fs=HT·fa.The coupling diagram is shown in Figure2.Figure 2
Coupled transfer diagram.


### 2.5. Laminate Calculation

The classical laminate theory has the following assumptions: (1) the bond layer between the laminates is very thin and firmly bonded, and no slip occurs between the laminates. (2) The laminates are thin plates, the thickness of the laminates is constant, and the single-layer laminates are analyzed according to the plane stress state. (3) The bending deformation of the laminate is in the range of small deflection. Before deformation, it is perpendicular to the straight line of the middle plane. After deformation, it remains straight and perpendicular to the middle plane, and the length of the straight line remains unchanged.Suppose that the displacement of any pointC in the distance plane z along the three coordinate axes x, y, and z direction is ux,y,z, vx,y,z, and wx,y,z, respectively, then the classical laminate theory can be obtained [39]:
(15)ux,y,z=u0x,y−z∂wx,y∂x,vx,y,z=v0x,y−z∂wx,y∂y,wx,y=w0x,y,where u0x,y, v0x,y, and w0x,y are the displacements of the middle plane of the laminates.Then, the strain-displacement geometric relationship is as follows:(16)εx=∂u∂x,γyz=∂w∂y+∂z∂x,εy=∂v∂y,γzx=∂w∂x+∂u∂z,εz=∂w∂z,γxy=∂u∂y+∂v∂x.It is calculated as follows:(17)εx=∂u0∂x−z∂2w∂x2=εx0+zkx,εy=∂v0∂y−z∂2w∂y2=εy0+zky,γxy=∂u0∂x+∂v0∂y−2z∂2w∂x∂y=γxy0+zkxy,where εx0, εy0, and γxy0 are the strains in the middle of the laminates, kx and ky are the bending rates of the middle surface of the laminate, and kxy is the distortion rate in the middle of the laminate.The above equation can be expressed as a matrix:(18)εxεyγxy=εx0εy0γxy0+zkxkykxy.In the Cartesian coordinate system, the off-axis stress-strain relationship of the single-layer plate is as follows:(19)σxσyτxy=Q¯11Q¯12Q¯16Q¯12Q¯22Q¯26Q¯16Q¯26Q¯66kεxεyγxy.Among them,(20)Q11=E11−V12V21,Q22=E21−V12V21,Q12=V21E21−V12V21=V12E11−V12V21,Q¯66=G12,where σx and σy are the off-axis normal stresses, τxy is the partial shear stress, εx and εy are the off-axis normal strains, γxy is the off-axis shear strain, E1 is the longitudinal modulus of elasticity, E2 is the transverse modulus of elasticity, G12 is shear elastic modulus, V12 is the longitudinal Poisson’s ratio, V21 is the transverse Poisson’s ratio, and Q¯ is the transformation matrix of reduced stiffness matrix Q, which is expressed as follows:
(21)Q¯=T−1QT−1T,where T is the transpose matrix:
(22)T=cos2θsin2θ2cosθsinθsin2θcos2θ−2cosθsinθ−cosθsinθcosθsinθcos2θ−sin2θk.The stress-strain relationship of thek-layer single layer of the laminated plate in the coordinate system XYZ can be obtained by simultaneous equations (18) and (19):
(23)σxσyτx=Q¯kεx0εy0γxy0+zkxkykxy.By integrating the stress of each single-layer plate, the internal force of the laminated plate can be expressed as follows:(24)NxNyNxy=∑k=1n∫zk−1zkσxσyτxykdz=∑k=1n∫zk−1zkQ¯kεx0εy0γxy0+zkxkykxydz,MxMyMxy=∑k=1n∫zk−1zkσxσyτxykzdz=∑k=1n∫zk−1zkQ¯kεx0εy0γxy0+zkxkykxyzdz,where n is the number of single-layer plates and zk−1 and zk are the z coordinates of the bottom and top surfaces of the single-layer plate of layer k, respectively.Finally, the expression of internal force of composite laminates is arranged into a matrix form:(25)NxNyNxyMxMyMxy=A11A12A16B11B12B16A12A16B11B12B16A22A26B12B22B26A26A66B16B26B66B12B16D11D12D16B22B26D12D22D26B26B66D16D26D66εx0εy0γxy0kxkykxy.Among them,(26)Aij=∑k=1nQ¯ijkzk−zk−1,Bij=12∑k=1nQ¯ijkzk2−zk−12,Bij=13∑k=1nQ¯ijkzk3−zk−13,where A11, A12, and A22 are the stiffness coefficients between tensile (compressive) force and middle plane tensile (compressive) strain. A66 is the stiffness coefficient between the shear force and the shear strain in the middle plane. A12 and A26 are the coupling stiffness coefficients between shear and tension. B11, B12, and B22 are the coupling stiffness coefficients between tensile and bending. B66 is the coupling stiffness coefficient between shear and torsion. B12 and B26 are the coupling stiffness coefficients between tension and torsion or between shear and bending. D11, D12, and D22 are the stiffness coefficients between bending and bending rate. D66 is the coupling stiffness coefficient between torsion and torsion rate. D12 and D26 are the coupling stiffness coefficients between torsion and bending.


### 2.6. Validation of Calculation Methods

The HIRENASD wing model was selected as the research object (the wing model is shown in Figure3 [40]), and the loose coupling unidirectional fluid-solid coupling theory was adopted as the calculation method [41]. The airfoil is supercritical BAC3-11/R ES/30/21, the reference area of the airfoil is 0.3926 m2, the reference length is 0.3445 m, and the aspect ratio is 12.08.The longitudinal deformations of point A near the leading edge of the wingtip at angles of attack of -1.5°, 0°, 1.5°, 3°, and 4.5°, with a Reynolds number of 7×106 and a flight speed of 0.8 Mach were calculated, respectively.Figure 3
HIRENASD wing model and parameters.
(a)
HIRENASD wing parameters(b)
The position of point A on the modelThe unstructured tetrahedral mesh was used to divide the flow field, and the mesh at the leading edge of the wing was encrypted. A total of 2,581,247 nodes and 1,869,810 grid cells were generated. The plane of the fuselage root is set as the symmetric boundary condition, the wing and fuselage are set as the material boundary condition, and the others are set as the far-field pressure boundary condition. In the structural calculation, the main wing finite element model was divided by unstructured tetrahedral mesh, with a total of 133,215 nodes and 75,790 mesh elements generated. The boundary condition was the constraint of the wing root fixed support. The overall meshing details are shown in Figure4.Figure 4
HIRENASD wing finite element model.
(a)
Finite element model of flow field(b)
Finite element model of structureAfter calculation, the numerical simulation results agree with the experimental ones. When the angle of attack is 4.5°, there is a big difference, the maximum deformation of the test value is 25.1433 mm, the maximum deformation of the calculated value is 22.751 mm, and the deviation between the two is about 10%.The reason for the small deviation in the overall calculation may be that the deformation marks of the wing structure in the test and the deformation calculation points of the numerical calculation model cannot coincide in a strict sense, so there is a certain numerical deviation. In addition, with the gradual increase of the angle of attack of the wing, the flow separation caused by the wing surface becomes more and more obvious, and the computational grid is difficult to accurately solve, so the data also produce a certain deviation. In conclusion, this method has a high simulation accuracy within a certain angle of attack calculation range and is suitable for the static aeroelastic analysis of high aspect ratio wing. Specific experimental values can be referred to References [40, 42]. Figure 5 shows the comparison between the calculated results and the simulation results.Figure 5
Comparison of calculated results and experimental results.


## 2.1. Calculation Method

Compared with the vortex lattice method in the common static aeroelastic analysis, the unidirectional fluid-solid coupling method based on loose coupling has the advantages of fast calculation speed and relatively high accuracy. In the numerical simulation, the Navier-Stokes equation is used as the governing equation to solve the aerodynamic load by the finite volume method. If the aerodynamic load converges, the interpolation method is used to insert it into the structural grid nodes. Finally, the finite element method and Newton-Raphson iterative method are used to solve the structural equation under aerodynamic force, and the deformation and stress-strain values of the wing structure with high aspect ratio are obtained. The specific implementation process is shown in Figure1.Figure 1
Calculation process.


## 2.2. Aerodynamic Calculation

Using the Euler description method, the differential form of three-dimensional conservation N-S equation in Cartesian coordinate system is as follows [30–33]:
(1)∂W∂t+∂f∂x+∂g∂y+∂q∂z=∂R∂x+∂S∂y+∂T∂z,where W is the conserved variable vector, i.e.,
(2)W=ρρuρvρwρE.f, g, and q are the convective flux vectors, i.e.,
(3)f=ρuρu2+pρuvρuwρuE+up,g=ρvρvuρv2+pρvwρvE+vp,q=ρwρwuρwvρw2+pρwE+wp.R, S, and T are viscous flux vectors, i.e.,
(4)R=0τxxτyxτzxφx,S=0τxyτyyτzyφy,T=0τxzτyzτzzφz.Among them,(5)φx=uτxx+vτxy+wτxz+k∂T∂x,φy=uτyx+vτyy+wτyz+k∂T∂yφz=uτzx+vτzy+wτzz+k∂T∂z.,For Newtonian fluid, the formula for calculating the viscous stressτij is as follows:
(6)τxx=2μux−23μux+vy+wz,τxy=τyx=μuy+vx,τyy=2μvy−23μux+vy+wz,τxz=τzx=μuz+wx,τzz=2μwz−23μux+vy+wz,τyz=τzy=μvz+wy.In order to make the N-S equation closed, the state equation is introduced as follows:(7)p=γ−1ρE−12u2+v2+w2,where p is the fluid pressure, ρ is the fluid density, E is the total fluid energy per unit mass, u, v, and w are the three velocity components of fluid velocity V in the Cartesian coordinate system, τij is the component of the viscous stress term τ¯¯, γ is the specific heat ratio, T is the fluid temperature, k is the fluid conductivity coefficient, and μ is the viscosity coefficient of the fluid.The viscosity coefficient of air is generally calculated according to Sutherland formula, and its expression is as follows:(8)μμ0=T0+CT+CTT01.5.In the equation,T0 and μ0 are the temperature and viscosity coefficient of the standard atmosphere at sea level, respectively. T0=288.15K and μ0=1.7894×10−5N·s/m2. The constant C=110.4k.The relationship between the heat conduction coefficient and viscosity coefficient of a complete gas is as follows:(9)k=cpμPr,where cp is the isobaric specific heat capacity, Pr is Prandtl number, and Pr of laminar air is 0.72.


## 2.3. Structural Calculation

In structural calculation, its equilibrium equation can be expressed as follows [33–37]:
(10)Mu¨+Cu̇+Ku=F.When the geometrically nonlinear problem of large deformation is involved, the stiffness matrixK can be expressed as follows [38]:
(11)K=Kinc+Ku−Ka,where M is the mass matrix,C is the damping matrix,K is the stiffness matrix,uis the displacement vector,u is the velocity vector,u¨is the acceleration vector, Ft is the force vector,Kinc is the main tangential stiffness matrix,Ku is the large displacement stiffness matrix, and Ka is the initial load matrix.


## 2.4. Coupled Data Transfer

In order to improve the computational accuracy of geometrical nonlinearity of structures, the interpolation accuracy of fluid-structure interface will play an important role in the calculation. The conserved interpolation method is used to satisfy the data exchange of the coupling surface. In the case of upper-solid coupling of coupling interface, it can be expressed as follows [39]:
(12)uf=Hus,δuf·fa=δusT·fs,whereuf is the pneumatic point displacement, us is the structural point displacement, H is the interpolation matrix, fa is the aerodynamic load, fs is the equivalent load of structure, and δuf and δus are the virtual displacements of aerodynamic points and structural points, respectively.IfH is independent of displacement, then
(13)δuf=H·δus.Then, the relationship between equivalent load of structure and aerodynamic load is obtained, i.e.,(14)fs=HT·fa.The coupling diagram is shown in Figure2.Figure 2
Coupled transfer diagram.


## 2.5. Laminate Calculation

The classical laminate theory has the following assumptions: (1) the bond layer between the laminates is very thin and firmly bonded, and no slip occurs between the laminates. (2) The laminates are thin plates, the thickness of the laminates is constant, and the single-layer laminates are analyzed according to the plane stress state. (3) The bending deformation of the laminate is in the range of small deflection. Before deformation, it is perpendicular to the straight line of the middle plane. After deformation, it remains straight and perpendicular to the middle plane, and the length of the straight line remains unchanged.Suppose that the displacement of any pointC in the distance plane z along the three coordinate axes x, y, and z direction is ux,y,z, vx,y,z, and wx,y,z, respectively, then the classical laminate theory can be obtained [39]:
(15)ux,y,z=u0x,y−z∂wx,y∂x,vx,y,z=v0x,y−z∂wx,y∂y,wx,y=w0x,y,where u0x,y, v0x,y, and w0x,y are the displacements of the middle plane of the laminates.Then, the strain-displacement geometric relationship is as follows:(16)εx=∂u∂x,γyz=∂w∂y+∂z∂x,εy=∂v∂y,γzx=∂w∂x+∂u∂z,εz=∂w∂z,γxy=∂u∂y+∂v∂x.It is calculated as follows:(17)εx=∂u0∂x−z∂2w∂x2=εx0+zkx,εy=∂v0∂y−z∂2w∂y2=εy0+zky,γxy=∂u0∂x+∂v0∂y−2z∂2w∂x∂y=γxy0+zkxy,where εx0, εy0, and γxy0 are the strains in the middle of the laminates, kx and ky are the bending rates of the middle surface of the laminate, and kxy is the distortion rate in the middle of the laminate.The above equation can be expressed as a matrix:(18)εxεyγxy=εx0εy0γxy0+zkxkykxy.In the Cartesian coordinate system, the off-axis stress-strain relationship of the single-layer plate is as follows:(19)σxσyτxy=Q¯11Q¯12Q¯16Q¯12Q¯22Q¯26Q¯16Q¯26Q¯66kεxεyγxy.Among them,(20)Q11=E11−V12V21,Q22=E21−V12V21,Q12=V21E21−V12V21=V12E11−V12V21,Q¯66=G12,where σx and σy are the off-axis normal stresses, τxy is the partial shear stress, εx and εy are the off-axis normal strains, γxy is the off-axis shear strain, E1 is the longitudinal modulus of elasticity, E2 is the transverse modulus of elasticity, G12 is shear elastic modulus, V12 is the longitudinal Poisson’s ratio, V21 is the transverse Poisson’s ratio, and Q¯ is the transformation matrix of reduced stiffness matrix Q, which is expressed as follows:
(21)Q¯=T−1QT−1T,where T is the transpose matrix:
(22)T=cos2θsin2θ2cosθsinθsin2θcos2θ−2cosθsinθ−cosθsinθcosθsinθcos2θ−sin2θk.The stress-strain relationship of thek-layer single layer of the laminated plate in the coordinate system XYZ can be obtained by simultaneous equations (18) and (19):
(23)σxσyτx=Q¯kεx0εy0γxy0+zkxkykxy.By integrating the stress of each single-layer plate, the internal force of the laminated plate can be expressed as follows:(24)NxNyNxy=∑k=1n∫zk−1zkσxσyτxykdz=∑k=1n∫zk−1zkQ¯kεx0εy0γxy0+zkxkykxydz,MxMyMxy=∑k=1n∫zk−1zkσxσyτxykzdz=∑k=1n∫zk−1zkQ¯kεx0εy0γxy0+zkxkykxyzdz,where n is the number of single-layer plates and zk−1 and zk are the z coordinates of the bottom and top surfaces of the single-layer plate of layer k, respectively.Finally, the expression of internal force of composite laminates is arranged into a matrix form:(25)NxNyNxyMxMyMxy=A11A12A16B11B12B16A12A16B11B12B16A22A26B12B22B26A26A66B16B26B66B12B16D11D12D16B22B26D12D22D26B26B66D16D26D66εx0εy0γxy0kxkykxy.Among them,(26)Aij=∑k=1nQ¯ijkzk−zk−1,Bij=12∑k=1nQ¯ijkzk2−zk−12,Bij=13∑k=1nQ¯ijkzk3−zk−13,where A11, A12, and A22 are the stiffness coefficients between tensile (compressive) force and middle plane tensile (compressive) strain. A66 is the stiffness coefficient between the shear force and the shear strain in the middle plane. A12 and A26 are the coupling stiffness coefficients between shear and tension. B11, B12, and B22 are the coupling stiffness coefficients between tensile and bending. B66 is the coupling stiffness coefficient between shear and torsion. B12 and B26 are the coupling stiffness coefficients between tension and torsion or between shear and bending. D11, D12, and D22 are the stiffness coefficients between bending and bending rate. D66 is the coupling stiffness coefficient between torsion and torsion rate. D12 and D26 are the coupling stiffness coefficients between torsion and bending.


## 2.6. Validation of Calculation Methods

The HIRENASD wing model was selected as the research object (the wing model is shown in Figure3 [40]), and the loose coupling unidirectional fluid-solid coupling theory was adopted as the calculation method [41]. The airfoil is supercritical BAC3-11/R ES/30/21, the reference area of the airfoil is 0.3926 m2, the reference length is 0.3445 m, and the aspect ratio is 12.08.The longitudinal deformations of point A near the leading edge of the wingtip at angles of attack of -1.5°, 0°, 1.5°, 3°, and 4.5°, with a Reynolds number of 7×106 and a flight speed of 0.8 Mach were calculated, respectively.Figure 3
HIRENASD wing model and parameters.
(a)
HIRENASD wing parameters(b)
The position of point A on the modelThe unstructured tetrahedral mesh was used to divide the flow field, and the mesh at the leading edge of the wing was encrypted. A total of 2,581,247 nodes and 1,869,810 grid cells were generated. The plane of the fuselage root is set as the symmetric boundary condition, the wing and fuselage are set as the material boundary condition, and the others are set as the far-field pressure boundary condition. In the structural calculation, the main wing finite element model was divided by unstructured tetrahedral mesh, with a total of 133,215 nodes and 75,790 mesh elements generated. The boundary condition was the constraint of the wing root fixed support. The overall meshing details are shown in Figure4.Figure 4
HIRENASD wing finite element model.
(a)
Finite element model of flow field(b)
Finite element model of structureAfter calculation, the numerical simulation results agree with the experimental ones. When the angle of attack is 4.5°, there is a big difference, the maximum deformation of the test value is 25.1433 mm, the maximum deformation of the calculated value is 22.751 mm, and the deviation between the two is about 10%.The reason for the small deviation in the overall calculation may be that the deformation marks of the wing structure in the test and the deformation calculation points of the numerical calculation model cannot coincide in a strict sense, so there is a certain numerical deviation. In addition, with the gradual increase of the angle of attack of the wing, the flow separation caused by the wing surface becomes more and more obvious, and the computational grid is difficult to accurately solve, so the data also produce a certain deviation. In conclusion, this method has a high simulation accuracy within a certain angle of attack calculation range and is suitable for the static aeroelastic analysis of high aspect ratio wing. Specific experimental values can be referred to References [40, 42]. Figure 5 shows the comparison between the calculated results and the simulation results.Figure 5
Comparison of calculated results and experimental results.


## 3. Computational Model

### 3.1. Geometric and Material Parameters

In this paper, the wing of Global Hawk UAV is taken as the research object, and the wing mounting angle is not considered temporarily. The wing structure is a classic double-beam, multiribbed, upper and lower skin structure. The airfoil is NACA63212 laminar flow airfoil with a flight altitude of 15,000 meters, angle of attack is 2°, and a flight speed of 0.6 Mach when cruising. The specific parameters of the wing are shown in Table1, and the structure and wing profile are shown in Figure 6.Table 1
Geometrical parameters and airfoil parameters.
ParametersValueHalf spanb/216 mAspect ratioλ20.65Taper ratioη3.4Leading chord sweepχΘ6°Figure 6
Wing structure and airfoil profile.
(a)
Structures(b)
NACA63212 airfoilThe mechanical characteristics of the wing with high aspect ratio are that the bending moment at the wing root is relatively large, and the wing surface is prone to combined deformation of bending and torsion. According to the above characteristics, the reasonable application of composite materials can significantly improve the performance of the structure. Based on the above viewpoints, the wing skin was laminated with composite materials. Among them, the wing SPAR and wing ribs are made of aluminum alloy, the yield limit is 280 MPa, the safety factor is 2, the Young’s modulus of elasticity is 71 GPa, and the Poisson’s ratio is 0.33. UT700 is selected as the laminating material, and the specific material parameters are shown in Table2.Table 2
UT700 material parameters.
Density (kg·m3)Elastic modulus (GPa)Shear modulus (GPa)Poisson’s ratioExEyEzGxyGyzGxzPRxyPRxyPRxy18001156.436.436660.280.340.28


### 3.2. Layout Design

In the layup area, layup area I was defined as 40% and half of the wingspan along the wing root as the starting point, layup area II was defined as 20% of the extension along the wing tip, and layup area III was defined as the rest of the wingtip. The reference direction is the normal direction of the wing root plane. The specific division of layering area is shown in Figure7.Figure 7
Definition of layup area and reference direction.In layup design, in order to study the influence of the ratio of unbalanced coefficient on the layup wing structure, under the condition of gradually adjusting the unbalanced coefficient, the maximum wing tip displacements of 0° and ±45° lamination under constant thickness lamination and variable thickness lamination were studied, respectively. The maximum wing tip displacements with varying thickness and angle and varying thickness and angle are studied by adding the layup angle of 90°.In the case of constant thickness layup, it is assumed that the initial thickness of all-aluminum wing skin is 8 mm and the thickness of laminated plate is 0.2 mm per layer. Then, the reference composite layup scheme is shown in Table3.Table 3
Equal thickness layup scheme.
Unbalanced coefficientLayerThicknessLayer order0.2408.0 mm45°/0°/45°/0°/–45°2/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°s0.3408.0 mm45°/0°/45°/0°/45°/–45°/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°s0.4408.0 mm45°/0°/45°/0°/45°2/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°s0.5408.0 mm45°/0°/45°/0°/45°2/0°2/45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°s0.6408.0 mm45°/0°/45°/0°/45°2/0°2/45°/0°2/45°/–45°/0°2/–45°/0°/–45°/0°/–45°s0.7408.0 mm45°/0°/45°/0°/45°2/0°2/45°/0°2/45°2/0°2/–45°/0°/–45°/0°/–45°s0.8408.0 mm45°/0°/45°/0°/45°2/0°2/45°/0°2/45°2/0°2/45°/0°/–45°/0°/–45°sIn the layer with variable thickness, the proportion of 0° in each area was changed, and the size of the control unbalanced coefficient remained basically unchanged. The thickness at the wing root is the largest and gradually decreases along the span direction to the wing tip. The specific layup scheme is shown in Table4.Table 4
Variable thickness layout scheme.
Unbalanced coefficientAreaLayerThicknessLayer order0.2I408.0 mm45°/0°/45°/0°/–45°2/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/–45°/0°/–45°2/0°/–45°/0°/–45°/0°/–45°/0°/–45°2s0.3I408.0 mm45°/0°/45°/0°/45°/–45°/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/–45°2/0°/–45°/0°/–45°/0°/–45°/0°/–45°2s0.4I408.0 mm45°/0°/45°/0°/45°/45°/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/0°/45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2s0.5I408.0 mm45°/0°/45°/0°/45°/45°/0°2/45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/0°/–45°/0°/–45°/0°/–45°/0°/–45°2s0.6I408.0 mm45°/0°/45°/0°/45°/45°/0°2/45°/0°2/45°/–45°/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/0°/45°/0°/–45°/0°/–45°/0°/–45°2s0.7I408.0 mm45°/0°/45°/0°/45°/45°/0°2/45°/0°2/45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/0°/45°/0°/45°/0°/–45°/0°/–45°2s0.8I408.0 mm45°/0°/45°/0°/45°/45°/0°2/45°/0°2/45°2/0°2/45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/0°/45°/0°/45°/0°/45°/0°/–45°2sAfter adding 90° layup angle on the basis of variable thickness layup, the proportion of 0° layup angle in the overall layup angle will be reduced accordingly. Controlling 0°and 90° in each layer of the same proportion, both are 50%. At the same time, the size of the unbalanced coefficient is basically unchanged. The specific layup scheme is shown in Table5.Table 5
Variable thickness and angle layup scheme.
Unbalanced coefficientAreaLayerThicknessLayer order0.2I408.0 mm45°/0°/45°/0°/–45°2/0°/90°/–45°/90°2/–45°2/90°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/–45°/90°/–45°/90°/–45°/90°/–45°/90°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/–45°/0°/–45°2/90°/–45°/90°/–45°/90°/–45°/0°/–45°2s0.3I408.0 mm45°/0°/45°/0°/45°/–45°/0°/90°/–45°/90°2/–45°2/90°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/90°/–45°/90°/–45°/90°/–45°/90°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/–45°2/90°/–45°/90°/–45°/90°/–45°/90°/–45°2s0.4I408.0 mm45°/0°/45°/0°/45°/45°/0°/90°/–45°/90°2/–45°2/90°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/90°/45°/90°/–45°/90°/–45°/90°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°/–45°/90°/–45°/90°/–45°/90°/–45°/0°/–45°2s0.5I408.0 mm45°/0°/45°/0°/45°2/0°/90°/45°/90°2/–45°2/90°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/90°/45°/90°/45°/90°/–45°/90°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/90°/–45°/90°/–45°/90°/–45°/0°/–45°2s0.6I408.0 mm45°/0°/45°/0°/45°2/0°/90°/45°/90°2/45°/–45°/90°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/90°/45°/90°/45°/90°/45°/90°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/90°/45°/90°/–45°/90°/–45°/0°/–45°2s0.7I408.0 mm45°/0°/45°/0°/45°2/0°/90°/45°/90°2/45°2/90°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/90°/45°/90°/45°/90°/45°/90°/45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/90°/45°/90°/45°/90°/–45°/0°/–45°2s0.8I408.0 mm45°/0°/45°/0°/45°2/0°/90°/45°/90°2/45°2/90°2/45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/90°/45°/90°/45°/90°/45°/90°/45°/0°/45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/90°/45°/90°/45°/90°/45°/0°/–45°2sFinally, under the condition of changing the unbalanced coefficient, it is assumed that the layup angles of 90° are fixed four. When four 90° I layer angle is located in the area, is the area I orientations of 0°, ±45°, 90°, the rest of the two area layer angle only 0° and ±45°. Then the 90° layup angle is added to region II; the regions I and III are kept at 0° and ±45°, respectively. Similarly, to guarantee the area I, II is only 0° and ±45°, and area III contains added 90° laid angle. This is to determine the effect of the layup angle on each wing segment. Take the case of unbalanced coefficient 0.2 as an example, the specific scheme is shown in Table6, and the rest of the cases are similar.Table 6
Variable thickness and angle layup scheme example.
Unbalanced coefficientAreaLayerThicknessLayer order0.2I408.0 mm45°/0°/45°/0°/–45°2/0°/90°/–45°/90°2/–45°2/90°/0°/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/–45°/0°/–45°2/0°/–45°/0°/–45°/0°/–45°/0°/–45°2s0.2I408.0 mm45°/0°/45°/0°/–45°2/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/–45°/90°/–45°/90°/–45°/90°/–45°/90°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/–45°/0°/–45°2/0°/–45°/0°/–45°/0°/–45°/0°/–45°2s0.2I408.0 mm45°/0°/45°/0°/–45°2/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/90°/–45°2/90°/–45°/90°/–45°/90°/–45°/0°/–45°2s


### 3.3. Calculation and Solution

When the aerodynamic load is calculated, the flow field is solved and calculated based on the Spalart-Allmaras turbulence model. In the spatial discrete term, the viscous flux vector is selected as the central format, and the relative flux vector is selected as the Roe-FDS format. In the time discrete term, the LU-SGS implicit time discrete method is used to advance the solution [43]. Unstructured grids are used in the flow field calculation domain and the wing, and the grids near the object surface are refined and encrypted. Combined with the actual calculation conditions and through the independent verification of grid and iteration steps, a total of 986,262 nodes and 717,173 grid cells were finally generated. The boundary conditions at the wing root are set as symmetric boundary conditions, the wing surface is set as object surface boundary conditions, and the rest are set as far-field boundary conditions of pressure. The grid division of the computational domain of the wing flow field is shown in Figure 8.Figure 8
Calculation domain and mesh of flow field.
(a)
Overall computing domain(b)
Grid near the wallConsidering the influence of geometric nonlinearity, the Newton-Raphson iterative method is used for numerical analysis of structural deformation [44]. The wing finite element structure model was divided by unstructured tetrahedral mesh, and the wing root was applied with fixed constraints. A total of 121,171 nodes and 86,675 mesh elements were obtained. Secondly, the aerodynamic loads solved in the above flow fields are interpolated to the nodes of the structural grid. The specific grid partitioning and aerodynamic interpolation are shown in Figure 9.Figure 9
Mesh partitioning and aerodynamic interpolation results.
(a)
The results of finite element grid differentiation(b)
Wing aerodynamic load interpolation


## 3.1. Geometric and Material Parameters

In this paper, the wing of Global Hawk UAV is taken as the research object, and the wing mounting angle is not considered temporarily. The wing structure is a classic double-beam, multiribbed, upper and lower skin structure. The airfoil is NACA63212 laminar flow airfoil with a flight altitude of 15,000 meters, angle of attack is 2°, and a flight speed of 0.6 Mach when cruising. The specific parameters of the wing are shown in Table1, and the structure and wing profile are shown in Figure 6.Table 1
Geometrical parameters and airfoil parameters.
ParametersValueHalf spanb/216 mAspect ratioλ20.65Taper ratioη3.4Leading chord sweepχΘ6°Figure 6
Wing structure and airfoil profile.
(a)
Structures(b)
NACA63212 airfoilThe mechanical characteristics of the wing with high aspect ratio are that the bending moment at the wing root is relatively large, and the wing surface is prone to combined deformation of bending and torsion. According to the above characteristics, the reasonable application of composite materials can significantly improve the performance of the structure. Based on the above viewpoints, the wing skin was laminated with composite materials. Among them, the wing SPAR and wing ribs are made of aluminum alloy, the yield limit is 280 MPa, the safety factor is 2, the Young’s modulus of elasticity is 71 GPa, and the Poisson’s ratio is 0.33. UT700 is selected as the laminating material, and the specific material parameters are shown in Table2.Table 2
UT700 material parameters.
Density (kg·m3)Elastic modulus (GPa)Shear modulus (GPa)Poisson’s ratioExEyEzGxyGyzGxzPRxyPRxyPRxy18001156.436.436660.280.340.28


## 3.2. Layout Design

In the layup area, layup area I was defined as 40% and half of the wingspan along the wing root as the starting point, layup area II was defined as 20% of the extension along the wing tip, and layup area III was defined as the rest of the wingtip. The reference direction is the normal direction of the wing root plane. The specific division of layering area is shown in Figure7.Figure 7
Definition of layup area and reference direction.In layup design, in order to study the influence of the ratio of unbalanced coefficient on the layup wing structure, under the condition of gradually adjusting the unbalanced coefficient, the maximum wing tip displacements of 0° and ±45° lamination under constant thickness lamination and variable thickness lamination were studied, respectively. The maximum wing tip displacements with varying thickness and angle and varying thickness and angle are studied by adding the layup angle of 90°.In the case of constant thickness layup, it is assumed that the initial thickness of all-aluminum wing skin is 8 mm and the thickness of laminated plate is 0.2 mm per layer. Then, the reference composite layup scheme is shown in Table3.Table 3
Equal thickness layup scheme.
Unbalanced coefficientLayerThicknessLayer order0.2408.0 mm45°/0°/45°/0°/–45°2/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°s0.3408.0 mm45°/0°/45°/0°/45°/–45°/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°s0.4408.0 mm45°/0°/45°/0°/45°2/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°s0.5408.0 mm45°/0°/45°/0°/45°2/0°2/45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°s0.6408.0 mm45°/0°/45°/0°/45°2/0°2/45°/0°2/45°/–45°/0°2/–45°/0°/–45°/0°/–45°s0.7408.0 mm45°/0°/45°/0°/45°2/0°2/45°/0°2/45°2/0°2/–45°/0°/–45°/0°/–45°s0.8408.0 mm45°/0°/45°/0°/45°2/0°2/45°/0°2/45°2/0°2/45°/0°/–45°/0°/–45°sIn the layer with variable thickness, the proportion of 0° in each area was changed, and the size of the control unbalanced coefficient remained basically unchanged. The thickness at the wing root is the largest and gradually decreases along the span direction to the wing tip. The specific layup scheme is shown in Table4.Table 4
Variable thickness layout scheme.
Unbalanced coefficientAreaLayerThicknessLayer order0.2I408.0 mm45°/0°/45°/0°/–45°2/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/–45°/0°/–45°2/0°/–45°/0°/–45°/0°/–45°/0°/–45°2s0.3I408.0 mm45°/0°/45°/0°/45°/–45°/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/–45°2/0°/–45°/0°/–45°/0°/–45°/0°/–45°2s0.4I408.0 mm45°/0°/45°/0°/45°/45°/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/0°/45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2s0.5I408.0 mm45°/0°/45°/0°/45°/45°/0°2/45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/0°/–45°/0°/–45°/0°/–45°/0°/–45°2s0.6I408.0 mm45°/0°/45°/0°/45°/45°/0°2/45°/0°2/45°/–45°/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/0°/45°/0°/–45°/0°/–45°/0°/–45°2s0.7I408.0 mm45°/0°/45°/0°/45°/45°/0°2/45°/0°2/45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/0°/45°/0°/45°/0°/–45°/0°/–45°2s0.8I408.0 mm45°/0°/45°/0°/45°/45°/0°2/45°/0°2/45°2/0°2/45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/0°/45°/0°/45°/0°/45°/0°/–45°2sAfter adding 90° layup angle on the basis of variable thickness layup, the proportion of 0° layup angle in the overall layup angle will be reduced accordingly. Controlling 0°and 90° in each layer of the same proportion, both are 50%. At the same time, the size of the unbalanced coefficient is basically unchanged. The specific layup scheme is shown in Table5.Table 5
Variable thickness and angle layup scheme.
Unbalanced coefficientAreaLayerThicknessLayer order0.2I408.0 mm45°/0°/45°/0°/–45°2/0°/90°/–45°/90°2/–45°2/90°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/–45°/90°/–45°/90°/–45°/90°/–45°/90°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/–45°/0°/–45°2/90°/–45°/90°/–45°/90°/–45°/0°/–45°2s0.3I408.0 mm45°/0°/45°/0°/45°/–45°/0°/90°/–45°/90°2/–45°2/90°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/90°/–45°/90°/–45°/90°/–45°/90°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/–45°2/90°/–45°/90°/–45°/90°/–45°/90°/–45°2s0.4I408.0 mm45°/0°/45°/0°/45°/45°/0°/90°/–45°/90°2/–45°2/90°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/90°/45°/90°/–45°/90°/–45°/90°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°/–45°/90°/–45°/90°/–45°/90°/–45°/0°/–45°2s0.5I408.0 mm45°/0°/45°/0°/45°2/0°/90°/45°/90°2/–45°2/90°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/90°/45°/90°/45°/90°/–45°/90°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/90°/–45°/90°/–45°/90°/–45°/0°/–45°2s0.6I408.0 mm45°/0°/45°/0°/45°2/0°/90°/45°/90°2/45°/–45°/90°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/90°/45°/90°/45°/90°/45°/90°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/90°/45°/90°/–45°/90°/–45°/0°/–45°2s0.7I408.0 mm45°/0°/45°/0°/45°2/0°/90°/45°/90°2/45°2/90°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/90°/45°/90°/45°/90°/45°/90°/45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/90°/45°/90°/45°/90°/–45°/0°/–45°2s0.8I408.0 mm45°/0°/45°/0°/45°2/0°/90°/45°/90°2/45°2/90°2/45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/45°/90°/45°/90°/45°/90°/45°/90°/45°/0°/45°/0°/–45°2sIII326.4 mm45°2/0°/45°/0°/45°2/90°/45°/90°/45°/90°/45°/0°/–45°2sFinally, under the condition of changing the unbalanced coefficient, it is assumed that the layup angles of 90° are fixed four. When four 90° I layer angle is located in the area, is the area I orientations of 0°, ±45°, 90°, the rest of the two area layer angle only 0° and ±45°. Then the 90° layup angle is added to region II; the regions I and III are kept at 0° and ±45°, respectively. Similarly, to guarantee the area I, II is only 0° and ±45°, and area III contains added 90° laid angle. This is to determine the effect of the layup angle on each wing segment. Take the case of unbalanced coefficient 0.2 as an example, the specific scheme is shown in Table6, and the rest of the cases are similar.Table 6
Variable thickness and angle layup scheme example.
Unbalanced coefficientAreaLayerThicknessLayer order0.2I408.0 mm45°/0°/45°/0°/–45°2/0°/90°/–45°/90°2/–45°2/90°/0°/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/–45°/0°/–45°2/0°/–45°/0°/–45°/0°/–45°/0°/–45°2s0.2I408.0 mm45°/0°/45°/0°/–45°2/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/–45°/90°/–45°/90°/–45°/90°/–45°/90°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/–45°/0°/–45°2/0°/–45°/0°/–45°/0°/–45°/0°/–45°2s0.2I408.0 mm45°/0°/45°/0°/–45°2/0°2/–45°/0°2/–45°2/0°2/–45°/0°/–45°/0°/–45°sII367.2 mm45°/0°/45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°/0°/–45°2sIII326.4 mm45°2/0°/45°/90°/–45°2/90°/–45°/90°/–45°/90°/–45°/0°/–45°2s


## 3.3. Calculation and Solution

When the aerodynamic load is calculated, the flow field is solved and calculated based on the Spalart-Allmaras turbulence model. In the spatial discrete term, the viscous flux vector is selected as the central format, and the relative flux vector is selected as the Roe-FDS format. In the time discrete term, the LU-SGS implicit time discrete method is used to advance the solution [43]. Unstructured grids are used in the flow field calculation domain and the wing, and the grids near the object surface are refined and encrypted. Combined with the actual calculation conditions and through the independent verification of grid and iteration steps, a total of 986,262 nodes and 717,173 grid cells were finally generated. The boundary conditions at the wing root are set as symmetric boundary conditions, the wing surface is set as object surface boundary conditions, and the rest are set as far-field boundary conditions of pressure. The grid division of the computational domain of the wing flow field is shown in Figure 8.Figure 8
Calculation domain and mesh of flow field.
(a)
Overall computing domain(b)
Grid near the wallConsidering the influence of geometric nonlinearity, the Newton-Raphson iterative method is used for numerical analysis of structural deformation [44]. The wing finite element structure model was divided by unstructured tetrahedral mesh, and the wing root was applied with fixed constraints. A total of 121,171 nodes and 86,675 mesh elements were obtained. Secondly, the aerodynamic loads solved in the above flow fields are interpolated to the nodes of the structural grid. The specific grid partitioning and aerodynamic interpolation are shown in Figure 9.Figure 9
Mesh partitioning and aerodynamic interpolation results.
(a)
The results of finite element grid differentiation(b)
Wing aerodynamic load interpolation


## 4. Result Analysis

### 4.1. The Relation between Unbalanced Coefficient and Thickness

Figure10 shows the variation relationship between the unbalanced coefficient and the wing tip displacement. With the increase of the unbalanced coefficient, the wingtip displacement first decreases and then increases, and presents a symmetric trend along the unbalanced coefficient of 0.5 (i.e., symmetric layup). When the unbalanced coefficient is 0.5, the structural deformation is the least. Secondly, in this example, the maximum wing tip displacement calculated in the case of constant thickness layup is similar to that in the case of variable thickness layup. The distributions of the longitudinal displacements and the broad-span displacements at the leading edge of the wing are basically the same, which is due to the small variation of the gradient with variable thickness. The distribution of chord displacement is different with the change of unbalanced coefficient.Figure 10
Deformation of airfoil with changing thickness layup.
(a)
Maximum wing displacement(b)
Variation of leading edge longitudinal displacement(c)
Variation of leading edge wingspan displacement(d)
Variation of leading edge chordal displacementWhen the unbalanced coefficient is 0.3, it shows a linear deformation relationship along the wingspan under constant thickness and variable thickness. When the unbalanced coefficient is 0.5 and 0.7, the chord deformation no longer obeys the linear relationship, and its change situation increases first and then decreases, and the peak value of growth occurs at about half of the half-length. With the increase of the unbalanced coefficient, the chordal displacement of the wingtip gradually decreases. Therefore, in general composite wing skin-laying, when only 0° and ±45° lamination angles are considered, symmetrical lamination should be used as far as possible. Moreover, the thickness of lamination varies slightly along the span direction, so variable thickness lamination should be used to reduce the weight of wing structure.


### 4.2. The Relation between Unbalanced Coefficient and Angle

The influence of unbalanced coefficient and lamination angle on structural deformation is shown in Figure11. The maximum displacement of the wing structure is significantly reduced by adding 90° layup angle. Similarly, after adding the layup angle of 90°, the overall deformation still obeys the law of symmetrical distribution. When the unbalanced coefficient is 0.5, the maximum displacement of the wing tip is the smallest, and the displacements under the other coefficients are symmetrically distributed with the unbalanced coefficient of 0.5.Taking the unbalanced coefficient 0.5 as an example, the influence of the change of the unbalanced coefficient on the wing displacement in all directions is shown in Figures 10(b)–10(d).Figure 11
Deformation of the wing with the addition of angle layering.
(a)
Maximum wing displacement(b)
Variation of leading edge longitudinal displacement(c)
Variation of leading edge wingspan displacement(d)
Variation of leading edge chordal displacementThe significant reduction of displacement is due to the fact that the original laminates laid at 0° and ±45° only provide directional tensile and shear strength, and the addition of additional 90° laminates along different directions can produce quasi-isotropic laminates. Compared with the laminates with 0° and ±45°, the multiangle composite laminates with 90° have stronger deformation-resistant and damage-resistant abilities. Therefore, the multiangle composite laminates have stronger resistance to deformation than the double-angle laminates.


### 4.3. The Relation between Unbalanced Coefficient and Layup Area

The size of unbalanced coefficient and the selection of layup angle have different effects on different layup areas. As shown in Figure12, taking the unbalanced coefficient of 0.5 as an example, the outermost segment of the wing has the largest deformation, while the innermost segment of the wing has the smallest deformation. The influence of the change of the unbalanced coefficient on the lamination area III is smaller and larger than that of the lamination area. Compared with the general wing, the high aspect ratio wing has greater flexibility and lower stiffness, so the influence of geometric nonlinearity cannot be ignored. The layup area III is located at the outer end of the wing, and the geometric deformation is relatively larger, so after the composite layup design, there is a more obvious deformation change.Figure 12
Influence of laminated area on wing deformation.
(a)
Maximum wing displacement(b)
Variation of leading edge longitudinal displacement(c)
Variation of leading edge wingspan displacement(d)
Variation of leading edge chordal displacementFigure13, respectively, shows the displacement changes in all directions in the layup area I, and tunnel with the change of the unbalanced coefficient. Figures 13(a)–13(c) describe the longitudinal displacement of the wing under the change of the unbalanced coefficient. The maximum vertical displacement in the laying area I is 1176.6 mm, and the maximum vertical displacement in the laying area is 1407.8 mm and 1454.4 mm, respectively. The difference between the two is small. Secondly, when the unbalanced coefficient is 0.8 or 0.2 in all the bedding areas, the displacement reaches the maximum. When the unbalanced coefficient is 0.5, its displacement reaches the minimum value. Figures 13(d)–13(f) describe the variation of displacement in the spreading direction of each laminated area. The maximum displacement in the spreading direction of the bedding area I is -53.9 mm, the maximum displacement in the spreading direction of the bedding area is -74.5 mm, and the maximum displacement in the spreading direction of the bedding area is -79.9 mm. The variation rule is basically the same as that in the longitudinal direction, but the spanwise displacement when the unbalanced coefficient is 0.8 is slightly larger than that when the unbalanced coefficient is 0.2. Figures 13(g)–13(i) show that the chord direction wing displacement changes with the change of the layup area. It can be seen that, compared with the longitudinal displacement and the spanwise displacement, the chord direction is greatly affected by the change of the unbalanced coefficient. When the unbalanced coefficient is 0.8, as the layup area of the wing segment changes, the chordal displacement change in the layup area II is smaller than that in the layup area; that is, the chordal displacement change in the layup coupling area in the middle of the wing segment is smaller than that in the coupling area of the outer wing segment and the inner wing segment. Similarly, the unbalanced coefficient increases in the three layup regions, the chordal displacement changes from linear to nonlinear along the span direction, and the chordal deformation gradually decreases along the wing span direction. In general, the change of unbalanced coefficient has a certain influence on the division of the wing layup area. The specific causes are still under further research and exploration.Figure 13
Displacement comparison in different directions in each laminated area.
(a)
Leading edge longitudinal displacement in area I(b)
Leading edge longitudinal displacement in area II(c)
Leading edge longitudinal displacement in area III(d)
Leading edge wingspan displacement in area I(e)
Leading edge wingspan displacement in area II(f)
Leading edge wingspan displacement in area III(g)
Leading edge chord displacement in area I(h)
Leading edge chord displacement in area II(i)
Leading edge chord displacement in area III


## 4.1. The Relation between Unbalanced Coefficient and Thickness

Figure10 shows the variation relationship between the unbalanced coefficient and the wing tip displacement. With the increase of the unbalanced coefficient, the wingtip displacement first decreases and then increases, and presents a symmetric trend along the unbalanced coefficient of 0.5 (i.e., symmetric layup). When the unbalanced coefficient is 0.5, the structural deformation is the least. Secondly, in this example, the maximum wing tip displacement calculated in the case of constant thickness layup is similar to that in the case of variable thickness layup. The distributions of the longitudinal displacements and the broad-span displacements at the leading edge of the wing are basically the same, which is due to the small variation of the gradient with variable thickness. The distribution of chord displacement is different with the change of unbalanced coefficient.Figure 10
Deformation of airfoil with changing thickness layup.
(a)
Maximum wing displacement(b)
Variation of leading edge longitudinal displacement(c)
Variation of leading edge wingspan displacement(d)
Variation of leading edge chordal displacementWhen the unbalanced coefficient is 0.3, it shows a linear deformation relationship along the wingspan under constant thickness and variable thickness. When the unbalanced coefficient is 0.5 and 0.7, the chord deformation no longer obeys the linear relationship, and its change situation increases first and then decreases, and the peak value of growth occurs at about half of the half-length. With the increase of the unbalanced coefficient, the chordal displacement of the wingtip gradually decreases. Therefore, in general composite wing skin-laying, when only 0° and ±45° lamination angles are considered, symmetrical lamination should be used as far as possible. Moreover, the thickness of lamination varies slightly along the span direction, so variable thickness lamination should be used to reduce the weight of wing structure.


## 4.2. The Relation between Unbalanced Coefficient and Angle

The influence of unbalanced coefficient and lamination angle on structural deformation is shown in Figure11. The maximum displacement of the wing structure is significantly reduced by adding 90° layup angle. Similarly, after adding the layup angle of 90°, the overall deformation still obeys the law of symmetrical distribution. When the unbalanced coefficient is 0.5, the maximum displacement of the wing tip is the smallest, and the displacements under the other coefficients are symmetrically distributed with the unbalanced coefficient of 0.5.Taking the unbalanced coefficient 0.5 as an example, the influence of the change of the unbalanced coefficient on the wing displacement in all directions is shown in Figures 10(b)–10(d).Figure 11
Deformation of the wing with the addition of angle layering.
(a)
Maximum wing displacement(b)
Variation of leading edge longitudinal displacement(c)
Variation of leading edge wingspan displacement(d)
Variation of leading edge chordal displacementThe significant reduction of displacement is due to the fact that the original laminates laid at 0° and ±45° only provide directional tensile and shear strength, and the addition of additional 90° laminates along different directions can produce quasi-isotropic laminates. Compared with the laminates with 0° and ±45°, the multiangle composite laminates with 90° have stronger deformation-resistant and damage-resistant abilities. Therefore, the multiangle composite laminates have stronger resistance to deformation than the double-angle laminates.


## 4.3. The Relation between Unbalanced Coefficient and Layup Area

The size of unbalanced coefficient and the selection of layup angle have different effects on different layup areas. As shown in Figure12, taking the unbalanced coefficient of 0.5 as an example, the outermost segment of the wing has the largest deformation, while the innermost segment of the wing has the smallest deformation. The influence of the change of the unbalanced coefficient on the lamination area III is smaller and larger than that of the lamination area. Compared with the general wing, the high aspect ratio wing has greater flexibility and lower stiffness, so the influence of geometric nonlinearity cannot be ignored. The layup area III is located at the outer end of the wing, and the geometric deformation is relatively larger, so after the composite layup design, there is a more obvious deformation change.Figure 12
Influence of laminated area on wing deformation.
(a)
Maximum wing displacement(b)
Variation of leading edge longitudinal displacement(c)
Variation of leading edge wingspan displacement(d)
Variation of leading edge chordal displacementFigure13, respectively, shows the displacement changes in all directions in the layup area I, and tunnel with the change of the unbalanced coefficient. Figures 13(a)–13(c) describe the longitudinal displacement of the wing under the change of the unbalanced coefficient. The maximum vertical displacement in the laying area I is 1176.6 mm, and the maximum vertical displacement in the laying area is 1407.8 mm and 1454.4 mm, respectively. The difference between the two is small. Secondly, when the unbalanced coefficient is 0.8 or 0.2 in all the bedding areas, the displacement reaches the maximum. When the unbalanced coefficient is 0.5, its displacement reaches the minimum value. Figures 13(d)–13(f) describe the variation of displacement in the spreading direction of each laminated area. The maximum displacement in the spreading direction of the bedding area I is -53.9 mm, the maximum displacement in the spreading direction of the bedding area is -74.5 mm, and the maximum displacement in the spreading direction of the bedding area is -79.9 mm. The variation rule is basically the same as that in the longitudinal direction, but the spanwise displacement when the unbalanced coefficient is 0.8 is slightly larger than that when the unbalanced coefficient is 0.2. Figures 13(g)–13(i) show that the chord direction wing displacement changes with the change of the layup area. It can be seen that, compared with the longitudinal displacement and the spanwise displacement, the chord direction is greatly affected by the change of the unbalanced coefficient. When the unbalanced coefficient is 0.8, as the layup area of the wing segment changes, the chordal displacement change in the layup area II is smaller than that in the layup area; that is, the chordal displacement change in the layup coupling area in the middle of the wing segment is smaller than that in the coupling area of the outer wing segment and the inner wing segment. Similarly, the unbalanced coefficient increases in the three layup regions, the chordal displacement changes from linear to nonlinear along the span direction, and the chordal deformation gradually decreases along the wing span direction. In general, the change of unbalanced coefficient has a certain influence on the division of the wing layup area. The specific causes are still under further research and exploration.Figure 13
Displacement comparison in different directions in each laminated area.
(a)
Leading edge longitudinal displacement in area I(b)
Leading edge longitudinal displacement in area II(c)
Leading edge longitudinal displacement in area III(d)
Leading edge wingspan displacement in area I(e)
Leading edge wingspan displacement in area II(f)
Leading edge wingspan displacement in area III(g)
Leading edge chord displacement in area I(h)
Leading edge chord displacement in area II(i)
Leading edge chord displacement in area III


## 5. Conclusion

In this paper, the unidirectional fluid-structure coupling calculation method based on loose coupling is used to study the influence of the change of unbalanced coefficient on the static aeroelastic characteristics of the wing with high aspect ratio under the condition of geometric nonlinearity, and the following conclusions are drawn:(1)
There is not only a single corresponding relationship between the maximum wing deformation and the change of unbalanced coefficient. However, with the increase of the unbalanced coefficient, it shows a trend of decreasing first and then increasing. The minimum deformation of the wing occurs when the unbalanced coefficient is 0.5 (i.e., the opposite lamination), and the deformation of the wing is symmetrically distributed on both sides when the unbalanced coefficient is 0.5(2)
When the layup angle only has 0° and ±45°, the wing deformation obtained by the two schemes is basically the same as that obtained by the same thickness layup method and the variable thickness layup method. So, in the design of wing skin layup, variable thickness layup should be adopted as far as possible, which is beneficial to the weight reduction of wing structure(3)
Compared with the double-ply angle-ply scheme, the mixed angle-ply scheme has a greater inhibition on the wing deformation. This is because the mechanical properties of laminates will be changed by multiangle laminating, which will change the bending and torsion properties of laminates(4)
The effect of layup angle on the layup area III (outer wing section) was greater than that of the other two areas. Different wing segments have different sensibilities to the 90° layup angle. The specific layup sequence, layup angle, and the layup of wing segment should be selected according to the mission requirements of aircraft design(5)
The change of unbalanced coefficient will significantly change the chord displacement of the wing. With the increase of unbalanced coefficient, the chord displacement along the spanwise direction changes from linear distribution to nonlinear distribution. The abrupt change is generally located at the average half span of the wing, and the chordal displacement decreases along the wing span



---

*Source: 3949078-2021-09-15.xml*