Numerical analysis of flow-induced nonlinear vibrations of an airfoil with three degrees of freedom

Computers & Fluids 49 (2011) 110–127

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Numerical analysis of flow-induced nonlinear vibrations of an airfoil with three degrees of freedom Miloslav Feistauer a,⇑, Jaromír Horácˇek b, Martin Ru˚zˇicˇka a, Petr Svácˇek c a

Charles University Prague, Faculty of Mathematics and Physics, Sokolovská 83, 186 75 Praha 8, Czech Republic Institute of Thermomechanics, Academy of Sciences of the Czech Republic, Dolejškova 5, 182 00 Praha 8, Czech Republic c Czech Technical University Prague, Faculty of Mechanical Engineering, Karlovo nám. 13, 121 35 Praha 2, Czech Republic b

a r t i c l e

i n f o

Article history: Received 26 April 2010 Received in revised form 5 May 2011 Accepted 6 May 2011 Available online 18 May 2011 Keywords: Aeroelasticity Navier–Stokes equations Arbitrary Lagrangian–Eulerian formulation Finite element method Stabilization for high Reynolds numbers Non-linear oscillations Flutter instability

a b s t r a c t The subject of the paper is the numerical simulation of the interaction of two-dimensional incompressible viscous flow and a vibrating airfoil with large amplitudes. A solid airfoil with three degrees of freedom performs rotation around an elastic axis, oscillations in the vertical direction and rotation of a flap. The numerical simulation consists of the finite element solution of the Navier–Stokes equations coupled with a system of nonlinear ordinary differential equations describing the airfoil motion. The time-dependent computational domain and a moving grid are treated by the arbitrary Lagrangian–Eulerian (ALE) formulation of the Navier–Stokes equations. High Reynolds numbers require the application of a suitable stabilization of the finite element discretization. With the aid of numerical experiments we analyze the convergence of the method, if the computational mesh is refined and the time step is successively decreased. Also the effect of the loose (weak) and strong coupling of the flow and structure problems is tested. The applicability of the method is demonstrated by the comparison with NASTRAN computations and the solution of the coupled fluid–structure problem with increasing far-field velocity up to and behind the flutter instability. The developed method can be used in theoretical prediction of dynamic behaviour of aeroelastic systems especially when the system stability has been lost. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction The interaction of fluids and structures plays an important role in a number of fields. Important applications can be found in technology (design of engineering systems, as aircrafts, blade machines, bridges and other structures in civil or mechanical engineering) and in biomechanics (e.g. elastic artery modelling for stent design, simulation of the flow in blood vessels and studying the vocal folds vibrations and vocal sound production). The interaction between flowing fluids and vibrating structures is the main subject of aeroelasticity, which studies the influence of aerodynamic and elastic forces on an elastic structure. The flow-induced vibrations may affect negatively the operation and stability of the systems. Therefore, one of the main goals of aeroelasticity is the prediction and cure of the aeroelastic instability. This discipline achieved many results, particularly from engineering point of view (see, e.g. the monographs [3,13,28]). From the point of view of pure mathematical theory of fluid– structure interaction, there are not too many works, due to a high ⇑ Corresponding author. E-mail addresses: [email protected] (M. Feistauer), [email protected] (J. Horácˇek), [email protected] (M. Ru˚zˇicˇka), [email protected] (P. Svácˇek). 0045-7930/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2011.05.004

mathematical complexity of the problem, caused by the timedependence of the domain occupied by the fluid and coupling of systems describing flow and elastic structures. We can mention, e.g. the work Neustupa [29] devoted to the solvability of the Navier–Stokes system in time-dependent domains and the works Granmond [17], Guidorzi et al. [18] and Hoffmann and Starovoitov [19] dealing with the interaction of viscous incompressible flow and elastic plates or solid bodies. The subject of our attention is the numerical analysis of the interaction of viscous flow with a vibrating airfoil. In many aerospace engineering problems, it is important to simulate flow induced airfoil vibrations in order to predict the bounds of stability, to determine the rise of instabilities leading to flutter and to analyze postcritical regimes. The fluid–structure interaction can lead to unacceptable vibration level due to the self-induced oscillations, and to fatigue failures [3]. The linearized aeroelastic equations allow the determination of the flutter boundary, whereas the nonlinear approach allows the determination of the character of the flutter boundary. The dynamic behavior of the structure at the flutter boundary can be either acceptable, when sustained vibration amplitudes are moderate, or catastrophic, when the amplitudes are increasing in time above safety limits (see [4]). The terminology of benign or catastrophic flutter is synonymous with that of stable and unstable

111

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

limit cycle oscillations (LCO) (see [12]), in Paidoussis [31,32] also referred as supercritical and subcritical Hopf bifurcation. The nonlinearities in the aeroelastic model can yield a moderate vibration level at the flutter boundary or they can contribute to the catastrophic type of flutter boundary. Conversion of the catastrophic type of flutter boundary into a benign requires consideration of the nonlinear aeroelastic equations and the inclusion of a feedback control as shown in Librescu and Marzocca [22]. Many papers published on flutter controlling methods are very closely related to the nonlinear aeroelasticity of airfoils studied here. Suppression of the occurrence of the flutter instability has a great importance. For example, Librescu et al. [23] considered the subcritical aeroelastic response of two degrees of freefom (2-DOF) airfoil in 2-D incompressible flow to external time-dependent excitations, and the flutter instability of actively controlled airfoils. Analytical analysis and numerical simulations of the aeroelastic response of 3-DOF airfoil–flap system subjected to time-dependent loads in an incompressible flow were addressed by Shim et al. [35] focusing on feedback control in order to suppress the flutter instability. The aerodynamic forces are derived from the Theodorsen equations. In Na et al. [27] the aeroelastic response and control of 3-DOF flapped-wing system in an incompressible fluid flow and exposed to external pressure pulse was studied. Structural nonlinearities can have significant effects on the aeroelastic response. Recently, Chung et al. [8] analyzed of LCO for a 2-DOF airfoil motion containing a hysteresis structural nonlinearity. The mechanism of limit cycle excitation is investigated for an aeroelastic system with structural nonlinearities by Dessi and Mastroddi [11]. Flutter boundaries and LCO with influence of structural nonlinearities like freeplay and hysteresis in stiffness was studied by Jones et al. [21]. Recently, a history of 2-D unsteady incompressible airfoil theory using potential flow model has been overviewed by Peters [33]. In fluid–structure interaction problems, input parameters can suffer from randomness. This subject is analyzed, for example, in Witteveen and Bijl [39], where an adaptive stochastic finite element approach for rigid-body fluid–structure interaction is developed. In all cited papers, only small amplitudes of vibrations are considered and no effects of large rotation amplitudes and nonlinear mass matrix are taken into account. An important item in the area of fluid–structure interaction is the state of art of the numerical simulation of coupled fluid–structure problems. There are several works concerned with a general framework of strategies applied to this subject. In Badia and Codina [1], an iterative coupling based on the domain decomposition framework is described. In [2], partitioned procedures for fluid– structure interaction coupling based on the Robin-type transmission conditions are designed. The paper Farhat et al. [14] presents a methodology for designing formally second-order time accurate loosely coupled fluid–structure interaction procedures. Various aspects of coupling fluid–structure algorithms are treated by Fernandez and Gerbeau [15], particularly with respect to modelling and simulation of cardiovascular systems. In our paper we are concerned with a two-dimensional viscous incompressible flow past a moving airfoil, which is considered as a solid flexibly supported body with three degrees of freedom, allowing its vertical and torsional oscillations and the rotation of a flap. The airfoil vibrations are described by a system of second-order nonlinear ordinary differential equations for the vertical displacement of the airfoil, the rotation angle of the whole airfoil around its elastic axis and the rotation angle of the flap around its elastic axis. This system allows the simulation of airfoil vibrations with large amplitudes and includes also the torsional stiffness modelled by nonlinear cubic spring terms. The numerical simulation consists of the stabilized finite element solution of the Navier–Stokes equa-

tions written in the ALE form coupled with the system of equations describing the airfoil motion and solved by the Runge–Kutta method. We give a detailed description of all ingredients of the coupled fluid–structure interaction problem. Namely, the time and space discretization including the finite element stabilization allowing the solution of flow problem with high Reynolds numbers is treated, the construction of the ALE mapping is described and the coupling algorithm involving the loose and strong coupling is formulated. The developed technique is analyzed and tested on the basis of numerical experiments. We are concerned with the convergence analysis and results sensitivity with respect to the the time and mesh size. Moreover, the accuracy of the method is tested for loosely coupled and strongly coupled fluid–structure strategies. The method is also tested on a problem for which the results computed by the MSC.NASTRAN program code are available. The comparison of our computations and the NASTRAN results shows good agreement. Numerical experiments demonstrate that the developed method is sufficiently accurate and robust and allows to resolve both subcritical and postcritical regimes of flow-induced airfoil vibrations. 2. Formulation of the problem The two-dimensional non-stationary flow of viscous, incompressible fluid is considered in the time interval [0, T], where T > 0. The symbol Xt will denote the computational domain occupied by the fluid at time t. The flow is characterized by the velocity u = u(x, t) and the kinematic pressure p = p(x, t), for x 2 Xt and t 2 [0, T]. The kinematic pressure is defined as P/q, where P is the pressure and q = const. > 0 is the density of the fluid. By m we denote the kinematic viscosity. The quantites u and p satisfy the Navier–Stokes equations and the continuity equation

@u þ ðu  rÞu  mMu þ rp ¼ 0; @t r  u ¼ 0 in Xt ; t 2 ð0; TÞ:

ð2:1Þ

By Xt and @ Xt we shall denote the closure and the boundary of the domain Xt. We shall assume that @ Xt ¼ CD [ CO [ CW t , where the sets CD, CO and CW t are mutually disjoint and boundary conditions of different types are prescribed there. The part CD represents the inlet, where the fluid flows into the domain Xt, or a fixed, impermeable wall, CO denotes the outlet, where the fluid leaves Xt and CW t is the moving airfoil boundary at time t. It consists of two parts – profile Pt and flap Ft: CW t ¼ P t [ F t . We assume that CD and CO are independent of time in contrast to CW t . The motion of the airfoil is described by functions a(t), b(t) and h(t), representing the rotation angle of the whole airfoil around an elastic axis EA, the rotation angle of the flap around an elastic axis EF and the vertical displacement of the whole airfoil, respectively. The shape of the domain Xt depends on the functions a(t), b(t) and h(t). We consider two situations, namely the flap attached to the main body of the airfoil without any gap and the flap separated from the main body of the airfoil by a narrow gap with width g. See Fig. 1. 2.1. ALE formulation of the Navier–Stokes equations In order to simulate flow in a moving domain Xt, we employ the arbitrary Lagrangian–Eulerian (ALE) method (cf. [30]), based on a regular one-to-one ALE mapping

At : X0 # Xt ;

X 2 X0 # xðX; tÞ ¼ At ðXÞ 2 Xt ;

t 2 ½0; T:

ð2:2Þ

112

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

EA

EA

EF

h

h

EF

EA

EF

EF

EA g

Fig. 1. Model scheme – airfoil with three degrees of freedom without a gap (left) and with a gap (right).

The reference domain X0 is identical with the domain occupied by the fluid at the initial time t = 0. The coordinates of points x 2 Xt are called spatial coordinates, the coordinates of points X 2 X0 are called ALE coordinates or reference coordinates. First, we define the domain velocity

~ wðX; tÞ ¼

@ xðX; tÞ: @t

ð2:3Þ

This velocity can be expressed in the spatial coordinates as

  ~ A1 wðx; tÞ ¼ w t ðxÞ; t :

ð2:4Þ

Let us consider a function f ¼ f ðx; tÞ; x 2 Xt ; t 2 ½0; T; f ðx; tÞ 2 R, where R is the set of real numbers. Let us set ~f ðX; tÞ ¼ f ðAt ðXÞ; tÞ. We define the ALE derivative of the function f by

D @~f f ðx; tÞ ¼ ðX; tÞ; @t Dt

sional moment of the force acting on the flap with respect to the flap axis EF, Dhh, Daa, Dbb are the coefficients of a structural damping, Sa, Ia and m denote the static moment of the whole airfoil around the elastic axis EA, the inertia moment of the whole airfoil around the elastic axis EA and the mass of the whole profile, respectively. The coefficient Sb is the static moment of the flap of the airfoil around the flap axis EF and Ib is the moment of inertia of the flap of the airfoil around the flap axis EF. Constants khh, kaa, kbb denote the ^aa ; k ^bb are spring stiffnesses of the flexible support of the airfoil, k the nonlinear torsional spring stiffnesses and dPF is the distance between the elastic axis EA and the flap axis EF. The interaction between the flow and the airfoil is given by the force component L2 and the moments Ma and Mb defined by

L2 ¼ lq

A

X ¼ A1 t ðxÞ:

Z

ð2:5Þ

D @f þ w  rf : f ¼ @t Dt

ð2:6Þ

Using relation (2.6), we obtain the Navier–Stokes equations in the ALE form A

D u þ ððu  wÞ  rÞu þ rp  mDu ¼ 0 Dt div u ¼ 0

in Xt :

ð2:7Þ

2.2. Equations for the moving airfoil The nonlinear equations describing the vibrations with large amplitudes of the airfoil given by the functions a, b and h read

€ þ ½ðS  S Þ cos a þ S cosða þ bÞa € mh a b b € cosða þ bÞ  ðSa  Sb Þa_ 2 sin a þSb b _ 2 sinða þ bÞ þ Dhh h_ þ khh h ¼ L2 ; Sb ða_ þ bÞ € ½ðSa  Sb Þ cos a þ Sb cosða þ bÞh   € þ ðIa  2dPF Sb Þ þ 2dPF Sb cos b a €  dPF Sb b_ 2 sin b þ½Ib þ dPF Sb cos bb

ð2:9Þ

2 X

Mb ¼ lq

Z X 2 F t i;j¼1

T ij nj ð1Þi ðx1þd1i  xEA 1þd1i ÞdS;

T ij nj ð1Þi ðx1þd1i  xEF 1þd1i ÞdS;

ð2:10Þ ð2:11Þ

where l is the depth of the segment of the airfoil, n = (n1, n2) is the outer unit normal to @ Xt on CW t , the symbol dij is the Kronecker symbol defined by dij = 1 for i = j and dij = 0 for i – j, x1 and x2 are the coordinates of the point on CWt ; xEA i ; i ¼ 1; 2, are the coordinates of the current location of the elastic axis EA and xEF i , i = 1, 2, are the coordinates of the current location of the flap elastic axis EF. The components of stress tensor are computed from the flow velocity and the pressure by the relation

  @ui @uj : T ij ¼ pdij þ m þ @xj @xi

ð2:12Þ

2.3. Initial and boundary conditions The Navier–Stokes system (2.7) is completed by the initial condition

ð2:8Þ

^aa a3 ¼ M a ; 2dPF Sb a_ b_ sin b þ Daa a_ þ kaa a þ k € þ ½I þ d S cos ba € S cosða þ bÞh b

Z

Pt [F t i;j¼1

A

b

T 2j nj dS;

Pt [F t j¼1

Ma ¼ lq

The application of the chain rule gives

2 X

PF b

^bb b3 ¼ Mb : € þ dPF Sb a_ 2 sin b þ Dbb b_ þ kbb b þ k þIb b For the derivation, see Horácˇek et al. [20]. The symbol L2 denotes the lift, i.e. the component of the force acting on the whole airfoil in the vertical direction, Ma is the torsional moment of the force acting on the whole airfoil with respect to the axis EA, Mb is the tor-

uðx; 0Þ ¼ u0 ðxÞ;

x 2 X0 ;

ð2:13Þ

and the following boundary conditions. On CD we prescribe the Dirichlet condition

ujCD ¼ uD :

ð2:14Þ

On the outlet CO we consider the so-called do-nothing boundary condition

Þ n þ m ðp  p

@u ¼ 0 on CO ; @n

ð2:15Þ

113

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

 is a reference value of the pressure and n denotes the unit where p outer normal to @ Xt. On CW t we consider the condition

ujCW ¼ wjCW : t

ð2:16Þ

t

að0Þ ¼ a0 ; a_ ð0Þ ¼ a1 ; hð0Þ ¼ h0 ;

W ¼ ½H1 ðXÞ2 ;

_ bð0Þ ¼ b1 ; _hð0Þ ¼ h ;

ð2:17Þ

3. Discretization of the flow problem

X ¼ fv 2 W; v jCD [CW ¼ 0g

ð3:4Þ

and the pressure space

We start from the time semidiscretization. We construct an equidistant partition of the time interval [0, T], formed by time instants 0 = t0 < t1 <    < T, tk = ks, where s > 0 is a time step. We use the approximation u(tn)  un, p(tn)  pn of the exact solution and w(tn)  wn of the domain velocity at time tn. First, we shall be concerned with the ALE derivative DAu/Dt. If X 2 X0 is a given point from the reference configuration, then for a given time level tn we can write Atn1 ðXÞ ¼ xn1 2 Xtn1 ; Atn ðXÞ ¼ xn 2 Xtn ; Atnþ1 ðXÞ ¼ xnþ1 2 Xtnþ1 . Using definition (2.5), we can approximate the ALE derivative of the velocity at time tn+1 and point xn+1 by the second-order two-step backward difference formula

~ nþ1 ðXÞ  4u ~ n ðXÞ þ u ~ n1 ðXÞ DA u 3u ðxnþ1 ; t nþ1 Þ  Dt 2s 3unþ1 ðxnþ1 Þ  4un ðxn Þ þ un1 ðxn1 Þ ¼ : 2s

ð3:1Þ

If we take into account that Atnþ1 ðA1 t i ðxi ÞÞ 2 Xtnþ1 and introduce the ^ i ¼ ui  Ati  A1 functions u t nþ1 ; i ¼ n; n  1, obtained by the transformation of ui to the domain Xtnþ1 , we get the implicit scheme for unknown functions unþ1 : Xtnþ1 # R2 and pnþ1 : Xtnþ1 # R:

þ ððunþ1  wnþ1 Þ  rÞunþ1

þ rpnþ1 ¼ 0;

Q ¼ L2 ðXÞ;

ð3:5Þ

where L2(X) is the Lebesgue space of square integrable functions over the domain X and H1(X) is the Sobolev space of functions square integrable together with their first order derivatives. Further, we introduce the notation

ða; bÞx ¼

Z x

abdx; ð; Þ ¼ ð; ÞX

ð3:6Þ

for the scalar product in L2(x) for some set x and define the forms

3 ðu; v Þ þ mðru; rv Þ þ ðððu  wnþ1 Þ  rÞu; v Þ 2s  ðp; r  v Þ þ ðr  u; qÞ; ð3:7Þ Z 1 ^n  u ^ n1 ; v Þ  v  ndS; ð4u p f ðVÞ ¼ 2s CO

aðU  ; U; VÞ ¼

U ¼ ðu; pÞ; V ¼ ðv ; qÞ; U  ¼ ðu ; p Þ:

3.1. Time discretization

mDu

ð3:3Þ

1

where a0, a1, b0, b1, h0, h1 are input parameters of the model. The interaction of a fluid and an airfoil consists in the solution of the flow problem (2.7), (2.13)–(2.16) coupled with the structural model (2.8) and (2.17). This coupled problem represents a strongly nonlinear dynamical system. In what follows we shall be concerned with the discretization of the problem and describe the algorithm for the numerical solution of the complete fluid–structure interaction problem. We proceed in such a way that first we describe separately the discretization of the flow problem and the structural problem in Sections 3 and 4, respectively. Then, in Section 5, the realization of the complete coupled problem will be treated.

^ n þu ^ n1 3unþ1 4u 2s nþ1

in X;

considered with the boundary conditions (2.14)–(2.16). We define the velocity function spaces

Moreover, we equip system (2.8) with the initial conditions

bð0Þ ¼ b0 ;

^n þ u ^ n1 3u  4u þ ððu  wÞ  rÞu  mDu þ rp ¼ 0; 2s r  u ¼ 0;

in Xtnþ1 :

ð3:2Þ

r  unþ1 ¼ 0; This system is equipped with the boundary conditions (2.14)– (2.16). 3.2. Space discretization of the flow problem For the space discretization we shall use the finite element method (FEM), which appears suitable for the solution of the flow problem because of its flexibility and accuracy. First, we reformulate the continuous problem (3.2) in a weak sense, which is a starting point for the application of the finite element method. Because of simplicity we shall use the notation X ¼ Xtnþ1 ; u ¼ unþ1 ; p ¼ pnþ1 ; w ¼ wnþ1 ; CW ¼ CW tnþ1 and write system (3.2) in the form

Now, if we multiply the first and second equation in (3.3) by arbitrary functions v 2 X and q 2 Q, respectively, sum them, integrate over X, transform the terms containing Du and rp with the aid of Green’s theorem and use the boundary condition (2.15), we arrive at the concept of a weak solution of problem (3.3), (2.14)–(2.16) as a couple U = (u, p) 2 X  Q satisfying the identity

aðU; U; VÞ ¼ f ðVÞ;

forall V ¼ ðv ; qÞ 2 X  Q ;

ð3:8Þ

and conditions (2.14) and (2.16). In order to apply the finite element method to the numerical solution, we assume that the domain X is a polygonal approximation of the computational domain at time tn+1. By T D we denote a triangulation of X formed by a finite number of closed triangles. The parameter D denotes the maximal size of elements K 2 T D . We assume that any two different triangles are either disjoint or intersect each other in a common face or in a common vertex (cf., e.g. [7]). We use the Taylor–Hood P2/P1 elements [37]. To this end we introduce the finite-dimensional spaces

Q D ¼ fq 2 CðXÞ; qjK 2 P1K ; 8K 2 T D g; W D ¼ fv 2 ½CðXÞ2 ; v jK 2 ½P2K 2 ; 8K 2 T D g;

ð3:9Þ

X D ¼ fv 2 W D ; v jCD [ CW ¼ 0g: Here the symbol P kK denotes the space of all polynomials on K of degree 6k. The couple (XD, QD) satisfies the Babuška–Brezzi condition (see, e.g. [5, 6, 38]), which is important for the stability of the finite element scheme. The domain velocity w at time tn+1 is approximated by a func^i ^i tion wD ¼ wnþ1 D , we use the approximations u  uD ; i ¼ n; n  1, ⁄ and approximate the forms a(U , U, V) and f(V) by

3 ðu; v Þ þ mðru; rv Þ þ ðððu  wD Þ  rÞu; v Þ 2s ð3:10Þ  ðp; r  v Þ þ ðr  u; qÞ; Z  1 ^ nD  u ^ n1 v  ndS; 4u fD ðVÞ ¼ p D ;v  2s CO

aD ðU  ; U; VÞ ¼

U ¼ ðu; pÞ; V ¼ ðv ; qÞ; U  ¼ ðu ; p Þ 2 W D  Q D :

114

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

The approximate solution is defined as a couple UD = (uD, pD) 2 WD  QD such that biuD satisfies approximately conditions (2.14) and (2.16) and the identity

aD ðU D ; U D ; V D Þ ¼ fD ðV D Þ for all V D ¼ ðv D ; qD Þ 2 X D  Q D :

ð3:11Þ

 3   rÞu þ rp; ðw   rÞv ; u  mDu þ ðw 2s K K2T D  X 1 ^ nD  u ^ n1  dK ð4u ; ð3:12Þ F D ðVÞ ¼ D Þ; ðw  rÞv 2s K K2T D X PD ðU; VÞ ¼ sK ðr  u; r  v ÞK ;

LD ðU  ; U; VÞ ¼

X



dK

K2T D

3.3. Stabilization of the finite element method For high Reynolds numbers approximate solutions can contain nonphysical spurious oscillations. In order to avoid them, it is necessary to apply a suitable stabilization of the finite element scheme. Here we use the streamline-diffusion and div–div stabilization based on the forms

 ¼ u  wD and dK, sK P 0 where U = (u, p) U⁄ = (u⁄, p) V = (v, q), w are suitable parameters. The solution of the stabilized discrete problem is UD = (uD, pD) 2 WD  QD such that uD satisfies approximately the boundary conditions (2.14) on CD and (2.16) on CW and

aD ðU D ; U D ; V D Þ þ LD ðU D ; U D ; V D Þ þ PD ðU D ; V D Þ

Table 1 Comparison of the eigenfrequences of the system in vacuo

NASTRAN System (2.8)

h-Bending f (Hz)

a-Torsion f (Hz)

b-Torsion f (Hz)

5.55 5.50

15.49 15.46

11.41 11.39

The couple (uD, pD) represents the approximate solution on the time level tn+1 defined in the domain Xtnþ1 . The parameters dK and sK are defined on the basis of results from Gelhard et al. [16] and Lube [26] and our numerical experiments and tests. We put

dK ¼ d Table 2 Comparison of the results computed by NASTRAN ([25]) and by the developed finite element method for the airfoil without gap (FEM0) and with gap (FEM1): (a) eigenfrequencies f for far-field airflow velocity 2 m/s, (b) critical flutter velocities vcr and flutter frequencies fcr.

NASTRAN FEM0 FEM1

h-Bending f (Hz)

a-Torsion

b-Torsion f (Hz)

vcr

f (Hz)

(m/s)

fcr (Hz)

5.49 5.35 5.37

15.44 15.38 15.40

11.41 11.31 11.32

11.32 8 < vcr < 10 11 < vcr < 12

14.87 15.15 14.89

ð3:13Þ

8 V D ¼ ðv D ; qD Þ 2 X D  Q D :

¼ fD ðV D Þ þ F D ðV D Þ;

hK  nðRw K Þ;  L1 ðKÞ 2kwk

ð3:14Þ

 L1 ðKÞ ¼ maxK jwj,  hK is the size of K measured in the direcwhere kwk  and tion of the mean value of w 

Rw K ¼

 L1 ðKÞ hK kwk 2m

;



nðRw K Þ ¼ min

 w  RK ;1 : 6

ð3:15Þ

The parameters sK are defined by

 L1 ðKÞ ; s 2 ð0; 1: sK ¼ s hK kwk

ð3:16Þ

In practical computations we use the values d⁄ = 0.025 and s⁄ = 1.

Fig. 2. Coarse mesh for NACA 0012 airfoil with a gap adapted at the airfoil and its details.

115

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

3.4. Treatment of the nonlinearity in the flow model

6

h [mm]

4

The nonlinear problem (3.13) is (on each time level) solved with the aid of the Oseen iterative process. Starting from an initial approximation U nþ1 D;0 at time tn+1 and assuming that the iteration U nþ1 has already been computed, we define U nþ1 D;m D;mþ1 ¼ ðuD;mþ1 ; pD;mþ1 Þ 2 W D  M D by

2 0 -2 -4 -6

0

0.2

0.4

0.6

0.8

1

nþ1 nþ1 nþ1 nþ1 aD ðU nþ1 D;m ; U D;mþ1 ; V D Þ þ LD ðU D;m ; U D;mþ1 ; V D Þ þ P D ðU D;mþ1 ; V D Þ

4

α [deg]

2 0 -2 -4 0

0.2

0.4

0.6

0.8

1

t [s] 10

β [deg]

5

pD;mþ1 ¼

nQ P j¼1

-5

ð3:17Þ

We obtain a sequence U nþ1 D;m , m = 0, 1, . . ., and assume that it conÞ and verges to the solution U nþ1 of Eq. (3.13). We set U 1D;0 ¼ ðu0D ; p D ^ n ^ n1 n for each time level tn+1, n P 1, we set U nþ1 D;0 ¼ ð2uD  uD ; pD Þ. In order to realize the solution of the Oseen problem (3.17) we X introduce the system of functions fwi gni¼1 as a basis of the space XD  nQ and the system of functions fqj gj¼1 , which forms a basis of the space QD. Let us assume that uD 2 W D is a function satisfying at time tn+1 the boundary conditions (2.14) and (2.16) at the nodes (i.e. the vertices and midpoints of sides of elements K 2 T D ) lying on CD and CW t , respectively, and vanishing at other nodes. The couple U nþ1 D;mþ1 can be expressed in the form

uD;mþ1 ¼ uD þ

0

8V D ¼ ðv D ; qD Þ 2 X D  Q D :

¼ fD ðV D Þ þ F D ðV D Þ;

t [s]

nX P j¼1

U j wj ; U j 2 R; j ¼ 1; . . . ; nX ; ð3:18Þ

Pj qj ; Pj 2 R; j ¼ 1; . . . ; nQ :

Substituting (3.18) into (3.17), we obtain the linear system of nX + nQ equations of the saddle-point type

-10 0

0.2

0.4

0.6

0.8

1

t [s]



A BT

BþC O



U P



 ¼

 F : G

ð3:19Þ

8 6 4 2 0 -2 -4 -6 -8

h [mm]

h [mm]

Fig. 3. Results on three successively refined meshes: airfoil with gap, far-field velocity 10 m/s, ..... – coarse mesh,    – refinement 1, —— – refinement 2.

0

0.2

0.4

0.6

0.8

1

1.2

3.5

4

4.5

5

4.5

5

1.5

α [deg]

1

4

α [deg]

3

t [s] t [s]

2 0

0.5 0 -0.5 -1

-2

-1.5

-4 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

3.5

1.2

t [s] Fig. 4. Results obtained with the use of three successively decreased time steps: Airfoil with gap, far-field velocity 10 m/s, ..... – s = 8sref,    – s = sref, —— – s = sref/ 2.

4

t [s]

β [deg]

20 15 10 5 0 -5 -10 -15 -20

3

1.2

t [s]

β [deg]

0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5

5 4 3 2 1 0 -1 -2 -3 -4 -5

3

3.5

4

4.5

5

t [s] Fig. 5. Comparison of the results obtained with the use of the loose coupling – full line, and the strong coupling – dotted line, in the case of the airfoil with gap and farfield velocity 10 m/s.

116

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

Here A is a nX  nX matrix, B and C are nX  nQ matrices and O is T

T _ _ ZðtÞ ¼ ðhðtÞ; a_ ðtÞ; bðtÞÞ ;

nX þnQ

nQ  nQ zero matrix. The vector ðU; PÞ 2 R represents the solu

  nX n  onQ nþ1 tion U D;mþ1 with respect to the basis S ¼ wi ;0 i¼1 ; 0;qj j¼1

of the space XD  QD via formulas (3.18). The solution of system (3.19) is realized by the direct multifrontal method by Davis [9] implemented by package solver UMFPACK [10]. The Oseen iterative process is terminated under the condition

0

khh

B K¼@ 0

0

0 kaa 0

0

f ¼ ðL2 ; M a ; M b ÞT ;

1

0

C 0 A; kbb

Dhh

B D¼@ 0 0

nþ1 þPD ðU nþ1 D ; V D Þ  fD ðV D Þ  F D ðV D Þj 6 #kU D kL2 ðXÞ

ð3:20Þ

8V D 2 S; where # > 0 is a parameter representing the required accuracy of the numerical process. In our computations, # = 0.0001 was used. If the current iteration U nþ1 D;mþ1 satisfies the relation (3.20), we put U nþ1 ¼ U nþ1 D D;mþ1 and finish the Oseen iterative process. As numerical experiments show, only a few iterations (3.17), usually 3–4, have to be computed on each time level. 4. Discretization of the structural problem

M11 ¼ m;

C 0 A; Dbb

Daa 0

M 21 ¼ M 12 ;

2

3

4

M 32 ¼ M 23 ;

M 23 ¼ Ib þ dPF Sb cos b;

Further, we introduce the following notation: O – 3  3 zero matrix, I – unit 3  3 matrix, 0 – 3-dimensional zero vector and g – the vector of nonlinearities:

0

1 _ 2 sinða þ bÞ ðSa  Sb Þa_ 2 sin a þ Sb ða_ þ bÞ B C ^ 3C _2 _ g¼B @ dPF Sb b sin b þ 2ðdPF Sb Þa_ b sin b  kaa a A: ^bb b3 dPF Sb a_ 2 sin b  k

Z_ ¼ hðt; ZÞ;

5

dα/dt [deg/s]

α [deg]

2 1 0 -1 -2 2

3

4

5

t [s]

ð4:6Þ

250 200 150 100 50 0 -50 -100 -150 -200 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5

250 200 150 100 50 0 -50 -100 -150 -200 -250 -300 -3

-2

-1

0

1

2

3

10

15

α [deg]

15

dβ/dt [deg/s]

β [deg]

10 5 0 -5 -10 -15

0

1

2

3

t [s]

4

ð4:5Þ

Then system (2.8) is equivalent to the first-order system

3

1

ð4:4Þ

M 33 ¼ Ib :

h [mm]

0

ð4:2Þ

M 12 ¼ ðSa  Sb Þ cos a þ Sb cosða þ bÞ;

M31 ¼ M13 ;

t [s]

-3

1

ð4:3Þ

M13 ¼ Sb cosða þ bÞ;

dh/dt [mm/s]

h [mm]

In order to solve Eq. (2.8) of motion describing the airfoil vibrations, we transform them to a first-order system. We introduce the following notation:

1

0

M ¼ ðM ij Þ3i;j¼1 ;

M22 ¼ Ia  2dPF Sb þ 2dPF Sb cos b;

0

0

where the components of the nonlinear mass matrix M ¼ MðZÞ are

nþ1 nþ1 nþ1 jaD ðU nþ1 D ; U D ; V D Þ þ LD ðU D ; U D ; V D Þ

5 4 3 2 1 0 -1 -2 -3 -4 -5 -6

ð4:1Þ

5

1000 800 600 400 200 0 -200 -400 -600 -800 -1000 -15

-10

-5

0

5

β [deg] Fig. 6. Airfoil with gap: Functions h, a, b and their phase diagrams for far-field velocity 4 m/s.

117

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

where h is the vector function defined by

hðt; ZÞ ¼

M1 ðZÞ O O

!

f ðtÞ O

I



 



  D O g Zþ : O K 0

ð4:7Þ

r½ðk þ lÞr  d þ r  ðlrdÞ ¼ 0 in X0 ;

This system is equipped by the initial condition prescribing the value Z(0) given by conditions (2.17). The initial value problem for system (4.6) is solved by the fourth-order Runge–Kutta method. In the step from tn to tn+1 one needs the evaluation of the values f ð^tÞ at discrete instants ^t 2 ½t n ; tnþ1 . They are obtained by a linear extrapolation from the interval [tn1, tn] to [tn, tn+1]. If the values f(tn) and f(tn+1) have already been approximated, then f ð^tÞ is computed by the linear interpolation in the interval [tn, tn+1].

In this section we shall describe the algorithm of the numerical realization of the complete fluid–structure interaction problem. 5.1. Construction of the ALE mapping for three degrees of freedom In papers Ru˚zˇicˇka et al. [34] and Svácˇek et al. [36] we were concerned with flow induced vibrations of an airfoil with two degrees of freedom. In this case, the construction of the ALE mapping was relatively simple, based on the description of the airfoil as a solid body. However, for an airfoil with three degrees of freedom, this

where d = (d1, d2) is a displacement defined in X0. The Lamé coefficients k and l are computed by

k¼

Ea ra ; ð1 þ ra Þð1  2Ea Þ

l¼

Ea ; 2 þ 2ra

ð5:2Þ

where Ea is an artificial Young modulus and ra is an artificial Poisson number. The boundary conditions for d are prescribed by

ð5:3Þ

and on CW 0 they are determined by the functions h(t), a(t), b(t):

d1 ¼ X 1 cos a  X 2 sin a; d2 ¼ X 1 sin a þ X 2 cos a þ h;

X ¼ ðX 1 ; X 2 Þ 2 P0 ;

d1 ¼ X 1 cosða þ bÞ  X 2 sinða þ bÞ þ dPF cos a; d2 ¼ X 1 sinða þ bÞ þ X 2 cosða þ bÞ þ dPF sin a þ h; for the flap of the airfoil.

200

1

2

3

4

150 100 50 0 -50 -100

5

-150 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5

t [s]

h [mm] 3

dα/dt [deg/s]

α [deg]

2 1 0 -1 -2 -3

0

1

2

3

4

5

t [s]

250 200 150 100 50 0 -50 -100 -150 -200 -250 -300 -3

-2

-1

0

1

2

3

10

15

α [deg]

15

dβ/dt [deg/s]

β [deg]

10 5 0 -5 -10 -15

0

1

2

3

t [s]

4

X ¼ ðX 1 ; X 2 Þ 2 F 0 ; ð5:5Þ

250

0

ð5:4Þ

for the main part of the airfoil and

dh/dt [mm/s]

h [mm]

ð5:1Þ

djCD [CO ¼ 0

5. The realization of the complete coupled problem

5 4 3 2 1 0 -1 -2 -3 -4 -5 -6

approach cannot be applied. Therefore, the ALE mapping is constructed with the use of the linear equations describing the deformation of elastic bodies:

5

1000 800 600 400 200 0 -200 -400 -600 -800 -1000 -15

-10

-5

0

5

β [deg] Fig. 7. Airfoil with gap: Functions h, a, b and their phase diagrams for far-field velocity 8 m/s.

118

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

The solution of system (5.1) gives us the ALE mapping of X0 onto Xt by the relation

At ðXÞ ¼ X þ dðXÞ;

X 2 X0 ;

ð5:6Þ

for each time t. System (5.1) is discretized by the conforming piecewise linear finite elements on the mesh T 0D used for computing the velocity and pressure fields in the begining of the computational process in the polygonal approximation X0D of the domain X0. We introduce the finite element spaces

X D ¼ fdD ¼ ðdD1 ; dD2 Þ 2 CðX0D Þ2 ; dDi jK 2 P1K 8K 2 T 0D ; i ¼ 1; 2g;ð5:7Þ V D ¼ fuD 2 X D ; uD ðQ Þ ¼ 0 for all vertices Q 2 @ X0 g; and the form

 BD ðdD ; uD Þ ¼ ðk þ lÞðr  dD Þ; ðr  uD Þ X 0D  þ lrdD ; ruD X :

ð5:8Þ

0D

Then the approximate solution of problem (5.1), (5.3)–(5.5) is defined as a function dD 2 X D satisfying the Dirichlet boundary conditions defined by (5.3)–(5.5) with the values of h, a, b at time tn+1 and considered at the vertices lying on @ X0 and the identity

BD ðdD ; uD Þ ¼ 0 8uD 2 V D :

The use of linear finite elements is sufficient, because we need only to know the movement of the points of the mesh. It is possible to choose the Lamé coefficients k and l as constants, but it is more suitable to define them by (5.2), where the parameters Ea and ra are piecewise constant on the mesh T 0D . We define them by

ra jK ¼ 0:25; Ea jK ¼

Atnþ1 D ðXÞ ¼ X þ dD ðXÞ;

X 2 X0D :

ð5:11Þ

The knowledge of the ALE mapping at the time instants tn1, tn, tn+1 allows us to approximate the domain velocity with the aid of the second-order backward difference formula

wnþ1 D ðxÞ ¼

1 3x  4Atn D ðA1 t nþ1 D ðxÞÞ þ Atn1 D ðAt nþ1 D ðxÞÞ

2s

ð5:9Þ

250 200

dh/dt [mm/s]

2 0 -2 -4 0

1

2

3

4

150 100 50 0 -50 -100

5

-150 -6

t [s]

-4

-2

0

2

4

6

1

2

3

10

15

h [mm]

3

dα/dt [deg/s]

α [deg]

2 1 0 -1 -2 -3

0

1

2

3

4

5

t [s]

250 200 150 100 50 0 -50 -100 -150 -200 -250 -300 -3

-2

-1

0

α [deg]

15

dβ/dt [deg/s]

β [deg]

10 5 0 -5 -10 -15

0

1

2

3

t [s]

4

;

x 2 Xtnþ1 D : ð5:12Þ

4

h [mm]

ð5:10Þ

where meas (K) denotes the area of an element K. The mesh around the airfoil is typically refined into smaller triangles. Since smaller triangles imply the larger Young modulus Ea in (5.10), the mesh around the airfoil moves with the airfoil and therefore the deformation of elements close to the airfoil is small. If the displacement dD is computed at time tn+1, then in view of (5.6), the approximation of the ALE mapping is obtained in the form

6

-6

1 ; measðKÞ

5

1200 1000 800 600 400 200 0 -200 -400 -600 -800 -1000 -15

-10

-5

0

5

β [deg] Fig. 8. Airfoil with gap: functions h, a, b and their phase diagrams for far-field velocity 10 m/s.

119

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

2. Extrapolate linearly L2, Ma and Mb from the interval [tn1, tn] to [tn, tn+1]. Set m := 0. 3. Prediction of h, a, b: Compute the displacement h and angles a and b at time tn+1 as the solution of system (4.6) by the Runge–Kutta method. Denote it by h⁄, a⁄, b⁄. 4. On the basis of h⁄, a⁄, b⁄ determine the position of the airfoil at time tn+1, the domain Xtnþ1 D , the ALE mapping Atnþ1 D and the domain velocity wnþ1 D . 5. Solve the nonlinear discrete stabilized problem (3.13) at time tn+1 by the Oseen iterations. 6. Correction of h, a, b: Compute L2, Ma and Mb from (2.9)–(2.11) at time tn+1 and interpolate L2, Ma and Mb on [tn, tn+1]. Compute h, a, b at time tn+1 from (4.6) by the Runge–Kutta method. 7. If jh⁄  hj + ja⁄  aj + jb⁄  bj P e and m < M, set h⁄ = h, a⁄ = a, b⁄ = b, m: = m + 1 and go to 4. Otherwise, n:=n + 1 and go to 2.

5.2. Computation of aerodynamic forces acting on the airfoil In the case when the flap is not separated from the main body of the airfoil, the aerodynamic forces L2, Ma, Mb at time tn+1 are computed from (2.9)–(2.11) by using the approximation of the stress tensor (2.12) known from the solution UD = (uD, pD) of the stabilized discrete flow problem (3.13) and extrapolated to the boundary. The integrals in (2.9)–(2.11) are computed with the aid of numerical quadratures. In the case, when the flap is separated from the main body of the airfoil, i.e. Pt \ Ft = ;, the force and moments can be computed on the basis of a weak formulation similarly as in Svácˇek et al. [36]. 5.3. Coupling procedure In the solution of the complete coupled fluid–structure interaction problem it is necessary to apply a suitable coupling procedure. See, e.g. Badia and Codina [1] for a general framework. Here we apply the following algorithm. 0. Prescribe e > 0 – the measure of accuracy in the coupling procedure, and an integer M P 0 – the maximal number of iterations in the coupling procedure. 1. Assume that the solution UD = (uD, pD) of the discrete flow problem (3.13) and the force L2 and torsional moments Ma and Mb computed from (2.9)–(2.11) are known at time levels tn1 and tn.

If M = 0, then we get a loose (weak) coupling of the flow and structural problems. With increasing M and decreasing e, the coupling becomes stronger. 6. Numerical tests In this section we shall carry out test computations in order to analyze the robustness, accuracy and reliability of the described method and to demonstrate its applicability. In computations, we use the dimensionless form of the Navier–Stokes equations. We introduce the reference velocity U,

6

dh/dt [mm/s]

h [mm]

4 2 0 -2 -4 -6

0

1

2

3

4

5

t [s]

300 250 200 150 100 50 0 -50 -100 -150 -6

-4

-2

0

2

4

6

4 3 2 1 0 -1 -2 -3 -4

300 200

dα/dt [deg/s]

α [deg]

h [mm]

100 0 -100 -200

0

1

2

3

4

5 -300 -4

t [s]

-3

-2

-1

0

1

2

3

4

α [deg] 1500

15

1000

dβ/dt [deg/s]

β [deg]

10 5 0 -5 -10 -15

500 0 -500 -1000

0

1

2

3

t [s]

4

5 -1500 -15

-10

-5

0

5

β [deg] Fig. 9. Airfoil with gap: functions h, a, b and their phase diagrams for far-field velocity 12 m/s.

10

15

120

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

which is in our case the far-field flow velocity, and the reference length L, which is the length of the airfoil. We define the reference time Tref = L/U and the reference pressure pref = U2. The scaled dimensionless quantities and the Reynolds number are defined by

|G(h)|

0.1 0.05 0

5

10

15

0

20

0

5

0

5

10

15

20

5

10

15

20

0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0

0

5

15

20

15

20

10

frequency [Hz]

|G(β)| 0

10

15

20

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

5

10

frequency [Hz]

Fig. 10. Airfoil with gap: spectral analysis for far-field velocities 4 m/s (left) and 8 m/s (right).

0.25

0.25 0.2

|G(h)|

|G(h)|

0.2 0.15 0.1 0.05 0

0.15 0.1 0.05

0

5

10

15

0

20

0

5

0

5

10

15

20

1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

5

frequency [Hz]

|G(β)|

|G(β)|

2 1.5 1 0.5 5

10

frequency [Hz]

15

20

10

15

20

frequency [Hz]

2.5

0

10

frequency [Hz]

|G(α)|

|G(α)|

frequency [Hz]

0

m

frequency [Hz]

frequency [Hz]

0.6 0.5 0.4 0.3 0.2 0.1 0

UL

0.05

|G(α)|

|G(α)|

Re ¼

0.1

frequency [Hz]

|G(β)|

p ; pref

0.15

frequency [Hz]

0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

p ¼

0.2

0.15

0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0

x x ¼ ; L

0.25

0.2

|G(h)|

u ; U

u ¼

15

20

:

ð6:1Þ

Applying these substitutions in (2.7), we get the dimensionless Navier–Stokes equations in the form

0.25

0

t ; T ref

t ¼

8 7 6 5 4 3 2 1 0

0

5

10

15

frequency [Hz]

Fig. 11. Airfoil with gap: spectral analysis for far-field velocities 10 m/s (left) and 12 m/s (right).

20

121

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127 DA Dt 

1 u þ ððu  w Þ  rÞu þ rp  Re Mu ¼ 0;

ð6:2Þ

r  u ¼ 0:

We performed computations for the airfoil configurations considered by Losı´k and Cˇecˇrdle in [24,25], where the authors computed the flutter characteristics of a wing profile model by the MSC.NASTRAN package, which is based on a linear description of the structure behaviour and does not allow the solution in the case of large amplitudes, when nonlinear effects are significant. The numerical simulation was carried out for the airfoil NACA 0012 of length 0.3 m. The axis EA is placed at 40% of the length of the whole airfoil measured from the leading edge and the axis EF is placed at 80% of the length of the whole airfoil. We consider the situations, when the flap is attached to the main body of the airfoil without any gap, and when the flap is separated from the main body of the airfoil by the gap of the width g = 0.012 m. The following structural parameters were chosen chosen according to [24,25]:

m ¼ 0:086622 kg; khh ¼ 105:109 N=m;

The damping coefficients Dhh, Daa, Dbb are assumed to be zero. The fluid density q = 1.225 kg/m3 and the kinematic viscosity m = 1.5  105 m2/s. In our test computations, the Reynolds number was in the range between 4  104 and 4  105. The initial conditions for the structural system in all computations were set to

_ hð0Þ ¼ 5 mm; hð0Þ ¼ 0 mm=s;  að0Þ ¼ 3 ; a_ ð0Þ ¼ 0  =s; _ bð0Þ ¼ 0  =s: bð0Þ ¼ 0  ;

The computational process started from the solution of the flow past a fixed airfoil at time t = 0.01 s. At time t = 0.0 s the airfoil was released and the computation of the interaction started. The frequency analysis of the numerically simulated displacements c = h, a, b was carried out with the aid of the Fourier transform

Gðfn Þ ¼

Sa ¼ 0:000779598 kg m; Sb ¼ 0 kg m; 2

T

cðtÞe2pifn t dt

ð6:4Þ

with fn = nDf 2 [0, 50], Df = 0.1 Hz, approximated by the rectangle formula

2

Ia ¼ 0:000487291 kg m ; Ib ¼ 0:0000341104 kg m ; ^bb ¼ 0 or 5 N m=rad3 ; ^aa ¼ 0 N m=rad3 ; k k

Gðfn Þ ¼

dPF ¼ 0:12 m; l ¼ 0:079 m:

N1 X

cðtk Þe2pifn tk Dt:

ð6:5Þ

k¼0

dh/dt [mm/s]

h [mm]

Z 0

kaa ¼ 3:69558 N m=rad; kbb ¼ 0:2 N m=rad;

5 4 3 2 1 0 -1 -2 -3 -4 -5

ð6:3Þ

0

1

2

3

4

5

t [s]

250 200 150 100 50 0 -50 -100 -150 -200

-5 -4 -3 -2 -1

0

1

2

3

4

5

h [mm]

3

dα/dt [deg/s]

α [deg]

2 1 0 -1 -2 -3

0

1

2

3

4

5

t [s]

250 200 150 100 50 0 -50 -100 -150 -200 -250 -300 -3

-2

-1

0

1

2

3

10

15

α [deg] 15

dβ/dt [deg/s]

β [deg]

10 5 0 -5 -10 -15

0

1

2

3

t [s]

4

5

1000 800 600 400 200 0 -200 -400 -600 -800 -1000 -15

-10

-5

0

5

β [deg] Fig. 12. Airfoil without gap: functions h, a, b and their phase diagrams for far-field velocity 2 m/s.

122

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

Here i is the imaginary unit, Dt = T/N and N is the number of time steps in the interval [0, T). The results of the frequency analysis are shown in graphs of the quantity

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jGðfn Þj ¼ R2 ðGðfn ÞÞ þ I2 ðGðfn ÞÞ: In what follows, we present the results of several test computations. In Sections 6.1 and 6.2 we use the structural parameters men^aa ¼ k ^bb ¼ 0. In Section 6.3 we consider this case tioned above with k as well as the case of nonzero nonlinear torsional stiffness.

6.2. Convergence analysis In this section we shall present the results of experimental convergence analysis. We solve the coupled fluid–structure interaction problem using successively refined space meshes and successively decreasing time steps. Moreover, we are concerned with the effect of the loose and strong coupling. 1) First we solve the coupled problem in the case of the airfoil with gap for the far-field velocity 10 m/s with the use of the small time step s = 0.0003 s on three non-nested meshes with 4203 vertices and 8256 elements (coarse mesh), 28,924 vertices and 57,364 elements (refinement 1) and 87,905 vertices and 175,032 elements (refinement 2). Fig. 2 shows the coarse mesh and its details. In Fig. 3 we show the graphs of h, a, b in dependence on time obtained on these meshes. One can observe the convergence of the results, when the mesh is successively refined. Since the results computed on both refined meshes are practically identical, the solution of other test problems presented in Section 6.3 was carried out with the aid of the mesh ‘‘refinement 1’’, in order to save computational time. 2) The second test analyzes experimentally the convergence for decreasing time step. In this case, the space mesh ‘‘refinement 2’’ was used (overkill in the space discretization) and the coupled problem with the far-field velocity 10 m/s was solved with three time steps s = 8sref, s = sref = 0.0003 s,

6.1. Comparison with NASTRAN computations

5 4 3 2 1 0 -1 -2 -3 -4 -5 -6

dh/dt [mm/s]

h [mm]

In Table 1 we compare the eigen frequences of the airfoil in vacuo, i.e. without flow past the airfoil, computed with the aid of the MSC.NASTRAN package and on the basis of the nonlinear system (2.8) with linear torsional stiffness and Fourier analysis. We see a good agreement. Further, in Table 2 we present the vibration eigenfrequences of the functions h, a, b, for the far-field flow velocity 2m/s, critical flutter velocity and frequency computed by the presented method, compared with the NASTRAN computations. FEM0 and FEM1 denote the results for the airfoil without gap and with gap, respectively. The shown results are again in good agreement. The comparison of the critical velocity computed by MSC.NASTRAN package with our results indicates that MSC.NASTRAN is more suitable for the analysis of an airfoil with gap than an airfoil without gap.

0

1

2

3

4

5

t [s]

250 200 150 100 50 0 -50 -100 -150 -200 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5

h [mm] 3

dα/dt [deg/s]

α [deg]

2 1 0 -1 -2 -3

0

1

2

3

4

5

t [s]

250 200 150 100 50 0 -50 -100 -150 -200 -250 -300 -3

-2

-1

0

1

2

3

10

15

α [deg] 15

dβ/dt [deg/s]

β [deg]

10 5 0 -5 -10 -15 0

1

2

3

t [s]

4

5

1000 800 600 400 200 0 -200 -400 -600 -800 -1000 -15

-10

-5

0

5

β [deg] Fig. 13. Airfoil without gap: functions h, a, b and their phase diagrams for far-field velocity 4 m/s.

123

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

s = sref/2. Fig. 4 shows the comparison of h, a, b in depen-

1) First we assume that the coefficients in the nonlinear tor^aa ¼ k ^bb ¼ 0. In Figs. 6–11, the graphs sional spring stiffness k of the functions h, a, b, their phase diagrams and spectral analysis are shown for far-field velocities 4, 8, 10, 12 m/s in the case of the airfoil with gap. For small far-field flow velocities the vibration amplitudes of vertical displacement h and the rotation angles a, b of the airfoil are decreasing in time (see Figs. 6 and 7). The system is stable and low level sustained vibrations are caused by vortex separation on the airfoil. The vibration frequencies of the airfoil (see Fig. 10) are close to the eigenfrequencies of the structure in vacuo and correspond to the eigenfrequencies computed by NASTRAN (see Tables 1 and 2). By increasing the far-field flow velocity, the second and third resonance frequencies are getting closer, which is a typical behaviour for pre-flutter regimes of the airfoil vibration, when the coupling between the rotation of the main part of the airfoil and the flap rotation becomes stronger (see Figs. 8, 9 and 11 for 10 m/s and 12 m/s). The vibration flutter regime at 12 m/s can be considered as a stable limit cycle oscillation (LCO) regime with large amplitudes of the flap vibrations around 14 deg for the rotation angle b, and the flutter frequency 14.89 Hz. In Fig. 11 one can see that both higher frequencies of the self-oscillating airfoil for the rotation angles a and b merge together. Significant vortices are

dence on time for these three time steps. On the basis of these results, the computations presented in Section 6.3 were computed with time step sref. 3) An important issue is the use of the loose or strong coupling. Fig. 5 shows the comparison between the time dependence of h, a, b computed by the algorithm presented in Section 5.3 with the aid of the loose coupling (M = 0) and the ‘‘strong’’ coupling with prescribed accuracy e = 105 and the maxinal number of iterations M = 50. In our numerical experiments the prescribed accuracy was achived after m 6 5 iterations in each time step. It is interesting that there is a very small difference between the rotation angles computed with the aid of the loose and strong coupling. On the other hand, the influence of the choice of the type of coupling on the behaviour of the displacement h is quite significant. This is the reason that in our numerical tests the strong coupling was used.

6.3. Airfoil vibrations analysis in dependence on the far-field velocity

5 4 3 2 1 0 -1 -2 -3 -4 -5 -6

250 200

dh/dt [mm/s]

h [mm]

In this section, we analyze the development of airfoil vibrations in dependence on the far-field velocity starting from damped regimes up to limit cycle oscillation (LCO) and flutter.

0

1

2

3

4

150 100 50 0 -50 -100

5

-150 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5

t [s]

h [mm] 3

dα/dt [deg/s]

α [deg]

2 1 0 -1 -2 -3

0

1

2

3

4

5

t [s]

250 200 150 100 50 0 -50 -100 -150 -200 -250 -300 -3

-2

-1

0

1

2

3

10

15

α [deg] 15

dβ/dt [deg/s]

β [deg]

10 5 0 -5 -10 -15

0

1

2

3

t [s]

4

5

1000 800 600 400 200 0 -200 -400 -600 -800 -1000 -15

-10

-5

0

5

β [deg] Fig. 14. Airfoil without gap: functions h, a, b and their phase diagrams for far-field velocity 8 m/s.

124

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

2) In Figs. 12–16 we present results obtained for the airfoil without the gap. The development of the amplitude of vibrations in dependence on the far-field velocity is qualitatively similar to the results obtained for the airfoil with the gap, if the far-field velocity is small, but the flutter is obtained faster, when the far-field velocity increases. It can be explained by an interference of the flow in the gap with the main flow. Fig. 17 shows the vibrations of the airfoil without gap for the far-field velocity 20 m/s and initial conditions for the airfoil position h(0) = 10 mm, a(0) = 5°, b(0) = 0°. This example, when the angle b of the flap rotation attains the values up to 25°, shows that the developed method is applicable to

generated on the flap surface and behind the profile (see Figs. 19 and 20). For the higher flow velocity (14 m/s), the LCO is unstable and the amplitude of the rotation angle b is increasing in time reaching higher values, up to about 18 deg for the flap motion during 0.7 s (see Fig. 18). The simulation stopped due to numerical problems caused by the large mesh deformation. These results are in agreement with the NASTRAN computations, according to which the system becomes unstable by flutter in torsion for the far-field flow velocity at 11.32 m/s and the flutter frequency 14.87 Hz (see Table 2).

0.25

0.25 0.2

|G(h)|

|G(h)|

0.2 0.15 0.1

0.15 0.1

0.05 0

0.05 0

5

10

15

0

20

0

5

10

0.14 0.12 0.1 0.08 0.06 0.04 0.02 0

15

20

0.2 0.15 0.1 0.05 0

5

10

15

0

20

0

5

10

frequency [Hz]

|G(β)|

|G(β)|

20

0.25

frequency [Hz] 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

15

frequency [Hz]

|G(α)|

|G(α)|

frequency [Hz]

0

5

10

15

20

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

5

10

frequency [Hz]

15

20

frequency [Hz]

Fig. 15. Airfoil without gap: spectral analysis for far-field velocities 2 m/s (left) and 4 m/s (right).

0.25

|G(α)|

0.15 0.1 0.05 0

0

5

10

15

20

0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

0

5

frequency [Hz]

|G(β)|

|G(h)|

0.2

1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

10

frequency [Hz]

5

10

15

20

frequency [Hz] Fig. 16. Airfoil without gap: spectral analysis for far-field velocity 8 m/s.

15

20

125

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

15

6

h [mm]

4 10

-h [mm]

5

2 0 -2 -4 -6

0

0.2

0.4

0.6

0

α [deg]

-5

-10

-15

0.8

1

1.2

0.8

1

1.2

1

1.2

t [s]

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

time [s]

4 3 2 1 0 -1 -2 -3 -4

0

0.2

0.4

0.6

t [s] 5 4

β [deg]

3 2

α [°]

1 0

20 15 10 5 0 -5 -10 -15 -20

0

0.2

0.4

0.6

0.8

t [s]

-1

Fig. 18. Airfoil with gap: comparison of airfoil vibrations for far-field velocity 14 m/ s in the case of linear spring stiffness (large amplitudes marked by full line) and nonlinear spring vibrations (smaller amplitudes marked by dotted line).

-2 -3 -4 -5

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

6

0.4

4

h [mm]

time [s] 30

2 0 -2 -4

20

-6 0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

0.6

0.8

t [s]

10

2

0

α [deg]

β [°]

3

-10

1 0 -1 -2 -3

-20

0

0.2

0.4

t [s] -30

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

15

0.4

10

Fig. 17. Airfoil without gap: functions h, a, b for far-field velocity 20 m/s.

β [deg]

time [s]

5 0 -5 -10

the simulation of flow induced postcritical airfoil vibrations with large amplitudes. Both in the case of the airfoil with gap and without gap, we can observe an interesting phenomenon that in the range of low far-field velocities with a small increase of its magnitude

-15

0

0.2

0.4

1

t [s] Fig. 19. Airfoil with gap: points on the graphs of the functions h, a, b corresponding to the velocity fields computed for the far-field velocity 12 m/s, shown in Fig. 20.

126

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

Fig. 20. Airfoil with gap: velocity distribution for the far-field velocity 12 m/s.

the vibrations amplitude is decreasing (compare Fig. 6 with Fig. 7 and Fig. 12 with Fig. 13). It is caused by aerodynamic damping forces. This phenomenon disappears with a sufficient increase of the far-field velocity. 3) Further, we investigated the influence of the nonlinear torsional spring stiffness. Fig. 18 show the comparison of vibrations of the airfoil with gap for linear flap spring stiffness ^bb ¼ 5 N m=rad3 in the and nonlinear flap spring stiffness k case of the far-field velocity 14 m/s. In the case of linear stiffness (similarly as in Fig. 17), the amplitude of vibrations becomes large with increasing time, which causes the collaps of the computational process due to a large deformation

of the computational mesh. The additional cubic nonlinear spring stiffness reduces the vibrations amplitude and changes the flutter type from a catastrophic to the benign. 7. Conclusion The paper is concerned with the development, analysis and applications of the numerical method for the simulation of flow induced vibrations of an airfoil with three degrees of freedom. The fluid flow is described by the 2D incompressible Navier Stokes equations in the ALE formulation allowing to take into account time dependence of the computational domain for large

M. Feistauer et al. / Computers & Fluids 49 (2011) 110–127

deformations of the structure in postcritical regimes, when the aeroelastic system is unstable by flutter or divergence. The high values of the vibration amplitudes encountered in such regimes (e.g., in rotation up to 25°) need to describe the dynamic behaviour of the studied airfoil by three nonlinear equations of motion with nonlinear mass and stiffness matrices. The coupled system of PDEs for the fluid flow and ODEs for the airfoil motion is solved by the FEM for the flow field around the oscillating profile and the Runge–Kutta method for the integration of equations for the profile motion in time domain. Good agreement of the results was obtained, both for the flutter frequency and the critical flutter flow velocity determined from the numerically simulated airfoil motions in time domain, with the MSC.NASTRAN computations using the frequency-modal analysis considering the linear model of vibrations and the doublet-lattice aerodynamic theory [24,25]. From the mathematical point of view, the paper presents an accurate and robust technique for the simulation of flow induced airfoil vibrations. It contains a detailed description of all important ingredients necessary for the realization of the general framework proposed, e.g. in [1] or [15]. Moreover, the paper contains an experimental convergence analysis with respect to the mesh refinement in space and time. Presented numerical results prove the applicability and robustness of the method. The most significant contribution of the paper from the point of view of computational aeroelasticity and fluid dynamics consists in developing the numerical method for the simulation of postcritical behaviour of the 3-DOF airfoil system after losing the aeroelastic stability. Estimation of vibrations with large amplitudes of the airfoil with flap could be used instead of expensive experiments on aeroelastic models performed in wind tunnels with a possible destruction of the models, or dangerous flutter tests carried out on airplane prototypes. The numerical simulations allow the analysis of number of cases in a relatively short computational time and by much more lower financial expenses than needed for flutter experiments. There are the following subjects for further work: - to include turbulence modelling in the numerical simulation of flow induced airfoil vibrations, - to extend the developed technique to the coupling of compressible flow with a vibrating airfoil. Acknowledgements This research was supported under the Research Plan MSM 0021620839 of the Ministry of Education of the Czech Republic. It was also partly supported under the Grants No. 201/08/0012 and (in the year 2011) P101/11/0207 of the Czech Science Foundation. References [1] Badia S, Codina R. On some fluid–structure iterative algorithms using pressure segregation methods. Application to aeroelesticity. Int J Numer Methods Eng 2007;72:46–71. [2] Badia S, Nobile F, Vergara C. Fluid–structure partioned procedures based on Robin transmission conditions. J Comput Phys 2008;227:7027–51. [3] Bisplinghoff RL, Ashley H, Halfman RL. Aeroelasticity. New York: Dover; 1996. [4] Bolotin VV. Nonconservative problems of the theory of elastic stability. New York: Macmillan; 1963. [5] Brezzi F, Falk RS. Stability of higher-order Hood–Taylor methods. SIAM J Numer Anal 1991;28:581–90.

127

[6] Brezzi F, Fortin M. Mixed and hybrid finite element method. Springer series in computational mathematics 15. New York: Springer-Verlag; 1991. [7] Ciarlet PG. The finite element method for elliptic problems. Amsterdam: North-Holland; 1979. [8] Chung KW, He YB, Lee BHK. Bifurcation analysis of a two-degree-of-freedom aeroelastic system with hysteresis structural nonlinearity by a perturbationincremental method. J Sound Vib 2009;320:163–83. [9] Davis TA. A column pre-ordering strategy for the unsymmetric-pattern multifrontal method. ACM Trans Math Soft 2004;30:165–95. [10] T.A.Davis, UMFPACK V4.0, University of Florida.. [11] Dessi D, Mastroddi F. A nonlinear analysis of stability and gust response of aeroelastic systme. J Fluid Struct 2008;24:436–45. [12] Dowell EH, editor. Aeroelasticity of plates and shells. Dordrecht: Kluwer; 1974. [13] Dowell EH. A modern course in aeroelasticity. Dodrecht: Kluwer; 1995. [14] Farhat C, van der Zee KG, Geuzaine P. Provably second-order time-accurate loosely-coupled solution algorithms for transient nonlinear computational aeroelasticity. Comput Methods Appl Mech Eng 2006;195:1973–2001. [15] Fernández MA, Gerbeau J-F. Algorithms for fluid–structure interaction problems. In: Formaggia L, Quarteroni A, Veneziani A, editors. Cardiovascular mathematics. Modeling and simulation of the circulatory sestem. SpringerVerlag Italia; 2009. p. 307–46. [16] Gelhard T, Lube G, Olshanskii MA, Starcke JH. Stabilized finite element schemes with LBB-stable elements for incompressible flows. J Comput Appl Math 2005;177:243–67. [17] Grandmont C. Existence of weak solutions for the unsteady interaction of a viscous fluid with an elastic plate. SIAM J Math Anal 2008;40:716–37. [18] Guidorzi M, Padula M, Plotnikov PI. Hopf solutions to a fluid–elastic interaction model. Math Models Methods Appl Sci 2008;18:215–69. [19] Hoffmann KH, Starovoitov VN. On a motion of a solid body in a viscous fluid. Two-dimensional case. Adv Math Sci Appl 1999;9:633–48. [20] Horácˇek J, Svácˇek P, Ru˚zˇicˇka M, Feistauer M. Contribution to finite element modelling of airfoil aeroelastic instabilities. Appl Comput Mech 2007;1:43–52. [21] Jones DP, Roberts I, Gaitonde AL. Identification of limit cycles for piecewise nonlinear aeroelastic systems. J Fluid Struct 2007;23:1012–28. [22] Librescu L, Marzocca P. Advances in the linear/nonlinear control of aeroelastic structural systems. Acta Mech 2005;178:147–86. [23] Librescu L, Marzocca P, Silva WA. Aeroelasticity of 2-D lifting surfaces with time-delayed feedback control. J Fluid Struct 2005;20:197–215. [24] Losı´k V, Cˇecˇrdle J. Aeroelastic analysis of a verifying model of an airplane construction with three degrees of freedom – part I. Research Report No. V1833/3210/05. Aircraft Research and Test Institute (ARTI), Prague-Letany, 2005 [in Czech]. ˇ ecˇrdle J. Flutter computation of 3-DOF wing profile model. Technical [25] Losı´k V, C Report P-PL-0061/07. Aircraft Research and Test Institute (ARTI), PragueLetany, 2007 [in Czech]. [26] Lube G. Stabilized Galerkin finite element methods for convection dominated and incompressible flow problems. Numer Anal Methods Model 1994;29: 85–104. [27] Na S, Librescu L, Kim M-H, Jeong I-J, Marzocca P. Robust aeroelastic control of flapped wing systems using a sliding mode observer. Aerosp Sci Technol 2006;10:120–6. [28] Naudasher E, Rockwell D. Flow-induced vibrations. Rotterdam: A.A. Balkema; 1994. [29] Neustupa J. Existence of a weak solution to the Navier–Stokes equation in a general time-varying domain by the Rothe method. Math Methods Appl Sci 2009;32:653–83. [30] Nomura T, Hughes TJR. An arbitrary Lagrangian–Eulerian finite element method for interaction of fluid and a rigid body. Comput Methods Appl Mech Eng 1992;95:115–38. [31] Paidoussis MP. Fluid–structure interactions. Slender structures and axial flow, vol. I. San Diego: Academic Press; 1998. [32] Paidoussis MP. Fluid–structure interactions. Slender structures and axial flow., vol. II. Academic Press, London, UK: Elsevier; 2004. [33] Peters DA. Two-dimensional incompressible unsteady airfoil theory – an overview. J Fluid Struct 2008;24:295–312. [34] Ru˚zˇicˇka M, Feistauer M, Horácˇek J, Svácˇek P. Interaction of incompressible flow and a moving airfoil. Electron Trans Numer Anal 2008;32:123–33. [35] Shim J-H, Na S, Chung C-H. Aeroelastic response of an airfoil–flap system exposed to time-dependent disturbances. KSME Int J 2004;18:560–72. [36] Svácˇek P, Feistauer M, Horácˇek J. Numerical simulation of flow induced airfoil vibrations with large amplitudes. J Fluid Struct 2007;23:391–411. [37] Taylor C, Hood P. A numerical solution of the Navier–Stokes equations using the finite element technique. Comput Fluids 1973;1:73–100. [38] Verfürth R. Error estimates for a mixed finite element approximation of the Stokes equations. RAIRO Anal Numer/Numer Anal 1984;18:175–82. [39] Witteveen JAS, Bijl H. An unsteady adaptive stochastic finite elements formulation for rigid-body fluid–structure interaction. Comput Struct 2008;86:2123–40.