# Static Aeroelastic Modeling and Rapid Analysis of Wings in Transonic Flow

Authors: Jingyuan Yang; Yilang Liu; Weiwei Zhang
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

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.


---

## Body

## 1. 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 three-dimensional configurations, this method is time-consuming and has high computational complexity.More and more researchers now used the surrogate method in their studies [9–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–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 two-dimensional 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.


## 2. Introduction of Method

### 2.1. CFD/CSD Coupling Method

The flow governing equations used to solve the aerodynamics can be written as(1)∂∂t∭ΩQdV+∬∂ΩFQ⋅ndS=∬∂ΩGQ⋅ndS,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, FQ is the vector of the inviscid fluxes, and GQ denotes viscous fluxes. More details of Q, FQ, and GQ 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:(2)Kξ=F,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.(3)f=FQ,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:(4)xnew=xrigid+∑i=1Nξiφi,where xrigid 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 xnew 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:
(5)xnewj+1=η⋅xnewj+1−η⋅xoldj,where xnewj+1 represents the coordinate matrix of the boundary grid nodes at the j + 1th iteration and xoldj represents the coordinate matrix of the boundary grid nodes at the jth iteration. xnewj 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.


### 2.2. 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 via the statistical method of stochastic process. The functional expression [28] can be given as
(6)yx=Pβ,x+zx=fTxβ+zx,where Pβ,x is the regression model and β is the regression parameters. zx is the nonparametric random function, and its statistical properties are written asMean value(7)Ezx=0,Variance(8)Varzx=σ2,Covariance(9)Covzxi,zxj=σ2Rijθ,xi,xj,where xi and xj are the design sites, and Rijθ,xi,xj is the function with parameter θ and represents the spatial relativity among the design sites. The spatial relativity between every two design sites is related to their spatial distance. Hence, it can be depicted by the following equation:
(10)Rijθ,xi,xj=∏k=1nRkθk,dk,where n is the number of design variables and dk is the distance between every two design sites. The concrete function is given as
(11)dk=xik−xjk,i,j=1,2,…,n,where xik and xjk 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σ2xnew∗, after the mathematical derivation, the predictor is computed as
(12)y^xnew∗=fxnew∗Tβ∗+rxnew∗Tγ∗,where β∗=PTR−1P−1PTR−1Y and Rγ∗=Y−Pβ∗, where P=p1,p2,…,pnT, 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 fxnew∗ and rxnew∗, the predicted response y^xnew∗ will soon be obtained when xnew∗ is given.We utilized the Gauss Function asR matrix:
(13)Rkθk,dk=exp−θkdk2,where θk=0.1.


## 2.1. CFD/CSD Coupling Method

The flow governing equations used to solve the aerodynamics can be written as(1)∂∂t∭ΩQdV+∬∂ΩFQ⋅ndS=∬∂ΩGQ⋅ndS,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, FQ is the vector of the inviscid fluxes, and GQ denotes viscous fluxes. More details of Q, FQ, and GQ 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:(2)Kξ=F,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.(3)f=FQ,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:(4)xnew=xrigid+∑i=1Nξiφi,where xrigid 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 xnew 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:
(5)xnewj+1=η⋅xnewj+1−η⋅xoldj,where xnewj+1 represents the coordinate matrix of the boundary grid nodes at the j + 1th iteration and xoldj represents the coordinate matrix of the boundary grid nodes at the jth iteration. xnewj 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.


## 2.2. 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 via the statistical method of stochastic process. The functional expression [28] can be given as
(6)yx=Pβ,x+zx=fTxβ+zx,where Pβ,x is the regression model and β is the regression parameters. zx is the nonparametric random function, and its statistical properties are written asMean value(7)Ezx=0,Variance(8)Varzx=σ2,Covariance(9)Covzxi,zxj=σ2Rijθ,xi,xj,where xi and xj are the design sites, and Rijθ,xi,xj is the function with parameter θ and represents the spatial relativity among the design sites. The spatial relativity between every two design sites is related to their spatial distance. Hence, it can be depicted by the following equation:
(10)Rijθ,xi,xj=∏k=1nRkθk,dk,where n is the number of design variables and dk is the distance between every two design sites. The concrete function is given as
(11)dk=xik−xjk,i,j=1,2,…,n,where xik and xjk 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σ2xnew∗, after the mathematical derivation, the predictor is computed as
(12)y^xnew∗=fxnew∗Tβ∗+rxnew∗Tγ∗,where β∗=PTR−1P−1PTR−1Y and Rγ∗=Y−Pβ∗, where P=p1,p2,…,pnT, 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 fxnew∗ and rxnew∗, the predicted response y^xnew∗ will soon be obtained when xnew∗ is given.We utilized the Gauss Function asR matrix:
(13)Rkθk,dk=exp−θkdk2,where θk=0.1.


## 3. The Static Aeroelastic Analysis Method

### 3.1. Sampling Method

As the model is expected to calculate the generalized force coefficients at the different conditions of Mach numbers (Ma), angles of attack (AOA), and dynamic pressures (Q), model samples need to be chosen at three steps.The first step is to choose some sets ofMa 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:
(14)Stotal=samplingMa,AOA×samplingQ×samplingξ,where Stotal represents the sum of training cases, samplingMa,AOA represents the sampling sets of Ma and AOA, samplingQ represents the sampling sets of Q, and samplingξ represents the sampling sets of generalized displacements.


### 3.2. Modeling Method

A model is required to calculate the generalized forces at differentMa, AOA, and deflections. Hence, besides a set of generalized displacements, Ma and 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
(15)f=fgfMa,AOA,ξ,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(16)CL,CDT=fdfMa,AOA,ξ,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(17)f0=fgf_0Ma,AOA,where f0 is the vector of the generalized force coefficients for the rigid wing. In routine, we named this model as rigid generalized force (RGF) model.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(18)CL0,CD0T=fdf_0Ma,AOA,where CL0 is the lift coefficient of the rigid wing and CD0 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 Figures1 and 2.Remark 2.
Figures1 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.Figure 1
The flow chart of analysis.Figure 2
The flow chart of predicting the aerodynamic force coefficients during the analysis.The yellow and blue charts in both Figures1 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 Figures1 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 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.


## 3.1. Sampling Method

As the model is expected to calculate the generalized force coefficients at the different conditions of Mach numbers (Ma), angles of attack (AOA), and dynamic pressures (Q), model samples need to be chosen at three steps.The first step is to choose some sets ofMa 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:
(14)Stotal=samplingMa,AOA×samplingQ×samplingξ,where Stotal represents the sum of training cases, samplingMa,AOA represents the sampling sets of Ma and AOA, samplingQ represents the sampling sets of Q, and samplingξ represents the sampling sets of generalized displacements.


## 3.2. Modeling Method

A model is required to calculate the generalized forces at differentMa, AOA, and deflections. Hence, besides a set of generalized displacements, Ma and 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
(15)f=fgfMa,AOA,ξ,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(16)CL,CDT=fdfMa,AOA,ξ,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(17)f0=fgf_0Ma,AOA,where f0 is the vector of the generalized force coefficients for the rigid wing. In routine, we named this model as rigid generalized force (RGF) model.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(18)CL0,CD0T=fdf_0Ma,AOA,where CL0 is the lift coefficient of the rigid wing and CD0 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 Figures1 and 2.Remark 2.
Figures1 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.Figure 1
The flow chart of analysis.Figure 2
The flow chart of predicting the aerodynamic force coefficients during the analysis.The yellow and blue charts in both Figures1 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 Figures1 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 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.


## 4. Computation and Results

### 4.1. 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.Figure 3
HIRENASD wing model planform.Table 1
The parameter list of the model.
Model parameters
Data
Reference area (m2)
0.3926
Mean aerodynamic chord (m)
0.3445
Semispan (m)
1.2857We 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 Figure4. The modal circular frequencies are shown in Table 2. And the span sections are illustrated in Table 3.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).
(a)
The 1st mode
(b)
The 2nd mode
(c)
The 3rd mode
(d)
The 4th modeTable 2
The modal circular frequencies (unit: rad/s).
ω1
ω2
ω3
ω4
ω5
ω6
ω7
ω8
166.76
540.47
1189.48
1714.44
2021.77
2830.61
3120.73
3910.70Table 3
The span sections (y/b).
Sec1
Sec2
Sec3
Sec4
Sec5
Sec6
Sec7
0.145
0.323
0.456
0.589
0.655
0.804
0.953


### 4.2. Mesh Information

The information of the hybrid mesh is given in Table4.Table 4
The grid information.
Grid information parameters
Data
Viscous wall spacing (m)
1.0 × 10−5
Growth rate of the boundary layer
1.3
Prism layer cells
20
Surface mesh size (elements)
47,176
Surface mesh size (nodes)
23,679
Grid size (elements)
1,484,903
Grid size (nodes)
577,111The details of the grid are shown in Figure5.Figure 5
The details of the grid (x: streamwise, y: spanwise, z: height).
(a)
The grid of the wing surface
(b)
The grid of the wing root
(c)
The grid of the wing tip
(d)
The plane of the volume grid aty/b = 0.99


### 4.3. Sampling Method

The 1st step: the ranges of theMa and AOA are 0.72 ~ 0.85 and 1° ~ 4°, respectively, and 6 sets of cases are chosen by LHS. The details are given in Table 5. The 2nd step: at each set of Ma and AOA, the static aeroelastic analyses are performed at 3 sets of Q, 20,000 Pa; 50,000 Pa; and 80,000 Pa, by the CFD/CSD method. The 3rd step: 10 sets of generalized displacements are sampled by LHS at each flow condition. All the mode generalized displacements are sampled in the range from ±10% at their equilibrium positions. Then, CFD method is used to calculate the corresponding generalized and aerodynamic force coefficients for each training case. In summary, there are 180 (6 × 3 × 10) training cases for the aerodynamic model.Remark 3.
For the initial model, we calculated the generalized force coefficients at each set of Mach number and angle of attack, which means that there are 6 training cases for the initial model.Table 5
The sampling data of Mach number and angle of attack.
Mach number
Angle of attack (°)
0.7699
2.047
0.7293
1.638
0.7876
3.024
0.8158
1.208
0.8385
3.553
0.7469
2.739


### 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 ofMa=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 [31] and the experimental data. The results are shown in Figure 6, and the deflection of the wing tip is shown in Figure 7.Figure 6
The deflection results of static aeroelasticity (x: streamwise, y: spanwise, z: height and blue: the rigid wing, white: the elastic wing).
(a)
The comparison between the rigid and elastic wings
(b)
The comparison between the rigid and elastic wing tipsFigure 7
The deflection of the wing tip.To validate the accuracy of the model, several cases of the static aeroelasticity for the HIRENASD model were analyzed at the conditions ofMa=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.Table 6
The comparison of the deflection results of the reference and experimental data.
Reference result
Computational result
Experimental data
LE_deflection (m)
TE_deflection (m)
LE_deflection (m)
TE_deflection (m)
Wing tip deflection (m)
0.01289
0.01374
0.01372
0.01467
0.01250Figure 8
The rigid steady and static aeroelastic surface pressure comparison at sections 1, 2, 5, and 7.At the condition ofMa=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 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.Table 7
The comparison of the displacement of the wing tip by the CFD/CSD method and modeling method.
Dynamic pressure (Pa)
LE_CFD (m)
LE_model (m)
LE_error
TE_CFD (m)
TE_model (m)
TE_error
20,000
0.006859
0.006994
+1.96%
0.007617
0.007609
−0.10%
40,055.4
0.01313
0.01335
+1.70%
0.01468
0.01471
+0.20%
50,000
0.01635
0.01655
+1.20%
0.01792
0.01809
+0.95%
60,000
0.01932
0.01949
+0.85%
0.02112
0.02114
+0.10%
80,000
0.02480
0.02497
+0.72%
0.02727
0.02725
−0.06%Table 8
The comparison of the aerodynamic force coefficient by the CFD/CSD method and modeling method.
Dynamic pressure (Pa)
CL (CFD)
CL (model)
CL (error)
CD (CFD)
CD (model)
CD (error)
20,000
0.387007
0.38864
+0.42%
0.023206
0.023548
+1.47%
40,055.4
0.378365
0.38004
+0.44%
0.022848
0.023065
+0.95%
50,000
0.374299
0.37583
+0.41%
0.022695
0.022903
+0.92%
60,000
0.370305
0.37172
+0.38%
0.022549
0.022765
+0.96%
80,000
0.362644
0.36397
+0.37%
0.022267
0.022462
+0.88%Figures9 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.Figure 9
The convergence history of lift coefficient and the 1st mode displacements by CFD/CSD method.Figure 10
The displacement of the wing tip (dashed linear line is showed as reference).Figure 11
The wing deflection along the span (dynamic pressure: 40055.5 Pa).Figure 12
The convergence history of lift coefficient and the 1st mode displacement by the model underMa=0.8, AOA=1.5°, and Q=40055.4 Pa condition.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.


## 4.1. 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.Figure 3
HIRENASD wing model planform.Table 1
The parameter list of the model.
Model parameters
Data
Reference area (m2)
0.3926
Mean aerodynamic chord (m)
0.3445
Semispan (m)
1.2857We 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 Figure4. The modal circular frequencies are shown in Table 2. And the span sections are illustrated in Table 3.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).
(a)
The 1st mode
(b)
The 2nd mode
(c)
The 3rd mode
(d)
The 4th modeTable 2
The modal circular frequencies (unit: rad/s).
ω1
ω2
ω3
ω4
ω5
ω6
ω7
ω8
166.76
540.47
1189.48
1714.44
2021.77
2830.61
3120.73
3910.70Table 3
The span sections (y/b).
Sec1
Sec2
Sec3
Sec4
Sec5
Sec6
Sec7
0.145
0.323
0.456
0.589
0.655
0.804
0.953


## 4.2. Mesh Information

The information of the hybrid mesh is given in Table4.Table 4
The grid information.
Grid information parameters
Data
Viscous wall spacing (m)
1.0 × 10−5
Growth rate of the boundary layer
1.3
Prism layer cells
20
Surface mesh size (elements)
47,176
Surface mesh size (nodes)
23,679
Grid size (elements)
1,484,903
Grid size (nodes)
577,111The details of the grid are shown in Figure5.Figure 5
The details of the grid (x: streamwise, y: spanwise, z: height).
(a)
The grid of the wing surface
(b)
The grid of the wing root
(c)
The grid of the wing tip
(d)
The plane of the volume grid aty/b = 0.99


## 4.3. Sampling Method

The 1st step: the ranges of theMa and AOA are 0.72 ~ 0.85 and 1° ~ 4°, respectively, and 6 sets of cases are chosen by LHS. The details are given in Table 5. The 2nd step: at each set of Ma and AOA, the static aeroelastic analyses are performed at 3 sets of Q, 20,000 Pa; 50,000 Pa; and 80,000 Pa, by the CFD/CSD method. The 3rd step: 10 sets of generalized displacements are sampled by LHS at each flow condition. All the mode generalized displacements are sampled in the range from ±10% at their equilibrium positions. Then, CFD method is used to calculate the corresponding generalized and aerodynamic force coefficients for each training case. In summary, there are 180 (6 × 3 × 10) training cases for the aerodynamic model.Remark 3.
For the initial model, we calculated the generalized force coefficients at each set of Mach number and angle of attack, which means that there are 6 training cases for the initial model.Table 5
The sampling data of Mach number and angle of attack.
Mach number
Angle of attack (°)
0.7699
2.047
0.7293
1.638
0.7876
3.024
0.8158
1.208
0.8385
3.553
0.7469
2.739


## 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 ofMa=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 [31] and the experimental data. The results are shown in Figure 6, and the deflection of the wing tip is shown in Figure 7.Figure 6
The deflection results of static aeroelasticity (x: streamwise, y: spanwise, z: height and blue: the rigid wing, white: the elastic wing).
(a)
The comparison between the rigid and elastic wings
(b)
The comparison between the rigid and elastic wing tipsFigure 7
The deflection of the wing tip.To validate the accuracy of the model, several cases of the static aeroelasticity for the HIRENASD model were analyzed at the conditions ofMa=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.Table 6
The comparison of the deflection results of the reference and experimental data.
Reference result
Computational result
Experimental data
LE_deflection (m)
TE_deflection (m)
LE_deflection (m)
TE_deflection (m)
Wing tip deflection (m)
0.01289
0.01374
0.01372
0.01467
0.01250Figure 8
The rigid steady and static aeroelastic surface pressure comparison at sections 1, 2, 5, and 7.At the condition ofMa=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 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.Table 7
The comparison of the displacement of the wing tip by the CFD/CSD method and modeling method.
Dynamic pressure (Pa)
LE_CFD (m)
LE_model (m)
LE_error
TE_CFD (m)
TE_model (m)
TE_error
20,000
0.006859
0.006994
+1.96%
0.007617
0.007609
−0.10%
40,055.4
0.01313
0.01335
+1.70%
0.01468
0.01471
+0.20%
50,000
0.01635
0.01655
+1.20%
0.01792
0.01809
+0.95%
60,000
0.01932
0.01949
+0.85%
0.02112
0.02114
+0.10%
80,000
0.02480
0.02497
+0.72%
0.02727
0.02725
−0.06%Table 8
The comparison of the aerodynamic force coefficient by the CFD/CSD method and modeling method.
Dynamic pressure (Pa)
CL (CFD)
CL (model)
CL (error)
CD (CFD)
CD (model)
CD (error)
20,000
0.387007
0.38864
+0.42%
0.023206
0.023548
+1.47%
40,055.4
0.378365
0.38004
+0.44%
0.022848
0.023065
+0.95%
50,000
0.374299
0.37583
+0.41%
0.022695
0.022903
+0.92%
60,000
0.370305
0.37172
+0.38%
0.022549
0.022765
+0.96%
80,000
0.362644
0.36397
+0.37%
0.022267
0.022462
+0.88%Figures9 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.Figure 9
The convergence history of lift coefficient and the 1st mode displacements by CFD/CSD method.Figure 10
The displacement of the wing tip (dashed linear line is showed as reference).Figure 11
The wing deflection along the span (dynamic pressure: 40055.5 Pa).Figure 12
The convergence history of lift coefficient and the 1st mode displacement by the model underMa=0.8, AOA=1.5°, and Q=40055.4 Pa condition.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.


## 5. 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.



---

*Source: 5421027-2018-05-13.xml*