Finite Element Error Analysis of a Variational Multiscale Method

Finite Element Error Analysis of a Variational Multiscale Method for the Navier-Stokes Equations Volker John∗ and Songul Kaya† Abstract The paper presents finite element error estimates for a variational multiscale method (VMS) for the incompressible Navier-Stokes equations. The constants in these estimates do not depend on the Reynolds number but on a reduced Reynolds number or on the mesh size of a coarse mesh. Keywords. Variational multiscale method, finite element method, error analysis, incompressible Navier-Stokes equations AMS subject classification. 65M15, 65M60, 76F65 1 Introduction The flow of an incompressible fluid is governed by the incompressible NavierStokes equations ut − 2ν · D(u) + (u · )u + p ·u u u(0, x) Ω = = = = f 0 0 u0 in in in in (0, T ] × Ω, [0, T ] × Ω, [0, T ] × ∂Ω, Ω, (1) p dx = 0, in (0, T ]. Here, Ω ⊂ Rd , d ∈ {2, 3}, is a bounded domain with boundary ∂Ω, [0, T ] a finite time interval, u(t, x) is the fluid velocity and p(t, x) the fluid pressure. Faculty of Mathematics, Otto-von-Guericke-University, Postfach 4120, 39016 Magdeburg, Germany; email:john@mathematik.uni-magdeburg.de † Department of Mathematics, University of Pittsburgh, Pittsburgh, PA,15260, U.S.A.; email: sokst20@pitt.edu; This work is partially supported by NSF grants DMS9972622, DMS20207627 and INT9814115. ∗ 1 The viscosity ν > 0, which is inverse proportional to the Reynolds number Re = O(ν −1 ), the body forces f (t, x) and the initial velocity field u0 (x) are given. The velocity deformation tensor is the symmetric part of the velocity gradient D(u) = ( u + uT )/2. Turbulent flows are characterized by a small viscosity (or a large Reynolds number). The characteristic feature of such flows is the richness of scales. There are very large flow structures but also rather small ones, e.g., imagine a tornado. A direct discretization of (1), e.g., by a Galerkin finite element method, would try to simulate the behavior of all flow structures. This requires, at least, that all scales can be resolved by the given mesh. However, this is in general, in particular in three dimensions, far beyond the capacity of present day computers. A practical approach, which is the main goal of Large Eddy Simulation (LES), is to simulate only the behavior of the large scales appropriately. Following Collis [5], the flow can be decomposed into three scales - the (resolved) large scales, - the (resolved) small scales, - the unresolved scales. An important aspect of turbulent flows is a mutual influence of all scales. Thus, simply neglecting the unresolved scales would in general result in a drastic change of the behavior of the resolved ones. The numerical simulation would predict a wrong behavior of them. Instead, one needs to model the influence of the unresolved scales onto the resolved ones by using a turbulence model. It still has to be clarified how the resolved scales are defined. In the classical LES approach [17, 12], they are defined as an average in space given by convolution with an appropriate filter function. A crucial assumption in the derivation of equations for the resolved scales is the commutation of convolution and differentiation. However, this property only holds in special cases. In particular, it does not hold in the case that Ω is a bounded domain. A so-called commutation error is committed if convolution and differentiation are commuted nevertheless. Recent analytical results of different commutation errors, [6, 3, 4], show that this error is in general large near the boundary of the domain and it cannot be neglected. This error is one reason why numerical simulations with LES models need in general a special treatment at the boundary (wall models, van Driest damping). In [11], Hughes and co-workers proposed to define the large scales in a different way, namely by a projection into appropriate subspaces. This approach is called Variational Multiscale Method (VMS). VMS might be a 2 remedy of the problems which LES encounters near the boundary. There have been several VMS proposed in the literature [8, 11, 13]. This paper presents a finite element error analysis of the VMS presented in [13]. This VMS has as parameters a large scale space LH consisting of symmetric tensors and a turbulent viscosity νT . The analysis will be carried out for the case that νT is a positive constant. Error estimates are proven for the difference of the velocity u of the Navier-Stokes equations (1) and the finite element velocity uh of the VMS. From the numerical point of view, the VMS introduces additional viscosity which acts directly, however, only on the (resolved) small scales of the discrete velocity. Because of this artificial viscosity, one can expect to obtain error estimates whose constants do not depend on the Reynolds number like for the Galerkin finite element discretization of the Navier-Stokes equations but on some reduced Reynolds number. This will be exactly the results proved in this paper. We will present two results of this kind. The paper is organized as follows. In Section 2, the VMS is introduced. Section 3 provides some auxiliary results which are needed in the proofs and gives an outline of the proofs. The error estimates with constants depending on a reduced Reynolds number are presented in Sections 4 and 5. 2 The variational multiscale method Standard notations for Lebesgue and Sobolev spaces are used throughout this paper, e.g., see Adams [1]. The inner product in (L2 (Ω))d , d ∈ N, is denoted by (·, ·). Generic constants which do not depend on the Reynolds number Re and the mesh width h are denoted by C. 1 v L2 and Q = Let V = (H0 (Ω))d equipped with the norm v V = 2 (Ω). A variational formulation of the Navier-Stokes equations (1) is as L0 follows: Find u : [0, T ] → V, p : (0, T ] → Q satisfying (ut , v) + (2νD(u), D(v)) + bs (u, u, v) − (p, (q, and u(0, x) = u0 (x) ∈ V . Here, bs (u, v, w) = 1 (((u · 2 )v, w) − ((u · )w, v)) · v) = (f , v) ∀ v ∈ V · u) = 0 ∀q∈Q (2) is the skew-symmetric form of the convective term. A Galerkin finite element discretization of (2) is unstable in the case of small viscosity (or high Reynolds number). The use of a turbulence model 3 becomes necessary which should model the action of the unresolved scales onto the resolved scales and which serves as stabilization in a numerical simulation. We consider an approach which was presented in [13], see also [15] for this approach in the context of scalar convection-diffusion equations. We will study the continuous-in-time finite element approach, i.e., only a semi-discretization in space is considered. Let V h ⊂ V and Qh ⊂ Q be conforming finite element spaces which fulfill the inf-sup condition inf sup · vh , q h vh L2 q h ≥ C > 0, L2 q h ∈Qh vh ∈V h (3) where C is independent of h. For the VMS, a large scale space LH ⊂ L = {L ∈ (L2 (Ω))d×d , L = LT } and a so-called turbulent viscosity νT = νT (t, x, uh , ph ) ≥ 0 are introduced. The semi-discrete problem reads as follows: Find uh : [0, T ] → V h , ph : (0, T ] → Qh and GH : [0, T ] → LH satisfying ((2ν + νT )D(uh ), D(vh )) +bs · vh ) − (νT GH , D(vh )) = (f , vh ) ∀ vh ∈ V h (q h , · uh ) = 0 ∀ q h ∈ Qh (GH − D(uh ), LH ) = 0 ∀ LH ∈ LH (4) h (0, x) = uh ∈ V h is a discretely divergence free approximation of u . and u 0 0 The parameters in (4) are the large scale space LH and the turbulent viscosity νT . The turbulent viscosity is added to all resolved scales in the viscous term (the second term in the first equation) and it is subtracted from the resolved large scales in the last term on the left hand side in this equation. The large scales are defined in the last equation of (4) by projection of the resolved scales onto LH . Altogether, νT acts directly only on the resolved small scales. This is the basic idea of a VMS, see [13] for a discussion of the connection of (4) to the VMS proposed in [11]. The implementation of the VMS (4) into a code for solving the Navier-Stokes equations and numerical results are presented in [13]. Consider now the limit cases for LH : - LH = D(V h ), that means, D(vh ) ∈ LH for all vh ∈ V h and if LH ∈ LH then exists a vh ∈ V h such that LH = D(vh ). In this case, D(uh ) ∈ LH , GH = D(uh ) and the turbulence model is subtracted for all scales. System (4) becomes a Galerkin finite element discretization of the Navier-Stokes equations. (uh , vh ) + t h , uh , vh ) − (ph , (u 4 - LH = {O}. Then GH = O and the last term on the left hand side of the first equation in (4) vanishes. In this case, the model νT acts on all resolved scales. We require LH ⊆ D(V h ). If this relation is violated, the turbulence model νT would be subtracted from scales where it was not added previously. Let PLH : L → LH , D(v) → PLH D(v) with (PLH D(v) − D(v), LH ) = 0 ∀ LH ∈ LH (5) denote the L2 -projection from L onto LH . Then GH = PLH D(uh ) in (4). Since PLH is an L2 -projection, it follows for v ∈ V and D(v) L2 > 0 νT (I − PLH )D(v) 2 L2 = = νT νT D(v) 1− 2 L2 − PLH D(v) 2 L2 2 L2 PLH D(v) 2 2 L D(v) 2 2 L 2 L2 . D(v) =: νadd (v) D(v) In addition, from 0 ≤ PLH D(v) L2 (6) ≤ D(v) L2 follows (7) 0 ≤ νadd (v) ≤ νT . Note, if v depends on t then νadd (v), too. From (7) follows νadd (v(t, ·)) ∈ L∞ (0, T ) if νT is bounded almost everywhere in the time interval. If D(v) L2 = 0 then v = 0 since v ∈ V . In this case, we set νadd (v) = 0. In this paper, we consider the case that νT is a non-negative constant. A straightforward calculation shows that (νT D(uh ), D(vh )) − (νT PLH D(uh ), D(vh )) = (νT (I − PLH )D(u ), (I − PLH )D(v )). Thus, System (4) can be reformulated as follows: Find uh V h , ph : (0, T ] → Qh satisfying : [0, T ] → h h (8) (uh , vh ) + (2νD(uh ), D(vh )) t +bs (uh , uh , vh ) − (ph , · vh ) +(νT (I − PLH )D(uh ), (I − PLH )D(vh )) = (f , vh ) ∀ vh ∈ V h (q h , · uh ) = 0 ∀ q h ∈ Qh . (9) h Let Vdiv = {vh ∈ V h : ( ·vh , q h ) = 0 ∀ q h ∈ Qh } the space of discretely divergence free functions. From the inf-sup condition (3) follows that this 5 h space is not empty. Then (9) is equivalent to: Find uh : [0, T ] → Vdiv such that (uh , vh ) + (2νD(uh ), D(vh )) + bs (uh , uh , vh ) t +(νT (I − PLH )D(uh ), (I − PLH )D(vh )) = (f h , vh ) (10) h for all vh ∈ Vdiv . The finite element spaces (V h , Qh ) contain all resolved scales. Let V H ∈ 1 (Ω))d be a finite element space such that LH = D(V H ). The space V H (H should be coarser than V h . The functions of V H may be defined on the same grid as the functions of V h with a lower piecewise polynomial degree or on a coarser grid. But no boundary conditions, like no-slip conditions, are incorporated in the definition of V H , i.e., in general V H ⊂ V h . The large eddies of a turbulent flow generally do not fulfill no-slip boundary conditions, e.g., the large eddies of a tornado move along the surface of the earth instead of sticking on the surface. The pair of spaces for the resolved large scales is given by (V H , QH ). Here, QH is chosen such that an inf-sup condition of type (3) is fulfilled for (V H , QH ). The large scales PH u of the velocity are defined by an elliptic projection into V H and the large scales PH p of the pressure by the L2 -projection into QH ; PH : (V, Q) → (V H , QH ) (D(u − PH u), D(vH )) = 0 (u − PH u, 1) = 0, (p − PH p, q H ) = 0 ∀ vH ∈ V H , ∀ q H ∈ QH . It was proven in [13] that for LH = D(V H ) holds PLH D(v) = D(PH v) ∀ v ∈ V. (11) That means, the definition of the large scales by projection and differentiation commute. This property does not hold in general for the classical LES. 3 Preliminaries and the outline of the proofs The following sections present a finite element error analysis for the space discrete velocity solution of (10). For simplicity, we assume that the characteristic length scale and velocity scale in the Navier-Stokes equations are chosen such that Re = ν −1 . Considering the limit cases of the choice of LH in Section 2, one has: 6 - LH = D(V h ): a finite element error estimate which constants depending on 2Re or (2ν)−1 , e.g., see Heywood and Rannacher [9, 10]; - LH = {O}: a finite element error estimate where the most constants, in particular the constant in the dominating exponential factor, depend on ReνT := (2ν +νT )−1 , e.g., see the analysis for the Smagorinsky turbulence model in [14, 12]. Since in our analysis, the finite element solution of the VMS is compared to the solution of the continuous Navier-Stokes equations (2), some constants may depend on 2Re instead on ReνT . The limit cases lead to the expectation that if LH is chosen in between them allows finite element error estimates with constants depending on a reduced Reynolds number Rered with ReνT < Rered < 2Re. (12) In following sections, such error estimates will be derived for the case that νT is constant. For the finite element error analysis, we need some assumptions on the regularity of the parameters of the Navier-Stokes equations and the solution. We assume that f ∈ (L2 (0, T ; L2 ))d , u0 ∈ V, (13) and that (2) possesses a solution (u, p) with u ∈ (L4 (0, T ; L2 ))d×d , ut ∈ (L2 (0, T ; H −1 ))d , p ∈ L4 (0, T ; L2 ). (14) Note, these assumptions imply that Serrin’s condition is fulfilled from what follows that the solution of (2) is unique, e.g., see Temam [19], Galdi [7] or Sohr [18]. For simplicity let f = f h . In addition, we assume that Ω has a polygonal (in 2d) or polyhedral (in 3d) boundary such that no boundary approximation in the application of the finite element method becomes necessary. Inequalities which will be used frequently are Young’s inequality t t−q/p q ab ≤ ap + b , p q a, b, p, q, t ∈ R, 1 1 + = 1, p q p, q ∈ (1, ∞), t > 0, (15) Poincar´’s inequality in V e v and Korn’s inequality in V v L2 L2 ≤C v L2 ∀v∈V (16) ≤ C D(v) 7 L2 ∀ v ∈ V. (17) The proof of the finite element error estimate uses an approach by Rannacher and Heywood [9, 10]. We will first give an outline: 1. Prove stability of u and uh , i.e., certain norms of u and uh are bounded a priori by the data of the problem: f , u0 , ν. 2. Derive an error equation by subtracting (10) from (2) for test functions h from Vdiv . Split the error into an approximation term η and a (finite element) remainder φh ˜ ˜ e = u − uh = (u − uh ) − (uh − uh ) =: η − φh (18) h ˜ where uh ∈ Vdiv is a projection of u which fulfills certain interpolation estimates. An example for such a projection is the Stokes projection. Then, take φh as test function in the error equation. 3. Estimate the right hand side of the error equation such that one obtains an inequality of the form d φh dt 2 L2 + g1 (t, φh ) ≤ g2 (t, η, u) + g3 (t, u) φh 2 L2 (19) where all functions are non-negative for almost all t ∈ [0, T ]. 4. Show that Gronwall’s lemma can be applied to (19), i.e., show that all functions in (19) belong to L1 (0, T ). Apply Gronwall’s lemma to get an estimate for φh . 5. Prove the error estimate for e by applying the triangle inequality to (18). Along these lines, two estimates with constants depending on a reduced Reynolds number will be proved. 4 First error estimate with constants depending on a reduced Reynolds number This error estimate uses the parameter νadd defined in (6). We will first prove the stability of u and uh . Lemma 4.1 The solution uh of the finite element problem (4) fulfills uh ∈ (L∞ (0, T ; L2 ))d and D(uh ) ∈ (L2 (0, T ; L2 ))d×d . The velocity solution of the continuous problem (2) fulfills u ∈ (L∞ (0, T ; L2 ))d and D(u) ∈ (L2 (0, T ; L2 ))d×d . 8 Proof: The proof for uh and u is very similar. We will show the result for uh . Set vh = uh in (10), use the skew symmetry of bs (·, ·, ·), (6), the standard estimate of the dual pairing, Korn’s inequality (17) and integrate over (0, t) with t ≤ T : t 1 h u (t) 2 2 + (2ν + νadd (uh (τ ))) D(uh )(τ ) 2 2 dτ L L 2 0 t 1 h 2 ≤ u f (τ ) H −1 uh (τ ) L2 dτ 2 + 2 0 L 0 t 1 h 2 C 2ν + νadd (uh (τ )) ≤ u0 L2 + f 2 2 (0,t;H −1 ) + D(uh )(τ ) L 2 ν 2 0 2 L2 dτ. Subtraction of the last term and the regularity (13) of the coefficients give D(uh ) ∈ (L2 (0, T ; L2 ))d×d . Taking then the supremum of t ∈ (0, T ) gives the statement uh ∈ (L∞ (0, T ; L2 ))d . The splitting of the error (18) is performed with the help of a projection h ˜ uh ∈ Vdiv of u. Let t ∈ [0, T ] be arbitrary. We require that this projection fulfills η ηt L2 L2 + h D(η) L2 L2 ≤ Chk ( u(t, ·) ≤ Ch ( u(t, ·) k Hk Hk + ν −1 p(t, ·) +ν −1 H k−1 ), H k−1 ), (20) (21) + h (D(η))t p(t, ·) where the constants depend only on Ω. Korn’s inequality (17), (20) with k = 1 and the regularity assumptions (14) imply η ∈ (L4 (0, T ; L2 ))d×d . (22) An example for an appropriate projection is the Stokes projection which is h ˜ the solution of: Find uh ∈ Vdiv such that ˜ (2νD(u(t, ·) − uh ), D(vh )) = (p(t, ·), h · vh ) ∀ vh ∈ Vdiv . Let u(t, ·) ∈ (H k (Ω))d , p(t, ·) ∈ H k−1 (Ω), k ≥ 1 and V h possess a (k − 1)th order approximation property, e.g., V h is the finite element space P k−1 on simplicial meshes or Qk−1 on quadrilateral/hexahedral meshes, then a simple scaling argument of Lemma 5.3. in [10] gives (20), (21). For t = 0, the pressure can be well defined, e.g., see [9, 19]. Now, Step 2 of the proof is carried out by a straightforward calculation. 9 One obtains 1d φh 2 2 + (2ν + νadd (φh )) D(φh ) 2 2 (23) L L 2 dt = (η t , φh ) + (2νD(η), D(φh )) + (νT (I − PLH )D(η), (I − PLH )D(φh )) +bs (u, u, φh ) − bs (uh , uh , φh ) − (νT (I − PLH )D(u), (I − PLH )D(φh )) −(p − λh , · φh ) with arbitrary λh ∈ Qh . To get an inequality of form (19), the terms on the right hand side of (23) have to be estimated. All bilinear terms are estimated essentially in the same way: using the Cauchy-Schwarz inequality (or the estimate for the dual pairing), Korn’s inequality (17) and Young’s inequality (15). In addition, (6) is used. One obtains (η t , φh ) ≤ ≤ ηt H −1 φh h L2 ≤ C ηt 2 L2 H −1 D(φh ) L2 2 H −1 , 2ν + νadd (φ ) D(φh ) 8 + C ηt 2ν + νadd (φh ) (2νD(η), D(φh )) ≤ 2ν D(η) L2 D(φh ) L2 ν ≤ D(φh ) 2 2 + 8ν D(η) 2 2 , L L 8 h h h h (p − λ , · φ ) ≤ p − λ L2 · φ L2 ≤ C p − λh ≤ 2ν + νadd (φ ) D(φh ) 8 h 2 L2 L2 D(φh ) L2 2 L2 , + C p − λh h 2ν + νadd (φ ) (νT (I − PLH )D(η), (I − PLH )D(φh )) νT ≤ (I − PLH )D(φh ) 2 2 + 4νT (I − PLH )D(η) L 16 νadd (φh ) = D(φh ) 2 2 + 4νadd (η) D(η) 2 2 , L L 16 h (νT (I − PLH )D(u), (I − PLH )D(φ )) ≤ νT (I − PLH )D(u) = νadd (φ ) D(φh ) 16 h 2 L2 L2 2 L2 (I − PLH )D(φh ) L2 2 L2 . + 4νT (I − PLH )D(u) The trilinear term is first decomposed into three terms. A direct calculation gives bs (u, u, φh ) − bs (uh , uh , φh ) = bs (η, u, φh ) − bs (φh , u, φh ) + bs (uh , η, φh ). 10 The terms on the right hand side are estimated separately using the inequality 1/2 1/2 bs (u, v, w) ≤ C u L2 D(u) L2 D(v) L2 D(w) L2 . (24) This estimate is well known. It can be derived by applying H¨lder’s ino equality, Sobolev imbeddings, interpolation theorems in Sobolev spaces, Poincar´’s and Korn’s inequality, e.g., see Layton and Tobiska [16]. One e obtains by applying (24) and Young’s inequality (15) bs (η, u, φh ) ≤ C η ≤ 1/2 L2 1/2 L2 D(η) h D(u) 2 L2 L2 D(φh ) L2 L2 2ν + νadd (φ ) D(φh ) 8 + C η 2ν + νadd (φh ) D(η) L2 D(u) 2 L2 , bs (φh , u, φh ) ≤ C φh ≤ 1/2 L2 D(u) L2 D(φh ) 2 L2 3/2 L2 2ν + νadd (φh ) D(φh ) 8 + C (2ν + νadd (φ h ))3 φh 2 L2 D(u) 4 L2 , bs (uh , η, φh ) ≤ C uh ≤ 1/2 L2 D(uh ) 1/2 L2 D(η) 2 L2 L2 D(φh ) L2 L2 2ν + νadd (φh ) D(φh ) 8 1d φh 2 dt + C uh 2ν + νadd (φh ) D(uh ) L2 D(η) 2 L2 . Collecting terms gives 2ν + νadd (φh ) D(φh ) 2 2 L 4 C ≤ η t 2 −1 + (8ν + 4νadd (η)) D(η) 2 2 H L 2ν + νadd (φh ) C + p − λh 2 2 + 4νT (I − PLH )D(u) 2 2 L L h 2ν + νadd (φ ) C + η L2 D(η) L2 D(u) 2 2 L 2ν + νadd (φh ) 2 L2 + + uh + L2 D(uh ) C L2 D(η) 2 L2 (2ν + νadd (φh ))3 D(u) 4 L2 φh 2 L2 . 11 We define Rered := (2ν + inf νadd (φh (t)))−1 . t∈(0,T ] (25) It follows that Rered is smaller or equal than 2Re. Using νadd (η) ≤ νT finishes Step 3 of the proof: d φh dt 2 L2 + Re−1 red D(φh ) 2 2 H −1 2 L2 2 L2 ≤ C Rered η t + (Re−1 + νT ) D(η) 2 L2 L2 + Rered p − λh 2 L2 +νT (I − PLH )D(u) +Rered η L2 D(η) 4 L2 D(u) 2 L2 + uh L2 D(uh ) L2 D(η) 2 L2 +CRered D(u) φh 2 L2 . (26) To perform Step 4 of the proof, the L1 (0, T )-regularity of the terms appearing in (26) has to be studied. Let t ∈ (0, T ] be arbitrary. We have by Poincar´’s inequality (16), Korn’s inequality (17), the Cauchy-Schwarz e inequality, (14) and (22) t η(τ ) 0 L2 t D(η)(τ ) D(η)(τ ) L2 2 L2 D(u)(τ ) D(u)(τ ) 2 L2 dτ 2 L2 dτ ≤ C 0 ≤ C D(η) 2 L4 (0,t;L2 ) D(u) 2 L4 (0,t;L2 ) < ∞. Similarly follows by H¨lders inequality, Lemma 4.1 and (22) o t uh (τ ) 0 L2 D(uh )(τ ) t L2 D(η)(τ ) L2 2 L2 dτ ≤ ≤ uh uh L∞ (0,t;L2 ) 0 L∞ (0,t;L2 ) D(uh )(τ ) D(uh ) 2 L2 D(η)(τ ) D(η) 2 L2 dτ L2 (0,t;L2 ) 2 L4 (0,t;L2 ) ≤ C Re1/2 uh 0 + Re3/2 f 2 L2 (0,t;H −1 ) D(η) 2 L4 (0,t;L2 ) < ∞. The L1 (0, T )-regularity of the other terms is a direct consequence of (14), (20), (21) and (22). The application of Gronwall’s inequality and the last step of the proof are straightforward. 12 h Theorem 4.2 Let (u, p) ∈ V × Q be the solution of (2) and let uh ∈ Vdiv be the solution of (10) where νT ≥ 0 is a constant. Let the regularity assumph ˜ ˜ tions (14) be fulfilled, uh be a projection of u into Vdiv such that η = u − uh fulfills (20) and (21). Let the reduced Reynolds number Rered be defined in (25). Then, the error u − uh satisfies for T ≥ 0 (u − uh )(T ) ≤ C 2 L2 + Re−1 red D(u − uh ) 2 η(T ) 2 L2 2 L2 (0,T ;L2 ) 2 L2 (0,T ;L2 ) inf h ˜ uh ∈L4 (0,T ;Vdiv ) λh ∈L2 (0,T ;Qh ) + Re−1 D(η) red (27) + exp CRe3 D(u) red 4 L4 (0,T ;L2 ) u0 − uh 0 2 L2 ˜ + u0 − uh (0) 2 L2 +Rered + D(η) ηt 2 L2 (0,T ;H −1 ) + p − λh 2 L4 (0,t;L2 ) 2 L2 (0,T ;L2 ) 2 L4 (0,t;L2 ) 2 L2 D(u) + Re1/2 uh 0 + Re3/2 f 2 L2 (0,t;H −1 ) D(η) 2 L4 (0,t;L2 ) 2 L2 (0,T ;L2 ) +(Re−1 + νT ) D(η) 2 L2 (0,T ;L2 ) + νT (I − PLH )D(u) . Remark 4.1 We consider the convergence of the right hand side of (27). The crucial term is the last one on this side since it does not possess a factor where the interpolation error η appears. This term tends to zero as the mesh width h → 0 if νT → 0 or if LH tends to D(V ). In both cases, the Galerkin finite element discretization of the Navier-Stokes equations is recovered asymptotically. Otherwise, in particular if νT and LH are fixed and h → 0, one cannot expect that the solution of the discrete system converges to the solution of the continuous problem. Let (u, p) ∈ (H k+1 (Ω))d × H k (Ω) for all times, k ≥ 1, let the velocity finite element space be of piecewise order k and let the pressure finite element space be of piecewise order k − 1. Then, the optimal order of convergence of the left hand side of (27) is hk where h is the mesh parameter connected with the space V h . Again, the crucial term on the right hand side is the last one. Its optimal order of convergence is H k . Here, H is the mesh parameter connected with LH . Depending on the ratio of h and H, the artificial viscosity νT can be chosen in such a way that the last term on the 13 right hand side of (27) has at least the same order of convergence like the left hand side of (27). For fixed h and νT → 0, the estimate in Theorem 4.2 tends to the estimate for the Galerkin discretization of the Navier-Stokes equations. Remark 4.2 There is no improvement in the constant in the exponential, i.e. Rered = 2Re, if there is a time t at which νadd (φh (t)) = 0. This is equivalent to PLH D(φh (t)) 2 2 = D(φh (t)) 2 2 or L L (I − PLH )D(uh ) = (I − PLH )D(˜ h ). u (28) ˜ That means, the resolved small scales of uh and uh are the same. However, this situation is unlikely for turbulent flows since these scales of uh are considerably influenced by the model which is used for the unresolvable ˜ small scales whereas the interpolation uh does not posses any information h is defined by the Stokes projection which is ˜ about this model, e.g., if u asymptotically optimal. In this case, (28) is only likely if there are only large scales in the flow, which is not the case in turbulent flows. From the mathematical point of view, the difficulty consists in the fact that the equations for laminar flows and turbulent flows are the same, namely the Navier-Stokes equations (1). Since the analysis is carried out for (1), it is not possible to distinguish between the two kinds of flows and the results must also hold for the case of laminar flows. For such flows, νadd (φh (t)) may vanish and the error estimates of [9, 10] are recovered in which the constants depend on Re. Remark 4.3 A finite element error estimate for the L2 (Ω)-error in the pressure can also be derived following Heywood and Rannacher [10], Section 7. Using the inf-sup condition (3) and the estimates for the Stokes projection (20) and (21), the pressure error can be estimated by approximation errors and the velocity error (u − uh )(t) L2 . Then, the result of Theorem 4.2 finishes the error estimate for the pressure. Since the analysis is lengthy and follows closely [10], we will not present it here. 5 Second error estimate with constants depending on a reduced Reynolds number This section will present an error estimate with a mesh-dependent reduced Reynolds number. This reduced Reynolds number will be considerably smaller than Re if νT >> ν and if the mesh width H connected with the 14 space LH is also much larger than ν. This is the typical situation in turbulent flow simulations. The starting point for this error estimate is the error equation 1d φh 2 2 + 2ν D(φh ) 2 2 + νT (I − PLH )D(φh ) 2 2 (29) L L L 2 dt = (η t , φh ) + (2νD(η), D(φh )) + (νT (I − PLH )D(η), (I − PLH )D(φh )) +bs (u, u, φh ) − bs (uh , uh , φh ) − (νT (I − PLH )D(u), (I − PLH )D(φh )) −(p − λh , · φh ) with arbitrary λh ∈ Qh which follows directly by subtracting (10) from (2) h and taking φh ∈ Vdiv as test function. We assume that the condition LH = D(V H ) holds and that the finite element spaces V h , V H rely on quasiuniform triangulations of Ω. The first assumption was needed to prove the commutation of differentiation and the definition of the large scales by projection. From the latter assumption follows that inverse estimates for finite element functions hold. Starting with (11) and applying the inverse estimate for V H gives PLH D(φh ) L2 = D(PH φh ) L2 ≤ CH −1 PH φh L2 , (30) where H is the mesh parameter connected with V H . We assume now that h the elliptic projection is L2 -stable for functions from Vdiv , i.e., there is a constant C such that PH φ h L2 ≤ C φh L2 h ∀ φh ∈ Vdiv . (31) Together with (30), this gives PLH D(φh ) L2 ≤ CH −1 φh L2 h ∀ φh ∈ Vdiv . (32) Assumption (31) is true, e.g., for quasiuniform meshes. From Babuˇka and s h h Osborn [2, eq. (6.6)] follows PH φ L2 ,h ≤ C φ L2 ,h for a mesh-dependent norm · L2 ,h . For finite element functions, this mesh-dependent norm is equivalent to · L2 , from which the L2 -stability of the elliptic projection is obtained immediately. The first step in the estimate of the terms on the right hand side of (29) is the same as in the derivation of the error estimate leading to Theorem 4.2. But then, the term D(φh ) L2 is estimated further using the triangle inequality and (32) D(φh ) L2 ≤ ≤ (I − PLH )D(φh ) (I − PLH )D(φ ) 15 h L2 L2 + PLH D(φh ) + CH −1 L2 φ h L2 . This gives, e.g., (η t , φh ) ≤ ηt H −1 (I − PLH )D(φh ) −1 −2 L2 + CH −1 η t + 2 L2 . ≤ C max{ν, νT } + H max{ν, νT } + (I − PLH )D(φh ) 64 η t 2 −1 H H −1 h 2 φ L2 φh L2 (33) The last term on the right hand side of this estimate has to be hidden on the left hand side of (29). If νT ≥ ν, then it is absorbed from the third term on the left hand side. In the case νT < ν, we obtain by using the triangle inequality and the H 1 (Ω)-stability of the elliptic projection with constant 1 max{ν, νT } (I − PLH )D(φh ) 64 2 L2 ≤ ≤ = 3ν D(φh ) 2 2 + PLH D(φh ) L 64 3ν D(φh ) 2 2 + D(φh ) 2 2 L L 64 3ν D(φh ) 2 2 . L 32 2 L2 In this case, this term is absorbed by the second term on the left hand side of (29). The trilinear term which determines the most important coefficient in the exponential factor of the final estimate can be estimated in the following way bs (φh , u, φh ) ≤ C φh 1/2 L2 D(u) L2 (I − PLH )D(φh ) L2 3/2 L2 + PLH D(φh ) 3/2 L2 −3 ≤ C νT D(u) 4 2 + H −3/2 D(u) L νT + (I − PLH )D(φh ) 2 2 L 64 φh 2 L2 or in the standard way bs (φh , u, φh ) ≤ ν D(φh ) 64 2 L2 + Cν −3 D(u) 4 L2 φh 2 L2 . For the final estimated of the VMS, we will take the estimate which gives the smaller constant. All other terms in (29) are estimated in the same fashion as (33). The conditions which allow the application of Gronwall’s lemma are checked in the same way as for the error estimate given in Theorem 4.2. The second error estimate is formulated in the following theorem. 16 h Theorem 5.1 Let (u, p) ∈ V × Q be the solution of (2) and let uh ∈ Vdiv be the solution of (10) where νT ≥ 0 is a constant. Let the regularity assumph ˜ ˜ tions (14) be fulfilled, uh be a projection of u into Vdiv such that η = u − uh H = D(V H ), let V h and V H be defined on quasifulfills (20) and (21). Let L uniform triangulations of Ω and let (31) be fulfilled. Then, the error u − uh satisfies for T ≥ 0 Re−1 (u − uh )(T ) 2 2 + D(u − uh ) L 16 νT + (I − PLH )D(u − uh ) 2 2 L 16 ≤ C inf h ˜ uh ∈L4 (0,T ;Vdiv ) λh ∈L2 (0,T ;Qh ) 2 L2 (0,T ;L2 ) η(T ) 2 L2 + Re−1 D(η) red 2 L2 (0,T ;L2 ) + exp 4T + C min ν −3 D(u) −3 νT D(u) 4 L4 (0,T ;L2 ) 4 L4 (0,T ;L2 ) , + H −3/2 D(u) 2 L2 L1 (0,T ;L2 ) × u0 − uh 0 2 L2 ˜ + u0 − uh (0) + D(η) + Rered D(u) ηt 2 L2 (0,T ;H −1 ) + p − λh 2 L2 (0,T ;L2 ) 2 L2 2 L4 (0,t;L2 ) 2 L4 (0,t;L2 ) 2 L4 (0,t;L2 ) 2 L2 (0,T ;L2 ) (34) + Re1/2 uh 0 + Re3/2 f 2 L2 (0,t;H −1 ) D(η) +(Re−1 + νT ) D(η) where 2 L2 (0,T ;L2 ) + νT (I − PLH )D(u) Rered = max{ν, νT }−1 + H −2 . (35) Remark 5.1 The viscosity ν is very small for turbulent flows. To have a stabilizing effect by using an artificial viscosity, one has to choose νT > ν, −1 in general even νT >> ν. Hence max{ν, νT }−1 = νT in (35). Then, the reduced Reynolds number (35) is not dominated by the mesh size if νT ≤ H 2 . (36) Making the same assumptions and considerations on the convergence of the individual terms in (34) as in Remark 4.1, one finds that the second term on 17 the left hand side of (34) behaves like hk , the third term on the left hand side like νT hk and the last term on the right hand side like νT H k . Given h and either H or νT , one can choose by equilibrating these orders of convergence an appropriate value for the remaining parameter, taking into account also restriction (36). Let ν be small, ν << h ≤ H, let h be fixed and νT → 0. Then, Rered is very close to ν −1 = Re and the constants in estimate (34) have the same dependency on Re as for the Galerkin finite element discretization of the Navier-Stokes equations. Acknowledgment. We would like to acknowledge W. Layton for fruitful discussion on the subject of this paper. References [1] R.A. Adams. Sobolev spaces. Academic Press, New York, 1975. [2] I. Babuˇka and J. Osborn. Analysis of finite element methods for second order s boundary value problems using mesh dependent norms. Numer. Math., 34:41 – 62, 1980. [3] L.C. Berselli, C.R. Grisanti, and V. John. On commutation errors in the derivation of the space averaged Navier-Stokes equations. Preprint 12, Universit` di a Pisa, Dipartimento di Matematica Applicata “U.Dini”, 2004. [4] L.C. Berselli and V. John. On the comparison of a commutation error and the Reynolds stress tensor for flows obeying a wall law. Preprint 18, Universit` di a Pisa, Dipartimento di Matematica Applicata “U.Dini”, 2004. [5] S.S. Collis. Monitoring unresolved scales in multiscale turbulence modeling. Physics of Fluids, 13:1800 – 1806, 2001. [6] A. Dunca, V. John, and W.J. Layton. The commutation error of the space averaged Navier-Stokes equations on a bounded domain. J. Math. Fluid Mech., 2003. accepted for publication. [7] G.P. Galdi. An introduction to the Navier-Stokes initial-boundary value problem. In G.P. Galdi, J.G. Heywood, and R. Rannacher, editors, Fundamental Directions in Mathematical Fluid Dynamics, pages 1 – 70. Birkh¨user, 2000. a [8] J.-L. Guermond. Stabilization of Galerkin approximations of transport equations by subgrid modeling. M2AN, 33:1293 – 1316, 1999. [9] J.G. Heywood and R. Rannacher. Finite element approximation of the nonstationary Navier-Stokes problem I: Regularity of solutions and second order error estimates for spatial discretization. SIAM J. Numer. Anal., 19:275 – 311, 1982. 18 [10] J.G. Heywood and R. Rannacher. Finite element approximation of the nonstationary Navier-Stokes problem III: Smoothing property and higher order estimates for spatial discretization. SIAM J. Num. Anal., 25:489 – 512, 1988. [11] T.J. Hughes, L. Mazzei, and K.E. Jansen. Large eddy simulation and the variational multiscale method. Comput. Visual. Sci., 3:47 – 59, 2000. [12] V. John. Large Eddy Simulation of Turbulent Incompressible Flows. Analytical and Numerical Results for a Class of LES Models, volume 34 of Lecture Notes in Computational Science and Engineering. Springer-Verlag Berlin, Heidelberg, New York, 2003. [13] V. John and S. Kaya. A finite element variational multiscale method for the Navier-Stokes equations. SIAM J. Sci. Comp., 2004. in press. [14] V. John and W.J. Layton. Analysis of numerical errors in large eddy simulation. SIAM J. Numer. Anal., 40:995 – 1020, 2002. [15] S. Kaya and W.J. Layton. Subgrid-scale eddy viscosity models are variational multiscale methods. Technical Report TR-MATH 03-05, University of Pittsburgh, 2003. [16] W. Layton and L. Tobiska. A two-level method with backtracking for the Navier-Stokes equations. SIAM J. Numer. Anal., 35(5):2035 – 2054, 1998. [17] P. Sagaut. Large Eddy Simulation for Incompressible Flows. Springer-Verlag, Berlin, Heidelberg New York, 2nd edition, 2003. [18] H. Sohr. The Navier-Stokes Equations, An Elementary Functional Analytic Approach. Birh¨user Advanced Texts. Birkh¨user Verlag Basel, Boston, Berlin, a a 2001. [19] R. Temam. Navier-Stokes Equations. Theory and Numerical Analysis, volume 2 of Studies in Mathematics and Its Applications. North-Holland Publishing Company, Amsterdam, New York, Oxford, 1977. 19