Printer Friendly

Interaction of incompressible flow and a moving airfoil.

1. Introduction. The interaction of fluid flow with vibrating bodies plays a significant role in many areas of engineering. We can mention, for example, development of airplanes or turbines and some problems from civil engineering; see, e.g., [6]. In the design of airplanes the point of interest is the analysis of deformations and vibrations of wings induced by flowing air. In this paper we are concerned with numerical solution of an aeroelastic problem of two dimensional viscous incompressible flow over an airfoil with two degrees of freedom in a wind tunnel. The airfoil is represented by a solid body, which can perform vertical and torsional vibrations. A mathematical model of the flow is formed by the system of twodimensional non-stationary Navier-Stokes equations and the continuity equation, equipped with initial and mixed boundary conditions.

Due to the moving airfoil, the computational domain is time-dependent. This requires the use of a suitable technique for the simulation on a moving computational grid. Here we apply the Arbitrary Lagrangian-Eulerian (ALE) method [11].

The flow problem is discretized by the stabilized finite element method, which leads to a large system of nonlinear algebraic equations. We overcome the nonlinearity by the use of the Oseen linearization, resulting in a sequence of linear systems of saddle-point type. They are solved by the direct solver UMFPACK [3]. The numerical results are compared with wind tunnel measurements.

2. Mathematical model. We consider two-dimensional non-stationary viscous incompressible flow past a vibrating airfoil inserted into a channel (wind tunnel) in a time interval [0, T], where T > 0. The symbol [[OMEGA].sub.t] denotes the computational domain occupied by the fluid at time t. The boundary [delta][[OMEGA].sub.t] of the domain [[OMEGA].sub.t] consists of mutually disjoint sets [[GAMMA].sub.D], [[GAMMA].sub.O] and [[GAMMA].sub.Wt] on which different types of boundary conditions are prescribed. By [[GAMMA].sub.D] we denote impermeable walls and the inlet, through which the fluid flows into the domain [[OMEGA].sub.t], [[GAMMA].sub.O] denotes the outlet, where the fluid flows out and [[GAMMA].sub.Wt] is the boundary of the profile at the time t. In contrast to [[GAMMA].sub.Wt], we assume that [[GAMMA].sub.D] and [[GAMMA].sub.O] are independent of time. The flow is characterized by the velocity field u = u(x, t), and the kinematic pressure p = p(x, t), for x [member of] [[OMEGA].sub.t] and t [member of] [0, T]. The kinematic pressure is defined as P/[rho], where P is the pressure and [rho] > 0 is the constant fluid density. The motion of the profile is described by functions [alpha](t), representing the rotation around the elastic axis TR, and H(t), denoting the vertical displacement; see Figure 2.1.


2.1. ALE formulation of the Navier-Stokes equations. The time dependent computational domain can be treated with the aid of a smooth, one-to-one ALE mapping [11]

(2.1) [A.sub.t] : [[OMEGA].sub.0] [??] [[OMEGA].sub.t], X [??] x (X,t) = [A.sub.t](X), t [member of] [0,T].

The coordinates of points x [member of] [[OMEGA].sub.t] are called spatial coordinates, the coordinates of points X [member of] [[OMEGA].sub.0] are called ALE coordinates or reference coordinates. The ALE mapping reflects the deformation of the computational domain; see Figure 2.2.

We define the domain velocity in the following way


This velocity can be expressed in spatial coordinates as

(2.3) w(x,t) = [??] ([A.sup.-1.sub.t](x),t).

Let us consider a function f = f (x, t), x [member of] [[OMEGA].sub.t], t [member of] [0, T], f (x, t) [member of] R, where R is the set of real numbers. Let f (X, t) = f ([A.sub.t] (X), t). We define the ALE derivative of the function f by



The application of the chain rule gives

(2.5) [D.sup.A]/[D.sub.t] f = [partial derivative]f/[partial derivative]f + w * [nabla] f.

By using this relation we can obtain the Navier-Stokes equations in the form


The symbol v denotes the kinematic viscosity of the fluid. We assume that v > 0 is constant.

2.2. Equations for the moving airfoil. The equations of airfoil motion are derived from the Lagrange equations for the generalized coordinates H and [alpha]. These equations have the form [6]

(2.7) Kd(t) + Dd(t) + Md(t) = f (t),

where the stiffness matrix K, the viscous damping D and the mass matrix M have the form


The force vector f and the generalized coordinates d read


The symbol F stands for the component of the aerodynamic force acting on the profile in the vertical direction, M is the torsional aerodynamic moment with respect to the elastic axis, [D.sub.HH], [D.sub.H[alpha]], [D.sub.[alpha]H], [D.sub.[alpha][alpha]], are the coefficients of structural damping, [S.sub.[alpha]], [I.sub.[alpha]], m and [k.sub.HH], [k.sub.[alpha][alpha]], [k.sub.H[alpha]], [k.sub.[alpha]H] denote the static moment around the elastic axis TR, the moment of inertia around TR, the mass of the profile and the stiffnesses of the profile elastic support. The force F and moment M acting on the airfoil are given by the relations



Here l is the airfoil depth in the direction orthogonal to the plane [x.sub.1], [x.sub.2], representing the length of a wing segment in consideration. Further, n is the unit outer normal to [partial derivative][[OMEGA].sub.t] on [[GAMMA].sub.Wt], [[delta].sub.ij] is the Kronecker symbol, i.e.,

[[delta].sub.ij] = 1 for i = j, and [[delta].sub.ij] = 0 for i [not equal to] j,

[x.sub.1], [x.sub.2] are the coordinates of points on [[GAMMA].sub.Wt], [x.sup.TR.sub.i], i = 1, 2, are the coordinates of the elastic axis [x.sup.TR], and


2.3. Initial and boundary conditions. The Navier-Stokes equations are completed by the initial condition

(2.11) u(x, 0) = [u.sub.0], x [member of] [[OMEGA].sub.0],

and the following boundary conditions. On [[GAMMA].sub.D] we prescribe the Dirichlet boundary condition

(2.12) u|[[GAMMA].sub.D] = [u.sub.D].

On the outlet [[GAMMA].sub.O] we consider the so-called do-nothing boundary condition

(2.13) -(p-[p.sub.ref]) n + v [partial derivative]u/[partial derivative]n = 0 on [GAMMA].sub.O],

where [p.sub.ref] is a given reference pressure. On [[GAMMA].sub.Wt] we consider the condition

(2.14) u|[GAMMA].sub.Wt] = w|[[GAMMA].sub.Wt].

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


where [[alpha].sub.0], [[alpha].sub.1], [H.sub.0], [H.sub.1] are input parameters of the model. The initial value problem (2.7), (2.15) is transformed to a problem for a first-order system and then discretized by the fourthorder Runge-Kutta method.

The interaction of a fluid and an airfoil consists in the solution of the flow problem (2.6), (2.11)-(2.14) coupled with the structural model (2.7) and (2.15). In what follows, we shall be concerned with the discretization of the flow problem and describe the algorithm for the numerical solution of the complete fluid-structure interaction problem.

3. Discretization of the flow problem.

3.1. Discretization in time. We construct an equidistant partition of the time interval [0, T], formed by time instants 0 = [t.sub.0] < [t.sub.1] < ... < T, [t.sub.k] = k[tau], where [tau] > 0 is a time step. We use the approximation u([t.sub.n]) [approximately equal to] [u.sup.n], p([t.sub.n]) [approximately equal to] [p.sup.n] of the exact solution and w([t.sub.n]) [approximately equal to] [w.sup.n] of the domain velocity at time [t.sub.n]. On each time level [t.sub.n+1], under the assumption that the domain [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] is already known, using the second-order backward difference formula, we obtain the problem to find functions [u.sup.n+1]: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] [??] [R.sup.2] and [p.sup.n+1]: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] [??] R, such that


This system is considered with the boundary conditions (2.12), (2.13), and (2.14). The symbols [[??].sup.n] and [[??].sup.n-1] mean the functions [u.sup.n] and [u.sup.n-1] transformed from the domain [[OMEGA]] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] to the domain [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] using the ALE mapping. This means that [[??].sup.i] = [u.sup.i] o [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

3.2. Discretization in space. The starting point for the approximate solution is the weak formulation of problem (3.1). Because of simplicity we shall use the notation [OMEGA] = [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] The appropriate function spaces are


We introduce the forms


where U = (u, p) [member of] W x M, [U.sup.*] = ([u.sup.*], p) [member of] W x M, V = (v, q) [member of] X x M, and (*,*) denotes the scalar product in the spaces [L.sup.2]([OMEGA]), [[[L.sup.2]([OMEGA])].sup.2], and [[[L.sup.2](Q)].sup.2x2].

The weak solution is defined as a couple U = (u, p), such that it satisfies the conditions

(3.4) U [member of] W x M, a(U,U,V) = f(V), [for all]V = (v,q) [member of] X x M,

and u satisfies the boundary conditions (2.12) and (2.14).

Now we define an approximate solution. The spaces W X, M are approximated by their finite-dimensional subspaces [W.sub.h], [X.sub.h], [M.sub.h], h [member of] (0, [h.sub.0]), [h.sub.0] > 0, where


The approximate solutions is defined as a couple Uh = (uh,Ph) E Wh x Mh, such that


and [u.sub.h] satisfies a suitable approximation of the boundary conditions (2.12) and (2.14).

In the construction of the spaces [W.sub.h], [M.sub.h] we assume that the domain [OMEGA] is a polygonal approximation of the computational domain at time [t.sub.n+1]. By [tau]h (h [member of] (0, [h.sub.0])) we denote a triangulation of [OMEGA] with standard properties from the finite element method; see, e.g., [2]. Then [W.sub.h] and [M.sub.h] are defined as continuous piecewise polynomial functions satisfying the Babuska-Brezzi condition; see [1]. Here we use the well-known Taylor-Hood [P.sup.2]/[P.sup.1] elements over a triangulation [[tau].sub.h] of [OMEGA]. This means that the velocity components are continuous in [??], quadratic on each element K [member of] [T.sub.h], and the pressure is continuous piecewise linear.

3.3. Stabilization of the finite element method. For high Reynolds numbers approximate solutions can contain nonphysical spurious oscillations. In order to avoid them, we shall apply the streamline-diffusion and div-div stabilization based on the forms




U = (u,p) [U.sup.*] = ([u.sup.*],p) V = (v,q),

[[delta].sub.K], [[tau].sub.K] [greater than or equal to] 0 are suitable parameters, [??] = [u.sup.*] - [w.sup.n+1] is the transport velocity, and [(*,*).sub.K] is the scalar product in the space [L.sup.2] (K) or [[[L.sup.2](K)].sup.2].

The solution of the stabilized discrete problem is [U.sub.h] = ([U.sub.h], Ph) [member of] [W.sub.h] x [M.sub.h], such that the component [u.sub.h] satisfies the boundary conditions (2.12) on [[GAMMA].sub.D] and (2.14) on [[GAMMA].sub.W], and


If we solve problem (3.8), we obtain an approximate solution at time [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] and [p.sup.n+1.sub.h] = ph defined in the domain [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

Now, we describe how to choose parameters [[delta].sub.K] and [[tau].sub.K]. We follow the works [7] and [10]. The magnitude of the velocity field varies in different subdomains of [OMEGA]. That is why we split the domain into two subdomains. The diffusion component dominates on the first subdomain and the convective component on the second. On both subdomains we choose these parameters in a different way. The parameter [[delta].sub.K] is based on the transport velocity [??] and the viscosity v. We put




is the so-called local Reynolds number and [h.sub.K] is the size of the element K measured in the direction of [??]. The function [xi](*) is non-decreasing in dependence on [R.sup.[??].sub.K] in such a way that for local convective dominance ([R.sup.[??].sub.K]) > 1) [xi] [right arrow] 1 and for local diffusion dominance (([R.sup.[??].sub.K] < 1) [xi] [right arrow] 0. The parameter [[delta].sup.*] [member of] (0,1] is chosen suitably. The function [xi](*) can be defined, e.g., by


The parameters [[tau].sub.K] are defined by


In practical computations we use the values [[delta].sup.*] = 0.025 and [[tau].sup.*] = 1.

4. Simulation of the flow induced airfoil vibrations. In the solution of the complete coupled fluid-structure interaction problem we apply the following algorithm:

1) Assume that the approximate solution U = ([u.sub.h],Ph) of the flow problem (3.8) at time levels [t.sub.n-1] and [t.sub.n] is known and the force F and torsional moment M are computed from (2.8) and (2.9).

2) Extrapolate F and M on the time interval [[t.sub.n], [t.sub.n+]].

3) Compute the displacement H and angle [alpha] at time [t.sub.n+1] as the solution of system (2.7).

4) Determine the position of the airfoil at time [t.sub.n+1], the domain [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] the ALE mapping [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] and the domain velocity [w.sup.n+1].

5) Solve the discrete stabilized problem (3.8) at time [t.sub.n+1].

6) Compute F and M from (2.8) and (2.9) at time [t.sub.n+1] and interpolate F and M on [[t.sub.n], [t.sub.n+1]].

7) Is a higher accuracy needed? YES: go to 3); NO: n := n + 1 and go to 2).

The nonlinear flow problem (3.8) is solved with the aid of the Oseen iterations


where [U.sup.(0).sub.h]) is defined on the basis of the approximate solution on the previous time level.

The ALE mapping is constructed in such a way that the reference domain [OMEGA].sub.0] is divided in three subdomains by two ellipses with center at the elastic axis of the airfoil. In the interior ellipse containing the airfoil the ALE mapping is defined as the rigid body motion, outside of the exterior ellipse the ALE mapping is identity. In the domain between both ellipses the ALE mapping is defined by interpolation.

The knowledge of the ALE mapping at time instants [t.sub.n-1], [t.sub.n], [t.sub.n+1] allows us to approximate the domain velocity with the aid of the second-order backward difference formula


5. Numerical solution. The described method was applied to the simulation of flow induced vibrations of a profile (shown in Figure 5.1) inserted into a channel (wind tunnel) in the case of the following data: v = 1.5 * [10.sup.-5] [m.sup.2]/s, [k.sub.HH] = 1711.6 N/m, [k.sub.[alpha][alpha]] = 4.5 Nm/rad, [k.sub.H[alpha]] = 0.0 N/rad, [k.sub.[alpha]H] = 0.0 N, m = 0.0821 kg, [S.sub.[alpha]], = -0.00013 kgm, [I.sub.[alpha]] = 0.000095 [kgm.sup.2], DHH = 5.0 Ns/m, [D.sub.[alpha][alpha]], = 0.003 Nms/rad, DHa = 0.0 Ns/rad, [D.sub.[alpha]H] = 0.0 Ns/m, l = 0.08 m, c = length of the profile chord = 0.12 m.



The computation was carried out for several values of the inlet velocity U. We define the corresponding Reynolds number by

(5.1) Re = Uc/v.

The triangulation of the domain is realized by the method and software ANGENER [4], [5], which can be used for the construction of an initial isotropic triangulation and also for an anisotropic adaptive mesh refinement; see Figure 5.1.

By the numerical solution of the complete problem we obtain the velocity and pressure fields and also the time development of the displacement H and the rotation angle a. From this information we derive frequency characteristics obtained with the aid of the Fourier transform


with g = H or g = [alpha] and [[??].sub.n] = n[DELTA][??] [member of] [0, 50], [DELTA][??] = 0.1 Hz, approximated by the rectangle formula


Here i is the imaginary unit and N is the number of time steps in the interval [0, T) with length [tau].

The main vibrational frequencies [f.sub.1] and [f.sub.2] are defined in our computations as maximum points of the function [absolute value of G] corresponding to g = H and g = [alpha], respectively.

The numerical simulation started by the computation of the flow velocity and pressure fields for a fixed profile in the disequilibrium position [[alpha].sub.0] = 2[degrees] and [H.sub.0] = 3.6 mm. After a short time the airfoil was released, i.e., the motion of the airfoil starts to behave according to system (2.7) with initial conditions formed by the above data [[alpha].sub.0], [H.sub.0], [[alpha].sub.1] = 0, and [H.sub.1] = 0; cf., (2.15).

In what follows we present the dynamic response [alpha](t), H(t) of the fluid-structure system in time domain for different inlet flow velocities.

Inlet velocity 0 [ms.sup.-1]--Re = 0. For [alpha] and H we obtained the signals shown in Figure 5.2. The frequency analysis gives main frequencies of the signals 19.5 Hz and 38.7 Hz. Inlet velocity 40 [ms.sup.-1]--Re = 320000. We obtained the results shown in Figure 5.3. The frequency analysis gives us main frequencies of the signals 19.9 Hz and 36.4 Hz.


Inlet velocity 60 [ms.sup.-1]--Re = 480000. We obtained the signals shown in Figure 5.4. The frequency analysis gives us main frequencies of the signals 19.8 Hz a 36.2 Hz.


Inlet velocity 80 [ms.sup.-1]--Re = 640000. The signals a and H are shown in Figure 5.5. The frequency analysis gives us only one main frequency of the signals 30.4 Hz. The second one is very hard to detect.

In all cases studied the vibration amplitudes decrease in time and the system is stable. Figure 5.6 shows the comparison of the computed main frequencies with wind tunnel experiments described in [8] and [9].

6. Conclusion. In this article we derived a procedure for obtaining the numerical solution of the interaction of a moving airfoil inserted in a channel with a running fluid. We used this approach for solving a particular problem, which was studied experimentally in a wind tunnel [8], [9]. The computational results show that the presented method is sufficiently robust and useful for a given type of problem. The computed and experimentally obtained main frequencies of the flow induced vibrations of the profile for several values of the inlet flow velocity are in good agreement; see Figure 5.6. A further goal will be the implementation of a turbulence model to the description of the flow and the validation of the method in situations with large displacements of the airfoil and a higher inlet flow velocity, when the system may loose aeroelastic stability.



Acknowledgments. This research was supported under the grant of the Grant Agency of the Academy of Sciences of the Czech Republic No. IAA200760613. It was also partly supported under the Research Plan MSM 0021620839 (M. Feistauer) of the Ministry of Ed ucation of the Czech Republic and the project No. 48607 of the Grant Agency of the Charles University in Prague (M. Ruzicka).

* Received November 27, 2007. Accepted for publication June 3, 2008. Published online on February 26, 2009. Recommended by O. Steinbach.


[1] F. BREZZI AND R. S. FALK, Stability of higher-order Hood-Taylor methods, SIAM J. Numer. Anal., 28 (1991), pp. 581-590.

[2] P. G. CIARLET, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1979.

[3] T. A. DAVIS, UMFPACK V4.0, University of Florida. Available at

[4] V. DOLEJSI, Anisotropic mesh adaptation technique for viscous flow simulation, East-West J. Numer. Math., 9 (2001), pp. 1-24.

[5] V. DOLEJSI, ANGENER V3.0, Faculty of Mathematics and Physics, Charles University Prague, Prague, 2000. Available at

[6] E. H. DOWELL, A Modern Course in Aeroelasticity, Kluwer Academic Publishers, Dodrecht, 1995.

[7] T. GELHARD, G. LUBE, M. A. OLSHANSKII, AND J.-H. STARCKE, Stabilized finite element schemes with LBB-stable elements for incompressible flows, J. Comput. Appl. Math., 177 (2005), pp. 243-267.

[8] J. HORACEK, J. KOZANEK, AND J. VESEL Y, Dynamic and stability properties of an aeroelastic model, in Engineering Mechanics 2005, V. Fuis, P. Krejcf, and T. Ndvrat, eds., May 9-12, 2005, Institute of Thermomechanics of the Academy of Sciences, Prague, 2005, pp. 121-122.

[9] J. HORACEK, M. LUXA, F. VANEK, J. VESEL Y, AND V. VLCEK, Design of an experimental set-up for the study of unsteady 2D aeroelastic phenomena by optical methods (in Czech), Research Report of the Institute of Thermomechanics of AS CR Prague, No. Z 1347/04, November 2004.

[10] G. LUBE, Stabilized Galerkin finite element methods for convection dominated and incompressible flow problems, Banach Center Pubt, 29 (1994), pp. 85-104.

[11] T. NOMURA AND T. J. R. HUGHES, An arbitrary Lagrangian-Eulerian finite element method for interaction of fluid and a rigid body, Comput. Methods Appl. Mech. Engrg., 95 (1992), pp. 115-138.

([dagger]) Charles University Prague, Faculty of Mathematics and Physics, Sokolovske 83,186 75 Praha 8, Czech Republic (,

([double dagger]) Institute of Thermomechanics of the Academy of Sciences of the Czech Republic, Dolejskova 5, 182 00 Praha 8, Czech Republic (

([section]) Czech Technical University Prague, Faculty of Mechanical Engineering, Karlovo n. 13, 121 35 Praha 2, Czech Republic (
COPYRIGHT 2008 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2008 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Ruzicka, Martin; Feistauer, Miloslav; Horacek, Jaromir; Svacek, Petr
Publication:Electronic Transactions on Numerical Analysis
Article Type:Report
Geographic Code:4EXCZ
Date:Jan 1, 2008
Previous Article:Local projection stabilization for incompressible flows: equal-order vs. inf-sup stable interpolation.
Next Article:On the evaluation of finite element sensitivities to nodal coordinates.

Related Articles
A Study of Airfoil Lift.
Physics, engineering & computer science.
Geometric theory of incompressible flows with applications to fluid dynamics.
Fluid dynamics; theoretical and computational approaches, 3d ed.
Convection--radiation interaction in buoyancy-induced channel flow.
Homogeneous turbulence dynamics.

Terms of use | Copyright © 2017 Farlex, Inc. | Feedback | For webmasters