On several aspects of finite element approximations of laminar/turbulent transition models

In this paper the modeling of the laminar/turbulent transition is considered. Particularly, the model based on an algebraic equation for the intermittency coefficient is used. The detailed mathematical description is given. The problem is numerically approximated by stabilized finite element method and the numerical results for flows over an airfoil are shown.


Introduction
In this paper, the numerical approximation of an incompressible transition flow is considered. Particularly, we are interested in the interaction of flowing fluids and vibrating structures. The main goal of aero-elasticity is the prediction of the bounds of the structure stability. In the considered problem, we are interested in post-critical behaviour of the aeroelastic system. In such a case the turbulence can play an important role. The mathematical simulation of fluid and structure interaction requires to consider 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. Considering the Reynolds averaged Navier-Stokes equations and a vibrating structure with large displacements, the change of the fluid domain cannot be neglected. The methods with moving meshes [1], [2] must be employed and the application of efficient and robust methods for the numerical solution is required.
The turbulence and turbulence modelling effects can be significant in aeroelastic computations. Srinivasan et al. [3] used 2D RANS equations model for approximation of flow around the oscillating airfoil NACA0015 in rotation together with five models of turbulence, and found, that a significant improvement was obtained when better turbulence model had been used. Poirel et al. [4] studied the low amplitude self-sustained pitch airfoil for oscillations in incompressible flow by 2D numerical simulations, where both laminar and URANS calculations using the SST k − ω model with a low-Reynolds-number correction have been performed. The paper by Wang and Zha [5] investigates the NLR7301 airfoil limit cycle oscillation (LCO) in transonic flow caused by the flow nelinearity of the fluid-structure interaction using detached eddy simulation (DES) of tura e-mail: petr.svacek@fs.cvut.cz bulence. The aeroelastic behaviour can be also influenced by the transition from laminar to turbulent flow. Nevertheless, the transition is a complicated process, where many sources needs to be counted for: particularly, it depends on the turbulence intensity of the free flow, pressure gradient, roughness of the surface, temperature of the fluid and curvature of the surface. The classical mathematical models of transition are based on integral method of boundary layer, e.g. [6] considered the mean flow in the transition zone as a linear combination of the laminar and turbulent boundary layer in proportions determined by the transitional intermittency. Similar approach was used also in [7]. Integral method of boundary layer was used for airfoil in [8], where a simple transition model based on e N method was included.
Recently, there are three main concepts how to include transition in the turbulence model. One is the application of low-Reynolds number turbulence models, but without any coupling to an intermittency equation it appears to be unreliable for predicting transition and they are not suitable for aerodynamic flows. The second approach is the e N method which uses the local, linear stability theory and the parallel flow assumption. For isolated airfoils, the e N method has been shown to produce very good transition predictions compared to wind tunnel measurements. The third approach to predicting transition is the use of experimental correlations, but its application in practical implementations is complicated. Here, we use a transition model based on solution of two transport equations for intermittency and the transition onset momentum-thickness Reynolds number, [9], [10]. This second transport equation is connected with the empirical correlation to the onset criteria in the intermittency equation and thus it allows to use in general geometries. The intermittency function is coupled with the k-ω turbulence model.
In the present paper we are concerned with a numerical simulation of 2D viscous incompressible turbulent 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. The main attention is paid to the modelling and approximation of the k−ω turbulence model ( [11], [12]), with a transition model included [9]. Particularly, the approximation of the RANS equations and the approximation of the turbulence and transition models by algebraic flux correction scheme are discussed.

Governing equations
Viscous incompressible flow is described by the velocity u = u(x, t) and the kinematic pressure p = p(x, t) depending on x ∈ Ω and t ∈ [0, T ]. The density of the fluid ρ is assumed to be constant. The character of the flow depends on the magnitude of the Reynolds number Re = U ∞ c/ν , where ν is the kinematic viscosity, U ∞ denotes the far field velocity and c is the length of the airfoil chord. The turbulent flow is characterized by the fact that the fluid velocity field varies significantly and irregularly both in position and in time. The turbulence is a complicated motion, which results from the nonlinear advection that creates interactions between different scales of motion, which are the principal current (or the large eddies) and the eddying, random and reverse fluctuations. There are several strategies for the modelling of turbulent flow. For main concepts see, [13], [14], [11]. One possibility is to use the Reynolds decomposition of the flow velocity u and the kinematic pressure p into mean and fluctuating parts.

Reynolds equations
In what follows we shall be concerned with the approximation of two-dimensional incompressible turbulent flow in time-dependent computational domain Ω t ⊂ R 2 , t ∈ [0, T ]. In order to treat the time dependence of the domain occupied by fluid, the Arbitrary Lagrangian-Eulerian (ALE) method will be used. We assume that there exists a smooth one-to-one transformation A t (ALE mapping) of a reference computational domain Ω 0 onto Ω t for any t ∈ [0, T ]. Further, by w D (x, t) the (ALE) domain velocity and by D A /Dt the ALE derivative shall be denoted, cf. [1], [15], [16]. The flow is modeled in Ω t by the incompressible RANS equations written in the ALE form: t)) denotes the mean part of the velocity vector, p = p(x, t) denotes the mean part of the kinematic pressure (i.e. pressure divided by the constant fluid density ρ), ν e f f = ν + ν T , ν denotes the kinematic viscosity, ν T is the turbulent viscosity, w = u − w D , w D is the domain velocity and S(v) = 1 2 ∇u + ∇ T u . System (1) is equipped with initial and boundary conditions, see [15].

γ − Re θ transition model
Here, the transition from laminar to turbulence regimes is modelled with the aid of a Local-Correlation-based Transport model, cf. [9]. The intermittency equation (in ALE form for incompressible flow) is formulated as follows where γ is the intermittency, and P γ and E γ are the transition source and destruction terms, respectively. The transition source term is defined as where S and Ω are the strain rate and vorticity magnitudes, respectively. Further, F onset = max(F onset2 − F onset3 , 0), where y is the nearest wall distance, Re V is the vorticity Reynolds number, Re θc is the critical Reynolds number and Re θt is the transition Reynolds number. Particularly, Re θc can be thought of as the location where turbulence starts to grow, while Re θt is the location where the velocity profile first starts to deviate from the purely laminar profile. The connection between the two must be obtained from an empirical correlation. Another empirical correlation needs to be obtained for F length , which controls the length of transition region. For details see [9] or [10]. Both correlations Re θc and F length are obtained as a function of Re θt . The following constants for the intermittency equation were used c e1 = 1, c a1 = 2, c e2 = 50, c a2 = 0.06, σ f = 1.

Empirical correlation
Here, the following empirical correlation was used for relation of the transition onset momentum thickness Reynolds number Re θt , see also [9]. The following parameters were used ∂u is the acceleration in the streamwise direction. Then the Re θt is defined for T u < 1.3 by and for T u > 1.3 by Here, the function F(λ θ ) is defined for λ θ ≤ 0 by 3 Numerical approximation

Time discretization
In order to discretize the problem, an equidistant partition 0 = t 0 < t 1 < · · · < T , t k = kΔt of the time interval [0, T ] with a constant time step Δt is considered. The velocity and pressure are approximated at each time level by u(t n ) ≈ u n , p(t n ) ≈ p n . Similarly k n , ω n and ν n T are approximations of k(t n ), ω(t n ) and ν T (t n ), respectively. Further w n D is the approximation of the domain velocity w D at time t n , and we set w n+1 = u n+1 − w n+1 D . For simplicity, the ALE (time) derivative is approximated by backward Euler formula where we use the notationû n is obtained by the ALE transformation of u n onto Ω := Ω t n+1 . Similarly, the ALE derivative of functions k and ω are approximated by wherek n andω n are transformations of functions k n , ω n onto domain Ω t n+1 . Similarly, the ALE time derivatives are approximated also for the intermittency and for transition momentum thickness Reynolds number, i.e.

Space discretization of RANS equations
For the appproximation of the incompressible RANS equations the finite element velocity/pressure pair, which satisfy the Babuška-Breezi condition (cf. [17]), is chosen. Due to the dominating convection requires we also introduce some additional stabilization. In practical computations we assume that the domain Ω is a polygonal approximation of the region occupied by the fluid at time t n+1 and T h is a regular triangulation in Ω. The fluid velocity and pressure are sought in the finite element space X h of continuous piecewise quadratic (vector) functions and the finite element space Q h of continuous piecewise linear functions, respectively.
The fully stabilized problem reads: Find U = (u, p) ∈ X h × Q h such that u satisfies approximately the Dirichlet boundary conditions and the identity (7) for all V = (ϕ, q) ∈ X Δ × Q h . Here, the Galerkin terms are defined by and the SUPG/PSPG stabilizing terms are defined by where ψ(U; V) = w n+1 · ∇ ϕ + ∇q and the div-div stabilizing term The parameters τ K , δ K are chosen according to [18].

Space discretization of the k-ω turbulence model
The numerical approximation of the time discretized equations of the turbulence model is realized by the application of the SUPG stabilized finite element method, where the nonlinear terms of k − ω equations (3) are linearized. In order to avoid non-physical undershoots/overshoots of the approximations of k := k n+1 and ω := ω n+1 , the additional nonlinear crosswind diffusion method can be applied, cf. [19], [20], [21], but there is no guarantee of no undershoots/overshoots. Particularly for the solution of intermittency equation as well as equation for k and ω, these undershoots/overshoots can lead to incorrect results. In order to obtain better model, we introduce the idea of algebraic flux corrected transport.

Algebraic Flux Corrected Transport
The finite element approximation of the two equation turbulence model requires the use of more sophisticated stabilization as e.g. crosswind diffusion. Here, the algebraic flux corrections scheme is applied [22]. The description of the main idea is shown for the finite element approximation of the time dependent continuity equation ∂u ∂t + ∇ · (uu) = 0, i.e.
where u = (u i ) is vector of the nodal values, M C is the consistent mass matrix and K is the discrete transport operator. In order to obtain scheme both without undershoots/overshoots as well as not too diffusive, we need to switch between linear "upwind-like" approximations and the original scheme.
In the finite element context the idea of algebraic flux corrections reads: replace the consistent mass matrix M C by lumped mass matrix M L , and add an artificial diffusion operator D to operator K to eliminate all negative off-diagonal coefficients of K. The linear local extremum diminishing scheme then reads M Lu = Lu, L = K + D. The artificial diffusion operator D can be rewritten as The original scheme can be then recovered M Lu = Lu − Du + (M L − M C )u, or component by component where m i are coefficients of the lumped mass matrix. In order to prevent the oscillations of the solution, the fluxes f i j are multiplied by suitable correction factors Inserting these fluxes into (9) we get the nonlinear combination of the low order scheme (α i j = 0) and the original higher order scheme (α i j = 1). More detailed description can be found in [22].

Numerical results. Conclusion.
First, the developed flux corrected transport scheme was tested for finite element implementation in 1d and 2d (see  Figure 2 shows the exact solution (left) and the numerical solution(right). Further, the numerical results for approximation of laminar/turbulent flow are presented for the airfoil configurations considered in [23], where the authors computed the stability bounds of a wing profile model based on a linear description of the structure behaviour.
The numerical simulation was carried out for the airfoil NACA 0012 of length 0.3m. 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 considered a gap between the main lifting surface and the flap varied from g = 0% to g = 7% of the flap chord. Figures 3 show examples of the computed functions h(t), α(t), β(t), the corresponding spectra and the phase diagrams for k − ω turbulence model and for close-to-critical far-field flow velocity U ∞ = 11 m/s. The flow velocity patterns can be found in figure 4.