Time-Domain Aeroelasticity Analysis by a Tightly Coupled Fluid-Structure Interaction Methodology

: A tightly coupled ﬂuid-structure interaction (FSI) methodology is developed for aeroelasticity analysis in the time domain. The preconditioned Navier–Stokes equations for all Mach numbers are employed and the structural equations are tightly coupled with the ﬂuid equations by discretizing their time derivative term in the same pseudo time-stepping method. A modiﬁed mesh deformation method based on reduced control points radial basis functions (RBF) is utilized, and a RBF based mapping algorithm is introduced for data exchange on the interaction interface. To evaluate the methodology, the ﬂutter boundary and the limit cycle oscillation of Isogai wing and the ﬂutter boundary of AGARD 445.6 wing are analyzed and validated. previous studies. The methodology predicts the ﬂutter frequency ratios well even in the supersonic region as well as in subsonic and transonic regimes, while simulation results computed by other three methods still have larger differences with the experiment. The simulation results show that the tightly coupled FSI methodology with preconditioned N-S equations and CSD models can predict the aeroelasticity in the time domain accurately, and can also reveal the nonlinear characteristics of the aeroelastic system. Contributions: Conceptualization, manuscript.


Introduction
Aeroelasticity phenomenon, which arises from the interaction of the aerodynamic, elastic and inertial forces, may have detrimental effect on the reliability, cost and safety of aircraft. However, the complicated interaction mechanism and multidisciplinary requirement of aeroelasticity problem make it difficult to be analyzed with efficiency and accuracy.
With the development of high performance computers, coupled fluid-structure interaction (FSI) methodology is made available for aeroelasticity precondition and efforts have been made to develop the precondition capability of these methods. Generally, coupled FSI methodology can be divided into three types: the fluid and structural models are fully coupled, loosely coupled or tightly coupled. In the fully coupled method, the fluid and structural equations are integrated in both Eulerian and Lagrangian systems and solved in time simultaneously. The fully coupled method has a requirement in limitations on grid size. Moreover, this leads to the matrices being orders of magnitudes stiffer for structure system than fluid system, which makes it virtually impossible to solve the equations for large-scale problems [1]. For the loosely and tightly coupled methods, they use the existing solvers for computational fluid dynamics (CFD) and computational structure dynamics (CSD) [2,3]. Because the aerodynamic and structure domains are discretized in a different manner, data communications at the fluid-structure interface are required. The difference between the loosely and tightly coupled methods in essence lies in the data communication times. The loosely coupled strategy transfers the aerodynamic forces and structural displacements between fluid and structure systems only once in a physical time step, so it is usually only first-order time accuracy. For tightly coupled methods, data communications are performed several times within each physical time step and second-order time accuracy can be obtained. The linear or nonlinear types of models can be coupled together in FSI methodology. Current industrial practice relies heavily on linear methods, which can lead to conservative design and flight envelope restrictions [4]. However, FSI methodology with linear aerodynamic model cannot reveal the nonlinear aeroelastic effects, such as the transonic flutter dip and the limit cycle oscillations.
During the analysis of aeroelasticity problem, flexible structure deforms under the influence of air-loads, so it's necessary to update the fluid volume mesh according to the structural displacement. In the past decades, several mesh deformation methods have been investigated and successfully applied to aeroelastic analysis. Recently, extensive attention has been paid to the mesh deformation method based on radial basis functions (RBF), which can provide high-quality mesh deformation and perform reliable capability for arbitrary topology. The efficiency of RBF method is related to the number of control points on the wall surfaces of fluid mesh, and a reduced control points method proposed by Rendall [5] greatly improves the efficiency, which makes it possible to deal with large meshes. Unfortunately, the gain in deformation efficiency may cause the problem of deformation error on deformed surfaces and inverted boundary layer meshes near wall.
Computational fluid dynamics (CFD) and computational structural dynamics (CSD) based FSI methodology become attractive for aeroelasticity analysis in the time domain. As CFD and computer technology progress, higher-order methods based on RANS coupled with CSD are preferable because they are able to model more accurately transonic, nonlinear, and viscous effects, compared with linear aerodynamic models. However, the traditional compressible RANS solver may lead to low precision and low efficiency solutions at low Mach number, on account of the overlarge numerical dissipation and condition number [6]. To overcome these problems, an effective method is to use the preconditioned method with multiplying the preconditioning matrix to the time-derivative term, which can reduce the disparity between the convective and acoustic wave velocities through the modification of the eigensystem of the governing equations. The preconditioned method has already been applied to the flow analysis [7][8][9], but it hasn't been coupled with CSD for FSI application yet. In the current study, a complete integrated CFD-CSD computational method is implemented for aeroelasticity analysis in the time domain. CFD solver adopts the preconditioned Navier-Stokes governing equations, which can provide accurate flow solutions for both compressible and incompressible flows, are solved to evaluate aerodynamic loads, especially the nonlinear loads in the transonic regime. The structural model is tightly coupled with the fluid model by a pseudo time-stepping temporal discretization method. A modified mesh deformation method based on reduced control points RBF is utilized to update the fluid volume mesh, and an RBF based mapping algorithm is employed for data exchange on the fluid-structure interface. To evaluate the validity of the developed methodology, the flutter boundary and the limit cycle oscillation of Isogai wing and the flutter boundary of AGARD 445.6 wing, both of which have typical nonlinear characteristics, are analyzed and validated by experimental data.

Preconditioned Governing Equations
The unsteady preconditioned Navier-Stokes governing equations are adopted to predict the aerodynamic loads during the analysis of aeroelasticity from incompressible to compressible flows [6]: whereŴ = (p u v w T) T /J is the primitive variable vectorQ = ρ ρu ρv ρw ρe T /J is the conservative variable vector. The preconditoner Γ has various forms and the type proposed by Weiss-Smith is introduced here: , and U r is the reference velocity given by: a is the sound speed, M ∞ is the free-stream Mach number, and k is a constant to avoid the singularity of eigenvector matrix near stagnation points. Let k be 0.5 here.
The preconditioned equations may degenerate into the traditional compressible ones for the compressible flows. For low Mach number flows, the preconditioned method relieves the disparity between the convective and acoustic wave velocities, which is beneficial for the efficiency and accuracy of the solution.

Temporal Discretization Method
The dual time-stepping method is introduced to solve the unsteady flow around the moving structure. The method iterates over the pseudo time step several times until it converges, then updates the flow field in the physical time step. Discretize the pseudotemporal term ∂Ŵ ∂τ with first-order forward difference algorithm and the real-temporal term ∂Q ∂t with second-order backward difference method, as a result, the governing equations can be written as: , M = ∂Q ∂Ŵ , m is the iteration number of pseudo time step, and n is the iteration number of the physical time step.
Implicit lower upper decomposition method, using the symmetric Gauss-Seidel method (LU-SGS) time marching scheme is employed to accelerate the convergence rate of the solution. Then update the flow field when the pseudo time-stepping iteration stops:

Isogai Wing
The test case focuses on the cross section of a sweptback wing and investigates its aeroelastic response to the unsteady aerodynamic loads. Isogai wing section model consists in a 2D elastically mounted airfoil, which moves in pitching about its elastic axis and plunging vertically.
The structural dynamic governing equations of the elastically mounted system can be described as follows: m where h and α are plunging and pitching displacements, respectively, K h and K α are mounted spring stiffnesses in two degrees of freedom, m is the mass of the wing per unit span length, S α and I α are the static moment and the inertia moment about the elastic axis, L and M are the aerodynamic lift and the moment about the elastic axis, respectively. The equations can be nondimensionalized with the half-chord length b and the free-stream sound speed α ∞ , and the resultant equations can be given as: where V f = U ∞ ω α b is the reduced velocity, r d = m/(ρπb 2 ), r ω = ω h ω a , C l and C m are lift coefficient and moment coefficient respectively. Using a matrix form, Equation (8) can be written as: [M] The nondimensional structural dynamic governing equations are solved as follow where p is subiteration step number of the pseudo time-stepping method.

AGARD 4445.6 Wing
Based on Rayleigh-Reiz method, the generalized structural dynamics equations be given in the form of second-order ordinary differential equations: mass matrix, damping matrix, and stiffness matrix of the structure, respectively. The g eralized force is determined as: where Φ is Ritz basis function satisfying the displacement boundary conditions having sufficient continuity, and p C is the aerodynamic load coefficient on the struct The generalized force is a key factor coupling the structural model and the fluid mo and it depends merely on the fluid model if the linear structural model is adopted.
When the generalized displacement is solved, the actual displacement of the str ture can be determined as: The nondimensional structural dynamic governing equations are solved as follows: where p is subiteration step number of the pseudo time-stepping method.

AGARD 4445.6 Wing
Based on Rayleigh-Reiz method, the generalized structural dynamics equations can be given in the form of second-order ordinary differential equations: are generalized mass matrix, damping matrix, and stiffness matrix of the structure, respectively. The generalized force is determined as: where Φ is Ritz basis function satisfying the displacement boundary conditions and having sufficient continuity, and C p is the aerodynamic load coefficient on the structure. The generalized force is a key factor coupling the structural model and the fluid model, and it depends merely on the fluid model if the linear structural model is adopted. When the generalized displacement is solved, the actual displacement of the structure can be determined as: Taking the generalized normal mode displacement as the Ritz basis function for the generalized structural dynamic equations and Equation (11) can be transformed into linear equations: ω and ς are modal frequency and damping vector respectively. To be tightly coupled with the fluid model, the linear Equation (14) is discretized into a pseudo time-stepping scheme: when φ = 0.5 and p → ∞ , Equation (15)

Mesh Deformation Method
For volume mesh deformation driven by surface motion in FSI problems, mesh deformation method based on radial basis functions (RBF) is a powerful solution. RBF mesh deformation method with reduced point selection proposed by Rendall can achieve high quality deformed mesh efficiently for arbitrary topology. However, the improvements in deformation efficiency cause the problems of deformation errors on deformed surfaces and inverted boundary layer meshes near the wall.
An improved option based on reduced point selection RBF method is proposed, details of the method can be found in the previous work [10]. The method selects two groups of control points from the object surfaces, and the one selected by the greedy method is used to roughly calculate the mesh position and the deformation error. Then surface points with deformation errors exceeding the tolerance are picked as the second group, which is utilized to interpolate the deformation errors to the mesh points nearby. Thus, the deformation errors on the surface are greatly reduced and the inverted boundary layer meshes near the wall are effectively avoided.

Fluid-Structure Coupling Procedure
In an aeroelasticity problem, the flexible structure vibrates in response to the unsteady aerodynamic force, and in turn, the vibration affects the flow field around the structure. By discretizing the fluid and structural dynamic governing equations with the same pseudo time-stepping method, a tightly coupled fluid structure interaction method is established to analyze the aeroelasticity problem in the time domain. Figure 1 depicts the solution procedure of the tightly coupled FSI methodology. At the beginning of the solution, the flow field and the structure properties are initialized with given parameters. Then in each physical time step, there are several sub-iteration times between the unsteady fluid dynamic equations and structural dynamic equations. Within each physical time step, the unsteady CFD solver predicts the aerodynamic loads, according to which the structural displacements are determined. At the end of each iteration, volume mesh is updated based on the structural displacements for the next physical time step until the stop criterion is satisfied.
The surface mesh points in fluid model and structural nodes usually don't coincide at the fluid-structure interface, hence a mapping algorithm is required in the coupled method. The algorithm transfers aerodynamic loads from the aerodynamic surface to the structural nodes, and the predicted displacement on the structural nodes are interpolated to the surface mesh points of the fluid mesh. Various algorithms, including finite-plate spline (FPS) [11], infinite-plate spline (IPS) [12], multiquadric-biharmonic (MQ) [13,14], inverse isoparametric mapping (IIM) [15], and non-uniform B-spline (NUBS) [16], have been developed for exchanging information on the interface of fluid and structure. The IPS method is extensively used in programs such as MSC/NASTRAN, ENS3DAE, CFL3DAE. An RBF-based mapping algorithm [17], which is a three-dimensional extension of the IPS method, is employed for interface mapping here: where the subscripts f and s stand for variables in fluid model and structural model, respectively, and the interpolation matrix [G] is determined as follows: 1 x s ns y s ns z s ns φ s 1 s ns · · · φ s ns s ns where φ is the basis function, n f is the total number of the surface mesh and ns is structure nodes number.  The surface mesh points in fluid model and structural nodes usually don't coincide at the fluid-structure interface, hence a mapping algorithm is required in the coupled method. The algorithm transfers aerodynamic loads from the aerodynamic surface to the structural nodes, and the predicted displacement on the structural nodes are interpolated to the surface mesh points of the fluid mesh. Various algorithms, including finite-plate spline (FPS) [11], infinite-plate spline (IPS) [12], multiquadric-biharmonic (MQ) [13,14], inverse isoparametric mapping (IIM) [15], and non-uniform B-spline (NUBS) [16], have been developed for exchanging information on the interface of fluid and structure. The IPS method is extensively used in programs such as MSC/NASTRAN, ENS3DAE, CFL3DAE. An RBF-based mapping algorithm [17], which is a three-dimensional extension of the IPS method, is employed for interface mapping here:

Results and Discussion
In the following FSI simulations, fluid dynamics computations are performed by an in-house structured Navier-Stokes solver developed in Institute of Mechanics, Chinese Academy of Sciences for aerodynamic applications in the field of aeronautics and astronautics [6]. The flow solver, designed for both 2D and 3D CFD simulations, is a finite volume structured flow solver which solves the Reynolds-averaged Navier-Stokes equations using a cell-centered approach. Tecplot is used as data post-processing software.
To evaluate the methodology, the flutter boundary and the limit cycle oscillation of Isogai wing and the flutter boundary of AGARD 445.6 wing are analyzed.

Flutter Boundary and Limit Cycle Oscillation of Isogai Wing
Firstly, the tightly coupled FSI methodology is applied to analyze the Isogai wing. The model simulates the bending and torsional motion of a cross section of a sweptback wing. It consists of two DOF, plunging and pitching for a NACA 64A010 airfoil. Two aeroelastic phenomena with typically nonlinear characteristics, including the transonic flutter boundary and the LCO response, are simulated to validate the method.
The test case A of Isogai wing is investigated, and the nondimensional structural parameters are: x α = 1.8, r ω = 1, r 2 α = 3.48, r d = 60 and the elastic axis located at the half-chord length. The responses to three different reduced velocities near the flutter boundary at the Mach number of 0.825 are calculated, and the corresponding displacement histories are plotted in Figure 2. The displacements in pitching and plunging coordinates oscillate with equi-amplitudes when V * = 0.63, the system reaches a neutrally stable state known as flutter and the corresponding reduced velocity is the flutter speed index. At a lower reduced velocity, V * = 0.60, the displacements decay with time and the system is stable for the condition. However, the displacements diverge when a higher reduced velocity is imposed, which may bring disastrous damage to the aeroelastic system. Similar simulations are implemented at different Mach numbers in transonic region, and the flutter boundary is obtained and shown in Figure 3. The speed index flutter boundary is S-shaped with a dip near the Mach number of 0.85, and the simulation result agrees well with other simulation results [18][19][20].     For the transonic flow, the shock wave appears in the flow field around the structure, which arouses nonlinear aerodynamic loads for the structure. As a response, the aeroelastic system performs a nonlinear behavior. At the Mach number of 0.85, LCOs are observed with V * = 0.80, a reduced velocity higher than the flutter speed index. During the simulation, the model is forced to oscillate sinusoidally with a given amplitude for a cycle, the plunging coordinate is activated since the second cycle and the aeroelastic response begins in the meantime. Figure 4 describes amplitudes of LCO responses to different initial movements, α 0 = 1 • , α 0 = 2 • and α 0 = 8 • . The pitching and plunging displacements diverge in the first several cycles, then come to equi-amplitude vibrations and the amplitudes are independent of the initial movements. It is found that the responses get to the LCO state more rapidly with the initial movement amplitude increasing. Figure 5 depicts the corresponding phase portraits of the two coordinates, and the responses to different initial movements in both Dofs reach the same orbit, respectively. For the transonic flow, the shock wave appears in the flow field around the structure, which arouses nonlinear aerodynamic loads for the structure. As a response, the aeroelastic system performs a nonlinear behavior. At the Mach number of 0.85, LCOs are observed with 0 80 V = * . , a reduced velocity higher than the flutter speed index. During the simulation, the model is forced to oscillate sinusoidally with a given amplitude for a cycle, the plunging coordinate is activated since the second cycle and the aeroelastic response begins in the meantime. Figure 4 describes amplitudes of LCO responses to different initial movements, 0 1 α =  , 0 2 α =  and 0 8 α =  . The pitching and plunging displacements diverge in the first several cycles, then come to equi-amplitude vibrations and the amplitudes are independent of the initial movements. It is found that the responses get to the LCO state more rapidly with the initial movement amplitude increasing. Figure 5 depicts the corresponding phase portraits of the two coordinates, and the responses to different initial movements in both Dofs reach the same orbit, respectively.      For the two-dimensional test case, the flutter boundary in the transonic region is predicted by the tightly coupled FSI method and it is similar to those predicted by other researchers. The simulation results reveal the nonlinear response to the periodic nonlinear aerodynamic loads. The coupled CFD-CSD method can estimate the flutter boundary accurately.

Flutter Boundary Analysis of AGARD 445.6 Wing
AGARD 445.6 wing is a standard aeroelastic configuration and its flutter responses are extensively simulated to evaluate and assess the FSI methodology. The geometry and experimental data are provided by the NASA Langley Research Center. The flutter response of the weakened structure model mounting on the wall is analyzed.
To evaluate the uncertainties of the CFD results brought by the mesh resolution, grid independence studies are conducted on five multiblock structural meshes, named by Mesh1 with about 7000 cells, Mesh2 with about 15,030 cells, Mesh3 with about 30,000 cells, Mesh4 with about 42,000 cells and Mesh5 with 78,340 cells. In the test case, the free-stream Mach number is 0.96 and the angle of attack is 2 • . The drag coefficients corresponding to the five meshes are compared in Figure 6. It can be seen that at first the drag coefficient shows a decreasing trend with grid number increasing, then it tends towards stability even the grid number increases further. Hence, Mesh4 is employed as the final grid.
A multiblock structural fluid mesh with 42,000 cells and 121 structure nodes are adopted in the simulation, the surface mesh and the structure nodes are plotted in Figure 7. The first four vibration modes, first bending, first torsion, second bending and second torsion, are utilized in the structure dynamic model. Further, the frequencies of these four modes are 9.6 Hz, 38.17 Hz, 48.35 Hz and 91.54 Hz, respectively. The modal shapes of the four modes are shown in Figure 8.
The cross-section of the AGARD 445.6 wing is the NACA 65A004 symmetric airfoil, and the responses are simulated at the angle of attack of 0 • with a small disturbance imposed to motivate the model at the beginning of the simulation. The dynamic aeroelastic responses at six different Mach numbers ranging from 0.499 to 1.141 are analyzed and corresponding generalized displacement responses on the flutter boundary are shown in Figure 9. On the flutter boundary of each case, all the generalized displacements perform the equi-amplitude oscillation with the same frequency, the aeroelastic system is neutral stable at the corresponding speed index.
independence studies are conducted on five multiblock structural meshes, named by Mesh1 with about 7000 cells, Mesh2 with about 15,030 cells, Mesh3 with about 30,000 cells, Mesh4 with about 42,000 cells and Mesh5 with 78,340 cells. In the test case, the free-stream Mach number is 0.96 and the angle of attack is 2°. The drag coefficients corresponding to the five meshes are compared in Figure 6. It can be seen that at first the drag coefficient shows a decreasing trend with grid number increasing, then it tends towards stability even the grid number increases further. Hence, Mesh4 is employed as the final grid. A multiblock structural fluid mesh with 42,000 cells and 121 structure nodes are adopted in the simulation, the surface mesh and the structure nodes are plotted in Figure  7. The first four vibration modes, first bending, first torsion, second bending and second torsion, are utilized in the structure dynamic model. Further, the frequencies of these four modes are 9.6 Hz, 38.17 Hz, 48.35 Hz and 91.54 Hz, respectively. The modal shapes of the four modes are shown in Figure 8.  The cross-section of the AGARD 445.6 wing is the NACA 65A004 symmetric airfoil,  The cross-section of the AGARD 445.6 wing is the NACA 65A004 symmetric airfoil, and the responses are simulated at the angle of attack of 0° with a small disturbance imposed to motivate the model at the beginning of the simulation. The dynamic aeroelastic The predicted flutter speed indexes and flutter frequency ratios are shown in Tables  1 and 2, respectively, together with the wing tunnel results. The results computed by the In transonic regime, M = 0.901 and M = 0.960, the results are reasonable with the relative error under 7%. For the supersonic cases, M = 1.072 and M = 1.141, however, coupled method overpredicts the flutter velocities by about 30% higher than the experimental data. The simulations predict the flutter frequency ratio in better accuracy compared with the flutter speed indexes, and the largest relative error is less than 10%. The simulation results are also compared with those predicted by other researchers. The comparisons of flutter speed indexes and the flutter frequency ratios are shown in Figures 10 and 11. As is noted that [20,21] also overpredicted the flutter speed indexes in the supersonic region, and the larger differences are observed in those predicted results. Compared to their results, our results show some improvement by performing preconditioned Navier-Stokes calculation. The amount of improvement is relatively large, especially in the prediction of flutter frequency compared to the differences between other computational results and the experiment. The flutter frequency ratios predicted in the paper are closer to the experiment compared with those of [20,21].    The results show that the tightly coupled preconditioned N-S equations with CSD methodology provides reasonable flutter speed indexes in the subsonic and the transonic regions, for the supersonic cases there is some improvement compared with other simulations. Note that the methodology predicts the flutter frequency ratios well even in the supersonic region, while other computational methods cannot still provide satisfying results in the supersonic cases.
Generally speaking, the flutter speed values are higher in the computational solutions than the experimental ones in the supersonic regime. Some researchers attribute this phenomenon to the inaccuracy of the experimental data [22,23]. It was believed that it was The results show that the tightly coupled preconditioned N-S equations with CSD methodology provides reasonable flutter speed indexes in the subsonic and the transonic regions, for the supersonic cases there is some improvement compared with other simulations. Note that the methodology predicts the flutter frequency ratios well even in the supersonic region, while other computational methods cannot still provide satisfying results in the supersonic cases.
Generally speaking, the flutter speed values are higher in the computational solutions than the experimental ones in the supersonic regime. Some researchers attribute this phenomenon to the inaccuracy of the experimental data [22,23]. It was believed that it was difficult to model correctly the physics processes occurring in the experiment at supersonic case for only minimal information on the experimental flow field was available. The differences between the predicted and measure results have also been noted in the previous studies, which is commonly believed that there exist measurement errors in the experiment of supersonic regime.

Conclusions
A tightly coupled FSI methodology based on computational fluid dynamics (CFD) and computational structural dynamics (CSD) is developed for aeroelasticity analysis in the time domain. Preconditioned N-S governing equations are used as flow solver and tightly coupled with structural model by the pseudo time-stepping temporal discretization method. A modified mesh deformation method on reduced control points RBF is utilized to update the fluid volume mesh, and an RBF based mapping algorithm is employed for data exchange on the interaction interface.
The flutter boundary of the Isogai wing Case A model is analyzed. An S-shaped flutter boundary with a dip near the Mach number of 0.85 is obtained, and the simulation results are in good agreement with those of other simulation results. The LCO response of the case is observed, with a higher reduced velocity index than the flutter speed index, the pitching and plunging displacements come to equi-amplitude oscillation finally and the vibrational amplitudes are independent of the initial movements. The responses of different initial movements in both two coordinates reach the same orbit respectively in phase portraits.
The flutter boundary of AGARD 445.6 wing is predicted. Results show that the methodology provides reasonable flutter speed indexes in the subsonic and the transonic regions, for the supersonic cases, the simulation results are overpredicted. However, some improvement is achieved by our method compared with other simulation results from previous studies. The methodology predicts the flutter frequency ratios well even in the supersonic region as well as in subsonic and transonic regimes, while simulation results computed by other three methods still have larger differences with the experiment.
The simulation results show that the tightly coupled FSI methodology with preconditioned N-S equations and CSD models can predict the aeroelasticity in the time domain accurately, and can also reveal the nonlinear characteristics of the aeroelastic system.