Finite element method application for turbulent and transitional flow

This paper is interested in numerical simulations of the interaction of the fluid flow with an airfoil. Particularly, the problem of the turbulent flow around the airfoil with elastic support is considered. The main attention is paid to the numerical approximation of the flow problem using the finite element approximations. The laminar turbulence transition of the flow on the surface airfoil is considered. The chois of the transition model is discussed. The transition model based on the two equation k−ω turbulence model is used. The structure motion is described with the aid of two degrees of freedom. The motion of the computational domain is treated with the aid of the arbitrary Lagrangian-Eulerian method. Numerical results are shown.


Introduction
Recently the mathematical modelling and numerical approximation play important role in the engineering practice in the civil, aerospace and mechanical engineering (see e.g.[1], [2]).The mutual interactions of air flow and vibrating structure is the main subject of aero-elasticity, where the flow-induced vibrations may affect negatively the operation and stability of aircrafts, bladed machines, bridges, and many other structures.Typically, the main goal of aeroelasticity is the prediction of the bounds of the structure stability, see, e.g. the monographs [1], [3].The mathematical simulation of fluid and structure interaction requires to consider the viscous, usually turbulent flow, changes of the flow domain in time, nonlinear behaviour of the elastic structure and to solve simultaneously the evolution systems for the fluid flow and for the oscillating structure.The changes fluid domain cannot be neglected and a moving mesh method ( [4], [5]) must be employed.
The subject of our attention is the numerical analysis of the interaction of viscous turbulent flow with a vibrating airfoil, where the influence of the transition on the aeroelastic response of the airfoil is studied.In order to include the transition model into aeroelastic simulations and avoid using the non-local operations, the transition model based on two transport equations for intermittency and the onset momentum-thickness Reynolds number were used, see

Mathematical description
The mathematical formulation of the problem consists of the flow model, the structure model and the interface conditions.

Fluid flow
The fluid flow is described in the two-dimensional time dependent computational domain Ω t ⊂ R 2 with the Lipschitz continuous boundary Here, Γ D denotes the inlet part of the computational domain, Γ O denotes the outlet part and Γ Wt denotes the surface of the airfoil at time t.The Arbitrary Lagrangian-Eulerian (ALE) method is used to treat the moving domain, see [10], [11].The differentiable ALE mapping A = A(ξ, t) = A t (ξ) is assumed to be defined for all t ∈ (0, T ) and ξ ∈ Ω re f 0 = Ω 0 .Then the domain velocity w D reads w D (x, t) = ∂ t A(ξ, t) for any x = A(ξ, t) and the derivative with respect to a fixed point of the reference configuration Ω re f 0 is called the ALE derivative denoted by D A /Dt, see [12].
The fluid motion in the domain Ω t is modelled using the Navier-Stokes system of equations in Ω t written in the ALE form (i = 1, 2) div u = 0, where u = (u 1 , u 2 ) is the fluid velocity vector, and the stress tensor σ with components σ i j is given by where S = S(u) = (∇u + ∇ T u)/2 is the symmetric part of the gradient with the components S i, j = ( ∂u i ∂x j + ∂u j ∂x i )/2, p is the kinematic pressure (i.e., the pressure divided by the constant fluid density ρ), ν e f f = ν + γν T is the effective kinematic viscosity, ν is the constant fluid kinematic viscosity, ν T is a turbulent viscosity obtained by an additional model, see [8], and γ is the intermittency coefficient.Further, the system (1) is equipped with an initial condition u(x, 0) = u 0 (x), x ∈ Ω 0 .and with the boundary conditions where α − = min(0, α) denotes the negative part of the number α ∈ R and n = (n 1 , n 2 ) denotes the unit outward normal to ∂Ω t .

Structure motion
The flow is interacting with a flexibly supported airfoil, which is allowed to be vertically displaced by h (downwards positive) and rotated by the angle α (clockwise positive).Its motion is then described by the nonlinear equations of motion where m is the mass of the airfoil, S α is the static moment around the elastic axis (EA), and I α is the inertia moment around EA, see [12].The parameters k h and k α denote the bending and torsional spring stiffness coefficients, respectively.

Interface conditions
The aerodynamical lift force L(t) is given by and the aerodynamical torsional moment M(t) (clockwise positive) is given by where l denotes the depth of the airfoil section, δ i j is the Kronecker's delta, is the position of EA of the airfoil at the time instant t.

Numerical approximation
In order to discretize the problem in time the equidistant partition t k = kΔt of the time interval I is considered with a time step Δt > 0. We denote the approximations u k ≈ u(•, t k ) and p k ≈ p(•, t k ).Moreover we approximate the domain velocity w D at time level t k by w k D and focus on the description of the discretization at an arbitrary fixed time instant t = t n+1 .For the sake of simplicity we shall omit the subscripts t or t n+1 in what follows.We shall consider all the function spaces X, W , Q defined for the time instant t = t n+1 on the domain Ω := Ω t n+1 .Then the ALE time derivative in the weak formulation of (1) is approximated at the time t = t n+1 by the second order backward difference formula, i.e.
where by ũi = u i • A t i • A −1 t n+1 the transformation of u i from Ω t i on Ω for i = n − 1 and i = n is denoted.
In order to spatially discretize the problem (1) by the finite element method, the spaces X, W and Q are approximated by finite element subspaces X Δ , W Δ and Q Δ , respectively.The Taylor-Hood family of finite element spaces defined over a triangulation T Δ of the computational domain Ω = Ω t n+1 are used, i.e. the continuous piecewise quadratic velocities and the continuous piecewise linear pressures are used.Moreover for the involved high Reynolds numbers the fully stabilized scheme is used, which consists of streamline-upwind/Petrov-Galerkin(SUPG) and pressurestabilizing/Petrov-Galerkin(PSPG) stabilization combined with the div-div stabilization, see [13].The non-linear problem is then solved using the Oseen linearization process, where the solution of each of the linear problems is performed using the efficient direct solver, see [14].  1) is replaced by the effective viscosity ν eff = ν + ν T , and the turbulent viscosity ν T is modelled using the turbulent kinetic energy k = k(x, t) and the turbulent specific dissipation rate ω = ω(x, t) satisfying for any t ∈ (0, T ) in Ω t equations
The effective intermittency is then modelled using the equation for the intermittency coefficient γ written in the ALE form where P γ and E γ are the transition source and destruction terms given by P γ = F length c a1 S γF onset (1 − c e1 γ), and , where S and Ω are the strain rate and vorticity magnitudes, respectively.Further, The following constants for the intermittency equation were used Further, Re θc is the critical Reynolds number given by an empirical correlation, and another empirical correlation is used for the function F length , which controls the length of the transition region.The correlations are based on newly defined transported unknown Re θt governed by the equation written in the ALE form where the source term , where c θt = 0.03, σ θt = 2, t ∞ = 500ν/U 2 is the time scale, U is the magnitude of velocity U = u 2 and the blending function F θt is defined as with δ = (375Ωyθ)/U, θ = Re θt ν/U, F wake = e −(Re ω /10 5 ) 2 , and Re ω = ωy 2 /ν.The right hand side includes the Reynolds number Re θt = θ t ν/U given by the empirical correlation, see [16].The equations ( 7) and (8) are equipped with the Dirichlet boundary conditions (γ = 1, Re θt obtained by the previously specified correlation) at the inlet and Neummann boundary conditions at the outlet and on the airfoil surface (Γ O ∪ Γ Wt ).For the turbulence and transition models the SUPG stabilization is used and the resulting system of linear equations is solved numerically using the direct solver.

Numerical results
In the paper we present the numerical simulation of the coupled problem.The inlet velocity is taken as U ∞ = 40m/s (post-critical regime), the structural parameters are taken same as in [12].In that case the divergence type of instability occurs (with further increase of the far field velocity the flutter instability should be observed).Here, the we investigated the response of the system depending on the different values of the initial values of α and h.The comparison of the aeroelastic response is shown in Figures 1-2, where the divergence type of instability can be observed.
The instability type changes with the further increase of the inlet velocity in Figure 3, where flutter type of instability was observed.

Conclusion
The present paper describes the solution of FSI problem of flow induced vibrations of an elastically supported airfoil in the post-critical regime.The problem was discretized in time and in space using the finite element method.The main attention was paid to the the description of the transition model, originally proposed by [6].The numerical results were shown.
[6], [7].In the present paper we are concerned with the numerical simulation of the aeroelastic problem of 2D viscous incompressible flow past a moving airfoil.The main attention is paid to the description of the application of the k −ω turbulence model (see [8], [9]) together with the transition model included, see [6].
a) u = u D on Γ D , DOI: 10.1051/ C Owned by the authors, published by EDP Sciences,

For
modelling of the turbulent flow with the transition the Menter's SST k − ω turbulence model, see [15] is used, with the γ − Re θt transition model introduced by [6].The viscosity coefficient ν in the equation (

Fig. 1 .
Fig.1.The aeroelastic response of the airfoil NACA 0012.The inlet velocity is taken as U ∞ = 40m/s and the instability of type divergence can be observed.The initial condition is α(0) = 0.03 o and h(0) = 0 mm.

Fig. 2 .Fig. 3 .
Fig.2.The aeroelastic response of the airfoil NACA 0012.The inlet velocity is taken as U ∞ = 40m/s and the instability of type divergence can be observed.The initial condition is α(0) = 0.03 o and h(0) = 0 mm.
4, y denotes the wall distance, and Re θt is the transition Reynolds number and