Simulation of Free Airfoil Vibrations in Incompressible Viscous Flow — Comparison of FEM and FVM

This paper deals with a numerical solution of the interaction of two-dimensional (2-D) incompressible viscous flow and a vibrating profile NACA 0012 with large amplitudes. The laminar flow is described by the Navier-Stokes equations in the arbitrary Lagrangian-Eulerian form. The profile with two degrees of freedom (2-DOF) can rotate around its elastic axis and oscillate in the vertical direction. Its motion is described by a nonlinear system of two ordinary differential equations. Deformations of the computational domain due to the profile motion are treated by the arbitrary Lagrangian-Eulerian method. The finite volume method and the finite element method are applied, and the numerical results are compared.


Introduction
Coupled problems describing the interactions of fluid flow with an elastic structure are of great importance in many engineering applications [23,9,1,24].Commercial codes are widely used (e.g.NASTRAN or ANSYS) in technical practice, and aeroelasticity computations are performed mostly in the linear domain, which enables the stability of the system to be determined.
Recently, the research has also focused on numerical modeling of nonlinear coupled problems, because nonlinear phenomena in post-critical states with large vibration amplitudes cannot be captured within linear analysis (see, e.g., [28]).Nonlinear phenomena are important namely for cases, when the structure loses its aeroelastic stability due to flutter or divergence.A nonlinear approach allows the character of the flutter boundary to be determined.The dynamic behavior of the structure at the stability boundary can be either acceptable, when the vibration amplitudes are moderate, or catastrophic, when the amplitudes increase in time above the safety limits.The terminology of benign or catastrophic instability is synonymous with that of stable and unstable limit cycle oscillation (LCO), [9], also referred to as supercritical and subcritical Hopf bifurcation [24].
The nonlinear aeroelasticity of airfoils is very closely related to flutter controlling methods.Librescu and Marzocca in [21] published a review of control techniques and presented in [22] the aeroelastic response of a 2-DOF airfoil in 2-D incompressible flow to external time-dependent excitations, and the flutter instability of actively controlled airfoils.Flutter boundaries and LCO of aeroelastic systems with structural nonlinearities for a 2-DOF airfoil in 2-D incompressible flow was studied by Jones et al. [18].The limit-cycle excitation mechanism was investigated for a similar aeroelastic system with structural nonlinearities by Dessi and Mastroddi in [7].Recently, Chung et al. in [4] analyzed LCO for a 2-DOF airfoil with hysteresis nonlinearity.
Here, two well-known numerical methods, the finite volume method (FVM, see [15,16]) and the finite element method (FEM, see [28]), were employed for a numerical solution of the interaction of 2-D incompressible viscous flow and the elastically supported profile NACA 0012.The mathematical model of the flow is represented by the incompressible Navier-Stokes equations.The profile motion is described in Section 2 by two nonlinear ordinary differential equations (ODEs) for rotation around an elastic axis and oscillation in the vertical direction.
The application of the FVM numerical scheme is described in Section 3. The dual-time stepping method is applied to the numerical solution of unsteady simulations and the Runge-Kutta method is used to solve ODEs numerically.The Arbitrary Lagrangian-Eulerian (ALE) method is employed to cope with strong distortions of the computational domain due to profile motion.The numerical scheme used for khh (u ,v ) T unsteady flow calculations satisfies the geometric conservation law (GCL), cf [19].The application of FEM is described in Section 4. Section 5 presents numerical results for both methods.

Equations of profile motion
A vibrating profile with 2-DOF can oscillate in the vertical direction and rotate around the elastic axis EA (see Figure 1).The motion is described by the following system of two nonlinear ordinary differential equations [17]: where h is vertical displacement of the elastic axis where d [m] is the airfoil depth, ρ [kg m −3 ] is the constant fluid density, U ∞ [m s −1 ] denotes the magnitude of the far field velocity, c [m] denotes the airfoil chord, n = (n 1 , n 2 ) is the unit inner normal to the profile surface Γ Wt (here, the concept of non-dimensional variables was used, cf.[11]).Furthermore, vector r ⊥ is given by r airfoil chord c), Γ Wt denotes the moving part of the boundary (i.e., the surface of the airfoil), and (x EA , y EA ) are the non-dimensional coordinates of the elastic axis (see Figure 2).The components of the (non-dimensional) stress tensor σ ij are given by where u = (u, v) is the non-dimensional fluid velocity vector, p is the non-dimensional pressure, Re is the Reynolds number defined as and ν [m 2 s −1 ] is the fluid kinematic viscosity.Equations ( 2), (3) together with the boundary conditions for the velocity on the moving part Γ Wt of the boundary represent the coupling of the fluid with the structure.
System of equations ( 1) is supplemented by the initial conditions prescribing the values h(0), ϕ(0), ḣ(0), φ(0).Furthermore, it is transformed to the system of first-order ordinary differential equations and solved numerically by the fourth-order multistage Runge-Kutta method.

Mathematical model of the flow
The flow of viscous incompressible fluid in the computational domain Ω t is described by the twodimensional Navier-Stokes equations written in the conservative form where W is the vector of conservative variables are the inviscid physical fluxes defined by and are the viscous physical fluxes defined by where the partial derivatives with respect to x, y and t are denoted by the subscripts x, y and t , respectively.Symbol D denotes the diagonal matrix D = diag(0, 1, 1), non-dimensional time is denoted by t, and the non-dimensional space coordinates are denoted by x, y.Symbols u = (u, v) T and p stand for the dimensionless velocity vector and the dimensionless pressure, respectively.The detailed relation of the dimensionless form of the governing equations and the relationship between dimensional and nondimensional variables can be found, e.g., in [11], [15].In order to time discretize (4), a partition 0 = is considered, and the (variable) time step ∆t n = t n+1 − t n is denoted by ∆t n .Further, W n denotes the approximation of the solution W at the time instant t n .Eqs. ( 4) are discretized in time, and the solution of the non-linear problem for each time step is performed with the use of the concept of dual time, cf.[12].We start with an approximation of the time derivative W t by the second order backward difference formula (BDF2) where ).The artificial compressibility method, see [3,16], together with a timemarching method is used for the steady state computations.This means that the derivatives D β W τ with respect to a fictitious dual time τ are added to the time discretized equations ( 4) with the matrix D β = diag(β, 1, 1), and β > 0 denotes the artificial compressibility constant.The solution of the nonlinear problem is then found as a solution of the steady-state problem in dual time τ , i.e.
where R(W ) is the unsteady residuum defined by and R(W ) is the steady residuum defined by The required real-time accurate solution at the time level n + 1 satisfies R(W n+1 ) = 0, and it is found by marching equation ( 6) to a steady state in the dual time τ .System (1) is equipped with boundary conditions on ∂Ω t .The Dirichlet boundary condition u = u ∞ = (U ∞ , 0) is prescribed on the inlet Γ D , whereas on the outlet Γ O a mean value of pressure p is specified.The non-slip boundary condition is used on walls, i.e. u = w on Γ Wt is prescribed, where w denotes the velocity of the moving wall.

Numerical scheme
The governing equations for the dual-time approach are represented by Eq. ( 6).The arbitrary Lagrangian-Eulerian method, satisfying the so-called geometrical conservation law, is used for application in the case of moving meshes see [19].The integral form of equation ( 4) in the ALE formulation is given by where T is the domain velocity, and C i = C i (t) is a control volume, which moves in time with the domain velocity w, see [19].
In what follows let us assume that Ω t is an polygonal approximation of the domain occupied by a fluid at time t, being discretized by a mesh consisting of a finite number of cells C j (t) satisfying Let us denote by N (i) the set of all neighbouring cells C j (t) of cell C i (t), i.e. the set of all cells whose intersection with C i (t) is a common part of their boundary denoted by Γ ij (t).Only the quadrilateral cells are considered in the computations.In this case, set Γ ij (t) is the common side of cells C i (t) and C j (t).The mean value of vector W over cell C i (t) at time instant t n is approximated by W n i , and symbol W n denotes the vector formed by the collection of all W n i , i ∈ J.In addition, the numerical approximation depends on the ALE mapping A t and on the domain velocity w(t), see also [19].In order to show this dependence we shall denote at time instant t n the vector of grid point positions by X n and the volume of cell C i (t) by Thus the resulting equation for the i-th finite volume cell reads with where F c i (W, X, w) and F v i (W, X) represent the numerical approximations of the convective and diffusive fluxes, respectively.
Further, the BDF2 formula is used for the approximation of the time derivative in (7) by the difference In order to solve the arising non-linear problem, an additional term related to the artificial time τ is added to the scheme, so that to find the solution on time level t n+1 the following problem written in the ALE formulation needs to be solved where and ) stands for the artificial dissipation term (see [15]).The approximation of the convective fluxes F c,n±1/2 i is given by where , ∆y m ij are the x-and y-coordinate differences of the segment Γ m ij .The approximation of the fluxes given by R i (W n+1 ) guarantees the satisfaction of GCL, see [19].The approximation of the viscous fluxes is evaluated on the mesh configuration at the time instant t n+1 and reads where are the approximations of the first partial derivatives of u and v with respect to x and y evaluated by the use of dual finite volume cells, see [15].
Finally, to find the steady state solution of (9) in the dual-time τ the four-stage Runge-Kutta scheme is used.Denoting 9) can be written in the form: and solved with the aid of Runge-Kutta method, cf.[15].

Finite element method 4.1 Mathematical model of the flow
For the application of the finite element method, the flow was described by the incompressible Navier-Stokes equations written in the ALE form where D A Dt is the ALE derivative, u denotes the velocity vector, p denotes the nondimensional pressure, w is the domain velocity known from the ALE method, and Re denotes the Reynolds number.Furthermore, equation ( 10) is equipped with the boundary conditions where n denotes the unit outward normal vector, u D is the Dirichlet boundary condition, and (α) − denotes the negative part of a real number α, see also [2], [14].Further, system (10) is equipped with the initial condition u(x, 0) = u 0 (x).

Finite element approximation
FEM is well known as a general discretization method for partial differential equations.However, straightforward application of FEM procedures often fails in the case of incompressible Navier-Stokes equations.The reason is that momentum equations are of advectiondiffusion type, with the dominating advection , the Galerkin FEM leads to unphysical solutions if the grid is not fine enough in regions of strong gradients (e.g. in the boundary layer).In order to obtain physically admissible correct solutions it is necessary to apply suitable mesh refinement (e.g. an anisotropically refined mesh, cf.[8]) combined with a stabilization technique, cf.[5], [27].In this work, FEM is stabilized with the aid of the streamline upwind/pressure stabilizing Petrov-Galerkin (SUPG/PSPG) method (the so called fully stabilized scheme, cf.[13]) modified for the application on moving domains (cf. [28]).
In order to discretize the problem (10), we approximate the time derivative by the second order backward difference formula, where u n+1 is the flow velocity at time t n+1 defined on the computational domain Ω n+1 , and u k is the transformation of the flow velocity at time t k defined on Ω k transformed onto Ω n+1 .Further, equation ( 10) is formulated weakly, and the solution is sought on a couple of finite element spaces W ∆ ⊂ H 1 (Ω n+1 ) and ) for approximation of the velocity components and pressure, respectively, and the subspace of the test functions is denoted by Let us mention that the finite element spaces should satisfy the Babuška-Brezzi (BB) condition (see e.g.[25], [26] or [29]).In practical computations we assume that the domain Ω = Ω n+1 is a polygonal approximation of the region occupied by the fluid at time t n+1 and the finite element spaces are defined over a triangulation T ∆ of domain Ω t as piecewise polynomial functions.In our computations, the wellknown Taylor-Hood P 2 /P 1 conforming finite element spaces are used for the velocity/pressure approximation.This means that the pressure approximation p = p ∆ and the velocity approximation u = u ∆ are linear and quadratic vector-valued functions on each element K ∈ T ∆ , respectively.The stabilized discrete problem at the time instant t = t n+1 reads: Find U = (u, p) ∈ W ∆ × Q ∆ , p := p n+1 , u := u n+1 , such that u satisfies approximately the conditions (11 a-b) and holds for all V = (v, q) ∈ X ∆ × Q ∆ .Here, the Galerkin terms are defined for any U = (u, p), V = (v, q), U * = (u * , p * ) by where w = u * − w n+1 , and the scalar product in L 2 (Ω) is denoted by (•, •) Ω .Further, the SUPG/PSPG stabilization terms are used in order to obtain a stable solution also for large Reynolds number values, where ψ(V ) = (w • ∇)v + ∇q, w stands for the transport velocity at time instant t n+1 , i.e. w = v * −w n+1 , and (•, •) K denotes the scalar product in L 2 (K).The term P(U, V ) is the additional grad-div stabilization defined by Here, the choice of the parameters δ K and τ K is carried out according to [13] or [27] on the basis of the local element length Furthermore, problem ( 12) is solved with the aid of Oseen linearizations, and the arising large system of linear equations is solved with the aid of a direct solver, e.g.UMFPACK (cf.[6]), where different stabilization procedures can easily be applied, even when anisotropically refined grids are employed.
The motion equations describing the motion of the flexibly supported airfoil are discretized with the aid of the 4th order Runge-Kutta method, and the coupled fluid-structure model is solved with the aid of a partitioned strongly coupled scheme.This means that for every time step the fluid flow and the structure motion are approximated repeatedly in order to converge to a solution which will satisfy all interface conditions. to 8 × 10 5 ).The initial conditions h(0) = −20 mm, ḣ(0) = 0, ϕ(0) = 6°, φ(0) = 0 were used.The vibrations are evidently damped by aerodynamic forces, and the decay of the oscillations is faster with the fluid velocity up to about U ∞ = 35 m/s (see Figures 3 -5).Further, the system loses its static stability at about U ∞ = 40 m/s by divergence.More severe instability occurs in the FEM simulation, see Figure 6, where the rotation angle reaches nearly 14°(at that moment the computation failed due to mesh distortion).The aeroelastic response for this case computed by FVM rather limits to smaller static deflections (see Figures 5 and 6).However, the FVM results for upstream velocity U ∞ = 42.5 m/s shown in Figure 7 apparently demonstrate unstable behaviour of the system.Here, the vibration amplitudes increase rapidly and reach values of about 9°for rotation and −70 mm for vertical displacement just before the computation collapses.Note that the critical flow velocity for the aeroelastic instability computed by NASTRAN was 37.7 m/s (see [28]).This corresponds to the results of the numerical simulations presented since the system was found stable for upstream flow velocity 35 m/s and unstable for upstream flow velocity 40 m/s in the case of FEM and 42.5 m/s in the case of FVM.

Numerical flow induced vibration results for both
Further, the Fourier transformations of rotation angle ϕ(t) and vertical displacement h(t) signals were computed.Figures 8 and 9 compare the frequency spectrum analysis of FVM and FEM results for the far-field flow velocities U ∞ = 5, 10, 15 and 20 m/s, respectively.The two dominant frequencies can be seen for all the considered far field velocities.The lower frequency f 1 ≈ 5.3 Hz refers to the vertical translation of the profile, and the higher frequency f 2 ≈ 13.6 Hz refers to the profile rotation around the elastic axis.When the flow velocity is increased, both resonances are more damped and the frequencies are getting closer (f 1 ≈ 5.5 Hz, f 2 ≈ 12.8 Hz), which is a typical phenomenon in flutter analyses.The spectra show very good agreement of the dominant frequencies between the FVM and FEM results for these flow velocities.The result of the Fourier transformations shown in Fig. 10 for U ∞ = 25 m/s is different, and it is difficult to identify the dominant frequencies precisely.This is above all due to much higher aerodynamical damping, which is also similar for the cases U ∞ = 35-45 m/s (the results are not shown here).

Conclusions
Aeroelastic instability due to divergence appeared prior to flutter in both numerical approaches.The results of the numerical simulations in the time domain computed by FVM are in good qualitative agreement with the FEM computations of the same aeroelastic problem, and both numerical methods estimate the critical flow velocity in good agreement with the NAS-TRAN computation using a classic linear approach.
The numerical results presented here, computed by FVM and FEM, are in a very good quantitative agreement for upstream velocities up to 35 m/s.Here, both computations fully and almost identically represent the aeroelastic behaviour of the system.FVM and FEM results differ markedly when the upstream flow velocity reaches 40 m/s.The system is unstable for upstream velocity close to this value.The artificial numerical dissipation implemented in the FVM scheme may be responsible for the excessively strong influence of the flow, and consequently for the behaviour of the system for upstream velocities close to 40 m/s.
The numerical scheme in FVM is implicit in real time but explicit in dual time, while the FEM method is fully implicit.This implies that the computations in FVM are more time consuming than in FEM.In the computations presented here, the mesh for FVM consisted of approximately 18 thousand of cells and 54 thousand unknowns, and the computations took approximately 50-80 seconds for a single time step (depending on the far field velocity), on a PC computer with an Intel Quad Core processor with 4 GB of memory.The computation of FEM was performed on a mesh with approximately 34 thousand of elements with approximately 150 thousand unknowns, and the computations took approximately 30-50 seconds for a single time step.In both cases the time consumed is also influenced by the number of iterations of the coupling algorithm between fluid and structure.However, the implemented numerical scheme for FVM does not require any big matrices to be stored, unlike the fully implicit numerical scheme in FEM.Hence, the computer memory requirements are much more moderate for FVM than for FEM (approximately 1.2 GB of memory).

Figure 1 :
Figure 1: Elastic support of the profile on translational and rotational springs.
(downwards positive) [m], ϕ is rotation angle around the elastic axis (clockwise positive) [rad], m is the mass of the profile [kg], S ϕ is the static moment around the elastic axis [kg m], k h is the bending stiffness [N/m], I ϕ is the inertia moment around the elastic axis [kg m 2 ], and k ϕ is the torsional stiffness [N m/rad].The aerodynamic lift force L [N] acting in the vertical direction (upwards positive) and the torsional moment M [N m] (clockwise positive) are defined as

Figure 10 :
Figure 10: Frequency spectrum analysis of vertical displacement h [mm] and angle of rotation ϕ [°] signals for U ∞ = 25 m/s.Solid line -FVM results, dashed line -FEM results.