Stabilized finite element methods in the inertial range: Monitoring artificial dissipation

c

TU Delft, The Netherlands, 2006

STABILIZED FINITE ELEMENT METHODS IN THE

INERTIAL RANGE: MONITORING ARTIFICIAL

DISSIPATION

Erik Burman∗

∗_{Institut d’Analyse et Calcul Scientifique (IACS/CMCS)} Station 8, Ecole Polytechnique Federale de Lausanne

CH–1015 Lausanne, Switzerland e-mail: erik.burman@epfl.ch web page: http://iacs.epfl.ch/ burman/

Key words: interior penalty, stabilization, finite element method

Abstract. In this note we will first discuss finite element discretizations of the incom-pressible Navier-Stokes equations with some artificial dissipation added in the form of either a turbulence model or some numerical stabilization term. We will then review some results on interior penalty stabilized finite element methods for incompressible flow. These results depend on the regularity of the solution and are typically relevant only when the mesh parameter h is sufficiently small.

In the case of high Reynolds number flow in three space dimensions it is unrealistic in most cases of industrial importance to resolve all scales and hence the computation must aim at the resolution only of the large eddies. In this case we propose to monitor the relative artificial dissipation and relate it to the power spectrum in the inertial range as a means of validation of large eddy simulations. Finally a numerical example is given on a two dimensional test case as an illustration to the discussion.

1 INTRODUCTION

The computation of flows at high Reynolds number remains a challenging problem. For industrial applications a full simulation of all the phenomena at hand must be expected to be out of reach for the next generations in spite of the rapid increase in computational power. The alternative is to attempt to compute only the large scales and model the effect of the unresolved fine scales on the large scales using a turbulence model, in a large eddy simulation (LES).

With the introduction of the variational multiscale method (VMS) by Hughes and coworkers [18] an attempt was made to apply ideas from stabilized finite element methods in an LES framework. The idea was to use the concept of subgrid viscosity stabilization introduced and analyzed by Guermond [15] for the computation of turbulent flows by adding the turbulent dissipation only to the finest scales.

Recently a number of new stabilized finite element methods have been proposed for Oseen’s equations or for the incompressible Navier-Stokes’ equations with the common feature that they may be interpreted as subgrid viscosity methods and hence may be cast in the variational multiscale framework. Artificial dissipation is added only on the finest scales. In some cases the motivation for the methods come from numerical analysis such as the orthogonal subscale method by Codina [12], the continuous interior penalty method of Burman and Hansbo [9, 10] or the local projection stabilization of Becker and Braack [2, 3] and in some the motivation is directly related to large eddy simulation, for instance in the subgrid viscosity method by John, Layton and Kaya [19, 25] or the variational multiscale method by Gravemeier, Ramm and Wall [14].

The main motivation of the present communication is to discuss, without claim to mathematical rigor, some connections between stabilized finite element methods and the physics of large eddy simulations. Particular focus will be given to the convergence rate of the artificial dissipation. We will consider convergence rates predicted by a priori error estimates for smooth solutions of the incompressible Navier-Stokes equations independent of the Reynolds number (Burman and Fern`andez [7]) and compare them with the conver-gence rates expected in case the stabilization term is considered a model for the energy content of the unresolved fine scales. In the latter case the convergence rate is obtained from the decay rate of the power spectrum in the inertial range.

2 THE NAVIER-STOKES EQUATIONS

In this paper we will mainly be concerned with the time-dependent incompressible Navier-Stokes with homogeneous boundary conditions

         ∂tu+ u · ∇u − ν∆u + ∇p = f in Ω × (0, T ), ∇ · u = 0 in Ω × (0, T ), u = 0 on ∂Ω × (0, T ), u(·, 0) = u0 in Ω. (1)

These equations describe the motion of a viscous incompressible fluid confined in Ω. In (1), ν > 0 corresponds to the kinematic fluid viscosity coefficient, f : Ω × (0, T ) −→ Rd

represents a given source term and u0 : Ω −→ Rd stands for the initial velocity.

The scalar product in L2_{(Ω) is denoted by (·, ·) and its norm by k · k. The scalar}

product on the boundary of Ω is denoted by h·, ·i . We denote by H1

0(Ω) and L20(Ω) closed

subspaces consisting respectively of functions in H1_{(Ω) with zero trace on ∂Ω, and of}

The weak formulation corresponding to equation (1) takes the form: Find (u, p) ∈ [H1 0(Ω)]d× L20(Ω),     

(∂tu, v) + c(u; u, v) + a(u, v) + b(p, v) = (f, v), a.e. in (0, T )

b(q, u) = 0, a.e. in (0, T ), u(0) = u0, a.e. in Ω,

(2)

for all (v, q) ∈ [H1

0(Ω)]d× L20(Ω), where

c(w; u, v)def= (w · ∇u, v), a(u, v)def= (ν∇u, ∇v), b(p, v)def= − (p, ∇ · v). 2.1 Finite element formulation

A general finite element formulation of problem (2) takes the form: Find (uh, ph) in

Wh := [Vh]d× Qh such that

(∂tuh, vh) + Auh; (uh, ph), (vh, qh) + Juh; (uh, ph), (vh, qh) = (f(t), vh), (3)

for all (vh, qh) ∈ Wh, equipped with the following initial condition

(uh(0), vh) = (u0, vh), ∀vh ∈ [Vh]d. (4)

The finite element spaces Vh and Qh may either be chosen so that the velocity/pressure

spaces satisfy the discrete inf-sup condition, or they may be chosen equal: Qh = Vh in

which case the stabilization term J must include a pressure stabilization term. 2.1.1 The standard Galerkin formulation

To formulate the standard Galerkin method we simply choose Awh; (uh, ph), (vh, qh)

def

= c(wh; uh, vh) + a(uh, vh)

+ b(ph, vh) − b(qh, uh).

(5) In this context J can be chosen either to be zero, if the intention is to resolve all scales so that the viscous dissipation dominates (|uh|h

ν < 1) or, in case a large eddy simulation

is aimed, it can be chosen as some turbulence model. For instance a Smagorinsky type model: Jwh; (uh, vh) def = (νT(wh)∇uh, ∇vh) with νT(wh) def = CSǫ2|∇Swh|

where ∇S denotes the symmetric part of the gradient tensor, CS is a parameter to choose

typically of the order 10−1 _{and ǫ is a characteristic lengthscale that usually is associated}

with the computational mesh.

Once again Wh must consist of an inf-sup stable velocity-pressure pair unless pressure

2.1.2 The interior penalty stabilized formulation

Another example of a finite element formulation for the Navier-Stokes equations on the form (2) is the following interior penalty formulation proposed in [7]. In this case we use spaces with equal order interpolation. Vhis the standard finite element space of continuous

piecewise polynomial functions of polynomial degree at most k and Qh = Vh. Boundary

conditions are not built into the velocity space but imposed weakly using a formulation due to Nitsche [27]. No turbulence model is added, but stability of the velocity/pressure coupling and stability at high Reynolds number are both assured by the addition of a term penalizing the jump of the gradient over element faces (denoted by J∇uhK). The

motivation for the weak boundary conditions and the penalization of the gradient jumps are:

1. For high Reynolds number flow with underresolved boundary layers the no-slip con-dition is automatically relaxed, making the simulation pass seamlessly from slip conditions to no-slip conditions as the layers are resolved (see also [13, 24, 6, 1, 8] for details).

2. For high Reynolds number flow the gradient jump stabilization can be shown to act only on the part of the discrete gradient that can not be represented by the finite element space, i.e. the part of the gradient in the unresolved scales.

The forms A and J of (3) are defined as follows. Awh; (uh, ph), (vh, qh) def = ch(wh; uh, vh) + ah(uh, vh) + bh(ph, vh) − bh(qh, uh), (6) ch(wh; uh, vh) def = c(wh; uh, vh) + 1 2(∇ · whuh, vh) − 1 2hwh· nuh, vhi , (7) ah(uh, vh) def = a(uh, vh) − hν∇uhn, vhi − huh, ν∇vhni +Dγν huh, vh E + huh· n, vh· ni , (8) bh(ph, vh) def = b(ph, vh) + hph, vh · ni , (9) Jwh; (uh, ph), (vh, qh) def = γS(jwuh(uh, vh) + j p wh(ph, qh)), (10) with ju wh(uh, vh) def =P K∈Th R ∂K\∂Ωh 2 K(|wh· n| + γCD)J∇uhK : J∇vhK ds, and jp wh(ph, qh) def =P K∈Th R ∂K\∂Ωmin(|wh| −1_,hK ν )h 2

KJ∇phK · J∇qhK ds. Here γS, and γCD are

more terms have to be evaluated and that the matrix stencil is extended because of the gradient jump terms. On the other hand this formulation allows for a fairly complete analysis as far as stability and convergence is concerned provided the exact solution is sufficiently regular. In particular the error estimates are independent of the viscosity and hence hold true also in the case of the incompressible Euler equations. A priori error estimates for this formulation was proved in [7] as given in the following theorem

Theorem 2.1 If we assume that u and p are sufficiently regular in Ω × [0, T ] and as-suming h > ν then the following error estimates hold

ku − uhk2L∞_{(0,T ;L}2_(Ω)) ≤ Ch2k+1 and Z T 0 kν12∇(u − u h)k2+ Juh; (uh, ph), (uh, ph) dt ≤ Ch2k+1

where C depends mainly on exp(k∇ukL∞_{(0,T ;L}∞_(Ω))T ) and Sobolev norms of the exact

solution u and p.

This means that if the solution is regular, in particular if the gradient of the exact solution is bounded the numerical method will converge. As an immediate consequence of the theorem we have

Z T

0

Juh; (uh, ph), (uh, ph) dt ≤ Ch2k+1

giving the rate of convergence of the artificial dissipation for laminar flows.

3 THE CURSE OF ASYMPTOTICS AND THE INERTIAL RANGE

Error estimates obtained in numerical analysis are typically of interest only in what is known as the asymptotic range. This means that the meshsize must be of the same order as the smallest scale of the exact solution. The constants in all these estimates depend on Sobolev norms of the exact solution. Since the exact solution is unknown, the estimate can only be of any use for error predictions in model cases where the solution is known to be smooth. In this section the aim is to make a thought experiment and compare the numerical analysis picture, with results from physics and propose a qualitative criterium for what behavior we should ask for from a turbulence model or a stabilized method. In the final section we will compare the prediction made by the error analysis with this quali-tative measure in a twodimensional numerical example which exhibits a simple multiscale structure and transition, but with a solution that is known to be smooth.

3.1 Constants depending on Sobolev norms

As already pointed out the weakness of the a priori analysis is the presence of constants depending on Sobolev norms of the exact solution. In particular kD2_uk

0, kDuk∞ (if

(ξ)

E

Log

Log

(ξ)

Integral

scale

Inertial

range

Dissipation

scale

ξ

Figure 1: In three dimensions energy is pushed to higher frequencies in turbulent flow, due to the vortex stretching mechanism giving rise to the well known energy cascade E(ξ) ∼ Cξ−

5

3.

The fact that a function has a certain regularity is equivalent to a certain decay rate of high frequency modes in Fourier space. If we let ˆu denote the Fourier transform of u this may be exemplified as

kD2uk2₀ < C ≡ Z

|ξ|4uˆ2dξ < C

where ˆu denotes the Fourier transform of u. Or in words: the Fourier transform of the function has to go to zero sufficiently fast as ξ → ∞. In one space dimension the H2

-regularity implies that there exists ξc such that

ˆ

u(ξ) ≤ Cξ−52−ε, with some ε > 0, for ξ > ξ

c

By definition in LES we have not resolved all scales, h−1 _{<< ξ}

c, so the numerical method

will never “see” this fast decrease (unless the flow is laminar), since the computational mesh can only represent the scales of u that are larger than h.

3.2 Different decay rates

In physics turbulence is often caracterized by the decay rate of its power spectrum. Most wellknown is the ξ−5

3 law for homogeneous turbulence derived by Kolmogorov [22],

• Homogeneous isotropic turbulence in three space dimensions (Kolmogorov, [22]): E(ξ) ∼ Cξ−53

• Driven homogeneous isotropic turbulence in two space dimensions (Kraichnan, [23]):

E(ξ) ∼ Cξ−53 (11)

• Decaying homogeneous isotropic turbulence in two space dimensions (Kraichnan, [23]):

E(ξ) ∼ Cξ−3 (12)

In large scale problems such as atmospheric turbulence decay rates may vary depending on the scales considered

• In the synoptic and subsynoptic range ( > 103 _{km): E(ξ) ∼ Cξ}−3

• In the mesoscales (< 103 _{km): E(ξ) ∼ Cξ}−5 3.

It seems reasonable to ask that a large eddy simulation of a complex flow respects the decay rate of the power spectrum of the inertial range of the problem at hand under refinement. In other words that the energy distribution between resolved and unresolved scales respects the physics.

3.3 Scale separation and the energy equality For the Navier-Stokes equations there holds

1 2ku(T )k 2₊ Z T 0 kν12∇uk2 dt = 1 2ku(0)k 2₊ Z T 0 (f , u) dt. (13)

Let πhu denote the L2-projection of the exact solution of the Navier-Stokes equations

projected onto the computational mesh,

(πhu− u, vh) = 0, ∀vh ∈ [Vh]d.

If we associate a wavenumber ξh with the mesh parameter h in such a way that

1 2kπhuk 2 ₌ Z ξh 0 E(ξ) dξ it then follows that the unresolved scales may be written

In the following we will assume that such a ξh exists and that it scales as ξh ∼ h−1. The

motivation for this is that in scaling theory the relation between eddy size and eddy wave number is ξλ ∼ λ−1and we expect the size of the smallest eddy resolved to be proportional

to the mesh size.

The aim now is to associate the unresolved scales to the energy spectrum in this ad hoc fashion and try to find a reasonable criterion for how stabilization terms should behave when ξh ∼ h−1 is in the inertial range.

3.4 The energy equality and the apriori error estimation picture

Let us first recall what we can hope for at best when using a (consistent) stabilized finite element method. Clearly by the consistency the artificial dissipation added must vanish under refinement, but what rate do we get from the a priori error estimate if we assume that the solution is regular? A relevant question is of course if it will vanish fast enough compared to the power spectrum in the inertial range.

Any Galerkin method (and most numerical methods) will satisfy 1 2kuh(T )k 2₊ Z T 0 n kν12∇u hk2+ S(uh, ph) o dt = 1 2kuh(0)k 2₊ Z T 0 (f , u) dt. (14) For some S(uh, ph) denoting the artificial dissipation added for the method to remain

stable. We will now consider two examples of artificial dissipation arising in stabilized finite element methods, first the interior penalty stabilization and then the classical SUPG stabilization

Example 1: interior penalty stabilization with piecewise linear elements: Z T 0 S(uh, ph) dt ≤ Ch3. By Theorem 2.1 we have Z T 0 S(uh, ph) dt = Z T 0 J[uh; (uh, ph), (uh, ph)] dt ≤ Ch3.

Example 2: SUPG, piecewise linear continuous elements in space and piecewise linear discontinuous elements in time [16]. In this case we obtain the rate

Z T

0

S(uh, ph) dt ≤ Ch

3 2

L2_{-norms are over the space time domain Ω × [0, T ])}

Z T

0

S(uh, ph) dt = (δ(∂tuh+ L(uh)Uh− f ), ((∂tuh+ L(uh)Uh)))

≤ kδ12(∂ tuh+ L(uh)Uh)k(kδ 1 2∂ t(uh− u)k + kδ 1 2(L(u h)Uh− L(u)U )k + kδ 1 2ν∆uk) ≤ Ch32

Here L(uh)Uh = (uh· ∇)uh+ ∇ph.

Notice that for the non-homogeneous equation, the convergence rate of the artificial dissipation of the SUPG-method is the square root of that of the interior penalty stabilized formulation. This is due to the fact that the SUPG method is not adjoint consistent. For the homogeneous equation on the other hand both methods have the same convergence rates.

3.5 The dissipation ratio

We now introduce the normalized artificial dissipation, D = RT 0 S(uh, ph) dt RT 0 kν 1 2∇u hk2 dt . (15)

Rewriting (14) using (15) we have 1 2kuh(T )k 2_{+ (1 + D)} Z T 0 kν12∇u hk2 dt = 1 2kuh(0)k 2₊ Z T 0 (f , uh) dt. (16)

For the sake of discussion assume f = 0 and consider the decomposition of u into the L2_{-projection π}

huand its orthogonal complement (I − πh)u we may then rewrite (13) as:

1 2kπhu(T )k 2₊1 2k(I − πh)u(T )k 2₊ Z T 0 kν12∇uk2 dt = 1 2ku(0)k 2_,

but (I − πh)u(T ) represents the unresolved scales and hence

1 2k(I − πh)u(T )k 2 ₌ Z ∞ ξh E(ξ) dξ by the definition of ξh, leading to

1 2 ku(0)k2_{− kπ} hu(T )k2 RT 0 kν 1 2∇uk2 dt = 1 + R∞ ξh E(ξ) dξ RT 0 kν 1 2∇uk2 dt

(ξ)

E

Log

Log

(ξ)

Integral

scale

Inertial

range

turbulent cascade into heat

Viscous dissipation

scale

Dissipation

(ξ)

ξ

ξ

Figure 2: Schematic representation of the inertial range for two-dimensional turbulence

We conclude that if kuhk ≈ kπhuk is to hold, that is if the energy content of the finite

element solution is to match that of the L2_{-projection of the exact solution, then we must}

have D ≈ R∞ ξh E(ξ) dξ RT 0 kν 1 2∇uk2 dt

It is interesting to note that for a standard Galerkin formulation D = 0 showing that no turbulence model at all can never be physical unless all the scales down to the viscous dissipation are resolved. The spurious oscillations observed in standard Galerkin approximations at high Reynolds number are a consequence of this.

3.6 Two-dimensional turbulence: a case study

SLIP CONDITION

SLIP CONDITION

PERIODIC BOUNDARY

PERIODIC BOUNDARY

σ

Figure 3: Schematic picture of the mixing layer problem

radial coordinate in two-space dimensions we have E(ξ) ∼ ξ|ˆu(ξ)|2_{. If ξ}

h ∼ h−1 is in the

inertial range where E(ξ) ∼ ξ−αE then, denoting by ξ

v the wave number at which viscous

dissipation sets in,

D ∼ C Rξv ξh ξ −αE dξ RT 0 kν 1 2∇uk2 dt ∼ Cξ −αE+1 h − ξ−α E+1 v RT 0 kν 1 2∇uk2 dt . Assuming ξ−αE+1 v negligeable we expect D ∼ hαE−1

It follows that if we look for the convergence rate αD of the dissipation ratio, D ∼ hαD,

then there should hold αD = αE − 1. We recall that for decaying homogeneous isotropic

turbulence in two dimensions there holds E(ξ) ∼ Cξ−3 _{and by the above argument we}

would therefore expect that αD = 2.

4 A NUMERICAL ILLUSTRATION: MIXING LAYER AT RE=10000

In this section we consider the numerical simulation of a mixing layer at Reynolds number 10000. The computational domain with flow directions and shear layer position indicated is presented in Figure 3. In this problem a shear layer instability creates four vortices that fuses into one big vortex in a two step transition process.

leading to a Reynolds number based on the layer width of Reσ = 10000. The initial data is taken as u₀ = 1 0  u∞ tanh  2y − 1 σ0  +  ∂yψ −∂xψ  .

The stream function ψ specifies the form of the perturbation. In our case four vortices given by ψ = cu∞exp  −(x2− 0.5) 2 σ2 0  cos(αωx)

where c is a parameter giving the strength of the perturbation. Here we have chosen c = 0.001, and αω = 8π. The very small viscosity makes the solution very sensitive

to under-resolution and over-dissipation in particular if one is to obtain a reasonable approximation of the transition sequence.

We have performed computations on three consecutive structured triangular meshes with equal order finite element spaces for velocity and pressure approximation and using the interior penalty formulation defined by equations (3), (6)–(10). The meshes were constructed using 81, 161 or 321 equidistributed degrees of freedom per side. The spaces were of polynomial order k = 1 (denoted P 1/P 1) or k = 2 (denoted P 2/P 2).

To avoid any influence from the time discretization we used a very small timestep, τ = 0.0015625 and time-stepped to the non-dimensional time T=200 (actual time scaled by σ0

u∞ =

1

28) using a backward differentiation scheme of second order (BDF2).

All the computations in this section has been performed using the finite element soft-ware FreeFem++ [17].

In all computations using affine elements the stabilization parameters were chosen as γS = 0.01, γCD = 0.05. In the case of quadratic elements γS was diminished by a factor

ten. This is approximately in accordance with the factor proposed in [11].

In Table 1 we give the convergence rates of the dissipation ratio for the P 1/P 1 dis-cretization to the left and the P 2/P 2 disdis-cretization to the right. Both computations yield a convergence order of approximately D ∼ h2 _{for the dissipation ratio (compared to the}

theoretical value of the error estimate for laminar flow of h3 _{for affine approximations and}

h5 _{for quadratic approximations), which corresponds well to the value α}

D = 2 predicted

by the argument in the previous section.

In the Figures 4–7 we have plotted the contour lines of the vorticity. Note that the transition points seem to be captured reasonably well on the 320 × 320 mesh with affine elements and already on the 40 × 40 mesh with quadratic elements. Comparing with Table 1 we see that this seems to coincide with a resolution for which D < 1 i.e. when the physical dissipation dominates the artificial dissipation. It is remarkable that this is the case already on the coarsest P 2/P 2 mesh having only 40 elements on the side of the domain. The P 1/P 1 computation on the 320 × 320 mesh (maxK∈ThReK ≈ 875 at t = 0)

has very similar large scale behavior as the P 2/P 2 computation on the 40 × 40 mesh (maxK∈ThReK ≈ 7000 at t = 0) and their dissipation ratios are similar. On the other

Table 1: Convergence of the dissipation ratio D for computations of the mixing layer using the P 1/P 1 interior penalty method (left) and the P 2/P 2 interior penalty method (right). The average number of (unpreconditionned) GMRES iterations per timestep is also presented

P1 el.

per side D O(h

α_{) GMRES} P2 el.

per side D O(h

α_{) GMRES}

80 5.6 - 38 40 0.38 - 93

160 1.4 2 66 80 0.1098 1.79 148

320 0.3 2.22 123 160 0.025 2.0 251

transition sequence on coarse meshes is better captured when using P 2/P 2 approximation. In fact in this case the large scale behavior on the 40 × 40 mesh and the 160 × 160 mesh is very similar (see Figures 6-7).

Qualitatively the affine and quadratic elements differ in the sense that the P 2/P 2 approximation captures the large scale features already on the coarsest mesh and then fill in fine scale features, whereas in the P 1/P 1 case the large and fine scale features seem to be resolved in parallel.

REFERENCES

[1] Y. Bazilevs and T.J.R Hughes. Weak imposition of Dirichlet boundary conditions in fluid mechanics. Technical Report ICES 05-25, The Institute For Computational Engineering And Sciences, University of Texas, 2005.

[2] R. Becker and M. Braack. A finite element pressure gradient stabilization for the Stokes equations based on local projections. Calcolo, 38(4):173–199, 2001.

[3] M. Braack and E. Burman. Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method. SIAM J. Numer. Anal. 43, no. 6, 2544–2566, 2006.

[4] M. Braack, E. Burman, V. John, and G. Lube. Stabilized finite element methods for the generalized Oseen problem. Technical Report NAM-Preprint, Institut f¨ur Numerische und Angewandte Mathematik, Goettingen, 2005.

[5] A. Brooks, T. J. R. Hughes. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. FENOMECH ’81, Part I (Stuttgart, 1981). Comput. Methods Appl. Mech. Engrg. 32 (1982), no. 1-3, 199–259.

[7] E. Burman and M. Fern`andez. A finite element method with edge oriented stabi-lization for the time-dependent Navier-Stokes equations: space discretization and convergence. Technical Report 11.2005, Ecole Polytechnique Federale de Lausanne, 2005.

[8] E. Burman, M.A. Fern´andez, and P. Hansbo. Continuous interior penalty finite element method for Oseen’s equations. SIAM J. Numer. Anal. to appear.

[9] E. Burman and P. Hansbo. Edge stabilization for Galerkin approximations of convection–diffusion–reaction problems. Comput. Methods Appl. Mech. Engrg., 193:1437–2453, 2004.

[10] E. Burman and P. Hansbo. Edge stabilization for the generalized Stokes problem: A continuous interior penalty method. Comput. Methods Appl. Mech. Engrg., 195:2393– 2410, 2006.

[11] E. Burman and A. Ern. Continuous interior penalty hp-finite element methods for transport operators. Math. Comp to appear (see also IACS Technical Report 02.2005, Ecole Polytechnique Federale de Lausanne) 2005.

[12] R. Codina. Stabilized finite element approximation of transient incompressible flows using orthogonal subscales. Comput. Methods Appl. Mech. Engrg., 191(39-40):4295– 4321, 2002.

[13] J. Freund and R. Stenberg. On weakly imposed boundary conditions for second order problems. In Proceedings of the Ninth Int. Conf. Finite Elements in Fluids, Venice 1995. M. Morandi Cecchi et al., Eds. pp. 327-336, 1995.

[14] V. Gravemeier, W. A. Wall, and E. Ramm. Large eddy simulation of turbulent incompressible flows by a three-level finite element method. Internat. J. Numer. Methods Fluids, 48(10):1067–1099, 2005.

[15] J.-L. Guermond. Stabilization of Galerkin approximations of transport equations by subgrid modeling. M2AN Math. Model. Numer. Anal. 33, no. 6, 1293–1316, 1999. [16] P. Hansbo and A. Szepessy. A velocity-pressure streamline diffusion finite element

method for the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Engrg., 84(2):175–192, 1990.

[17] F. Hecht and O. Pironneau. FreeFEM++ Manual. Laboratoire Jacques Louis Lions, 2005.

[19] V. John and S. Kaya. A finite element variational multiscale method for the Navier-Stokes equations. SIAM J. Sci. Comput. 26, no. 5, 1485–1503, 2005.

[20] V. John. An assessment of two models for the subgrid scale tensor in the rational LES model. J. Comput. Appl. Math., 173(1):57–80, 2005.

[21] C. Johnson, U. N¨avert. An analysis of some finite element methods for advection-diffusion problems. in Analytical and numerical approaches to asymptotic problems in analysis (Proc. Conf., Univ. Nijmegen, Nijmegen, 1980), pp. 99–116, North-Holland Math. Stud., 47, North-Holland, Amsterdam-New York, 1981.

[22] A. N. Kolmogorov Local structure of turbulence in an incompressible viscous fluid at very high Reynolds numbers. Dokl. Akad. Nauk SSSR 30, 299, 1941.

[23] R. H. Kraichnan. Inertial ranges in two-dimensional turbulence. Phys. Fluids, 10:1417–1423, 1967.

[24] W. Layton. Weak imposition of “no-slip” conditions in finite element methods. Com-put. Math. Appl., 38(5-6):129–142, 1999.

[25] W. Layton. A connection between subgrid scale eddy viscosity and mixed methods. Appl. Math. Comput., 133(1):147–157, 2002.

[26] M. Lesieur, C. Staquet, P. Le Roy, and P. Comte. The mixing layer and its coherence examined from the point of view of two-dimensional turbulence. J. Fluid Mech., 192:511–534, 1988.