Effect of ply angle on nonlinear static aeroelasticity of high-aspect-ratio composite wing

2.1. Computation methodology

The unidirectional fluid-solid coupling calculation method was adopted, and the N-S equation was used as the governing equation to calculate the flow field to obtain the aerodynamic load. After the flow field calculation was converged, the aerodynamic load was loaded onto the structural grid node through the numerical interpolation technique. The finite element method was used to solve the wing structure, and the deformation and stress of the wing under the aerodynamic load were obtained. The coupling calculation process is shown in Fig. 1.

2.1.1. Flow field solution technique

During the static aeroelastic solution process, the deformation speed is slow. The influence of acceleration and deformation speed on aerodynamics can be ignored. The integral form of the N-S equation has been used as the governing equation, and the expression in the Cartesian coordinate system is [13]:

1

$\frac{\partial }{\partial t}\underset{\Omega }{\iiint }\stackrel{-}{Q}d\mathrm{\Omega }+\underset{\partial \Omega }{\iint }\stackrel{-}{\stackrel{-}{F}}\cdot d\overrightarrow{S}=\frac{1}{R_e}\underset{\partial \Omega }{\iint }{\stackrel{-}{\stackrel{-}{F}}}^V\cdot d\overrightarrow{S},$

$\stackrel{-}{Q}=\left[\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho w\\ \rho E\end{array}\right],\overline{\overline{F}}=\left[\begin{array}{ccc}\rho u& \rho v& \rho w\\ \rho u^2+p& \rho vu& \rho wu\\ \rho uv& \rho v^2+p& \rho wv\\ \rho uw& \rho vw& \rho w^2+p\\ \rho uH& \rho vH& \rho wH\end{array}\right],{\stackrel{-}{\stackrel{-}{F}}}^V=\left[\begin{array}{ccc}0& 0& 0\\ {\tau }_{xx}& {\tau }_{yx}& {\tau }_{zx}\\ {\tau }_{xy}& {\tau }_{yy}& {\tau }_{zy}\\ {\tau }_{xz}& {\tau }_{yz}& {\tau }_{zz}\\ {\phi }_x& {\phi }_y& {\phi }_z\end{array}\right],$

where, $\stackrel{-}{Q}$ is the flow field conservation variable term. $\stackrel{-}{\stackrel{-}{F}}$ and ${\stackrel{-}{\stackrel{-}{F}}}^V$ are the flux term and the dissipation flux term, respectively. $\mathrm{\Omega }$ is the limit of the control volume. $\partial \mathrm{\Omega }$ is the control unit. $d\mathrm{\Omega }$ and $dS$ are the volume element and the area element, respectively. After introducing the Stokes hypothesis and the gas state equation $p=\rho RT$, the system of equations becomes closed. The Spalart-Allmaras turbulence model [14] was selected to simulate the flow field. The convection term was discretized using the second-order accuracy Roe-FDS [15], and the viscous terms were discretized by central differential interpolation [16].

Fig. 1Unidirectional fluid-solid coupling calculation process

For the steady flow, the discrete form of the governing Eq. (1) is:

2

$\frac{d}{dt}\left({\overrightarrow{Q}}_{i,j,k}\right)+{\overrightarrow{R}}_{i,j,k}=0,\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\overrightarrow{R}}_{i,j,k}=\frac{1}{V_{{\mathrm{}}_{i,j,k}}}\left({\overrightarrow{F}}_{{\mathrm{}}_{i,j,k}}-\frac{1}{R_e}{\overrightarrow{F}}_{{\mathrm{}}_{i,j,k}}^V\right),$

where, ${\overrightarrow{R}}_{i,j,k}$ is the residual value of the semi-discrete equation. ${\overrightarrow{F}}_{{\mathrm{}}_{i,j,k}}$ and ${\overrightarrow{F}}_{{\mathrm{}}_{i,j,k}}^V$ are the non-viscous flux term and the viscous flux term, respectively. The Eq. (2) can be solved by the LU-SGS implicit time-discrete method with first-order precision [17].

2.2. Mathematical model for optimization of ply angle

The ply angle $\theta$ was taken as the design variable. The wing beam and rib stress satisfied the strength condition, which was taken as the constraint. The composite skin of the wing satisfies the Tsai-Wu failure condition was taken as the constraint similarly. The maximum deformation of the wing is less than 15 % of the half-span, which was taken as the optimization objective. The ply angle for the composite skin of the high-aspect-ratio composite wing was optimized by the Screening method [18]. Its mathematical expression is as shown in Eq. (5):

where ${\theta }_j=\left[-\mathrm{45,0},\mathrm{45,90}\right](j=\mathrm{1,2},\cdot \cdot \cdot i)$, and $Max$$def$ is the maximum deformation. $Max$$str$ is the maximum wing beam and rib stress. $Tsai$ - $Wu$ is the maximum Tsai-Wu failure factor of the composite wing skin.

2.3. Computation methodology verification

Taking the HIRENASD wing model as the research object, the accuracy of the proposed static aeroelastic calculation method was verified. The geometry of the HIRENASD wing model [19] is shown in Fig. 2. The reference area is 0.3926 m²and the reference length is 0.3445 m [20]. Fig. 3 is the position of point $A$ to compare with experimental data.

Fig. 2The geometry of HIRENASD wing model

Fig. 3The position of point A compared with the experimental data

In the static aeroelastic calculation, the flow field adopted a hybrid mesh, and the mesh near the object surface was refined. In order to make more realistic fluid simulation, the inflation layer was inserted around the wing. A five-layer hexahedron or wedge mesh was generated around the wing, and then the rest tetrahedral unstructured meshes were made. The number of generated nodes was 793809, and the number of grids was 3308576. The grid near the object surface is shown in Fig. 4. The plane of the fuselage root was set to a symmetrical boundary condition, the wing and the fuselage were set to the object boundary condition, and the other ones were set to the pressure far away from the field boundary condition. The HIRENASD wing finite element structure model was meshed by unstructural tetrahedrons, with a total of 90963 nodes and 51747 grid elements. The finite element model of the wing was constrained by the fixed support of the wing root, as shown in Fig. 5.

Fig. 4Grid near the surface of the HIRENASD model

Fig. 5HIRENASD wing structure finite element mesh

The calculated Mach number was 0.8, and the Reynolds number was 7×10⁶. The angle of attack were –1.5°, 0°, 1.5°, 3° and 4.5°, respectively. Fig. 6 shows a comparison of the deformation of the point $A$ of the wing in the $Z$ direction at different angles of attack between calculation and experiment [21]. It can be seen that the calculated values agree well with the experimental values. It is shown that the static aeroelasticity numerical simulation method can be used for the analysis of static aeroelastic problems.

Fig. 6Deformation of the point A of the HIRENASD wing in the Z direction at different angles of attack

2.4. Composite wing ply setting and meshing

The high-aspect-ratio wing was adopted for analysis. The wing airfoil is the ONERA D section. The wing is a swept wing with zero upper and reverse angles. The geometric model and dimensions of the wing are shown in Fig. 7. The reference area is 3.675 m², and the reference length is 0.6125 m. The material of the wing beam and rib is titanium alloy. The yield strength is 930 MPa, and the safety factor is 1.5. The front beam has a thickness of 30 mm. The back beam has a thickness of 25 mm. The rib has a thickness of 20 mm. The material of the wing skin is carbon fiber composite, and the material properties are shown in Table 1. The area division of the wing skin ply is shown in Fig. 8. The wing skin is symmetrically ply. The number of plies at the wing root is the largest, and the number of plies is gradually reduced along the span, as shown in Table 2.

Fig. 7High-aspect-ratio composite wing model geometry

Fig. 8Area division of high-aspect-ratio composite wing skin ply

The structure of the wing model mainly includes skins, spars, and ribs, all of which were meshed by shell 181 shell elements. The number of generated nodes was 66504, and the number of grids was 23775. The finite element model of the wing is constrained by a fixed support of the wing root, as shown in Fig. 9.

Fig. 9Finite element model of the high-aspect-ratio composite wing

a) Boundary condition setting of finite element model

b) Internal structure of finite element model

Fig. 10 shows the flow field calculation domain and the local mesh division near the wing. The unstructured tetrahedral mesh was used for mesh division, and the number of grid cells was 2.8×10⁶. The wing root was set as a symmetrical boundary condition, the wing was set to the object boundary condition, and the rest ones were set as the pressure far away from the field boundary condition.

Fig. 10Calculation domain and mesh of flow field

a) Overall computing domain

b) Grid near the wall

3.1. Analysis of elastic deformation of high-aspect-ratio composite wing before optimization

The problem of nonlinear static aeroelastic deformation of the composite wing was solved under the condition that the Mach number was 0.7, the Reynolds number was 9.98×10⁶, and the angle of attack was 4°. Fig. 11 shows the deformation of the composite wing in the lift direction. Fig. 12 depicts the stress distribution of the wing beam and rib. Fig. 13 is the Tsai-Wu failure factor distribution of the wing composite skin. It can be seen from the figure that the maximum wing deformation is 0.55989 m. The maximum torsion angle of the wing is –3.215°. The maximum stress of the wing beam and rib is 237.21 MPa. The maximum Tsai-Wu failure factor of the wing skin is 0.84031. Although the structural strength condition is satisfied, the torsional deformation of the wing in the lift direction is too large, which reduces the aerodynamic characteristics of the high-aspect-ratio wing.

Fig. 11Deformation of the wing along the Z axis

Fig. 12Stress distribution of the wing beam and rib

3.2. Optimization result analysis

In order to find the optimum ply angle, the torsional deformation of the high-aspect-ratio wing was minimized. The ply angle of the wing skin was optimized. In the optimization process, the wing skin was asymmetrically ply. Thus, 16 independent variables were set. The wing deformation along the lift direction was taken as the optimization target. The structure strength was taken as the constraint. The screening method was taken as the optimization method. The total number of samples was defined as 100, and these 100 samples were analyzed. Finally, the three candidate points with the smallest deformation were obtained as shown in Table 3. Under the condition that the structural strength was satisfied, the deformation of the candidate point 1 was the smallest, so the candidate point 1 was selected as the final optimization result. The final results of the ply angle optimization of the wing skin are shown in Table 4.

Fig. 13Tsai-Wu failure factor distribution of the wing composite skin

Fig. 14 shows the wing deformation along the lift direction after optimization. Fig. 15 shows the stress distribution of the wing beam and rib after the optimization. Fig. 16 shows the Tsai-Wu failure factor distribution of the wing composite after the optimization. It can be seen from the figure that the maximum wing deformation is 0.34079 m. The maximum torsion angle of the wing is –2.576°. The maximum stress of the wing beam and rib is 144.59 MPa. The maximum Tsai-Wu failure factor of the wing skin is 0.44432. It can be seen that the optimized wing still meets the strength conditions.

Fig. 14Deformation of the wing along the Z axis

Fig. 15Stress distribution of the wing spar and rib

Fig. 16Tsai-Wu failure factor distribution of wing composite skin

Table 5 shows the comparison of optimization results. It can be seen from the table that the deformation in the lift direction and torsion angle of the wing are reduced: the maximum deformation is reduced by 39.1 %, and the maximum torsion angle is reduced by 19.8 %. The maximum stress of the wing beam and rib is reduced by 39.0 %, and the maximum Tsai-Wu failure factor of the wing skin is reduced by 47.1 %. This indicates that the torsional deformation of the wing is significantly reduced under the constraint of strength.