Acta Polytechnica Vol. 52 No. 6/2012

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

Petr Sváček1, Jaromír Horáček2, Radek Honzátko3, Karel Kozel1

1Faculty of Mechanical Engineering, Czech Technical University in Prague, Prague, Czech Republic
2Institute of Thermomechanics, Academy of Sciences of the Czech Republic, Prague, Czech Republic
3Faculty of Production Technology and Management, Jan Evangelisty Purkyně University in Ústí nad labem, Ústí nad
labem, Czech Republic.

Corresponding author: Petr.Svacek@fs.cvut.cz

Abstract
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.

Keywords: laminar flow, finite volume method, finite element method, arbitrary Lagrangian-Eulerian method, nonlinear
aeroelasticity.

1 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]. Com-
mercial codes are widely used (e.g. NASTRAN or
ANSYS) in technical practice, and aeroelasticity com-
putations are performed mostly in the linear domain,
which enables the stability of the system to be deter-
mined.
Recently, the research has also focused on numeri-

cal 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 mod-
erate, 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 oscilla-
tion (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 tech-
niques 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 instabil-
ity of actively controlled airfoils. Flutter boundaries
and LCO of aeroelastic systems with structural non-
linearities 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 ele-
ment method (FEM, see [28]), were employed for a nu-
merical solution of the interaction of 2-D incompress-
ible viscous flow and the elastically supported profile
NACA 0012. The mathematical model of the flow is
represented by the incompressible Navier-Stokes equa-
tions. 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 de-
scribed in Section 3. The dual-time stepping method
is applied to the numerical solution of unsteady sim-
ulations 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

105



Acta Polytechnica Vol. 52 No. 6/2012

khh

(u  ,v )T8 8

kϕϕ

EA
CG

L(t)M(t)

h

ϕ

Figure 1: Elastic support of the profile on transla-
tional and rotational springs.

unsteady flow calculations satisfies the geometric con-
servation law (GCL), cf [19]. The application of FEM
is described in Section 4. Section 5 presents numerical
results for both methods.

2 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]:

mḧ + Sϕϕ̈ cos ϕ−Sϕϕ̇2 sin ϕ + khh = −L,
Sϕḧ cos ϕ + Iϕϕ̈ + kϕϕ = M, (1)

where h is vertical displacement of the elastic axis
(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], kh is the bending stiffness [N/m], Iϕ
is the inertia moment around the elastic axis [kg m2],
and kϕ is the torsional stiffness [N m/rad].
The aerodynamic lift force L [N] acting in the ver-

tical direction (upwards positive) and the torsional
moment M [N m] (clockwise positive) are defined as

L = −dρU2∞c
∫

ΓWt

2∑
j=1

σ2jnjds, (2)

M = dρU2∞c
2
∫

ΓWt

2∑
i,j=1

σijnjr
⊥
i ds,

where d [m] is the airfoil depth, ρ [kg m−3] is the con-
stant fluid density, U∞ [m s−1] denotes the magnitude
of the far field velocity, c [m] denotes the airfoil chord,
n = (n1,n2) 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⊥ = (r⊥1 ,r⊥2 ) =
(−(y − yEA),x − xEA), where x,y are the non-
dimensional coordinates of a point on the profile sur-
face (i.e., the physical coordinates divided by the

r
d

n

r
EA

Figure 2: Airfoil segment.

airfoil chord c), ΓWt denotes the moving part of
the boundary (i.e., the surface of the airfoil), and
(xEA,yEA) are the non-dimensional coordinates of
the elastic axis (see Figure 2). The components of
the (non-dimensional) stress tensor σij are given by

σ11 =
(
−p +

2
Re

∂u

∂x

)
σ12 = σ21 =

1
Re

(
∂u

∂y
+
∂v

∂x

)
(3)

σ22 = ρ
(
−p +

2
Re

∂v

∂y

)
,

where u = (u,v) is the non-dimensional fluid velocity
vector, p is the non-dimensional pressure, Re is the
Reynolds number defined as

Re =
U∞c

ν
,

and ν [m2 s−1] is the fluid kinematic viscosity.
Equations (2), (3) together with the boundary con-

ditions 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.

3 Finite volume method

3.1 Mathematical model of the flow
The flow of viscous incompressible fluid in the
computational domain Ωt is described by the two-
dimensional Navier-Stokes equations written in the
conservative form

(DW)t + Fcx + G
c
y =

1
Re

(Fvx + G
v
y), (4)

106



Acta Polytechnica Vol. 52 No. 6/2012

where W is the vector of conservative variables

W = (p,u,v)T ,

Fc = Fc(W), Gc = Gc(W) are the inviscid physical
fluxes defined by

Fc = (u,u2 + p,uv)T , Gc = (v,uv,v2 + p)T ,

and Fv = Fv(W), Gv = Gv(W) are the viscous
physical fluxes defined by

Fv = (0,ux,vx)T , Gv = (0,uy,vy)T ,

where the partial derivatives with respect to x,y
and t are denoted by the subscripts x,y and t , re-
spectively. 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 dimen-
sionless pressure, respectively. The detailed relation
of the dimensionless form of the governing equations
and the relationship between dimensional and non-
dimensional variables can be found, e.g., in [11], [15].
In order to time discretize (4), a partition 0 =

t0 < t1 < · · · < tm = T of the time interval [0,T]
is considered, and the (variable) time step ∆tn =
tn+1 − tn is denoted by ∆tn. Further, Wn denotes
the approximation of the solution W at the time
instant tn. 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 Wt by the second order backward difference
formula (BDF2)

∂W

∂t
(tn+1) ≈

αn0 W
n+1 + αn1 Wn + αn2 Wn−1

∆t
, (5)

where αn0 = 2 − 1/(1 + θ), αn1 = −(1 + θ), αn2 =
θ2/(1 + θ), θ = (∆tn)/(∆tn−1). The artificial com-
pressibility method, see [3, 16], together with a time-
marching method is used for the steady state com-
putations. 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 non-
linear problem is then found as a solution of the
steady-state problem in dual time τ, i.e.

DβWτ = −R̃(W), (6)

where R̃(W) is the unsteady residuum defined by

R̃(W) = D
αn0 W + αn1 Wn + αn2 Wn−1

∆t
+ R(W),

and R(W) is the steady residuum defined by

R(W) = (Fc −
1
Re

Fv)x + (Gc −
1
Re

Gv)y.

The required real-time accurate solution at the time
level n + 1 satisfies R̃(Wn+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.

3.2 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

∂

∂t

∫
Ci(t)

DW dxdy +
∫
∂Ci

(F̃c dy − G̃c dx)

−
∫
∂Ci

1
Re

(Fv dy −Gv dx) = 0, (7)

where F̃c = F̃c(W,w1) = Fc(W) − w1DW, G̃c =
G̃c(W,w2) = Gc(W) −w2DW, w = (w1,w2)T is the
domain velocity, and Ci = Ci(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 Cj(t) satisfying

Ωt =
⋃
j∈J

Cj(t).

Let us denote by N(i) the set of all neighbouring
cells Cj(t) of cell Ci(t), i.e. the set of all cells whose
intersection with Ci(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 Ci(t) and Cj(t).
The mean value of vector W over cell Ci(t) at time
instant tn is approximated by Wni , and symbol Wn
denotes the vector formed by the collection of all
Wni , i ∈ J. In addition, the numerical approximation
depends on the ALE mapping At and on the domain
velocity w(t), see also [19]. In order to show this
dependence we shall denote at time instant tn the
vector of grid point positions by Xn and the volume
of cell Ci(t) by |Cni | = |Ci(tn)|. Thus the resulting
equation for the i-th finite volume cell reads

D
∂
(∣∣Ci(t)∣∣Wi)

∂t
+ Ri(W, t) = 0, (8)

107



Acta Polytechnica Vol. 52 No. 6/2012

with

Ri(W, t) = Fci (W(t),X(t), w(t)) −F
v
i (W(t),X(t)),

where Fci (W,X, w) and Fvi (W,X) represent the nu-
merical approximations of the convective and diffusive
fluxes, respectively.

Further, the BDF2 formula is used for the approxi-
mation of the time derivative in (7) by the difference

δn+1t Wi =
1

∆tn
(
αn0 |C

n+1
i |W

n+1
i + α

n
1 |C

n
i |W

n
i

+ αn2 |C
n−1
i |W

n−1
i

)
.

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 tn+1 the following problem written in the ALE
formulation needs to be solved

Dβ
∣∣Cn+1i ∣∣∂Wn+1i∂τ

+ Dδn+1t Wi = −R
n+1
i (W

n+1), (9)

where

Rn+1i (W
n+1) = ω1F

c,n+1/2
i (W

n+1)

+ ω2F
c,n−1/2
i (W

n+1) −Fv,n+1i (W
n+1)

+ ADn+1i (W
n+1)

and ω1 = αn+1, ω2 = −αn−1/θ. The symbol
ADn+1i (Wn+1) stands for the artificial dissipation
term (see [15]). The approximation of the convec-
tive fluxes Fc,n±1/2i is given by

Fc,mi (W
n+1) =

∑
j∈N(i)

(
F̃c(Wn+1ij ,w

m
1 )∆y

m
ij

− G̃c(Wn+1ij ,w
m
2 )∆x

m
ij

)
,

where Wn+1ij = (W
n+1
i + W

n+1
j )/2, X

n+1/2 =
(Xn+1 + Xn)/2, wn+1/2 = (Xn+1 − Xn)/∆tn and
∆xmij , ∆ymij are the x- and y-coordinate differences
of the segment Γmij . The approximation of the fluxes
given by Ri(Wn+1) guarantees the satisfaction of
GCL, see [19]. The approximation of the viscous
fluxes is evaluated on the mesh configuration at the
time instant tn+1 and reads

Fv,n+1i (W
n+1)

=
1
Re

∑
j∈N(i)

F̂v|n+1ij ∆y
n+1
ij − Ĝ

v|n+1ij ∆x
n+1
ij ,

where

F̂v|n+1ij =
(
0, ûx|n+1ij , v̂x|

n+1
ij

)T
,

Ĝv|n+1ij =
(
0, ûy|n+1ij , v̂y|

n+1
ij

)T
,

and ûx|n+1ij , ûy|
n+1
ij , v̂x|

n+1
ij , v̂y|

n+1
ij are the approxi-

mations 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

R̃n+1i (W
n+1) = Dδn+1t W

n+1
i + Ri(W

n+1),

equation (9) can be written in the form:

∂Wn+1i
∂τ

= −D−1β
1

|Cn+1i |
R̃i(Wn+1),

and solved with the aid of Runge-Kutta method,
cf. [15].

4 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

DAu
Dt

−
1
Re
4u + (u − w) ·∇u + ∇p = 0, (10)

div u = 0,

in Ωt ⊂ R2, where D
A

Dt
is the ALE derivative,

u denotes the velocity vector, p denotes the non-
dimensional 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

a) u(x,t) = uD(x) for x ∈ ΓD,
b) u(x,t) = w(x,t) for x ∈ ΓWt, (11)

c) −pn +
1
2
(u · n)−u + ν

∂u
∂n

= 0, on ΓO,

where n denotes the unit outward normal vector, uD
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) = u0(x).

4.2 Finite element approximation
FEM is well known as a general discretization method
for partial differential equations. However, straightfor-
ward application of FEM procedures often fails in the
case of incompressible Navier-Stokes equations. The
reason is that momentum equations are of advection-
diffusion type, with the dominating advection , the
Galerkin FEM leads to unphysical solutions if the
grid is not fine enough in regions of strong gradients

108



Acta Polytechnica Vol. 52 No. 6/2012

(e.g. in the boundary layer). In order to obtain phys-
ically 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 stabi-
lizing 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 approxi-
mate the time derivative by the second order backward
difference formula,

DAu
Dt

(x,t) ≈
3un+1 − 4ûn + ûn−1

2∆t
,

where un+1 is the flow velocity at time tn+1 defined
on the computational domain Ωn+1, and ûk is the
transformation of the flow velocity at time tk 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∆ ⊂ H1(Ωn+1) and
Q∆ ⊂ L2(Ωn+1) for approximation of the velocity
components and pressure, respectively, and the sub-
space of the test functions is denoted by X∆ ⊂ W∆.
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 as-
sume that the domain Ω = Ωn+1 is a polygonal ap-
proximation of the region occupied by the fluid at
time tn+1 and the finite element spaces are defined
over a triangulation T∆ of domain Ωt as piecewise
polynomial functions. In our computations, the well-
known Taylor-Hood P2/P1 conforming finite element
spaces are used for the velocity/pressure approxima-
tion. 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 = tn+1 reads: Find U = (u,p) ∈ W∆ × Q∆, p :=
pn+1, u := un+1, such that u satisfies approximately
the conditions (11 a-b) and

a(U; U,V ) + L(U; U,V ) + P(U,V )
= f(V ) + F(U; V ), (12)

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

a(U∗; U,V ) =
3

2∆t
(u, v)Ω +

1
Re

(∇u,∇v)Ω
+ (w ·∇u, v)Ω − (p,∇· v)Ω + (∇· u,q)Ω, (13)

f(u, v) =
1

2∆t
(4ûn − ûn−1, v)Ω,

where w = u∗ − wn+1, and the scalar product in
L2(Ω) is denoted by (·, ·)Ω. Further, the SUPG/PSPG

stabilization terms are used in order to obtain a stable
solution also for large Reynolds number values,

L(U∗; U,V ) =
∑
K∈T∆

δK

(3u
2τ
−

1
Re
4u

+ (w ·∇)u + ∇p,ψ(V )
)
K
,

F(U∗; V ) =
∑
K∈T∆

δK

(4ûn − ûn−1
2τ

,ψ(V )
)
K
,

where ψ(V ) = (w ·∇)v + ∇q, w stands for the trans-
port velocity at time instant tn+1, i.e. w = v∗−wn+1,
and (·, ·)K denotes the scalar product in L2(K). The
term P(U,V ) is the additional grad-div stabilization
defined by

P(U,V ) =
∑
K∈T∆

τK(∇· u,∇· v)K.

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 hK

δK = δ∗hK2, and τK = 1.

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 stabi-
lization 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 cou-
pled 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 struc-
ture motion are approximated repeatedly in order to
converge to a solution which will satisfy all interface
conditions.

5 Numerical results
Numerical flow induced vibration results for both
FVM and FEM are presented for the NACA 0012
profile, taking into account the following parameters
of the flowing air and the profile: m = 8.66×10−2 kg,
Sϕ = −7.797 × 10−4 kg m, Iϕ = 4.87 × 10−4 kg m2,
kh = 105.1 N/m, kϕ = 3.696 N m/rad, d = 0.05 m,
c = 0.3 m, ρ = 1.225 kg/m3, ν = 1.5 × 10−5 m/s2.
The elastic axis EA is located at 40 % of the profile.

Figures 3–6 compare the profile motion numerically
simulated by FVM and FEM. The angle of rotation
ϕ [°] and the vertical displacement h [mm] of the profile
in dependence on time t [s] are presented for upstream
flow velocities U∞ = {5, 10, 15, 20, 25, 35, 40}m/s (the
resulting Reynolds numbers are in the range from 105

109



Acta Polytechnica Vol. 52 No. 6/2012

to 8 × 105). The initial conditions h(0) = −20 mm,
ḣ(0) = 0, ϕ(0) = 6°, ϕ̇(0) = 0 were used. The vibra-
tions 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 numeri-
cal 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.
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 f1 ≈ 5.3 Hz refers to the vertical
translation of the profile, and the higher frequency
f2 ≈ 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 (f1 ≈ 5.5 Hz, f2 ≈ 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).

6 Conclusions

Aeroelastic instability due to divergence appeared
prior to flutter in both numerical approaches. The re-
sults 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 repre-
sent the aeroelastic behaviour of the system. FVM
and FEM results differ markedly when the upstream
flow velocity reaches 40 m/s. The system is unsta-
ble 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 com-
puter 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 (approxi-
mately 1.2 GB of memory).

Acknowledgements
This work has been supported by the Research Plan
of the Ministry of Education of the Czech Republic
No 6840770003 and by grants of the Czech Science
Foundation No. P101/11/0207, P101/12/1271.

References appear after the figures on page 115.

110



Acta Polytechnica Vol. 52 No. 6/2012

-20
-15
-10
-5
 0
 5

 10
 15
 20

 0  0.5  1  1.5  2

h
[m

m
]

t[s]

(a) Vertical displacement (U∞ = 5 m/s)

-6

-4

-2

 0

 2

 4

 6

 0  0.5  1  1.5  2

ϕ
[d

e
g

]

t[s]

(b) Angle of rotation (U∞ = 5 m/s)

-20
-15
-10

-5
 0
 5

 10
 15
 20

 0  0.5  1  1.5  2

h
[m

m
]

t[s]

(c) Vertical displacement (U∞ = 10 m/s)

-6

-4

-2

 0

 2

 4

 6

 0  0.5  1  1.5  2

ϕ
[d

e
g

]

t[s]

(d) Angle of rotation (U∞ = 10 m/s)

Figure 3: Vertical displacement h [mm] and angle of rotation ϕ [°] in dependence on time t [s] for U∞ = 5
(Re = 105, above) and 10 m/s (Re = 2 × 105, below). Solid line — FVM results, dashed line — FEM results.

-20
-15
-10

-5
 0
 5

 10
 15
 20

 0  0.2  0.4  0.6  0.8  1

h
[m

m
]

t[s]

(a) Vertical displacement (U∞ = 15 m/s)

-6

-4

-2

 0

 2

 4

 6

 0  0.2  0.4  0.6  0.8  1

ϕ
[d

e
g

]

t[s]

(b) Angle of rotation (U∞ = 15 m/s)

-20
-15
-10
-5
 0
 5

 10
 15
 20

 0  0.2  0.4  0.6  0.8  1

h
[m

m
]

t[s]

(c) Vertical displacement (U∞ = 20 m/s)

-6

-4

-2

 0

 2

 4

 6

 0  0.2  0.4  0.6  0.8  1

ϕ
[d

e
g

]

t[s]

(d) Angle of rotation (U∞ = 20 m/s)

Figure 4: Vertical displacement h [mm] and angle of rotation ϕ [°] in dependence on time t [s] for U∞ = 15
(Re = 3 × 105, above) and 20 m/s (Re = 4 × 105, below). Solid line — FVM results, dashed line — FEM results.

111



Acta Polytechnica Vol. 52 No. 6/2012

-25
-20
-15
-10

-5
 0
 5

 10
 15
 20

 0  0.2  0.4  0.6  0.8  1

h
[m

m
]

t[s]

(a) Vertical displacement (U∞ = 25 m/s)

-4

-2

 0

 2

 4

 6

 0  0.2  0.4  0.6  0.8  1

ϕ
[d

e
g

]

t[s]

(b) Angle of rotation (U∞ = 25 m/s)

-30
-25
-20
-15
-10

-5
 0
 5

 0  0.2  0.4  0.6  0.8  1

h
[m

m
]

t[s]

(c) Vertical displacement (U∞ = 35 m/s)

-2
-1
 0
 1
 2
 3
 4
 5
 6

 0  0.2  0.4  0.6  0.8  1

ϕ
[d

e
g

]

t[s]

(d) Angle of rotation (U∞ = 35 m/s)

Figure 5: Vertical displacement h [mm] and angle of rotation ϕ [°] in dependence on time t [s] for U∞ = 25
(Re = 5 × 105, above) and 35 m/s (Re = 7 × 105, below). Solid line — FVM results, dashed line — FEM results.

-45

-40

-35

-30

-25

-20

 0  0.2  0.4  0.6  0.8  1

h
[m

m
]

t[s]

(a) FVM

-180
-160
-140
-120
-100

-80
-60
-40
-20

 0

 0  0.2  0.4  0.6  0.8  1

h
[m

m
]

t[s]

(b) FEM

 2

 3

 4

 5

 6

 0  0.2  0.4  0.6  0.8  1

ϕ
[d

e
g

]

t[s]

(c) FVM

-2
 0
 2
 4
 6
 8

 10
 12
 14

 0  0.2  0.4  0.6  0.8  1

ϕ
[d

e
g

]

t[s]

(d) FEM

Figure 6: Vertical displacement h [mm] (above) and angle of rotation ϕ [°] (below) in dependence on time t [s]
for U∞ = 40 m/s (Re = 8 × 105).

112



Acta Polytechnica Vol. 52 No. 6/2012

t[s]

h
[m

m
]

0 0.2 0.4 0.6 0.8 1
-70
-60
-50
-40
-30
-20
-10

0
10

(a) Angle of rotation

t[s]

[d
e
g

]

0 0.2 0.4 0.6 0.8 1
-4
-2
0
2
4
6
8

10

ϕ

(b) Vertical displacement

Figure 7: Vertical displacement h [mm] and angle of rotation ϕ [°] in dependence on time t [s] for U∞ = 42.5 m/s
(Re = 8.5 × 105) computed using FVM.

 0
 0.5

 1
 1.5

 2
 2.5

 0  5  10  15  20

|G
(h

)|

frequency [Hz]
(a) Angle of rotation

 0

 0.15

 0.3

 0.45

 0  5  10  15  20

|G
(ϕ

)|

frequency [Hz]
(b) Vertical displacement

 0
 0.5

 1
 1.5

 2
 2.5

 0  5  10  15  20

|G
(h

)|

frequency [Hz]
(c) Angle of rotation

 0
 0.05

 0.1
 0.15

 0.2
 0.25

 0.3
 0.35

 0.4

 0  5  10  15  20

|G
(ϕ

)|

frequency [Hz]
(d) Vertical displacement

Figure 8: Frequency spectrum analysis of vertical displacement h [mm] and angle of rotation ϕ [°] signals for
U∞ = 5 (above) and 10 m/s (below). Solid line — FVM results, dashed line — FEM results.

113



Acta Polytechnica Vol. 52 No. 6/2012

 0
 0.5

 1
 1.5

 2
 2.5

 0  5  10  15  20

|G
(h

)|

frequency [Hz]
(a) Angle of rotation

 0

 0.15

 0.3

 0.45

 0  5  10  15  20

|G
(ϕ

)|

frequency [Hz]
(b) Vertical displacement

 0
 0.2
 0.4
 0.6
 0.8

 1
 1.2
 1.4
 1.6
 1.8

 0  5  10  15  20

|G
(h

)|

frequency [Hz]
(c) Angle of rotation

 0
 0.05

 0.1
 0.15

 0.2
 0.25

 0.3
 0.35

 0.4
 0.45

 0  5  10  15  20

|G
(ϕ

)|

frequency [Hz]
(d) Vertical displacement

Figure 9: Frequency spectrum analysis of vertical displacement h [mm] and angle of rotation ϕ [°] signals for
U∞ = 15 (above) and 20 m/s (below). Solid line — FVM results, dashed line — FEM results.

 0

 0.5

 1

 1.5

 2

 0  5  10  15  20

|G
(h

)|

frequency [Hz]
(a) Angle of rotation

 0

 0.15

 0.3

 0.45

 0  5  10  15  20

|G
(ϕ

)|

frequency [Hz]
(b) Vertical displacement

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.

114



Acta Polytechnica Vol. 52 No. 6/2012

References
[1] Bisplinghoff R. L., Ashley H., Halfman R. L.:

Aeroelasticity. New York : Dover, 1996.
[2] Bruneau Ch.-H., Fabrie P.: Effective downstream

boundary conditions for incompressible Navier-
Stokes equations, Int. J. Numer. Methods Fluids,
19, 8, 1994, 693–705.

[3] Chorin A. J.: A numerical method for solving
incompressible viscous flow problems, J. Comput.
Phys., 2, 1, 1967, 12–26.

[4] Chung K. W., He Y. B., Lee B. H. K: Bifurcation
analysis of a two-degree-of-freedom aeroelastic
system with hysteresis structural nonlinearity by
a perturbation-incremental method, J. Sound
Vib., 320, 2009, 163–183.

[5] Codina R.: Stabilization of incompressibility and
convection through orthogonal sub-scales in finite
element methods, Comput. Methods Appl. Mech.
Eng., 190, 2000, 1579–1599.

[6] Davis T. A., Duff I. S.: A combined
unifrontal/multifrontal method for unsymmet-
ric sparse matrices, ACM Trans. Math. Softw.,
25, 1999, 1–19.

[7] Dessi D., Mastroddi F.: A nonlinear analysis of
stability and gust response of aeroelastic systems,
J. Fluids Struct., 24, 2008, 436–445.

[8] Dolejší V.: Anisotropic mesh adaptation tech-
nique for viscous flow simulation, East-West J.
Numer. Math., 9, 2001, 1–24.

[9] Dowell E. H.: A modern course in aeroelasticity.
Dodrecht : Kluwer Academic Publisher, 1995.

[10] Feistauer, M., 1993. Mathematical Methods in
Fluid Dynamics. Longman Scientific & Technical,
Harlow.

[11] Feistauer, M., Felcman, J., Straškraba, J., 2003.
Mathematical and Computational Methods for
Compressible Flow. Clarendon Press, Oxford.

[12] Gaitonde A. L.: A dual-time method for two-
dimensional unsteady incompressible flow calcu-
lations, Int. J. Numer. Methods Eng., 41, 1998,
1153–1166.

[13] Gelhard T., Lube G., Olshanskii M. A., Star-
cke J.-H.: Stabilized finite element schemes with
LBB-stable elements for incompressible flows, J.
Comput. Appl. Math., 177, 2005, 243–267.

[14] Heywood J. G., Rannacher R., Turek S.: Artifi-
cial boundaries and flux and pressure conditions
for the incompressible Navier-Stokes equations,
Int. J. Numer. Math. Fluids, 22, 1992, 325–352.

[15] Honzátko R.: Numerical Simulations of Incom-
pressible Flows with Dynamical and Aeroelastic
Effects. Prague : PhD dissertation, Czech Tech-
nical University in Prague, Faculty of Nuclear
Sciences and Physical Engineering, 2007.

[16] Honzátko R., Horáček J., Kozel K.: Numerical
solution of flow induced vibrating of a profile. In:
Numerical Mathematics and Advanced Applica-
tions, ENUMATH 2007, Heidelberg : Springer,
2008, 547–554.

[17] Horáček J.: Nonlinear formulation of oscillations
of a profile for aero-hydroelastic computations,
In Dynamics of Machines, Prague : Institute of
Thermomechanics AS CR , 2003, 51–56.

[18] Jones D. P., Roberts I., Gaitonde A. L.: Iden-
tification of limit cycles for piecewise nonlinear
aeroelastic systems, J. Fluids Struct., 23, 2007,
1012–1028.

[19] Koobus B., Farhat Ch.: Second-order time-
accurate and geometrically conservative implicit
schemes for flow computations on unstructured
dynamic meshes, Comput. Methods Appl. Mech.
Eng., 170, 1999, 103–129.

[20] Lesoinne M., Farhat Ch.: Geometric conservation
laws for flow problems with moving boundaries
and deformable meshes, and their impact on
aeroelastic computations, Comput. Methods Appl.
Mech. Eng., 134, 1996, 71–90.

[21] Librescu L., Marzocca P.: Advances in the lin-
ear/nonlinear control of aeroelastic structural
systems, Acta Mech., 178, 2005, 147–186.

[22] Librescu L., Marzocca P., Silva W. A.: Aeroe-
lasticity of 2-d lifting surfaces with time-delayed
feedback control. J. Fluids Struct., 20, 2, 2005,
197–215.

[23] Naudasher E., Rockwell D.: Flow-induced vibra-
tions. Rotterdam : A. A. Balkema, 1994.

[24] Paidoussis M. P.: Slender Structures and Axial
Flow, volume 1 of Fluid-Structure Interactions.
San Diego : Academic Press, 1998.

[25] Raviart P.-A., Girault V.: Finite Element Meth-
ods for the Navier-Stokes Equations. Berlin :
Springer, 1986.

[26] R. L. Sani and P. M. Gresho.: Incompressible
Flow and the Finite Element Method. Chichester :
Wiley, 2000.

[27] Sváček P., Feistauer, M.: Application of a Sta-
bilized FEM to Problems of Aeroelasticity, In
Numerical Mathematics and Advanced Applica-
tion. Berlin: Springer, 2004, 796–805.

[28] Sváček P., Feistauer M., Horáček J.: Numerical
simulation of flow induced airfoil vibrations with
large amplitudes, J. Fluids Struct., 23, 3, 2007,
391–411.

[29] Verfürth R.: Error estimates for mixed finite
element approximation of the Stokes equations,
R.A.I.R.O. Analyse numérique/Numerical anal-
ysis,18, 1984, 175–182.

115