Static Aeroelastic Modeling and Rapid Analysis of Wings in Transonic Flow

A fast static aeroelastic analysis method, coupling with the modal method and Kriging surrogate model, is proposed in this paper. The deflection of the wing is described by the modal method, and the Kriging surrogate model is utilized to model the generalized forces under different deformations, angles of attack, and Mach numbers in order to replace the CFD solver. We analyzed the static aeroelasticity of HIRENASD wing in transonic flow field by coupling with the generalized force model by the static equilibrium equation. The results were compared with those of the experimental data and the references, and the comparison shows that the method is useful for the small deflections. After enough training cases are finished, the high-accuracy aerodynamic force coefficients and wing deflection will be obtained rapidly, which will only take several seconds. This method is more time saving than the CFD/CSD method, when it comes to a large quantity of the static aeroelastic analyses. Hence, it has good perspective for engineering applications during the aircraft design period.


Introduction
Most of the modern aircrafts have high-aspect-ratio swept wing [1] and use composite material [2], which makes the static aeroelasticity become more and more severe.For instance, the wing on the Boeing 787 Dreamliner will nominally deflect 10 feet at cruise [3].For the high-aspect-ratio swept wing, it will give rise to the wing deflection and torsion because of the static aeroelasticity problems, which will cause the angle of attack at the local section becomes smaller, and changes the distribution of surface pressure [3].For the transonic aircraft, the decrease of angle of attack at the local section will also influence the intensity of shock wave.Therefore, it is crucial to analyze the static aeroelasticity for modern aircrafts.
By CFD/CSD method, loosely coupling is usually utilized to increase the computational efficiency.Hence, we can analyze the structure and the flow field separately.When it comes to structural analysis, there are two methods: one is linear analysis method.The wing deflection can be depicted by the modal method, under the assumption of small deflection.The other is nonlinear method.Finite element analysis is usually applied.Mian et al. [4] used this method to analyze the nonlinear deflection of high-aspect-ratio wing.In addition, both the multibody method [5] and nonlinear aeroelastic scaling method [6] are also used.In this paper, we used the computational structural dynamics (CSD) to analyze the structural deflection.On the other hand, when it comes to aerodynamic analysis, there are also two methods, linear and nonlinear methods, which are similar to structural analysis methods.Traditional aerodynamic linear theory is usually applied in the linear analysis method.For certain conditions, linear theory works well.However, the conditions are limited to the subsonic and supersonic flows [7].When it refers to nonlinear aerodynamic calculation in transonic flow, computational fluid dynamics (CFD) is usually applied at present.When it comes to the static aeroelastic analysis, the steady CFD solution of the rigid wing is calculated firstly.Next, we get the generalized force and calculate the deflection of wing to obtain the new boundary by coupling with the static equilibrium equation.After that, we change the wall boundary by grid deformation method [8] and use CFD to calculate the aerodynamic force on the deformed wing.The results can be obtained by repeating the above steps in order, until the deflection converges or diverges.However, when it comes to the complex threedimensional configurations, this method is time-consuming and has high computational complexity.
More and more researchers now used the surrogate method in their studies [9][10][11][12], as the problems mentioned above exist in the fluid-structure coupling problems and optimum design.Surrogate methods can fit the nonlinear multiple-input/output function accurately [13,14] and has high computational efficiency.Due to these advantages, more and more researchers utilized these methods to model the nonlinear unsteady aerodynamics [15][16][17][18][19]. Lindhorst et al. [20] combined the parameter reduction via proper orthogonal decomposition and system identification methods to model nonlinear unsteady two-dimensional aerodynamics.And the model can accurately predict the static and transient response of the airfoil.Furthermore, they demonstrated the application on a three-dimensional case [21], the high-Reynolds-number aerostructural dynamics (HIRENASD), and the model can capture the influences of nonlinear aerodynamic effects on the forces.Moreover, the model can be used in both static and transient aeroelastic investigations at a fixed Mach number [22].Kou and Zhang [23] applied radial basis function neural network to model twodimensional nonlinear aerodynamics.And the approach can capture both linear and nonlinear characteristic.
This paper proposed an approach to model the nonlinear steady aerodynamics in transonic flow.The wing deflection is described by the modal method under linear structural assumption, while the aerodynamic force of different deformed wing is obtained by the surrogate model.And the deflection is calculated by the generalized forces coupling with the static equilibrium equation.The approach is used to perform the static aeroelastic analysis of HIRENASD wing.The results were compared with those obtained by CFD/CSD method, which shows the validity and the accuracy of the approach.The model has acceptable accuracy for engineering.This approach is more efficient than the CFD/CSD method in the acceptable accuracy, when it comes to plenty of analysis.

Introduction of Method
2.1.CFD/CSD Coupling Method.The flow governing equations used to solve the aerodynamics can be written as where Ω is the control volume, ∂Ω is the boundary of the control volume, n is the outer normal vector of the control volume boundary, V denotes the volume of the element, and S denotes the surface area of each surface.(When it comes to two dimension, V denotes the surface area, and S denotes the length.)Q is the vector of conservative variables, F Q is the vector of the inviscid fluxes, and G Q denotes viscous fluxes.More details of Q, F Q , and G Q are shown in [24].
The CFD solver based on the steady Reynolds-Averaged Navier-Stoke (RANS) [25] equations has the ability to simulate the flow with viscous effects.An unstructured RANS solver based on the finite volume is used in this paper.The Spalart-Allmaras (S-A) turbulence model [26] works well in describing the viscosity in the transonic flow, in which shock wave exists.Thus, this model is used in all cases in this study.
Under the assumption of linear small deflection, the modal method is utilized to describe the wing deflection.Since there is no need to take the structural inertial force into consideration, the deflection is computed by the static equilibrium equation: where ξ is the generalized displacement vector, K is the generalized stiffness matrix, and F is the generalized force vector.
To get rid of the dynamic pressure effect, the generalized forces are nondimensionalized by dynamic pressure, which are called as generalized force coefficients.
where f represents the generalized force coefficient vector and Q represents the dynamic pressure.
The new structural boundary can be depicted by the following equation: where x rigid is the coordinate matrix of the surface grid nodes for the rigid wing, ξ i is the vector of modal coordinates, φ i is the modal matrix of the surface mesh for the ith mode, and x new is the coordinate matrix of the grid nodes.
In aeroelastic analysis, the relaxation factor η is often used to enhance the computational stability, since large deformation may result in negative volume of the grid.And the new structural boundary would be obtained in the following formulation: where x j+1 new represents the coordinate matrix of the boundary grid nodes at the j + 1th iteration and x j old represents the coordinate matrix of the boundary grid nodes at the jth iteration.x j new represents the coordinate matrix of the grid nodes, which is obtained from the new wing deflection.It is calculated according to the generalized forces of the jth iteration j = 0, 1, 2, … and η ∈ 0, 1 0 .
The grid is deformed according to the new boundary by spring analogy method.In addition, loosely coupling is usually utilized to increase the computational efficiency.
The process of using CFD/CSD method has already been depicted in Introduction.

Kriging Model.
The Kriging surrogate model [27] is a kind of model aimed at minimizing variance and constructing an unbiased estimation of the spatial distribution data 2 International Journal of Aerospace Engineering via the statistical method of stochastic process.The functional expression [28] can be given as where P β, x is the regression model and β is the regression parameters.z x is the nonparametric random function, and its statistical properties are written as Mean value where x i and x j are the design sites, and R ij θ, x i , x j is the function with parameter θ and represents the spatial relativity among the design sites.The spatial relativity between 3 International Journal of Aerospace Engineering every two design sites is related to their spatial distance.Hence, it can be depicted by the following equation: where n is the number of design variables and d k is the distance between every two design sites.The concrete function is given as where x k i and x k j are the coordinate values of the ith and jth design sites in the kth direction and θ k is the constant parameter of the function in the kth direction.
Aimed to minimize σ 2 x * new , after the mathematical derivation, the predictor is computed as where , where P = p 1 , p 2 , … , p n T , which represents the vector of regression values for the design sites, and Y is the response array of the design sites.Since β * and γ * are related to the design sites instead of the predictor sites and the predictor sites are only related to f x * new and r x * new , the predicted response ŷ x * new will soon be obtained when x * new is given.We utilized the Gauss Function as R matrix: where θ k = 0 1.The first step is to choose some sets of Ma and AOA from a certain range by Latin hypercube sampling (LHS) [29].The nonlinear effect in transonic flow occurs when either Ma or AOA changes.For example, the aerodynamic force would change nonlinearly even if either AOA or Ma increases linearly due to the effect of the shock wave.Hence, the nonlinear effect of the Ma and AOA should be taken into consideration.The next step is to choose several sets of dynamic pressures at each set of Ma and AOA.Then, the static aeroelasticity will be analyzed at each condition to obtain the corresponding equilibrium position.The third step is to choose several sets of generalized displacements by LHS from a certain range close to each equilibrium position, in order to ensure that the numerical value of each mode generalized displacement is limited to a certain range.The obtained generalized displacements need to be able to describe the real wing deformations and enable the deformations to change in a certain range, since the model needs to predict the aerodynamic forces for the wing of different deformations.This sampling method can satisfy the requirements, so we gained the generalized displacements in this way.Finally, CFD method is used to calculate the corresponding generalized force coefficients and aerodynamic force coefficients at each training case.Hence, the sum of training cases can be calculated by the following function:

The Static Aeroelastic Analysis Method
where S total represents the sum of training cases, sampling Ma,AOA represents the sampling sets of Ma and AOA, sampling Q represents the sampling sets of Q, and sampling ξ represents the sampling sets of generalized displacements.AOA will also be input to the model.And the output is a set of the corresponding generalized force coefficients.The relation between the inputs and the output is given as

Modeling
where f is the vector of generalized force coefficients and ξ is the vector of generalized displacements.To get the generalized forces, we need to multiply the f by the dynamic pressure, for the coefficients are nondimensionalized by the dynamic pressure.The forces would be used to calculate the deformation of wing, which will put forward the procedure.This model is utilized to replace the CFD flow solver and will be called many times during the static aeroelastic analysis.We named this model as elastic generalized force (EGF) model for convenience.
In addition, we produced another model to predict the lift and drag coefficients of the deformed wing.The concrete function can be written as where CL represents the lift coefficient and CD is the drag coefficient.Different from the output of EGF, the outputs of this model are lift and drag coefficients, while the inputs are the same as those of EGF.In order to distinguish the models, we called this model as elastic aerodynamic force (EAF) model.
Remark 1.A model is required to calculate the generalized force coefficients of the rigid wing, since the corresponding generalized forces are the initial forces to calculate the wing deflection.The inputs of this model are the Mach number and angle of attack, since the forces are only related to flow condition parameters.And the output is the same as that of EGF.The function can be given as where f 0 is the vector of the generalized force coefficients for  Similar to EAF, we also produced another model to predict the aerodynamic force coefficient of the rigid wing.The inputs are same as those of RGF, while the outputs are lift and drag coefficients.Therefore, the function of this model can be depicted as where CL 0 is the lift coefficient of the rigid wing and CD 0 is the drag coefficient.This model will also be called only once.This model was named as rigid aerodynamic force (RAF) model in routine.
The flow charts of the analysis are shown in Figures 1  and 2.
Remark 2. Figures 1 and 2 are shown to illustrate the analysis procedure.The middle chart is the procedure of CFD/CSD method.At the beginning of the analysis, both the aerodynamic force of the rigid wing and the natural modes of the structure are needed and are not related to each other.Hence, the flow field and the structural analysis can be divided and computed separately.The 1st iteration of the flow field analysis consists of the red blocks ① and ②.Then, the flow field and structural analysis would be performed alternately until the wing deformation converges or diverges, as shown in the circle of the middle chart.Each iteration of flow field consists of the red blocks ③ and ②.However, the flow field analysis, shown by the red blocks, would cost a lot of time at each iteration.Therefore, the surrogate model is used to replace them.
The yellow and blue charts in both Figures 1 and 2 are the forms to obtain the training cases of the four models, which should be finished before the analysis.And the training cases of these four models, RGF (the yellow charts in Figure 1), EGF (the blue charts in Figure 1), RAF (the yellow charts in Figure 2), and EAF (the blue charts in Figure 2), would be obtained, and the models would be produced.This step would cost the most time of this method, since a number of training cases need to be calculated by the CFD method.
Then, the models (the green and purple charts), instead of CFD (the red blocks), would be used during the analysis.
The RGF and RAF (the green ones in Figures 1 and 2) would be only used once at the 1st iteration to obtain the corresponding forces, since they are produced for the rigid wing.In contrast, the EGF and EAF (the purple ones in Figures 1  and 2) would be used at each iteration except the 1st iteration.The outputs of the green and purple charts in Figure 1 would be used to calculate the wing deformation, while those in Figure 2 would be used to calculate the concerned lift and drag coefficients, nothing to do with the wing deformation calculation.After the model is produced, we can use it to analyze a static aeroelasticity analysis.The CFD solver will be replaced, so the mesh deformation used in CFD/CSD method will be also not needed for the proposed analysis method.Therefore, it could be finished in several seconds.
Furthermore, we would like to compare the characteristics between the different models.The outputs of RGF and EGF will be transferred to calculate the wing deflection, while RAF and EAF are produced to predict the aerodynamic force coefficients that we paid attention to.In addition, RAF and RGF are used to calculate the force coefficients of the rigid wing, while EAF and EGF are used to calculate those of the  7 International Journal of Aerospace Engineering deformed wing.Moreover, RGF and RAF will only be used once at the beginning of the progress, while EGF and EAF will be called many times during the process.
Last but not least, we would like to explain the limits of this method.This approach is a method based on the training cases to replace the CFD solver and shorten the computational time.Hence, the accuracy of this method depends on the number and accuracy of the training cases.On the other hand, this method is under the assumption of small deflection.When it comes to large deflection, the modal method is not suitable since it is based on the principle of linear superposition, so the proposed modeling method is not applicable.

Model Introduction.
The HIRENASD model is used to perform the static aeroelastic analysis.The model is downloaded from the AePW website [30].The model planform is shown in Figure 3 (copied from [30]) and the main parameter list is shown in Table 1.

International Journal of Aerospace Engineering
We chose the primary 10-order modes to describe the wing deflection.The 3rd and 7th mode are shear modes, which do not influence the static aeroelastic analysis.Hence, we used 8 modes in the modal method to describe the wing deflection.The natural mode shape of the first 4 modes we used is shown in Figure 4.The modal circular frequencies are shown in Table 2.And the span sections are illustrated in Table 3.

Mesh Information.
The information of the hybrid mesh is given in Table 4.
The details of the grid are shown in Figure 5.    [31] and the experimental data.The results are shown in Figure 6, and the deflection of the wing tip is shown in Figure 7.
To validate the accuracy of the model, several cases of the static aeroelasticity for the HIRENASD model were analyzed at the conditions of Ma = 0 8, AOA = 1 5 °, and dynamic pressure Q = 20,000 Pa; 40,055 4 Pa; 60,000 Pa; and 80,000 Pa.The leading edge and trailing edge deflections at the wing tip (99% semispan) have been measured, and the results will be shown below.
The code has been verified by many aeroelastic problems [32,33].Table 6 and Figure 8 illustrate that the computational results show a good agreement with those of the reference and experimental data.The results of grid A for HIRENASD, the coarsest grid in [31], were chosen as the reference results.
At the condition of Ma = 0 8 and AOA = 1 5 °, the deflection of wing tip is a little nonlinear along with the increase of Q, as illustrated in Figure 10.In addition, the results of the proposed modeling method match well with those obtained by the CFD/CSD method, which is also illustrated in Figure 11.The largest relative errors of the leading edge and trailing edge deflection are 1.96% and 0.95%, as illustrated in Table 7.And Table 8 shows that the EAF model works well on predicting the aerodynamic force coefficients of the deformed wings.The largest relative errors of lift and drag coefficients are +0.44% and +1.47%.Furthermore, the relative  10 International Journal of Aerospace Engineering errors of drag coefficients are obviously greater than those of lift coefficients, which reveals that it is more difficult to predict the drag coefficient than to predict the lift coefficient.Figures 9 and 12 illustrate that it converges at the 4th iteration by the model, while it converges at the 13th iteration by CFD/CSD method.The relaxation factor of CFD/CSD method was chosen as 0.3, since large deflection may give rise to the negative volume of the deformed mesh.However, after the model is produced, the mesh deformation used in CFD/CSD method is not needed anymore.The new positions of the wall boundary can be simply obtained by the superposition of the structural modes instead of moving all the computational grids.Therefore, the problem of the negative volume will not appear, and the relaxation factor of this method can be chosen as 1.0, which will accelerate the convergence significantly.
Furthermore, the modeling method improves the computational efficiency obviously, when it comes to a large quantity of the static aeroelastic analyses.To calculate a case of static aeroelastic at a certain condition, it takes about 19.6 hours to converge by CFD/CSD method via parallel computation with 6 cores.For the modeling method, under equal computational cost, it takes about 240 hours to finish training the cases, which costs such a long time as finishing the static aeroelastic analysis 12.2 times.After the model is produced, it takes only about 5 seconds to converge by the modeling method.Hence, for the static aeroelastic analyses at 12 or even less conditions, this method may not be as efficient compared with the CFD/CSD method.However, when it comes to doing the analyses at 13 or even more conditions, this method will become more and more time saving.
This approach has a good perspective of the engineering applications during the period of the aircraft design.For instance, using the CFD/CSD method would cost 1960 hours to analyze 100 cases in the flight envelope.However, it will take about 240.14 hours by the proposed modeling method, which can increase the computational efficiency for about 8.16 times.

Conclusions
This paper proposed a fast method to perform the static aeroelastic analysis of the wing in transonic flow.The developed model is aimed at replacing the CFD solver and can calculate the generalized force coefficients and aerodynamic force coefficients at different Mach numbers and different angles of attack.It is coupled with the static equilibrium equation to perform the static aeroelastic analysis.The model HIRENASD was analyzed by this method, and the results were compared with those obtained by the CFD/CSD method.The largest relative errors of the leading edge and trailing edge deflection are 1.96% and 0.95%, which shows that the model is accurate for the, respectively, small deflections.If there are enough training cases, the model can achieve the accurate wing static deflection and aerodynamic force coefficients rapidly.It can reduce the analysis time from 19.6 hours to 5 seconds at a certain condition.However, training cases will cost a lot of time.Hence, when it refers to a large quantity of the static aeroelastic analyses, this approach is much more efficient than the CFD/CSD method.For these advantages, the method has a good perspective for engineering applications and can save the computational resource and shorten the design period.

Figure 1 :Figure 2 :
Figure 1: The flow chart of analysis.

Figure 4 :
Figure 4: The shapes of the first four modes used in the modal method (the contour plot of the height (unit: m); x: streamwise, y: spanwise, z: height).
The comparison between the rigid and elastic wings z y x (b) The comparison between the rigid and elastic wing tips

4. 4 .
Computational Results.First, the static aeroelastic analysis of the HIRENASD model is performed by CFD/ CSD method.The solution is obtained at the conditions of Ma = 0 8, AOA = 1 5 °, Reynolds number Re = 7 million, and dynamic pressure Q = 40,055 4 Pa.And we compared the results with those of

Figure 12 :
Figure 12: The convergence history of lift coefficient and the 1st mode displacement by the model under Ma = 0 8, AOA = 1 5 °, and Q = 40055 4 Pa condition.

Table 1 :
The parameter list of the model.

Table 4 :
The grid information.

Table 5 :
The sampling data of Mach number and angle of attack.
the rigid wing.In routine, we named this model as rigid generalized force (RGF) model.

Table 6 :
The comparison of the deflection results of the reference and experimental data.

Table 7 :
The comparison of the displacement of the wing tip by the CFD/CSD method and modeling method.

Table 8 :
The comparison of the aerodynamic force coefficient by the CFD/CSD method and modeling method.