Nonlinear Panel Flutter Analysis Based on an Improved CFD/CSD Coupled Procedure

Xiaomin An,and Min Xu

1 Introduction

Aeroelasticity is an interesting subject of the interaction between flexible structure and the surrounding fluid.Nonlinearities involved with aeroelasticity arise from two aspects:the structural and aerodynamic points of view.Such as large shock motions, flow separation and geometrical nonlinearity in structure.At present,most of the nonlinear aeroealstic analysis depends on computation of unsteady aerodynamic forces and structural linear mode superposition method,which considers only the nonlinearity in fluid regime.However,the structural geometrical nonlinearity plays a crucial role and the linear theory is not able to determine accurately the response in the case of large deformation for large flexible structure.There are strong requirements in discussing approaches for solving two nonlinearities involved in aeroelasticity simultaneously and studying computational interfaces adequate to handle the nonlinear interdisciplinary interactions.Studying this phenomenon,well validated computational fluid dynamics(CFD)methods should be coupled with computational structural dynamics(CSD)in an accurate way[Baum,Luol,Mestreaul,Sharov,Lohner,Prlessone and Charman(2001);Patil and Hodges(2000)],and this will provide a path for developing a comprehensive prediction capability for simulating complex,nonlinear aerodynamics and structural dynamics[Gordnier and Fithen(2003)].

Dowell[Dowell(1970);Dowell,Thomas,and Hall(2001)]classified the nonlinear aeroelasticity into four different categories and studied limit cycle oscillations(LCOs)of the panel in supersonic flow,as well as proposed reduced order aerodynamic models-proper orthogonal decomposition and the harmonic balance technique.The latter method has been developed and applied by several researchers,such as Hall[Hall,Thomas,and Clark(2002)],Dai[Dai,Schnoor,and Atluri(2012)]and Liu[Liu,Thomas,and Dowell(2007)].Gordnier[Gordnier and Fithen(2003);Gordnier and Visbal(2004)]proposed a computational aeroelastic model consisted of a finite element solution of the nonlinear von Karman plate equations coupled to both an Euler and Navier-Stokes aerodynamic models,and performed the model on a three dimensional panel and a flexible delta wing.Mei[Mei(1977)]applied the finite element method and third-order piston theory to study the supersonic limit cycle oscillations of two-dimensional panels.Strganac[Strganac,Cizmas,and Nichkawde(2005)]studied bifurcation of vehicles with moderate-to-high deformations and presented the nonlinear structural equations of motion for cantilevered wing con figurations that possess in-plane,out-of-plane,and torsional couplings,and outlined a Navier-Stokes based unsteady aerodynamic model that addressed viscous and compressible flows for transonic flows with shock/boundary layer interaction.Patil and Hodges[Patil and Hodges(2000);Patil,Hodges,and Cesnik(2001)]presented a mixed variational formulation based on the exact intrinsic equations for the dynamics of beams in moving frames and studied finite-state airloads for deformable airfoils on fixed and rotating wings,and performed the theory on High Altitude Long Endurance(HALE)aircraft aeroelastic computation.Attar[Attar and Gordnier(2005)]studied the flutter and LCO behavior of a cropped delta wing using a newly developed computational aeroelastic solver which included Euler finite difference solver and nonlinear structural model.The structural solver for geometric nonlinearities was modeled by co-rotational formulation,which was ap-plied in aeroelastic analysis for the first time.

Most of flight vehicles’structure are generally modeled by beam or shell finite element in engineering analysis.Although many algorithms of large rotation solution have been developed as Total Lagrangian(TL)and Updated Lagrangian(UL),most of these early works were surprisingly sophisticated and worth improvement.The co-rotational approach,viewed as an effective way to describe nonlinear structural kinematics analysis,has generated an increasing amount of interest in the last decade[Cris field(1996)].The main idea is that the motion of the element can be split into a rigid part and a pure deformational response.Assuming the pure deformational part measured in a rotated reference frame to be small enough,a linear finite element theory can be used in the analysis[Battini(2002);Battini(2007);Pacoste(1998)].This leads to very simple expressions for the local internal force vector and tangent stiffness matrix.Cris field[Cris field(1996)]developed a unified co-rotational framework for solids,shells and beams,and applied it to facet shell elements and three-dimensional beams.Battini[Battini(2002)]implemented co-rotational beam elements and branch-switching procedures in order to analyze elastic and plastic instability problems,and proposed a new numerical method based on a minimal augmentation procedure for the direct computation of elastic critical points.Battini[Battini(2007)]also did some work to modifications in triangular shell elements.Pacoste[Pacoste(1998)]investigated the formulation of co-rotational fl at facet triangular elements for the numerical analysis of instability phenomena in shell structures considered three types of local formulations.Cai[Cai,Paik and Atluri(2009)]and Zhu[Zhu,Cai,Paik and Atluri(2010)]proposed an updated Lagrangian corotational frame for analyzing large rotations,and large deformation of plates and shells.

For structural dynamic response solution,Galvaneito and Cris field[Galvaneito and Cris field(1996)]presented an exact energy-conserving procedure for the implicit nonlinear dynamic analysis of planar beam structures.Chimakurth[Chimakurth and Cesnik(2011)]utilized a nonlinear finite element solver based on a co-rotational approach and simulated flapping wing structure under large displacements and rotations.Relvas[Relvas and Suleman(2007);Relvas and Suleman(2006)]applied the finite element co-rotational theory to model geometrically nonlinear structures and utilized a dynamic co-rotational energy-conserving algorithm to solve the nonlinear structural response.The unsteady Euler solver and vortexring method were used to model aerodynamic loads for a two-dimensional airfoil and a plate under low speed air fl ow respectively.

In subsystem,the CFD grid is different from the CSD grid on account of the different reference frame between Euler and Lagrange.In aeroelastic analysis,there are two kinds of data transfer on common boundary:one is interpolating the displace-ments calculated by the CSD solver to CFD grid,and the other is transforming the aerodynamic loads computed by CFD code on CSD points.At present,the most popular interpolation methods[Goura and Badcock(2001);Kim and Kim(2005);Sadeghi,Liu,and Lai(2004)]can be categorized as two groups:(1)Surface fitting method,which utilizes a global function with weighting coefficients,such as In finite-Plate Spline(IPS),this kind of method has dominantly used for displacement transfer;(2)Surface tracking method,which utilizes the shape functions of finite elements by projecting fluid nodes to nearest structural elements,for example Conservation Volume Transformation(CVT).Currently,surface tracking method is mostly used for load transfer.

However,the above popular methods:IPS and CVT cannot keep energy conservation on the boundary between CFD and CSD module.The boundary element method(BEM),which was generally used for solid mechanics analysis,was firstly proposed by Chen and Jadic[Chen and Jadic(1998)]to generate the universal spline matrix for CFD/CSD data transfer.Then,BEM scheme was developed to deal with complex con figuration,including structural discontinuity like control surfaces by Chen and Gao[Chen and Gao(2001)].The approach was developed to handle large and complex structures and transform the displacement vectors at a solid boundary to the interior of the field grid as mesh deformation algorithm by Lai[Lai,Tsai,and Lum(2002);Lai,Tsai,and Liu(2003)].It is a pity that the BEM technique has not been employed in aeroelastic simulation in any literatures.

So far,time marching algorithm of CFD/CSD coupled method can be classified broadly under four major categories[Kamakoti and Shyy(2004);Farhat and Van(2006)](as seen in Fig.1):

(i)Fully coupled method can solve the two solvers simultaneously and accurately,however,it is highly nonlinear and limited to two-dimensional problems.

(ii)Loosely coupled method has advantages of simplify explicit/implicit treatment,sub-iteration and software modularity,but it is only first order time accuracy,therefore,the time-step requires small enough.

(iii)Tightly coupled method is developed to eliminate the time-lag,and it can achieve second order time accuracy,however,the computational cost increase greatly.

(iv)The fourth CFD/CSD coupled procedure in Fig.1,named as second order loosely coupled method,was proposed by Farhat[Farhat and Van(2006);Geuzaine,Van,and Farhat(2004)]and proved to reach second order time accuracy if the solver of the subsystem was selected appropriately.It has the same advantages of the loosely coupled method.

Figure 1:Description of existing CFD/CSD coupled procedures.(a)fully coupled method,(b)loosely coupled method,(c)tightly coupled method,(d)second order loosely coupled method

In the present article,the CSD solver will use a predictor-corrector program with an approximate energy conservation algorithm to solve the nonlinear structure response.The CFD solver will utilize an AUSMpw+ flux splitting scheme,an implicit dual-time stepping technology and the geometric conservation law to solve Euler and Reynolds averaged Navier-Stokes equation.We use the BEM method which is most robust and accurate for energy mapping conservation by combining of the load and displacement transfer,and we will apply Farhat’s second order loosely coupled procedure for connecting the CFD program with CSD solver.This improved program is used to analyze the fl utter characteristics of two-dimensional and three-dimensional panels for supersonic,transonic and subsonic Mach numbers,as well performed on the solution of nonlinear elastic response of AGARD 445.6 wing.

2 Unsteady Aerodynamic Load Model

Euler equation and Reynolds averaged Navier-Stokes equation are applied to establish the model for solving aerodynamics,taking Euler equation as an example:

Wherecis sound velocity,pis fl ow pressure andwhereMis Mach number.

In order to solve unsteady flow,a sub-iterative course of the dual-time technology is introduced to eliminate the time marching error.Suppose that fluid grid moves with the body flexibly,and then the introduced dual-time term can be written as:

Where V is the velocity of the moving grid.

Note that the above program also can be used to solve Reynolds averaged Navier-Stokes equation.

Figure 2:Panel geometry

3 Geometrically Nonlinear Model of Panel

A flexible panel of lengtha,widthb,thicknessh,and mass densityρsshown in Fig.2 is investigated in the present work.The air flowing above the panel at Mach numberM∞and densityρ∞is in the positivexdirection.Two cases are studied:two-dimensional case asa/b=0 and three-dimensional case asa/b=1.The panel’s characteristics are:h/a=0.002,mass ratioμs=ρ∞a/ρsh=0.1,Poisson’s rationν=0.3 and modulus of elasticityEs=7.0×1010.The panel is pinned at both edges in the fl ow direction.

3.1 Two-Dimensional Panel Model

For two-dimensional casea/b=0,the panel can be considered as planar beam structure,which can be described by co-rotational frame as shown in Fig.3[Cris field(1996);Battini(2002)].

The vectors of global and local displacements of the element are defined bydanddlas

Letc=cos βands=sin β,through differentiation of Eq.9,and utilize geometry connection,then transformation matrix T between local and global displacements can be obtained by

Figure 3:Co-rotational frame of two-dimensional panel

Figure 4:Co-rotational frame of three-dimensional panel

Then utilize their differentiations,the second term in equation 13 can be obtained and the tangent stiffness equation can be finally written as

3.2 Three-Dimensional Panel Model

For three-dimensional casea/b=1,the panel can be modeled by triangular shell structure as described in Fig.4.

The origin is taken at the geometric center pointC,and an orthogonal matrix to specify the orientation of the local frame is defined as

Whereeidenote the three axes vectors.The rigid translation and rotation are described by the displacementsof nodeCand orthogonal matrixrespectively.The nodal local pure deformational values are defined byˉuiandas shown in Fig.4.

Then,the local displacements can be derived as

The rotational values of angular variable and rotational vector can be related as

The local displacements and internal forces can be written as

In local coordinate system,the above internal forces and displacements can be expressed in linear relation as

Where Klis the linear triangular shell element local stiffness matrix.The global displacements and forces are de fi ned as

Through differentiation of Eq.23,and utilize geometry connection,then transformation matrix between local and global frame can be obtained by

The above variables are computed by the defined spatial angular variation,it has to transform the rotational values to rotational vector.New definitions of the displacements and internal forces are given as

4 An Approximate Energy Conservation Algorithm

Newmark integral approach,an unconditional stable and predominant procedure to determine the structure dynamic response at each time step for linear problems,may encounter"lock"or"blow up"solution for nonlinear case.It has been verified by Galvanetto and Cris field[Galvaneito and Cris field(1996)],and an energy conservation algorithm was proposed for the implicit nonlinear dynamic analysis of planar beam structures in literature[Galvaneito and Cris field(1996)],in which the time stepwas used to achieve energy conservation duringnton+1 step.However,it may appear more complex and more data have to be stored.Based on this algorithm,an approximate energy conservation algorithm is developed,which retains the stability and accuracy,and simplifies the procedure without any modification.

Assume that there is no damp in the structural system,the equilibrium equation atn+1 step can be written as

Where M andFs,n+1are the mass matrix and the external force as unsteady aerodynamic load respectively.By introducing midpoint rule and the variables of step as

Letgn+1/2→0,then the energy change betweennandn+1step?E=gn+1/2?un+1/2could be conservative.Differentiate the Eq.30 and the new global tangent stiffness matrix can be obtained as

It is notable that Eq.32 can be used for both beam and shell structure.

4.1 Predictor Step

The initial values of displacement increment are determined by

4.2 Corrector Step

Secondly,the displacements and velocities can be updated byas

The above two equations 35 and 36 should be solved repeatedly in the corrector process until a convergence criterion withis satis fi ed.

5 Data Transfer Via Boundary Element Method

If the displacementsxon the entire boundary(CFD surface grid)are known,a transformation from boundary to interior displacementsuof the internal source(CSD grid)can then be performed by a transformation matrix B as

Where H and G are formed from the kernel integral of the displacement and traction,and the subscriptaandsstand for CFD surface grid and CSD grid values respectively.Via the minimum strain energyrequirement[Chen and Jadic(1998),the spline matrix is determined by

By embracing all CFD and CSD gird points,an intermediate BEM model is constructed as a third platform[Chen and Gao(2001)].Firstly,a universal spline matrix D can be generated to transfer the information of the CSD grid to the intermediate BEM model by Eq.38,then a BEM matrix B is constructed to transfer the information on the intermediate BEM model to the CFD grid by Eq.37.Finally the mapping matrix between two different grid systems can be obtained

Thus,the structural displacements and aerodynamic loads can be transferred by

Whereuandxare the displacements of structural and CFD surface grid respectively.FsandFadenote the vector of structural nodal forces and aerodynamic loads respectively.Thus,the requirements of the conservation of work between two transformations can be satisfied as

6 CFD/CSD Coupled Program

Based on the Farhat’s second order loosely coupled procedure[Farhat and Van(2006)],the two solvers are connected as following six steps:

(i)Predict the structural displacements at time-stepbynstep

(ii)Transfer the predicted motionun+1/2to the fluid system asxn+1/2by Eq.40.

(iii)Update the position of the fl uid grid by TFI technique and compute the new control finite volume by Eq.8.

(iv)Solve Eq.7 to obtain the loads on aerodynamic surface.

(v)Convert the aerodynamic loadsFa,n+1/2into structure element to get equivalent loadsFs,n+1/2by Eq.40.

(vi)Advance Eq.31 and approach the motion atn+1 step by the predictor corrector program.

7 Numerical Results

7.1 CFD Solver for AGARD CT5 Unsteady Test

Two cases of AGARD CT5 unsteady test are introduced to validate the present CFD solver.Unsteady motion of the NACA0012 airfoil is defined as

Figure 5:(a)Lift coefficient vs.angle of attack in case 1,(b)moment coefficient vs.angle of attack in case 1

Figure 6:(a)Lift coefficient vs.angle of attack in case 2,(b)moment coefficient vs.angle of attack in case 2

Flow conditions are set as follows:(1)Case1:M∞=0.6,αm=3.16o,α0=4.59oandk=0.0811.The axis position is 0.273 chord.(2)Case2:M∞=0.755,αm=0.016o,α0=2.51oandk=0.0814.The axis position is 0.25 chord.

In both two cases,the CFD gird is constructed by O grid consisting of 121×80 points and Euler equation is solved with the total time step numbers in one period is 80,and the number of sub-iteration within one time step is 10.The computed results of the lift and moment coefficient in Fig.5 and 6 are the first 7 periods.The results are also compared with experimental values.The computation demonstrates reasonable agreement with the experimental data for both two cases.

Figure 7:(a)Cantilever beam subjected to an end moment,(b)cantilever beam nonlinear static analysis results

7.2 Cantilever Beam Nonlinear Static Analysis

The cantilever beam subjected to an end moment(shown in Fig.7(a))which has an analytical solution is solved by present finite co-rotational formulation.

The beam’s parameters are defined as:l=20m,EIz=10N·mandEA=10N.The end momentMis set to 0.2π,0.4π,0.6π,0.8πandπrespectively.While the exactly analytical solution is expressed asθ=(MI)/(EIz).The results presented in Fig.7(b)and Tab.1 agree well with the analytical solution.

Table 1:Comparison of the tip displacements of cantilever beam between co-rotational formulation and exactly analytical solution

7.3 Hinged cylindrical Segment Nonlinear Dynamic Simulation

Figure 8:Geometry of a hinged cylindrical segment

The nonlinear problems of a hinged cylindrical segment with a total thicknesst=0.00635mpresented in Fig.8 are analyzed.The total length and circumferential length areLa=Lb=0.508m,and curvature of the circular edges isR=2.54m.The material properties areE=3.105×109N·mandμ=0.3.The hinged cylindrical segment is taken under the constant central point loadP=500N.The time step is set to?t=0.0001sand?t=0.001srespectively.Then the dynamic responses at point A are computed by the linear Newmark method,nonlinear Newmark method and the improved predictor-corrector procedure with the approximate energy conservation algorithm.

The time histories of the normal displacement and total energy change via various schemes are shown in Fig.9.It can be seen that the linear results are obviously different from the other two nonlinear results as the nonlinearity has much affection in the structure.The nonlinear Newmark method has the similar response to the improved scheme in original time steps.However,with time boosting,nonlinear Newmark method tends to augment after some time steps.Also,the nonlinear Newmark scheme has wide variations in total energy changes while the improved procedure tends to zero,which keeps the energy conservation.It is evident that the conventional nonlinear Newmark integration has a limited stability in the nonlinear case.

When the time step is made larger with?t=0.001s,as shown in Fig.10,the linear results also present error curve,and the two nonlinear results show obvious difference.The amplitude of the nonlinear Newmark method is increased markedly and ends up with high energy state,about 8×104J,which leads to a lack of convergence and provides improper information for designers.Clearly,improved predictor-corrector procedure is keeping stable versus time steps,and can be set to larger time steps and gives satisfactory results.

7.4 Nonlinear Panel Flutter Analysis

(1)Two-dimensional panel fl utter results

Figure 9:(a)Time-histories of normal displacement at?t=0.0001s,(b)time histories of total energy change at?t=0.0001s

Figure 10:(a)Time-histories of normal displacement at?t=0.001s,(b)time-histories of total energy change at?t=0.001s

The finite element model composes of 21 Bernoulli beam elements.The CFD grid is constructed by H grid consisting of 161×41 points,and Euler equation is solved with the number of sub-iteration within one time step is 10.Mach numbers are chosen as:supersonic caseM∞=1.2,M∞=1.414,M∞=1.8 andM∞=2.0;transonic and subsonic caseM∞=1.02,M∞=0.95,M∞=0.9 andM∞=0.7.The fl utter computation are initiated as

Fig.11 and Fig.12 show the time-histories of the oscillation and phase portrait ofunder dynamic pressure=100 forM∞=1.2 case atx/a=0.75.It is clear that a representative limited cycle oscillation appears when the nonlinearity is considered.Fig.13 shows the limit cycle oscillation amplitude values of Dowell(linearized potential theory),Gordnier(Euler equation)and present work for a supersonic Mach number atx/a=0.75.WhenM∞=1.2,the three amplitude values are compared well with each other though Dowell used linearized potential theory.ForM∞=1.414 case,the results have good agreement with Gordnier’s work;these two solutions have lower values than Dowell’s results.For a higher MachnumberM∞=1.8,there are small difference between the present work and Gordnier’s solution.WhenM∞＜1.0,static equilibrium positions are obtained under a broad range of dynamic pressureλ.The results are compared well with Gordnier’work atx/a=0.5 as shown in Fig.14.In the present work,only one static equilibrium values are obtained for a same initial conditionFig.15 shows the fl utter boundary.Note that it will be singular whenM∞=1.0 in Dowell’solution.

(2)Three-dimensional panel fl utter results

A finite element structural model is constructed 400 triangular shell elements and 1323 degrees of freedom.The CFD grid is constructed by H grid consisting of 121×121×31.In this case,Euler equation is still solved with the number of sub-iteration within one time step is still 10.Mach numbers are chosen as:supersonic caseM∞=1.2,M∞=1.414 andM∞=1.6.Fig.16 shows that the time-histories of the oscillation underλ=300 for case atx/a=0.75.Note that the peak values have a small variety in every cycle.Fig.17 shows the pressure on section under four different positions.The curves for the two peak position are nearly symmetric.Fig.18 shows the limit cycle oscillation amplitude values of Dowell,Gordnier and present work for three-dimensional panel atx/a=0.75.ForM∞=1.2 case,the three amplitude values are also compared well with each other.While for higher Mach numbers:M∞=1.414 andM∞=1.6,there are evident difference between the present work and Dowell’s solution as we apply nonlinear aerodynamic solver.

7.5 Nonlinear Aeroelastic Response of AGARD 445.6 Wing

Figure 11:Time-histories of the oscillation at x/a=0.75

Figure 12:Phase portrait of w/h vs.˙w/h at x/a=0.75

Figure 13:Limit cycle oscillation amplitude vs.dynamic pressure for 2D panel at x/a=0.75

Figure 14:Static deflection vs.dynamic pressure at x/a=0.5

Figure 15:Dynamic pressure vs.Mach number

Figure 16:Time-histories of the oscillation at x/a=0.75

Figure 17:Pressure of panel under four position at x/a=0.75

Figure 18:Limit cycle oscillation amplitude vs.dynamic pressure for 3D panel at x/a=0.75

Figure 19:(a)Sketch of AGARD 445.6 wing con figuration,(b)CFD model with intermediate BEM model

The simulation of the nonlinear aeroelastic response of the AGARD 445.6 wing(constructed by NACA 65A004 airfoil)is also considered to test the prediction capability for more complex con figuration of the developed program.The material properties of the flexible wing areE1=3.15Gpa,E2=0.42Gpa,μ=0.3,G=0.44Gpaandρ=382kg/m3.Fig.19(a)shows the con figuration of the wing.Fig.19(b)shows the CFD surface model surrounded by intermediate BEM model.A finite element structural model is constructed 200 triangular composite shell elements and 363 degrees of freedom.The CFD grid is constructed by C grid consisting of 159×45×45,the surface mesh contains 2940 points.In this case,Reynolds averaged Navier-Stokes equation is considered to instead of Euler equation,and B-L turbulence model is used.The number of sub-iteration within one time step is 10.

The free stream condition areM∞=1.05,α=5oandq∞=8171kg/?m·s2?,and the time step is set to?t=0.0002s.The aeroelastic responses of the wing under large aerodynamic loads are computed by the coupled system.In structural solution,the developed predictor-corrector procedure with approximate energy conservation algorithm,linear Newmark method and linear mode superposition with first four modes are considered.Fig.20 shows the displacement curves of normal direction and phase portrait normal displacement vs.velocity at the point on trailing edge tip by three structural dynamic algorithms.It can be seen that both two linear algorithms present the similar augment amplitude whereas the developed nonlinear approach tends to the same variation after some time steps.A representative LCO appears when geometrically nonlinear effect is considered.It means that the nonlinearity can alleviate the divergence of the structural response.

Figure 20:(a)Time-histories of normal displacement on trailing edge tip,(b)phase portrait normal displacement vs.velocity on trailing edge tip

Fig.21 shows the histories of lift and drag coefficient of the wing.It can be seen that the lift coefficient curve of the improved predictor-corrector nonlinear algorithm has a similar LCO phenomenon to the above structural response in Fig.20,while the linear results tend to diverge.From Fig.21,the drag coefficient of the present method shows a cycle oscillation characteristic,however,the two linear methods show markedly different trends with time marching.It indicates that the error caused by a linear approach in structural solver will induce a much discrepancy in fluid solution,which has a sensitive effect for fluid convergence.

Figure 21:(a)Time-histories of lift coefficient of the wing,(b)time-histories of drag coefficient of the wing

Fig.22 shows normal displacements and loads distribution of the wing at timet=0.11875svia three methods.It is obvious that the wing has small displacements in X and Y direction on the tip when considering geometrical nonlinearity.When considering the deformation and loads distribution in aeroelastic analysis,the linear mode superposition method may be of no avail,and the linear Newmark will be useless for a nonlinear solution under large load and deformation.

Figure 22:(a)Normal displacements distribution of the wing at t=0.11875s,(b)Normal loads distribution of the wing at t=0.11875s

8 Conclusions

In this paper,we design an improved CFD/CSD coupled procedure for solving both aerodynamic and structural nonlinearity involving four parts:

(i)CFD solver-AUSMpw+ flux splitting scheme with dual-time implicit technology and geometric conservation law are introduced.It can solve Euler and Reynolds averaged Navier-Stokes equation and obtain unsteady aerodynamic load.

(ii)CSD solver-a predictor-corrector program with an approximate energy conservation algorithm based on co-rotational theory is developed.It can treat with geometrically nonlinear dynamic response of the beam and shell structure.

(iii)Data transfer between two different solvers is performed by the BEM method,which keeps energy conservation accurately and combines of the load and displacement transformation together.

(iv)A second order loosely coupled method is applied to manage time marching algorithm between CFD and CSD solver.

It is evident that the predictor-corrector procedure with an approximate energy conservation algorithm based on co-rotational frame in our CSD solver can deal with geometrically nonlinear dynamic structure effectively.Combined with finite volume method and dual-time technology,and through CFD/CSD coupled algorithm,the improved procedure can capture the representative LCO behavior for the solution of two-dimensional,three-dimensional panel and a complex wing,and it will provide a path for developing a comprehensive prediction capability for simulating complex,nonlinear aerodynamics and structural dynamics.

Acknowledgement:The first author acknowledges the financial support provided by the National Natural Science Foundation of China(Grant No.11212165).

Attar P.;Gordnier R.(2005): Aeroelastic prediction of the limit cycle oscillations of a cropped delta Wing.AIAA 2005-1915.

Battini J.(2002):Co-rotational beam elements in instability problems,Technical Reports from Royal Institute of Technology Department of Mechanics SE-100 44 Stockholm,Sweden.

Battini J.(2007): A modified corotational framework for triangular shell elements.Computer methods in applied mechanics and engineering,vol.196,pp.1905–1914.

Baum,J.;Luol H.;Mestreaul E.;Sharov D;Lohner R.;Pelessone D.;Charman C.(2001): Recent developments of a coupled CFD/CSD methodology.A01-31097.

Cai Y.;Paik J.;Atluri S.(2009): Large deformation analyses of space-frame structures,with members of arbitrary cross-section,using explicit tangent stiffness matrices,based on a von karman type nonlinear theory in rotated reference frames.CMES-Computer Modeling in Engineering and Sciences,vol.53,no.3,pp.117–145.

Chen P.C.;Jadic L.(1998): Interfacing fluid and strcutural models via an innovative structural boundary element method.AIAA Journal,vol.36,no.2,pp.282–287.

Chen P.C.;Gao X.W.(2001): A multi-block boundary element method for CFD/CSD grid interfacing.AIAA 2001-0715.

Chimakurth S.;Cesnik C.(2011):Flapping-Wing Structural Dynamics Formulation Based on a Corotational Shell Finite Element.AIAA Journal,vol.49,no.1,pp.128–142.

Cris fi eld M.(1996): A unified co-rotational framework for Solids,shells and beams.Journal of Solids structures,vol.33,no.20-22,pp.2969–2992.

Dai H.;Schnoor M.;Atluri S.(2012): A Simple Collocation Scheme for Obtaining the Periodic Solutions of the Duffing Equation,and its Equivalence to the High Dimensional Harmonic Balance Method:Subharmonic Oscillations.CMES:Computer Modeling in Engineering and Sciences,vol.84,no.5,pp.459–497.

Dowell E.(1970):Panel Flutter:A Review of the Aeroelastic Stability of Plates and Shells.AIAA Journal,vol.8,no.3,pp.385–399.

Dowell E.;Thomas J.;Hall K.(2001):Transonic limit cycle oscillation analysis using reduced order aerodynamic models.AIAA 2001-1212.

Farhat C.;Van Z.K.(2006):Provably Second-order Time-Accurate Llooselycoupled Solution Algorithms for Transient Nonlinear Computational Aeroelasticity.Comput.Methods Appl.Mech.Engrg,vol.195,pp.1973–2001.

Galvaneito U.;Cris field M.(1996):An energy-conserving co-rotational procedure for the dynamics of planar beam structures.International journal for numerical methods in engineering,vol.39,pp.2265–2282.

Geuzaine P.;Van Z.K.;Farhat C.(2004): Second-Order Time-Accurate Loosely-Coupled Solution Algorithms for Nonlinear FSI Problems,European Congress on Computational Methods in Applied Sciences and Engineering.

Gordnier R.;Fithen R.(2003): Coupling of a nonlinear finite element structural method with a Navier-Stokes solvers.American Institute of Aeronautics and Astronautics Computers and Structures,vol.81,pp.75–89.

Gordnier R.;Visbal M.(2004): Computation of the aeroelastic response of a flexible delta wing at high angles of attack.Journal of Fluids and Structures,vol.19,pp.785–800.

Goura G.;Badcock K.(2001): A data exchange method for Fluid-Structure Interaction problems.The Aeronautical Journal,vol.4,pp.215–221.

Hall K.;Thomas J.;Clark W.(2002):Computation of unsteady nonlinear flows in cascades using a harmonic balance technique.AIAA Journal,vol.40,pp.879–886.

Kamakoti R.;Shyy W.(2004):Fluid-Cstructure interaction for aeroelastic applications.Progress in Aerospace Sciences,vol.40,pp.535–558.

Kim Y.;Kim J.(2005):A new hybrid interpolation method using surface tracking, fitting and smoothing function applied for aeroelasticity.AIAA 2005-2347.

Lai K.;Tsai H.;Lum K.(2002):CFD and CSD interaction algorithm for large and complex con fi gurations.AIAA 2002-2715.

Lai K.;Tsai H.;Liu F.(2003):Application of spline matrix for mesh deformation with dynamic multi-block grids.AIAA 2003-3514.

Liu L.;Dowell E.;Thomas J.(2007): A high dimensional harmonic balance approach for an aeroelastic airfoil with cubic restoring forces.Journal of Fluids and Structures,vol.23,no.7,pp.351–363.

Mei C.(1977): A Finite-Element Approach for Nonlinear Panel Flutter.AIAA Journal,vol.15,pp.1107–1110.

Pacoste C.(1998):Co-rotational fl at facet triangular elements for shell instability analysis.Computer methods in applied mechanics and engineering,vol.196,pp.75–110.

Patil M.;Hodges H.(2000):On the importance of aerodynamic and structural geometrical nonlinearities in aeroelastic behavior of high-aspect-ration Wings.AIAA 2000-1448.

Patil M.;Hodges D.;Cesnik C.(2001): Nonlinear aeroelasticity and flight dynamics of high-altitude long-endurance aircraft.Journal of Aircraft,vol.38,no.1,pp.88–94.

Relvas A.;Suleman A.(2007): Application of the corotational structural kinematics and Euler fl ow to two-dimensional nonlinear aeroelasticity.Computers and Structures,vol.85,no.1,pp.1372–1381.

Relvas A.;Suleman A.(2006):Fluid-Cstructure interaction modeling of nonlinear aeroelastic structures using the finite element corotational theory.Journal of Fluids and Structures,vol.22,pp.59–75.

Sadeghi M.;Liu F.;Lai L.(2004):Application of three-dimensional interfaces for data transfer in aeroelastic computations.AIAA 2004-5376.

Strganac T.;Cizmas P.;Nichkawde C.(2005): Aeroelastic analysis for future air vehicle concepts using a fully nonlinear methodology.AIAA 2005-2171.

Wong A.;Tsai H.(2000):Unsteady flow calculations with a multi-block moving mesh algorithm,AIAA 00-1002.

Yao W.G.;Xu M.(2008):Modified AUSMpw+Scheme and its Application.The 7th International Conference on System Simulation and Scientific Computing.

Zhu H.;Cai Y.;Paik J.;Atluri S.(2010):Locking-free thick-thin rod/beam element based on a von Karman type nonlinear theory in rotated reference frames for large deformation analyses of space-frame structures.CMES-Computer Modeling in Engineering and Sciences,vol.57,no.2,pp.175–204.