Numerical Simulation Research on The Transonic Aeroelasticity of a High-aspect-ratio Wing

2.1 CFD method

3D time-correlated conservative compressible RANS equations are used as governing equations. In general curvilinear coordinate system (ξ, η, ζ), the dimensionless form is

$\frac{\partial \boldsymbol{Q}}{\partial t}+\frac{\partial\left(\boldsymbol{F}-\boldsymbol{F}_{v}\right)}{\partial \xi}+\frac{\partial\left(\boldsymbol{G}-\boldsymbol{G}_{v}\right)}{\partial \eta}+\frac{\partial\left(\boldsymbol{H}-\boldsymbol{H}_{v}\right)}{\partial \zeta}=0$       (1)

Where, t is the time, Q is the conservation variable and F, G and H are inviscid flux terms. F_v, G_v and H_v represent the viscous flux terms.

The calculation is based on the finite volume method and the convection term utilizes Roe’s upwind flux-difference-splitting method for discretization. Also, a Venkat limiter is introduced to suppress the non-physical oscillation near shock waves and Van Leer’s Monotone Upstream-Centered Schemes for Conservation Laws (MUSCL) are used for the interpolation of the conservation variables of the flow field. In addition, viscous terms are discretized using a 2-order central difference scheme and an implicit LU-SGS (Lower-Upper Symmetry Gauss-Seide) method is utilized for time difference. Moreover, the one-equation SA (Spalart-Allmaras) model is chosen as the turbulence model, the object wall adopts an adiabatic and non-slip boundary condition, the far field is a pressure far field using non-reflecting boundary condition and the multi-grid technique is applied to accelerate the convergence of the calculation.

2.2 CSD method

For linear elastic problems, the matrix of flexibility influence coefficients can be used to solve the structural statics equation, which avoids solving complicated structural modal motion equations. However, structural elastic deformation can actually be obtained by several matrix operations, whose calculation time is shorter and processing is simpler and more suitable for studies of static aeroelastic problems. Generally, the elastic deformation of the wing is only considered in the process of static aeroelastic calculation of large aircraft and some other components such as fuselage, engine or hangings are assumed as rigid parts. The wing’s matrix of flexibility influence coefficients can be obtained by a structural finite element modeling calculation using MSC, NASTRAN or some other commercial software. Besides, it also can be directly achieved by test measurements. If aerodynamic loads, engine thrust and the mass force of three directions acting on the wing structure are known, the elastic deformation displacement of structural points can be obtained based on the structural equation whose form is

$\mathrm{u}_{\mathrm{s}}=\boldsymbol{C} \boldsymbol{F}_{s}$   (2)

Here, u_s is the deformation displacement vector of structure points, C is flexibility matrix of the structure points, and F_S is mass force of aerodynamic loads, engine thrust and mass force acting on the structure points.

2.3 Interpolation method between flow field and structure

CFD calculation of conventional static aeroelastic issues is based on the Euler coordinate system and CSD is described based on the Lagrange coordinate system. Therefore, the division of the CFD grid and construction of the structure model are independent, so data exchange needs to be introduced between the structure and the flow field. One method is to transform the displacement of structure points to that of grid points; a second method is to transform the aerodynamic loads of CFD grid points to an equivalent force acting on structure points. As for static aeroelastic calculation, its deformation of three-dimensional space is small, so it is more suitable to use a three-dimensional TPS method in consideration of the memory occupancy rate, calculating its efficiency and precision of interpolation [15].

In the TPS method, the structure displacement interpolation matrix H can be obtained through the coordinates of the CFD grids and structure points of the flexibility matrix. Define that u_s and u_a are the deformation shift vector of the structure points and deformation shift vector of CFD grid points respectively, then

$\boldsymbol{u}_{\mathrm{a}}=\boldsymbol{H} \boldsymbol{u}_{\mathrm{s}}$    (3)

According to the virtual work principle, the interpolation matrix between the aerodynamic loads and the displacement interpolation matrix between structural grids and CFD grid points are transposed, which means

$\boldsymbol{F}_{\mathrm{s}}=\boldsymbol{H}^{\mathrm{T}} \boldsymbol{F}_{\mathrm{a}}$     (4)

Therefore, the deformation shift vector u_a can be achieved through the equations (2) ~ (4) if aerodynamic loads of CFD surface grids F_a are known

$\boldsymbol{u}_{\mathrm{a}}=\boldsymbol{H} \boldsymbol{C} \boldsymbol{F}_{\mathrm{s}}=\boldsymbol{H} \boldsymbol{C} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{F}_{\mathrm{a}}$   (5)

2.4 Dynamic mesh generation method

It is apparent that previous CFD grids are no longer applicable because of the deformation of the object surface, so new grids must be generated. There are two methods to generate new grids. One is through regenerating the grid directly, whose computational efficiency is quite low and is likely to change the mesh topology and the number of grid cells. The other method is through adding an amount of correction to the previous mesh. Due to the static aeroelastic calculation being based on an assumption that the structural deformation is small and the deformation of the structure in the local area is usually not severe, the second method is feasible and applied extensively, and its computational efficiency is relatively higher. This paper integrates the RBF and TFI methods to generate dynamic structure grids.

The RBF method is a volume interpolation based on spline function and also can be regarded as the expansion in three dimensions of the surface spline function interpolation method [16]. Its interpolation formula is

$f(\boldsymbol{r})=\sum_{i=1}^{n} a_{i} \varphi\left(\left\|\boldsymbol{r}-\boldsymbol{r}_{i}\right\|\right)+\psi(\boldsymbol{r})$      (6)

Here, $\boldsymbol{r}_{i}=\left(x_{i}, y_{i}, z_{i}\right)$ are points whose displacements are known and their number is n. In addition, φ is the primary function of space distance $\left\|r-r_{i}\right\|$ . In this paper, $\varphi\left(\left\|\boldsymbol{r}-\boldsymbol{r}_{i}\right\|\right)=\left\|\boldsymbol{r}-\boldsymbol{r}_{i}\right\|^{3}$ , $\psi=b_{0}+b_{1} x+b_{2} y+b_{3} z$ . The coefficients of the interpolation equation can be obtained through displacement d_i of r_iand equilibrium condition

$f\left(r_{i}\right)=d_{i}$

$\sum_{i=1}^{n} a_{i}=\sum_{i=1}^{n} a_{i} x_{i}=\sum_{i=1}^{n} a_{i} y_{i}=\sum_{i=1}^{n} a_{i} z_{i}$      (7)

Their matrix forms are

$\left[\begin{array}{cc}\boldsymbol{\Phi} & \boldsymbol{P} \\ \boldsymbol{P}^{\mathrm{T}} & 0\end{array}\right]\left[\begin{array}{l}\boldsymbol{a} \\ \boldsymbol{b}\end{array}\right]=\left[\begin{array}{l}\boldsymbol{d} \\ 0\end{array}\right]$

$\boldsymbol{\Phi}=\left[\begin{array}{cccc}\varphi_{11} & \varphi_{12} & \mathrm{L} & \varphi_{1 n} \\ \varphi_{21} & \varphi_{22} & \mathrm{L} & \varphi_{2 n} \\ \mathrm{M} & \mathrm{M} & \mathrm{O} & \mathrm{M} \\ \varphi_{n 1} & \varphi_{n 2} & \mathrm{L} & \varphi_{n n}\end{array}\right]$      (8)

Where,

$P=\left[\begin{array}{llll}1 & x_{1} & y_{1} & z_{1} \\ 1 & x_{2} & y_{2} & z_{2} \\ \mathrm{M} & \mathrm{M} & \mathrm{M} & \mathrm{M} \\ 1 & x_{n} & y_{n} & z_{n}\end{array}\right]$ ,

$\boldsymbol{a}=\left[\begin{array}{c}a_{1} \\ a_{2} \\ \mathrm{M} \\ a_{n}\end{array}\right]$ ,$\boldsymbol{b}=\left[\begin{array}{l}b_{0} \\ b_{1} \\ b_{2} \\ b_{3}\end{array}\right]$ and $\boldsymbol{d}=\left[\begin{array}{c}d_{1} \\ d_{2} \\ \mathrm{M} \\ d_{n}\end{array}\right]$,

$\varphi_{i j}=\varphi\left(\left\|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\right\|\right)$ .

After obtaining the interpolation function coefficients, displacement d_inter of any interpolation point r_intercan be achieved: $\boldsymbol{d}_{\mathrm{int} e r}=f\left(\boldsymbol{r}_{\mathrm{int} e r}\right)$ .

Theoretically, the displacement of the whole space grid can be got directly by interpolation based on the primary function method, if the displacement of the surface grids is known and the displacement of the far field grid points is set to zero. However, the dimensions of matrixΦ in equation (8) will be too many and the calculation efficiency of grid deformation is quite low when the geometric configuration is complex and the number of surface grids is large. In order to improve the computational efficiency, this paper utilizes the greedy algorithm and only selects that part of the object surfaces as well as far-filed grids as points of the radial basis function whose coordinates are known in the process of interpolation using RBF, which can improve computational efficiency.

Based on the above methods, the displacement of the internal grid points of all grid blocks’ lines can be achieved through interpolation using the RBF method for dynamic grids, the displacement of internal surface grids of all grid blocks can be obtained by the TFI method [17] based on block surfaces’ lines and internal grid points in all blocks can be obtained by the TFI method based on the blocks’surfaces.

In order to improve the aerodynamic efficiency of large aircraft, their cruise Ma is usually in the transonic range. So the viscous effect cannot be ignored due to shock and boundary layer interaction, even if the angle of attack is very small, and an N-S equation should be used. For not losing generality, the static aeroelastic effects of a wing-body model when Ma=0.78, were calculated as an example based on the methods provided in this paper.

4.1 Geometrical deformation characteristics

Figure 4 illustrates changes of Δz and Δε (torsion angle induced by static aeroelasticity) along span-wise direction when Ma=0.78, q=35kPa and α=-4°~8°. It can be seen from

figure 4 that the static aeroelasticity results in upward bending deflection deformation of the wing and negative elastic torsion angle in the downwind wing sections when the angle of attack is positive. In addition, it can be found that deformation of the inner wing is smaller, and deformation becomes larger away from the location where the wing is folded and peaks at the wing tip, which is consistent with the distribution characteristics of wing stiffness. As for swept wings, these kinds of deformation characteristics will reduce the local angle of attack of wing along the span-wise direction, affect the load distribution, and change their aerodynamic characteristics, which are so-called static aeroelastic effects. Taking the calculation result in this paper when Ma=0.78, q=35kPa, α=2° as an example, bending deformation of the wing tip is nearly 2% of the wing span and the elastic torsion angle and the change of α of wing tip sections can be about -2° and -4° respectively.

4a.jpg

(a) Wing bending deformation

4b.jpg

(b) Wing elastic torsion deformation

Figure4. Wing geometrical deformation characteristics

4.2 Surface pressure distribution characteristics

Figures 5 and 6 demonstrate the surface pressure distribution contours of rigid and elastic models as well as pressure coefficients of different span-wise sections separately when Ma=0.78, q=35kPa and α=2°. As can be seen from figure 6, static aeroelasticity has a small impact on the surface pressure of the wing root, and decreases the negative pressure range of the upper outboard wing obviously, as well as affecting the pressure distribution characteristics of the wing tip to a large extent; this is consistent with wing geometric deformation characteristics. Figure 6 shows that the differences of pressure distribution among different wing sections are not very large but the bias becomes increasingly larger for the wing sections closer to the wing tip. This is mainly because the aerodynamic loads of a swept wing result in upward bending deflection deformation of the wing, induce decrease of the downwind local angle of attack due to the negative elastic torsion angle, and reduce the suck peak of wing leading edge when the angle of attack of the wing is small and the upper wing does not appear to have a large area of flow separation. Taking the calculation result in this paper when Ma=0.78, q=35kPa, α=2° as an example, C_L of an elastic wing is reduced by more than 15% compared with that of a rigid wing. Therefore, the impact of static aeroelasticity

on aerodynamic loads cannot be ignored for a high-aspect-ratio swept wing in a transonic speed range and must be underlined in the aircraft design stage.

5a.jpg

(a) Rigid model                    

5b.jpg

 (b) Elastic model

Figure 5. Effects of static aeroelasticity on wing surface pressure

6a.jpg

(a) y/b=20%

6b.jpg

(b) y/b=60%

6c.jpg

(c) y/b=95%

Figure 6. Comparison of pressure coefficients at different wing sections

4.3 Effect on aerodynamic characteristics and its correction

The dynamic pressure has a significant impact on the static aeroelastic effect of aircrafts. Therefore, the discipline as to how dynamic pressure affects the static aeroelastic correction factor of aerodynamic coefficients or derivatives needs to be studied.

Define  

K=C_{_F}/C_{_R}      (9)

Where, C_{_F} and C_{_R} represent aerodynamic coefficients or derivatives of an elastic wing and a rigid wing respectively. K is the correction factor of static aeroelasticity.

Define

K₁=C_Lα-F/C_Lα-R 

K₂= $C_{m C_{L_{-}} F}$ / $C_{m C_{L_{-}} R}$ 

K₃=K_α=2-F/ K_α=2-R         (10)

Here, C_Lα, $C_{m C_{L}}$ and K_α=2 represent the slope of lift coefficient , longitudinal static stability margin and the lift-to-drag ratio of α=2^o.

Figure 7 shows typical static aeroelastic correction factors of an elastic wing under different dynamic pressures when Ma=0.78. It can be seen that dynamic pressure has an impact on lift and drag characteristics as well as the static longitudinal stability of the elastic wing to different extents. Specifically, the slope of the lift coefficient K₁ has basically a linear downward trend, whose value is approximately 0.94 when q=35kPa at cruise height. However, this kind of linear discipline will change after the dynamic pressure reaches one particular large value. As can be seen from figure 7(a), the curve of K₁~q begins to upturn at q=65kPa, which probably results from geometrical nonlinearity of the wing elastic deformation when the dynamic pressure is relatively high. As for K₂, the longitudinal static stability margin, it generally has a nonlinear downward trend with an increase of dynamic pressure. However, it is noticed from figure 7(b) that the characteristics of K₂~q began to reverse when q>65kPa. Turning to K₃, the lift-to-drag ratio of α=2^o, static aeroelasticity will increase the lift-to-drag ratio as a whole, although it does decrease the lift due to the fact that it also reduces the drag induced by the lift. Figure 7 (c) demonstrates that linear characteristics can be found in K₃ when dynamic pressure is low, but the nonlinearity is quite obvious after q>45kPa.

7a.jpg

(a) K₁~q

7b.jpg

(b) K₂~q

7c.jpg

(c) K₃~q

Figure7.  Effects of static aeroelasticity on wing aerodynamic characteristics

In view of the above analysis, the probable nonlinearity of the static aeroelastic correction factor in high dynamic pressure should be considered and dynamic pressure range cannot be too small in calculation or experiment, in order to correct the static aeroelasticity of the high aspect ratio wing in transonic speed. Moreover, even if the dynamic pressure is low, its interval should not be too sparse in calculation, since some static aeroelastic correction factors of aerodynamic coefficients or derivatives may experience a nonlinear variation versus dynamic pressure. Therefore, considering both the reliability and economy of static aeroelasticity correction in a moderate or small angle of attack for large aircraft, a relatively reasonable way is by using the highly precise CFD/CSD method; however, the calculation method should be verified and corrected by a static aeroelastic wind tunnel test.