Advances in Numerical Analysis
Volume 2009 (2009), Article ID 494829, 20 pages
doi:10.1155/2009/494829
Research Article

A Subgrid Model for the Time-Dependent Navier-Stokes Equations

Yan Zhang and Yinnian He

Faculty of Science, Xi'an Jiaotong University, Shanxi 710049, China

Received 22 November 2008; Revised 5 May 2009; Accepted 25 June 2009

Academic Editor: Weimin Han

Copyright © 2009 Yan Zhang and Yinnian He. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

We propose a stabilized subgrid finite-element method for the two-dimensional (2D) nonstationary incompressible Naver-Stokes equation (NSE). This method yields a subgrid eddy viscosity which does not act on the large flow structures. The proposed eddy viscous term is constructed by a fluctuation operator based on an 𝐿 2 -projection. The fluctuation operator can be implemented by the 𝐿 2 -projection from high-order interpolation finite-element spaces to the low-order interpolation finite-element spaces. In this paper, 𝑃 2 / 𝑃 1 mixed finite-element spaces are adopted to implement the calculation and the analysis. The error analysis is given based on some regular assumptions. Finally, in the part of numerical tests, the numerical computations show that the numerical results agree with theoretical analysis very well. Meanwhile, the numerical investigations demonstrate that the proposed subgrid is very effective for high Reynolds number fluid flows and superior to other referred subgrid models.

1. Introduction

In this paper, we focus on formulating a subgrid eddy viscosity method for the nonstationary incompressible Navier-Stokes equations. For the subgrid method, we must admit that there exists a scale separation between large and small scales and this model can be viewed as a viscous correction for large-scale fluid flows. Meanwhile, for laminal fluid flows, the added subgrid viscosity term should not affect the large-scale structures of fluid flow field and should tend to vanish. This kind of subgrid method is just a flexible and effective method for high Reynolds number fluid flows.

It is well known that for most problems of fluid flows, the numerical algorithms capturing all scales of fluid flows are impossible and there exist several scales that span from the large-scales to the small Kolmogorov scales which can hardly be resolved by state-of-the-art computers for most engineering problems very efficiently. For the convection dominating fluid flows, we need to consider the dispersive effects of unresolved scales on resolved scales. The eddy viscosity models are often utilized to solve the convection-dominating or high Reynolds number NSE by engineers, which have been achieved many successes in engineering applications [1]. These kinds of models are firstly proposed by Frisch and Orszag [2], developed by Iliescu and Layton [3], and introduced a dissipative mechanism by Smagorinsky [4]. At present, these models have been further improved by various numerical methods [5–7]. In some recent models, the scale-separation subgrid terms are constructed by two different finite-element spaces based on a two-hierarchy mesh structure. Recently, Hughes et al. have proposed a variational multscale method (VMM) in which the diffusion acts only on the finest resolved scales. Generally, there exist different ways to define coarse and small scales under VMM frameworks [8]. According to the idea of VMM, the referred subgrid methods are variational multiscale methods.

In this paper, we will implement a novel way to circumvent the dispersive effects from small resolve scales. The subgrid term is established in a complement space of a low-order interpolation finite-element space in a high-order interpolation finite-element space. The complement space is established by an 𝐿 2 -projection decomposition of flow strain tensors. This method is easy to be implemented on a one-level mesh. This subgrid term does not act on the large-scale flow structures. For laminar low-Reynolds number flows, the action of this subgrid term tends to vanish. Numerical investigations demonstrate that the proposed model is effective and flexible for fluid flows of high Reynolds numbers ( R e = 1 0 5 and R e = 1 0 6 ).

The outline of the paper is organized as follows. In the following section, we introduce the nonstationary incompressible Navier-Stokes equations and the corresponding function settings. In Section 3, the subgrid viscous term is introduced into NSEs and the standard Galerkin discretization of the Navier-Stokes problem is given. In Section 4, we show the results of the error estimates. Some numerical results are presented in Section 5, which show the correctness and efficiency of the methods. Finally, we give some conclusions.

2. Navier-Stokes Equations and Functional Settings

Let Ω ⊂ ℝ 2 be a polygonal domain with Lipschitz continuous boundary Γ = 𝜕 Ω . We consider the time-dependent Navier-Stokes equations for incompressible flow as follows: 𝑢 𝑡 ] [ ] [ ] − 𝜈 Δ 𝑢 + ∇ 𝑝 + ( 𝑢 ⋅ ∇ ) 𝑢 = 𝑓 , i n ( 0 , 𝑇 × Ω , d i v 𝑢 = 0 , i n 0 , 𝑇 × Ω , 𝑢 = 0 , i n 0 , 𝑇 × Γ , 𝑢 ( 0 , 𝑥 ) = 𝑢 0 , i n Ω , Ω ] , 𝑝 𝑑 𝑥 = 0 , i n ( 0 , 𝑇 ( 2 . 1 ) where [ 0 , 𝑇 ] a finite-time interval, 𝑢 = ( 𝑢 1 , 𝑢 2 ) represents the velocity vector, 𝑝 is the pressure, 𝑓 is the body force, and 𝜈 > 0 is the viscosity.

We introduce the following functional settings: 𝑋 ∶ = 𝐻 1 0 ( Ω ) 2 , 𝑉 ∶ = { 𝑣 ∈ 𝑋 , d i v 𝑣 = 0 } , 𝑄 ∶ = 𝐿 2 0 ( Ω ) = 𝑞 ∈ 𝐿 2 ( Ω ) ∶ Ω . 𝑞 d 𝑥 = 0 ( 2 . 2 ) We denote by ( ⋅ , ⋅ ) and ‖ ⋅ ‖ 0 the inner product and norm in 𝐿 2 ( Ω ) or 𝐿 2 ( Ω ) 2 . The space 𝐻 𝑘 ( Ω ) or 𝐻 𝑘 ( Ω ) 2 denotes the standard Sobolev spaces within norm ‖ ⋅ ‖ 𝑘 and seminorm | ⋅ | 𝑘 . The space 𝐻 1 0 ( Ω ) or 𝐻 1 0 ( Ω ) 2 is equipped with the following scalar product and norm: ( ( 𝑢 , 𝑣 ) ) = ( ∇ 𝑢 , ∇ 𝑣 ) , | 𝑢 | 1 = ‖ ∇ 𝑢 ‖ 0 = ( ( 𝑢 , 𝑢 ) ) 1 / 2 . ( 2 . 3 ) For convenience, we introduce the following bilinear form 𝑎 ( ⋅ , ⋅ ) on 𝑋 × 𝑋 : 𝑎 ( 𝑢 , 𝑢 ) = ( ( 𝑢 , 𝑢 ) ) , ∀ 𝑢 ∈ 𝑋 , ( 2 . 4 ) and 𝑑 ( ⋅ , ⋅ ) on 𝑋 × 𝑄 defined by 𝑑 ( 𝑣 , 𝑞 ) = − ( 𝑣 , ∇ 𝑞 ) = ( 𝑞 , d i v 𝑣 ) , ∀ 𝑣 ∈ 𝑋 , 𝑞 ∈ 𝑄 . ( 2 . 5 ) The trilinear term is defined by 1 𝑏 ( 𝑢 ; 𝑣 , 𝑤 ) = ( ( 𝑢 ⋅ ∇ ) 𝑣 , 𝑤 ) + 2 = 1 ( ( d i v 𝑢 ) 𝑣 , 𝑤 ) 2 1 ( ( 𝑢 ⋅ ∇ ) 𝑣 , 𝑤 ) − 2 ( ( 𝑢 ⋅ ∇ ) 𝑤 , 𝑣 ) , ∀ 𝑢 , 𝑣 , 𝑤 ∈ 𝑋 , ( 2 . 6 ) which is the skew-symmetric form of the convective term. It is easy to gain 𝑏 ( 𝑢 ; 𝑣 , 𝑤 ) = − 𝑏 ( 𝑢 ; 𝑤 , 𝑣 ) . ( 2 . 7 ) Also, we have the following estimates [9]: | | 𝑏 | | ( 𝑢 , 𝑣 , 𝑤 ) ≤ 𝐶 ‖ 𝑢 ‖ 0 1 / 2 ‖ ∇ 𝑢 ‖ 0 1 / 2 ‖ ∇ 𝑣 ‖ 0 ‖ ∇ 𝑤 ‖ 0 , ∀ 𝑢 , 𝑣 , 𝑤 ∈ 𝑋 , ( 2 . 8 ) where 𝐶 is a positive constant depending only on the domain Ω , which stands for different values at different occurrences. The weak formulation of (2.1) reads as follows.

Find ( 𝑢 , 𝑝 ) ∈ ( 𝑋 , 𝑄 ) such that 𝑢 𝑡 , 𝑣 + 𝜈 𝑎 ( 𝑢 , 𝑣 ) + 𝑏 ( 𝑢 ; 𝑢 , 𝑣 ) − 𝑑 ( 𝑣 , 𝑝 ) = ( 𝑓 , 𝑣 ) , ∀ 𝑣 ∈ 𝑋 , 𝑑 ( 𝑢 , 𝑞 ) = 0 , ∀ 𝑞 ∈ 𝑄 , ( 2 . 9 ) and 𝑢 ( 0 , 𝑥 ) = 𝑢 0 ( 𝑥 ) ∈ 𝑋 .

For the finite-element analysis, we need some regularity assumptions of the Navier-Stokes equations and the solution.

Theorem 2.1 ([7, 10]). One assumes that 𝐿 𝑓 ∈ 2 0 , 𝑇 ; 𝐿 2 2 , 𝑢 0 ∈ 𝑋 , ( 2 . 1 0 ) and that (2.9) possesses a solution ( 𝑢 , 𝑝 ) with 𝐿 ∇ 𝑢 ∈ 4 0 , 𝑇 ; 𝐿 2 2 × 2 , 𝑢 𝑡 ∈ 𝐿 2 0 , 𝑇 ; 𝐻 − 1 2 , 𝑝 ∈ 𝐿 4 0 , 𝑇 ; 𝐿 2 . ( 2 . 1 1 )

Note. These assumptions imply that the solution of (2.9) is unique. For simplicity, let 𝑓 = 𝑓 ℎ . In addition, we assume that Ω has a polygonal boundary such that no boundary approximation in the application of the finite-element method becomes necessary.

3. Discretization of Navier-Stokes Equations and Subgrid Model

We give a family 𝜏 ℎ , which is a partition of Ω into triangles or quadrilaterals, assumed to be regular in the usual sense [11]. The diameter of the cell 𝐾 is denoted by ℎ 𝐾 . The mesh parameter ℎ describes the maximum diameter of the cells 𝐾 ∈ 𝜏 ℎ .

We introduce the finite-dimensional subspace 𝑋 ℎ and 𝑄 ℎ : 𝑋 ℎ 𝑣 ∶ = ℎ ∈ 𝐶 0 ( Ω ) 2 ∩ 𝑋 ∶ 𝑣 ℎ | | 𝑘 ∈ 𝑃 2 ( 𝐾 ) 2 , ∀ 𝐾 ∈ 𝜏 ℎ , 𝑄 ℎ 𝑞 ∶ = ℎ ∈ 𝑄 ∶ 𝑞 ℎ | | 𝑘 ∈ 𝑃 1 ( 𝐾 ) , ∀ 𝐾 ∈ 𝜏 ℎ . ( 3 . 1 ) We define the space of discretely divergence free functions denoted by 𝑉 ℎ : 𝑉 ℎ 𝑣 ∶ = ℎ ∈ 𝑋 ℎ 𝑣 ∶ 𝑑 ℎ , 𝑞 ℎ = 0 , ∀ 𝑞 ℎ ∈ 𝑄 ℎ . ( 3 . 2 ) Let ( 𝑢 , 𝑝 ) ∈ ( 𝐻 3 ( Ω ) 2 ∩ 𝑋 , 𝐻 2 ( Ω ) ∩ 𝑄 ) , 𝑡 ∈ [ 0 , 𝑇 ] be arbitrary, let ̃ 𝑢 ℎ ∈ 𝑉 ℎ be a projection of 𝑢 , and let the following approximation properties hold [12]: ‖ 𝜂 ‖ 0 + ℎ ‖ ∇ 𝜂 ‖ 0 ≤ 𝐶 ℎ 3 ‖ 𝑢 ‖ 3 + 𝜈 − 1 ‖ 𝑝 ‖ 2 , ‖ ‖ 𝜂 𝑡 ‖ ‖ 0 ‖ ‖ + ℎ ∇ 𝜂 𝑡 ‖ ‖ 0 ≤ 𝐶 ℎ 3 ‖ 𝑢 ‖ 3 + 𝜈 − 1 ‖ 𝑝 ‖ 2 , ( 3 . 3 ) here, 𝜂 = 𝑢 − ̃ 𝑢 ℎ , the constants depend only on Ω . The regularity assumptions (2.11) imply that 𝐿 ∇ 𝜂 ∈ 4 0 , 𝑇 ; 𝐿 2 2 × 2 . ( 3 . 4 ) Meanwhile, the velocity-pressure pair in ( 𝑋 ℎ , 𝑄 ℎ ) satisfies the following discrete inf-sup condition [13]: i n f 𝑞 ℎ ∈ 𝑄 ℎ s u p 𝑣 ℎ ∈ 𝑋 ℎ 𝑑 𝑣 ℎ , 𝑞 ℎ ‖ ‖ 𝑞 ℎ ‖ ‖ 0 | | 𝑣 ℎ | | 1 ≥ 𝛽 > 0 . ( 3 . 5 )

Remark 3.1 ([14, 15]). Let Π ∶ 𝑄 → 𝑅 0 be the standard 𝐿 2 -projection with the following properties: 𝑞 , 𝑞 ℎ = Π 𝑞 , 𝑞 ℎ , ∀ 𝑞 ∈ 𝑄 , 𝑞 ℎ ∈ 𝑅 0 , ‖ Π 𝑞 ‖ 0 ≤ 𝐶 ‖ 𝑞 ‖ 0 , ∀ 𝑞 ∈ 𝑄 , ‖ 𝑞 − Π 𝑞 ‖ 0 ≤ 𝐶 ℎ ‖ 𝑞 ‖ 1 , ∀ 𝑞 ∈ 𝐻 1 ( Ω ) ∩ 𝑄 , ( 3 . 6 ) where 𝑅 0 = { 𝑞 ℎ ∈ 𝑄 ∶ 𝑞 ℎ | 𝐾 i s c o n s t a n t , f o r a l l 𝐾 ∈ 𝜏 ℎ } . For a tensor 𝐓 , Π 𝐓 denotes that Π acts on each component of this tensor.

We know that for high Reynolds number fluid flows, when the fluid convection dominates fluid flow fields, under the finite-resolution of meshes, the flow becomes very instable. When the mesh scales cannot resolve the smallest scale in fluid flows, we must add some term into Navier-Stokes equations to smear out the effect from the unresolve scales. Here, we chose the following subgrid stabilization term to control the effect from the unresolved scales: 𝑀 𝑢 ℎ , 𝑣 ℎ = 𝛼 ( 𝐼 − Π ) ∇ 𝑢 ℎ , ( 𝐼 − Π ) ∇ 𝑣 ℎ , ( 3 . 7 ) where 𝛼 is the user-selected stabilization parameter and typically, 𝛼 = ℎ 𝑠 ( 𝑠 is a real number). The analogous stabilization is used to circumvent the pressure stabilization LBB condition for Stokes problems [16]. Since Π is an 𝐿 2 -projection, it follows for 𝑣 ℎ ∈ 𝑋 ℎ and ‖ ∇ 𝑣 ℎ ‖ 0 > 0 that 𝛼 ‖ ‖ ( 𝐼 − Π ) ∇ 𝑣 ℎ ‖ ‖ 2 0 ‖ ‖ = 𝛼 ∇ 𝑣 ℎ ‖ ‖ 2 0 − ‖ ‖ Π ∇ 𝑣 ℎ ‖ ‖ 2 0 ‖ ‖ = 𝛼 1 − Π ∇ 𝑣 ℎ ‖ ‖ 2 0 ‖ ‖ ∇ 𝑣 ℎ ‖ ‖ 2 0 ‖ ‖ ∇ 𝑣 ℎ ‖ ‖ 2 0 = ∶ 𝛼 a d d 𝑣 ℎ ‖ ‖ ∇ 𝑣 ℎ ‖ ‖ 2 0 . ( 3 . 8 ) In addition, from 0 ≤ ‖ Π ∇ 𝑣 ℎ ‖ 0 ≤ ‖ ∇ 𝑣 ℎ ‖ 0 , it follows that 0 ≤ 𝛼 a d d 𝑣 ℎ ≤ 𝛼 . ( 3 . 9 ) Note. If 𝑣 ℎ depends on 𝑡 , then 𝛼 a d d ( 𝑣 ℎ ) too. From (3.9) it follows that 𝛼 a d d ( 𝑣 ℎ ( 𝑡 , ⋅ ) ) ∈ 𝐿 ∞ ( 0 , 𝑇 ) if 𝛼 is bounded almost everywhere in the time interval. If ‖ ∇ 𝑣 ℎ ‖ 0 = 0 , then 𝑣 ℎ = 0 since 𝑣 ℎ ∈ 𝑋 ℎ . In this case, we set 𝛼 a d d ( 𝑣 ℎ ) = 0 . The analogous formula was proposed to define a reduced Reynolds number for a variational multiscale method of the Navier-Stokes equations [7].

Using the stabilization term 𝑀 ( 𝑢 ℎ , 𝑣 ℎ ) , we give the following stabilization finite-element discretization form of the variational form (2.9).

Find ( 𝑢 ℎ , 𝑞 ℎ ) ∈ ( 𝑋 ℎ , 𝑄 ℎ ) satisfying 𝑢 ℎ , 𝑡 , 𝑣 ℎ 𝑢 + 𝜈 𝑎 ℎ , 𝑣 ℎ 𝑢 + 𝑏 ℎ ; 𝑢 ℎ , 𝑣 ℎ 𝑣 − 𝑑 ℎ , 𝑝 ℎ 𝑢 + 𝑀 ℎ , 𝑣 ℎ = 𝑓 , 𝑣 ℎ , ∀ 𝑣 ℎ ∈ 𝑋 ℎ , 𝑑 𝑢 ℎ , 𝑞 ℎ = 0 , ∀ 𝑞 ℎ ∈ 𝑄 ℎ . ( 3 . 1 0 ) Under the inf-sup condition (3.5), formulation (3.10) is equivalent to the following problem.

Find 𝑢 ℎ ∈ 𝑉 ℎ such that 𝑢 ℎ , 𝑡 , 𝑣 ℎ 𝑢 + 𝜈 𝑎 ℎ , 𝑣 ℎ 𝑢 + 𝑏 ℎ , 𝑢 ℎ , 𝑣 ℎ 𝑢 + 𝑀 ℎ , 𝑣 ℎ = 𝑓 , 𝑣 ℎ , ∀ 𝑣 ℎ ∈ 𝑉 ℎ . ( 3 . 1 1 )

4. Error Analysis

The proof of the finite-element error estimate uses an approach by John and Kaya [7] and Heywood and Rannacher [12, 14]. The theoretical results of error analysis are classical [7]. We first show an outline [7].

(1)Prove stability of 𝑢 and 𝑢 ℎ , that is, certain norms of 𝑢 and 𝑢 ℎ are bounded by the data of the problem: 𝑓 , 𝑢 0 , 𝜈 .(2)Derive an error equation by subtracting (3.11) from (2.9) for test functions from 𝑉 ℎ . Split the error 𝑒 into an approximation term 𝜂 and a remainder 𝜙 ℎ : 𝑒 = 𝑢 − 𝑢 ℎ = 𝑢 − ̃ 𝑢 ℎ − 𝑢 ℎ − ̃ 𝑢 ℎ = ∶ 𝜂 − 𝜙 ℎ , ( 4 . 1 ) where ̃ 𝑢 ℎ ∈ 𝑉 ℎ is a projection of 𝑢 which fulfills the estimates (3.3). Then take 𝜙 ℎ ∈ 𝑉 ℎ as test function in the error equation.(3)Estimate the right hand side of the error equation which has the following form 𝑑 ‖ ‖ 𝜙 𝑑 𝑡 ℎ ‖ ‖ 2 0 + 𝑓 1 𝑡 , 𝜙 ℎ ≤ 𝑓 2 𝑡 , 𝜙 ℎ + 𝑓 3 𝑡 , 𝜙 ℎ ‖ ‖ 𝜙 ℎ ‖ ‖ 2 0 , ( 4 . 2 ) where all functions are nonnegative for all 𝑡 ∈ [ 0 , 𝑇 ] .(4)Apply Gronwall's lemma to (4.2), that is, derive all the functions in (4.2) belong to 𝐿 1 ( 0 , 𝑇 ) and get an estimate of 𝜙 ℎ .(5)Derive the error estimate of 𝑒 by applying the triangle inequality to (4.1).

Along these lines, the estimate will be proved [7]. This error estimate uses the parameter 𝛼 a d d defined in (3.8). The proving method is based on the classical scheme [7]. However, we still give the details of theoretical analysis for completeness. Firstly, we present the stability of 𝑢 and 𝑢 ℎ .

Lemma 4.1. The solution 𝑢 ℎ of the finite-element problem (3.10) fulfills 𝑢 ℎ ∈ ( 𝐿 ∞ ( 0 , 𝑇 ; 𝐿 2 ) ) 2 and ∇ 𝑢 ℎ ∈ ( 𝐿 2 ( 0 , 𝑇 ; 𝐿 2 ) ) 2 × 2 . The velocity solution of the continuous problem (2.9) fulfills 𝑢 ∈ ( 𝐿 ∞ ( 0 , 𝑇 ; 𝐿 2 ) ) 2 and ∇ 𝑢 ∈ ( 𝐿 2 ( 0 , 𝑇 ; 𝐿 2 ) ) 2 × 2 .

Proof. The proof for the 𝑢 ℎ and 𝑢 is very similar. We will show the result for 𝑢 ℎ . Setting 𝑣 ℎ = 𝑢 ℎ in (3.11), use the skew-symmetric form of 𝑏 ( ⋅ , ⋅ , ⋅ ) , (3.8), and integrate over ( 0 , 𝑡 ) with 𝑡 ≤ 𝑇 , we get 1 2 ‖ ‖ 𝑢 ℎ ( ‖ ‖ 𝑡 ) 2 0 + 𝑡 0 𝜈 + 𝛼 a d d 𝑢 ℎ ( ‖ ‖ 𝜏 ) ∇ 𝑢 ℎ ( ‖ ‖ 𝜏 ) 2 0 ≤ 1 𝑑 𝜏 2 ‖ ‖ 𝑢 0 , ℎ ‖ ‖ 2 0 + 𝑡 0 ‖ 𝑓 ( 𝜏 ) ‖ 𝐻 − 1 ‖ ‖ ∇ 𝑢 ℎ ‖ ‖ ( 𝜏 ) 0 ≤ 1 𝑑 𝜏 2 ‖ ‖ 𝑢 0 , ℎ ‖ ‖ 2 0 + 𝐶 𝜈 ‖ 𝑓 ‖ 2 𝐿 2 ( 0 , 𝑡 ; 𝐻 − 1 ) + 𝑡 0 𝜈 + 𝛼 a d d 𝑢 ℎ ( 𝜏 ) 2 ‖ ‖ ∇ 𝑢 ℎ ‖ ‖ ( 𝜏 ) 2 0 𝑑 𝜏 . ( 4 . 3 ) Here 𝑢 0 , ℎ is the value of 𝑢 ℎ at 𝑡 = 0 .

Subtraction of the last term and the regularity (2.10) gives ∇ 𝑢 ℎ ∈ ( 𝐿 2 ( 0 , 𝑇 ; 𝐿 2 ) ) 2 × 2 . Taking the supremum of 𝑡 ∈ ( 0 , 𝑇 ) gives the statement 𝑢 ℎ ∈ ( 𝐿 ∞ ( 0 , 𝑇 ; 𝐿 2 ) ) 2 .

Now, Step (2) of the proof is carried out by a strightforward calculation. We obtain 1 2 𝑑 ‖ ‖ 𝜙 𝑑 𝑡 ℎ ‖ ‖ 2 0 + 𝜈 + 𝛼 a d d 𝜙 ℎ ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 2 0 = 𝜂 𝑡 , 𝜙 ℎ + 𝜈 ∇ 𝜂 , ∇ 𝜙 ℎ + 𝛼 ( 𝐼 − Π ) ∇ 𝜂 , ( 𝐼 − Π ) ∇ 𝜙 ℎ + 𝑏 𝑢 , 𝑢 , 𝜙 ℎ 𝑢 − 𝑏 ℎ , 𝑢 ℎ , 𝜙 ℎ − 𝛼 ( 𝐼 − Π ) ∇ 𝑢 , ( 𝐼 − Π ) ∇ 𝜙 ℎ − 𝑝 − 𝜆 ℎ , ∇ ⋅ 𝜙 ℎ ( 4 . 4 ) with arbitrary 𝜆 ℎ ∈ 𝑄 ℎ .

To get the form of (4.2), the terms on the right-hand side of (4.4) have to be estimated. Using the Cauchy-Schwarz inequality, Young's inequality, and (3.8), we get 𝜂 𝑡 , 𝜙 ℎ ≤ ‖ ‖ 𝜂 𝑡 ‖ ‖ 𝐻 − 1 ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 0 ≤ 𝜈 + 𝛼 a d d 𝜙 ℎ 8 ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 2 0 + 𝐶 𝜈 + 𝛼 a d d 𝜙 ℎ ‖ ‖ 𝜂 𝑡 ‖ ‖ 2 𝐻 − 1 , 𝜈 ∇ 𝜂 , ∇ 𝜙 ℎ ≤ 𝜈 ‖ ∇ 𝜂 ‖ 0 ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 0 ≤ 𝜈 8 ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 2 0 + 2 𝜈 ‖ ∇ 𝜂 ‖ 2 0 , 𝑝 − 𝜆 ℎ , ∇ ⋅ 𝜙 ℎ ≤ ‖ 𝑝 − 𝜆 ℎ ‖ 0 ‖ ∇ ⋅ 𝜙 ℎ ‖ 0 ≤ 𝜈 + 𝛼 a d d 𝜙 ℎ 8 ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 2 0 + 𝐶 𝜈 + 𝛼 a d d 𝜙 ℎ ‖ ‖ 𝑝 − 𝜆 ℎ ‖ ‖ 2 0 , 𝛼 ( 𝐼 − Π ) ∇ 𝜂 , ( 𝐼 − Π ) ∇ 𝜙 ℎ ≤ 𝛼 ‖ ‖ 1 6 ( 𝐼 − Π ) ∇ 𝜙 ℎ ‖ ‖ 2 0 ‖ ‖ + 4 𝛼 ( 𝐼 − Π ) ∇ 𝜂 2 0 = 𝛼 a d d 𝜙 ℎ ‖ ‖ 1 6 ∇ 𝜙 ℎ ‖ ‖ 2 0 + 4 𝛼 a d d ( 𝜂 ) ‖ ∇ 𝜂 ‖ 2 0 , 𝛼 ( 𝐼 − Π ) ∇ 𝑢 , ( 𝐼 − Π ) ∇ 𝜙 ℎ ≤ 𝛼 ‖ ( 𝐼 − Π ) ∇ 𝑢 ‖ 0 ‖ ‖ ( 𝐼 − Π ) ∇ 𝜙 ℎ ‖ ‖ 0 ≤ 𝛼 ‖ ‖ 1 6 ( 𝐼 − Π ) ∇ 𝜙 ℎ ‖ ‖ 2 0 + 4 𝛼 ‖ ( 𝐼 − Π ) ∇ 𝑢 ‖ 2 0 = 𝛼 𝑎 𝑑 𝑑 𝜙 ℎ ‖ ‖ 1 6 ∇ 𝜙 ℎ ‖ ‖ 2 0 + 4 𝛼 ‖ ( 𝐼 − Π ) ∇ 𝑢 ‖ 2 0 . ( 4 . 5 ) The trilinear term is decomposed into three terms as follows: 𝑏 𝑢 , 𝑢 , 𝜙 ℎ 𝑢 − 𝑏 ℎ , 𝑢 ℎ , 𝜙 ℎ = 𝑏 𝜂 , 𝑢 , 𝜙 ℎ 𝜙 − 𝑏 ℎ , 𝑢 , 𝜙 ℎ 𝑢 + 𝑏 ℎ , 𝜂 , 𝜙 ℎ . ( 4 . 6 ) Using (2.8) and Young's inequality, we have 𝑏 𝜂 , 𝑢 , 𝜙 ℎ ≤ 𝐶 ‖ 𝜂 ‖ 0 1 / 2 ‖ ∇ 𝜂 ‖ 0 1 / 2 ‖ ∇ 𝑢 ‖ 0 ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 0 ≤ 𝜈 + 𝛼 a d d 𝜙 ℎ 8 ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 2 0 + 𝐶 𝜈 + 𝛼 a d d 𝜙 ℎ ‖ 𝜂 ‖ 0 ‖ ∇ 𝜂 ‖ 0 ‖ ∇ 𝑢 ‖ 2 0 , 𝑏 𝜙 ℎ , 𝑢 , 𝜙 ℎ ‖ ‖ 𝜙 ≤ 𝐶 ℎ ‖ ‖ 0 1 / 2 ‖ ∇ 𝑢 ‖ 0 ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 0 3 / 2 ≤ 𝜈 + 𝛼 a d d 𝜙 ℎ 8 ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 2 0 + 𝐶 𝜈 + 𝛼 a d d 𝜙 ℎ 3 ‖ ‖ 𝜙 ℎ ‖ ‖ 2 0 ‖ ∇ 𝑢 ‖ 4 0 , 𝑏 𝑢 ℎ , 𝜂 , 𝜙 ℎ ‖ ‖ 𝑢 ≤ 𝐶 ℎ ‖ ‖ 0 1 / 2 ‖ ‖ ∇ 𝑢 ℎ ‖ ‖ 0 1 / 2 ‖ ∇ 𝜂 ‖ 0 ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 0 ≤ 𝜈 + 𝛼 a d d 𝜙 ℎ 8 ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 2 0 + 𝐶 𝜈 + 𝛼 a d d 𝜙 ℎ ‖ ‖ 𝑢 ℎ ‖ ‖ 0 ‖ ‖ ∇ 𝑢 ℎ ‖ ‖ 0 ‖ ∇ 𝜂 ‖ 2 0 . ( 4 . 7 ) Collecting all the above terms, we obtain 1 2 𝑑 ‖ ‖ 𝜙 𝑑 𝑡 ℎ ‖ ‖ 2 0 + 𝜈 + 𝛼 a d d 𝜙 ℎ 4 ‖ ‖ ∇ 𝜙 ℎ ‖ ‖ 2 0 ≤ 𝐶 𝜈 + 𝛼 a d d 𝜙 ℎ ‖ ‖ 𝜂 𝑡 ‖ ‖ 2 𝐻 − 1 + 4 𝜈 + 𝛼 a d d ( 𝜂 ) ‖ ∇ 𝜂 ‖ 2 0 + 𝐶 𝜈 + 𝛼 a d d 𝜙 ℎ ‖ ‖ 𝑝 − 𝜆 ℎ ‖ ‖ 2 0 + 4 𝛼 ‖ ( 𝐼 − Π ) ∇ 𝑢 ‖ 2 0 + 𝐶 𝜈 + 𝛼 a d d 𝜙 ℎ ‖ 𝜂 ‖ 0 ‖ ∇ 𝜂 ‖ 0 ‖ ∇ 𝑢 ‖ 2 0 + ‖ ‖ 𝑢 ℎ ‖ ‖ 0 ‖ ‖ ∇ 𝑢 ℎ ‖ ‖ 0 ‖ ∇ 𝜂 ‖ 2 0 + 𝐶 𝜈 + 𝛼 a d d 𝜙 ℎ 3 ‖ ∇ 𝑢 ‖ 4 0 ‖ ‖ 𝜙 ℎ ‖ ‖ 2 0 . ( 4 . 8 ) We define R e r e d ∶ = 𝜈 + i n f ] 𝑡 ∈ ( 0 , 𝑇 𝛼 a d d 𝜙 ℎ ( 𝑡 ) − 1 . ( 4 . 9 ) It follows that R e r e d is smaller or equal than Re. Using (3.9) finishes Step (3) of the proof: 𝑑 ‖ ‖ 𝜙 𝑑 𝑡 ℎ ‖ ‖ 2 0 + R e − 1 r e d 2 ‖ ∇ 𝜙 ℎ ‖ 2 0 ≤ 𝐶 R e r e d ‖ 𝜂 𝑡 ‖ 2 𝐻 − 1 + R e − 1 + 𝛼 ‖ ∇ 𝜂 ‖ 2 0 + R e r e d ‖ 𝑝 − 𝜆 ℎ ‖ 2 0 + 𝛼 ‖ ( 𝐼 − Π ) ∇ 𝑢 ‖ 2 0 + R e r e d ‖ 𝜂 ‖ 0 ‖ ∇ 𝜂 ‖ 0 ‖ ∇ 𝑢 ‖ 2 0 + ‖ 𝑢 ℎ ‖ 0 ‖ ∇ 𝑢 ℎ ‖ 0 ‖ ∇ 𝜂 ‖ 2 0 + 𝐶 R e 3 r e d ‖ ∇ 𝑢 ‖ 4 0 ‖ 𝜙 ℎ ‖ 2 0 . ( 4 . 1 0 ) To perform Step (4) of the proof, we have to study the 𝐿 1 ( 0 , 𝑇 ) -regularity of the terms in (4.10). Let 𝑡 ∈ ( 0 , 𝑇 ] be arbitrary. Using Poincaré's inequality, the Cauchy-Schwarz inequality, then (2.11) and (3.4), we get 𝑡 0 ‖ 𝜂 ( 𝜏 ) ‖ 0 ‖ ∇ 𝜂 ( 𝜏 ) ‖ 0 ‖ ∇ 𝑢 ( 𝜏 ) ‖ 2 0 𝑑 𝜏 ≤ 𝐶 𝑡 0 ‖ ∇ 𝜂 ( 𝜏 ) ‖ 2 0 ‖ ∇ 𝑢 ( 𝜏 ) ‖ 2 0 𝑑 𝜏 ≤ 𝐶 ‖ ∇ 𝜂 ‖ 2 𝐿 4 ( 0 , 𝑡 ; 𝐿 2 ) ‖ ∇ 𝑢 ‖ 2 𝐿 4 ( 0 , 𝑡 ; 𝐿 2 ) < ∞ . ( 4 . 1 1 ) Similarly by Hölders inequality, Lemma 4.1, and (3.4), we have 𝑡 0 ‖ 𝑢 ℎ ( 𝜏 ) ‖ 0 ‖ ∇ 𝑢 ℎ ( 𝜏 ) ‖ 0 ‖ ∇ 𝜂 ( 𝜏 ) ‖ 2 0 𝑑 𝜏 ≤ ‖ 𝑢 ℎ ‖ 𝐿 ∞ ( 0 , 𝑡 ; 𝐿 2 ) 𝑡 0 ‖ ∇ 𝑢 ℎ ( 𝜏 ) ‖ 0 ‖ ∇ 𝜂 ( 𝜏 ) ‖ 2 0 𝑑 𝜏 ≤ ‖ 𝑢 ℎ ‖ 𝐿 ∞ ( 0 , 𝑡 ; 𝐿 2 ) ‖ ∇ 𝑢 ℎ ‖ 2 𝐿 2 ( 0 , 𝑡 ; 𝐿 2 ) ‖ ∇ 𝜂 ‖ 2 𝐿 4 ( 0 , 𝑡 ; 𝐿 2 ) ≤ 𝐶 R e 1 / 2 ‖ 𝑢 0 , ℎ ‖ 2 0 + R e 3 / 2 ‖ 𝑓 ‖ 2 𝐿 2 0 , 𝑡 ; 𝐻 − 1 ‖ ∇ 𝜂 ‖ 2 𝐿 4 0 , 𝑡 ; 𝐿 2 < ∞ . ( 4 . 1 2 ) The 𝐿 1 ( 0 , 𝑇 ) -regularity of other terms is a direct consequence of (2.11), (3.3), and (3.4). The application of Gronwall's inequality and the last step of the proof give the following theorem.

Theorem 4.2. Let ( 𝑢 , 𝑝 ) ∈ 𝑋 × 𝑄 be the solution of (2.9) and let 𝑢 ℎ ∈ 𝑉 ℎ be the solution of (3.11). Let the regularity assumptions (2.11) be fulfilled, let ̃ 𝑢 ℎ be a projection of 𝑢 into 𝑉 ℎ such that 𝜂 fulfills (3.3). Let the reduced Reynolds number R e r e d be defined in (4.9). Then, the error 𝑒 = 𝑢 − 𝑢 ℎ satisfies, for 𝑇 ≥ 0 , ‖ ( 𝑢 − 𝑢 ℎ ) ( 𝑇 ) ‖ 2 0 + R e − 1 r e d 2 ‖ ∇ ( 𝑢 − 𝑢 ℎ ) ‖ 2 𝐿 2 ( 0 , 𝑇 ; 𝐿 2 ) ≤ 𝐶 i n f ̃ 𝑢 ℎ ∈ 𝐿 4 ( 0 , 𝑇 ; 𝑉 ℎ ) , 𝜆 ℎ ∈ 𝐿 2 ( 0 , 𝑇 ; 𝑄 ℎ ) × ‖ ( 𝑢 − ̃ 𝑢 ℎ ) ( 𝑇 ) ‖ 2 0 + R e − 1 r e d ‖ ∇ ( 𝑢 − ̃ 𝑢 ℎ ) ‖ 2 𝐿 2 ( 0 , 𝑇 ; 𝐿 2 ) + e x p 𝐶 R e 3 r e d ‖ ∇ 𝑢 ‖ 4 𝐿 4 ( 0 , 𝑇 ; 𝐿 2 ) × ‖ 𝑢 0 − 𝑢 0 , ℎ ‖ 2 0 + ‖ 𝑢 0 − ̃ 𝑢 0 , ℎ ‖ 2 0 + R e r e d ‖ ( 𝑢 − ̃ 𝑢 ℎ ) 𝑡 ‖ 2 𝐿 2 ( 0 , 𝑇 ; 𝐻 − 1 ) + ‖ 𝑝 − 𝜆 ℎ ‖ 2 𝐿 2 ( 0 , 𝑇 ; 𝐿 2 ) + ‖ ∇ ( 𝑢 − ̃ 𝑢 ℎ ) ‖ 2 𝐿 4 ( 0 , 𝑇 ; 𝐿 2 ) ‖ ∇ 𝑢 ‖ 2 𝐿 4 ( 0 , 𝑇 ; 𝐿 2 ) + R e 1 / 2 ‖ 𝑢 0 , ℎ ‖ 2 0 + R e 3 / 2 ‖ 𝑓 ‖ 2 𝐿 2 ( 0 , 𝑇 ; 𝐻 − 1 ) ‖ ∇ ( 𝑢 − ̃ 𝑢 ℎ ) ‖ 2 𝐿 4 ( 0 , 𝑇 ; 𝐿 2 ) + R e − 1 + 𝛼 ‖ ∇ ( 𝑢 − ̃ 𝑢 ℎ ) ‖ 2 𝐿 2 ( 0 , 𝑇 ; 𝐿 2 ) + 𝛼 ‖ ( 𝐼 − Π ) ∇ 𝑢 ‖ 2 𝐿 2 ( 0 , 𝑇 ; 𝐿 2 ) . ( 4 . 1 3 )

Corollary 4.3. From Theorem 4.2, the approximation result (3.3) and Remark 3.1, one has ‖ ( 𝑢 − 𝑢 ℎ ) ( 𝑇 ) ‖ 2 0 + R e − 1 r e d 2 ‖ ∇ ( 𝑢 − 𝑢 ℎ ) ‖ 2 𝐿 2 ( 0 , 𝑇 ; 𝐿 2 ) ℎ ≤ 𝐶 4 + 𝛼 ℎ 2 . ( 4 . 1 4 ) In particular, ‖ ( 𝑢 − 𝑢 ℎ ) ( 𝑇 ) ‖ 0 ≤ 𝐶 ℎ 2 , i f 𝛼 = ℎ 2 . ( 4 . 1 5 )

Remark 4.4. A finite-element error estimate for the 𝐿 2 ( Ω ) -error in the pressure can also be derived following Heywood and Rannacher [12]. Using the inf-sup condition (3.5) and the estimates of (3.3), the pressure error can be estimated by approximation errors and the velocity error ‖ ( 𝑢 − 𝑢 ℎ ) ( 𝑡 ) ‖ 0 . Theorem 4.2 finishes the error estimate for the pressure. Since the analysis is lengthy and follows closely ([12]), we will not present it here.

5. Numerical Tests

Firstly, we give the algorithm used to deal with the nonlinear term and the subgrid eddy viscosity term. Given ( 𝑢 ℎ 𝑛 − 1 , 𝑝 ℎ 𝑛 − 1 ) ∈ ( 𝑋 ℎ , 𝑄 ℎ ) , we find ( 𝑢 𝑛 ℎ , 𝑝 𝑛 ℎ ) ∈ ( 𝑋 ℎ , 𝑄 ℎ ) satisfying 𝑑 𝑡 𝑢 𝑛 ℎ , 𝑣 ℎ 𝑢 + 𝜈 𝑎 𝑛 ℎ , 𝑣 ℎ 𝑣 − 𝑑 ℎ , 𝑝 𝑛 ℎ 𝑢 + 𝑑 𝑛 ℎ , 𝑞 ℎ ( + 𝛼 𝐼 − Π ) ∇ 𝑢 𝑛 ℎ , ( 𝐼 − Π ) ∇ 𝑣 ℎ 𝑢 + 𝑏 𝑛 ℎ , 𝑢 𝑛 ℎ , 𝑣 ℎ = 𝑓 , 𝑣 ℎ , ( 5 . 1 ) where 𝑑 𝑡 𝑢 𝑛 ℎ = ( 1 / Δ 𝑡 ) ( 𝑢 𝑛 ℎ − 𝑢 ℎ 𝑛 − 1 ) and Δ 𝑡 is the time interval. Also, 𝑢 𝑛 ℎ is the finite-element approximate solutions of 𝑢 ℎ ( 𝑡 𝑛 = 𝑛 Δ 𝑡 ) ( 1 < 𝑛 < 𝑁 ) . The fixed point iteration is adopted to solve (5.1) based on the Newtonian method [17] and the iteration tolerance is set to equal 1 0 − 8 .

For calculating the subgrid term 𝛼 ( ( 𝐼 − Π ) ∇ 𝑢 𝑛 ℎ , ( 𝐼 − Π ) ∇ 𝑣 ℎ ) , the subgrid term is rewritten as follows: 𝛼 ( 𝐼 − Π ) ∇ 𝑢 𝑛 ℎ , ( 𝐼 − Π ) ∇ 𝑣 ℎ = 𝛼 ∇ 𝑢 𝑛 ℎ , ∇ 𝑣 ℎ − 𝛼 Π ∇ 𝑢 𝑛 ℎ , ∇ 𝑣 ℎ . ( 5 . 2 ) In order to reduce the computational work, the following linear time extrapolation method is used to calculate Π ∇ 𝑢 𝑛 ℎ : Π ∇ 𝑢 𝑛 ℎ = 2 Π ∇ 𝑢 ℎ 𝑛 − 1 − Π ∇ 𝑢 ℎ 𝑛 − 2 + 𝒪 𝛿 𝑡 2 . ( 5 . 3 )

5.1. Validate Convergence Rate by Taylor Vortex

It is essential to investigate the subgrid model (3.7) for low viscosity fluid flow and validate the flexibility and convergence orders of this model. So, we need to choose a true solution. We consider (2.1) on the domain Ω = [ 0 , 1 ] × [ 0 , 1 ] , with a body force obtained such that the following true solution is given by 𝑢 = ( 𝑢 1 , 𝑢 2 ) : 𝑢 1 ( 𝑥 , 𝑦 , 𝑡 ) = − c o s ( 𝑛 𝜋 𝑥 ) s i n ( 𝑛 𝜋 𝑦 ) e x p − 2 𝑛 2 𝜋 2 , 𝑢 𝜈 𝑡 2 ( 𝑥 , 𝑦 , 𝑡 ) = s i n ( 𝑛 𝜋 𝑥 ) c o s ( 𝑛 𝜋 𝑦 ) e x p − 2 𝑛 2 𝜋 2 , 1 𝜈 𝑡 𝑝 ( 𝑥 , 𝑦 , 𝑡 ) = − 4 ( c o s ( 2 𝑛 𝜋 𝑥 ) + c o s ( 2 𝑛 𝜋 𝑦 ) ) e x p − 2 𝑛 2 𝜋 2 . 𝜈 𝑡 ( 5 . 4 ) The viscosity is 𝜈 = 0 . 0 0 1 and the corresponding Reynolds number is 1 0 0 0 . The mesh scales we choose are ℎ = 1 / 𝑚 ( 𝑚 = 8 , 1 6 , 3 2 , 4 0 , 4