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

Authors: Yan Zhang; Yinnian He
Journal: Advances in Numerical Analysis
 (2009)

---

## 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 anL2-projection. The fluctuation operator can be implemented by the
L2-projection from high-order interpolation finite-element spaces to the low-order
interpolation finite-element spaces. In this paper, P2/P1 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.


---

## Body

## 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 anL2-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 (Re=105 and Re=106).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 Section3, 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:
(2.1)ut-νΔu+∇p+(u·∇)u=f,in(0,T]×Ω,divu=0,in[0,T]×Ω,u=0,in[0,T]×Γ,u(0,x)=u0,inΩ,∫Ωpdx=0,in(0,T],
where [0,T] a finite-time interval, u=(u1,u2) represents the velocity vector, p is the pressure, f is the body force, and ν>0 is the viscosity.We introduce the following functional settings:(2.2)X:=H01(Ω)2,V:={v∈X,divv=0},Q:=L02(Ω)={q∈L2(Ω):∫Ωqdx=0}.
We denote by (·,·) and ∥·∥0 the inner product and norm in L2(Ω) or L2(Ω)2. The space Hk(Ω) or Hk(Ω)2 denotes the standard Sobolev spaces within norm ∥·∥k and seminorm |·|k. The space H01(Ω) or H01(Ω)2 is equipped with the following scalar product and norm:
(2.3)((u,v))=(∇u,∇v),|u|1=∥∇u∥0=((u,u))1/2.
For convenience, we introduce the following bilinear form a(·,·) on X×X:
(2.4)a(u,u)=((u,u)),∀u∈X,
and d(·,·) on X×Q defined by
(2.5)d(v,q)=-(v,∇q)=(q,divv),∀v∈X,q∈Q.
The trilinear term is defined by
(2.6)b(u;v,w)=((u·∇)v,w)+12((divu)v,w)=12((u·∇)v,w)-12((u·∇)w,v),∀u,v,w∈X,
which is the skew-symmetric form of the convective term. It is easy to gain
(2.7)b(u;v,w)=-b(u;w,v).
Also, we have the following estimates [9]:
(2.8)|b(u,v,w)|≤C∥u∥01/2∥∇u∥01/2∥∇v∥0∥∇w∥0,∀u,v,w∈X,
where C 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(u,p)∈(X,Q) such that
(2.9)(ut,v)+νa(u,v)+b(u;u,v)-d(v,p)=(f,v),∀v∈X,d(u,q)=0,∀q∈Q,
and u(0,x)=u0(x)∈X.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.10)f∈(L2(0,T;L2))2,u0∈X,
and that (2.9) possesses a solution (u,p) with
(2.11)∇u∈(L4(0,T;L2))2×2,ut∈(L2(0,T;H-1))2,p∈L4(0,T;L2).Note. These assumptions imply that the solution of (2.9) is unique. For simplicity, let f=fh. 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τh, which is a partition of Ω into triangles or quadrilaterals, assumed to be regular in the usual sense [11]. The diameter of the cell K is denoted by hK. The mesh parameter h describes the maximum diameter of the cells K∈τh.We introduce the finite-dimensional subspaceXh and Qh:
(3.1)Xh:={vh∈(C0(Ω))2∩X:vh|k∈P2(K)2,∀K∈τh},Qh:={qh∈Q:qh|k∈P1(K),∀K∈τh}.
We define the space of discretely divergence free functions denoted by Vh:
(3.2)Vh:={vh∈Xh:d(vh,qh)=0,∀qh∈Qh}.
Let (u,p)∈(H3(Ω)2∩X,H2(Ω)∩Q), t∈[0,T] be arbitrary, let ũh∈Vh be a projection of u, and let the following approximation properties hold [12]:
(3.3)∥η∥0+h∥∇η∥0≤Ch3(∥u∥3+ν-1∥p∥2),∥ηt∥0+h∥∇ηt∥0≤Ch3(∥u∥3+ν-1∥p∥2),
here, η=u-ũh, the constants depend only on Ω. The regularity assumptions (2.11) imply that
(3.4)∇η∈(L4(0,T;L2))2×2.
Meanwhile, the velocity-pressure pair in (Xh,Qh) satisfies the following discrete inf-sup condition [13]:
(3.5)infqh∈Qhsupvh∈Xhd(vh,qh)∥qh∥0|vh|1≥β>0.Remark 3.1 ([14, 15]).
LetΠ:Q→R0 be the standard L2-projection with the following properties:
(3.6)(q,qh)=(Πq,qh),∀q∈Q,qh∈R0,∥Πq∥0≤C∥q∥0,∀q∈Q,∥q-Πq∥0≤Ch∥q∥1,∀q∈H1(Ω)∩Q,
where R0={qh∈Q:qh|Kisconstant,forallK∈τh}. For a tensor T, ΠT 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)M(uh,vh)=α((I-Π)∇uh,(I-Π)∇vh),
where α is the user-selected stabilization parameter and typically, α=hs (s is a real number). The analogous stabilization is used to circumvent the pressure stabilization LBB condition for Stokes problems [16]. Since Π is an L2-projection, it follows for vh∈Xh and ∥∇vh∥0>0 that
(3.8)α∥(I-Π)∇vh∥02=α(∥∇vh∥02-∥Π∇vh∥02)=α(1-∥Π∇vh∥02∥∇vh∥02)∥∇vh∥02=:αadd(vh)∥∇vh∥02.
In addition, from 0≤∥Π∇vh∥0≤∥∇vh∥0, it follows that
(3.9)0≤αadd(vh)≤α.
Note. If vh depends on t, then αadd(vh) too. From (3.9) it follows that αadd(vh(t,·))∈L∞(0,T) if α is bounded almost everywhere in the time interval. If ∥∇vh∥0=0, then vh=0 since vh∈Xh. In this case, we set αadd(vh)=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 termM(uh,vh), we give the following stabilization finite-element discretization form of the variational form (2.9).Find(uh,qh)∈(Xh,Qh) satisfying
(3.10)(uh,t,vh)+νa(uh,vh)+b(uh;uh,vh)-d(vh,ph)+M(uh,vh)=(f,vh),∀vh∈Xh,d(uh,qh)=0,∀qh∈Qh.
Under the inf-sup condition (3.5), formulation (3.10) is equivalent to the following problem.Finduh∈Vh such that
(3.11)(uh,t,vh)+νa(uh,vh)+b(uh,uh,vh)+M(uh,vh)=(f,vh),∀vh∈Vh.


## 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 ofu and uh, that is, certain norms of u and uh are bounded by the data of the problem: f,u0,ν.(2)
Derive an error equation by subtracting (3.11) from (2.9) for test functions from Vh. Split the error e into an approximation term η and a remainder ϕh:
(4.1)e=u-uh=(u-ũh)-(uh-ũh)=:η-ϕh,
where ũh∈Vh is a projection of u which fulfills the estimates (3.3). Then take ϕh∈Vh as test function in the error equation.(3)
Estimate the right hand side of the error equation which has the following form(4.2)ddt∥ϕh∥02+f1(t,ϕh)≤f2(t,ϕh)+f3(t,ϕh)∥ϕh∥02,
where all functions are nonnegative for all t∈[0,T].(4)
Apply Gronwall's lemma to (4.2), that is, derive all the functions in (4.2) belong to L1(0,T) and get an estimate of ϕh.(5)
Derive the error estimate ofe by applying the triangle inequality to (4.1).Along these lines, the estimate will be proved [7]. This error estimate uses the parameter αadd 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 u and uh.Lemma 4.1.
The solutionuh of the finite-element problem (3.10) fulfills uh∈(L∞(0,T;L2))2 and ∇uh∈(L2(0,T;L2))2×2. The velocity solution of the continuous problem (2.9) fulfills u∈(L∞(0,T;L2))2 and ∇u∈(L2(0,T;L2))2×2.Proof.
The proof for theuh and u is very similar. We will show the result for uh. Setting vh=uh in (3.11), use the skew-symmetric form of b(·,·,·), (3.8), and integrate over (0,t) with t≤T, we get
(4.3)12∥uh(t)∥02+∫0t(ν+αadd(uh(τ)))∥∇uh(τ)∥02dτ≤12∥u0,h∥02+∫0t∥f(τ)∥H-1∥∇uh(τ)∥0dτ≤12∥u0,h∥02+Cν∥f∥L2(0,t;H-1)2+∫0tν+αadd(uh(τ))2∥∇uh(τ)∥02dτ.
Here u0,h is the value of uh at t=0.Subtraction of the last term and the regularity (2.10) gives ∇uh∈(L2(0,T;L2))2×2. Taking the supremum of t∈(0,T) gives the statement uh∈(L∞(0,T;L2))2.Now, Step (2) of the proof is carried out by a strightforward calculation. We obtain(4.4)12ddt∥ϕh∥02+(ν+αadd(ϕh))∥∇ϕh∥02=(ηt,ϕh)+ν(∇η,∇ϕh)+α((I-Π)∇η,(I-Π)∇ϕh)+b(u,u,ϕh)-b(uh,uh,ϕh)-α((I-Π)∇u,(I-Π)∇ϕh)-(p-λh,∇·ϕh)
with arbitrary λh∈Qh.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
(4.5)(ηt,ϕh)≤∥ηt∥H-1∥∇ϕh∥0≤ν+αadd(ϕh)8∥∇ϕh∥02+Cν+αadd(ϕh)∥ηt∥H-12,ν(∇η,∇ϕh)≤ν∥∇η∥0∥∇ϕh∥0≤ν8∥∇ϕh∥02+2ν∥∇η∥02,(p-λh,∇·ϕh)≤∥p-λh∥0∥∇·ϕh∥0≤ν+αadd(ϕh)8∥∇ϕh∥02+Cν+αadd(ϕh)∥p-λh∥02,α((I-Π)∇η,(I-Π)∇ϕh)≤α16∥(I-Π)∇ϕh∥02+4α∥(I-Π)∇η∥02=αadd(ϕh)16∥∇ϕh∥02+4αadd(η)∥∇η∥02,α((I-Π)∇u,(I-Π)∇ϕh)≤α∥(I-Π)∇u∥0∥(I-Π)∇ϕh∥0≤α16∥(I-Π)∇ϕh∥02+4α∥(I-Π)∇u∥02=αadd(ϕh)16∥∇ϕh∥02+4α∥(I-Π)∇u∥02.
The trilinear term is decomposed into three terms as follows:
(4.6)b(u,u,ϕh)-b(uh,uh,ϕh)=b(η,u,ϕh)-b(ϕh,u,ϕh)+b(uh,η,ϕh).
Using (2.8) and Young's inequality, we have
(4.7)b(η,u,ϕh)≤C∥η∥01/2∥∇η∥01/2∥∇u∥0∥∇ϕh∥0≤ν+αadd(ϕh)8∥∇ϕh∥02+Cν+αadd(ϕh)∥η∥0∥∇η∥0∥∇u∥02,b(ϕh,u,ϕh)≤C∥ϕh∥01/2∥∇u∥0∥∇ϕh∥03/2≤ν+αadd(ϕh)8∥∇ϕh∥02+C(ν+αadd(ϕh))3∥ϕh∥02∥∇u∥04,b(uh,η,ϕh)≤C∥uh∥01/2∥∇uh∥01/2∥∇η∥0∥∇ϕh∥0≤ν+αadd(ϕh)8∥∇ϕh∥02+Cν+αadd(ϕh)∥uh∥0∥∇uh∥0∥∇η∥02.
Collecting all the above terms, we obtain
(4.8)12ddt∥ϕh∥02+ν+αadd(ϕh)4∥∇ϕh∥02≤[Cν+αadd(ϕh)∥ηt∥H-12+4(ν+αadd(η))∥∇η∥02+Cν+αadd(ϕh)∥p-λh∥02+4α∥(I-Π)∇u∥02+Cν+αadd(ϕh)(∥η∥0∥∇η∥0∥∇u∥02+∥uh∥0∥∇uh∥0∥∇η∥02)]+[C(ν+αadd(ϕh))3∥∇u∥04]∥ϕh∥02.
We define
(4.9)Rered:=(ν+inft∈(0,T]αadd(ϕh(t)))-1.
It follows that Rered is smaller or equal than Re. Using (3.9) finishes Step (3) of the proof:
(4.10)ddt∥ϕh∥02+Rered-12∥∇ϕh∥02≤C[Rered∥ηt∥H-12+(Re-1+α)∥∇η∥02+Rered∥p-λh∥02+α∥(I-Π)∇u∥02+Rered(∥η∥0∥∇η∥0∥∇u∥02+∥uh∥0∥∇uh∥0∥∇η∥02)]+CRered3∥∇u∥04∥ϕh∥02.
To perform Step (4) of the proof, we have to study the L1(0,T)-regularity of the terms in (4.10). Let t∈(0,T] be arbitrary. Using Poincaré's inequality, the Cauchy-Schwarz inequality, then (2.11) and (3.4), we get
(4.11)∫0t∥η(τ)∥0∥∇η(τ)∥0∥∇u(τ)∥02dτ≤C∫0t∥∇η(τ)∥02∥∇u(τ)∥02dτ≤C∥∇η∥L4(0,t;L2)2∥∇u∥L4(0,t;L2)2<∞.
Similarly by Hölders inequality, Lemma 4.1, and (3.4), we have
(4.12)∫0t∥uh(τ)∥0∥∇uh(τ)∥0∥∇η(τ)∥02dτ≤∥uh∥L∞(0,t;L2)∫0t∥∇uh(τ)∥0∥∇η(τ)∥02dτ≤∥uh∥L∞(0,t;L2)∥∇uh∥L2(0,t;L2)2∥∇η∥L4(0,t;L2)2≤C(Re1/2∥u0,h∥02+Re3/2∥f∥L2(0,t;H-1)2)∥∇η∥L4(0,t;L2)2<∞.
The L1(0,T)-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(u,p)∈X×Q be the solution of (2.9) and let uh∈Vh be the solution of (3.11). Let the regularity assumptions (2.11) be fulfilled, let ũh be a projection of u into Vh such that η fulfills (3.3). Let the reduced Reynolds number Rered be defined in (4.9). Then, the error e=u-uh satisfies, for T≥0,
(4.13)∥(u-uh)(T)∥02+Rered-12∥∇(u-uh)∥L2(0,T;L2)2≤Cinfũh∈L4(0,T;Vh),λh∈L2(0,T;Qh)×{∥(u-ũh)(T)∥02+Rered-1∥∇(u-ũh)∥L2(0,T;L2)2+exp(CRered3∥∇u∥L4(0,T;L2)4)×[∥u0-u0,h∥02+∥u0-ũ0,h∥02+Rered[∥(u-ũh)t∥L2(0,T;H-1)2+∥p-λh∥L2(0,T;L2)2+∥∇(u-ũh)∥L4(0,T;L2)2∥∇u∥L4(0,T;L2)2+(Re1/2∥u0,h∥02+Re3/2∥f∥L2(0,T;H-1)2)∥∇(u-ũh)∥L4(0,T;L2)2]+(Re-1+α)∥∇(u-ũh)∥L2(0,T;L2)2+α∥(I-Π)∇u∥L2(0,T;L2)2]}.Corollary 4.3.
From Theorem4.2, the approximation result (3.3) and Remark 3.1, one has
(4.14)∥(u-uh)(T)∥02+Rered-12∥∇(u-uh)∥L2(0,T;L2)2≤C(h4+αh2).
In particular,
(4.15)∥(u-uh)(T)∥0≤Ch2,ifα=h2.Remark 4.4.
A finite-element error estimate for theL2(Ω)-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 ∥(u-uh)(t)∥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(uhn-1,phn-1)∈(Xh,Qh), we find (uhn,phn)∈(Xh,Qh) satisfying
(5.1)(dtuhn,vh)+νa(uhn,vh)-d(vh,phn)+d(uhn,qh)+α((I-Π)∇uhn,(I-Π)∇vh)+b(uhn,uhn,vh)=(f,vh),
where dtuhn=(1/Δt)(uhn-uhn-1) and Δt is the time interval. Also, uhn is the finite-element approximate solutions of uh(tn=nΔt)(1<n<N). The fixed point iteration is adopted to solve (5.1) based on the Newtonian method [17] and the iteration tolerance is set to equal 10-8.For calculating the subgrid termα((I-Π)∇uhn,(I-Π)∇vh), the subgrid term is rewritten as follows:
(5.2)α((I-Π)∇uhn,(I-Π)∇vh)=α(∇uhn,∇vh)-α(Π∇uhn,∇vh).
In order to reduce the computational work, the following linear time extrapolation method is used to calculate Π∇uhn:
(5.3)Π∇uhn=2Π∇uhn-1-Π∇uhn-2+𝒪(δt2).


### 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 u=(u1,u2):
(5.4)u1(x,y,t)=-cos(nπx)sin(nπy)exp(-2n2π2νt),u2(x,y,t)=sin(nπx)cos(nπy)exp(-2n2π2νt),p(x,y,t)=-14(cos(2nπx)+cos(2nπy))exp(-2n2π2νt).
The viscosity is ν=0.001 and the corresponding Reynolds number is 1000. The mesh scales we choose are h=1/m (m=8,16,32,40,48,56,64,72). Let n=4 and α=mean(h)2, where mean(h) denotes the averaging mesh scale. The time step δt=0.001  and the final time t=2. In Table 1, the relative errors and absolute errors are given by different mesh scale h. In Figures 1 and 2, L2 and H1 convergence orders are given by two log-log plots. From these two figures, L2 and H1 convergence orders are equal to 2.31194 and 1.34385, respectively. The calculated convergence orders coincide with the theoretical analysis in Corollary 4.3. The expected convergence orders for the velocity in L2 and H1 are the second-order and the first-order, respectively. The computational convergence orders are a little higher than the expected orders. These results are mainly attributed to a specific test case which has a good regularity of the solutions.Table 1
L2 and H1 relative errors and absolute errors of numerical velocity uh.
h∥u-uh∥0/∥u∥0∥u-uh∥0|u-uh|1/|u|1|u-uh|10.1250.015750.011110.098390.617640.06250.00430.003040.060390.379140.031259.0178E-46.37147E-40.026020.163380.0255.03846E-43.55991E-40.018060.113380.020833.03229E-42.14246E-40.012810.080440.017861.95654E-41.38239E-40.009260.058140.015631.37093E-49.68631E-50.006810.042750.013891.06511E-47.52554E-50.005090.03198Figure 1
L2 Convergence rate of uh by a log-log plot.Figure 2
H1 Convergence rate of uh by a log-log plot.


### 5.2. Lid-Driven Cavity Flows

In this part, the comparisons among three kinds of subgrid models (Current model, Guermond-Marra-Quartapelle (GMQ) model [5], and Kaya-Layton-Riviere (KLR) model [6, 7]) are investigated for lid-driven cavity flows. These investigations will address the actions of subgrid models on large-scale flow structures. Now, we introduce the GMQ model and KLR model. The GMQ model is based on the two-level Lagrange finite-element setting, which is defined by [5]
(5.5)SGMQ(uhH,vhH)=cb∑K∈τk|K|∫K∇uhH·∇vhH,
where |K| is the measure of K, uhH=(1-PH)uh, and vhH=(1-PH)vh (see  [5] for the details of the functional settings and the definition of the operator PH). The KLR model is also based on two-level finite-element space, which is defined by [6, 7]
(5.6)SKLR(gH,vh)=α((∇uh,∇vh)-(gH,∇vh)),
where gH is the L2 projection of ∇uh(5.7)(gH-∇uh,lH)=0.
However, α is a user-selected stabilization parameter and typically α=𝒪(h) (see [6, 7] for the details of the functional settings). In the numerical computations, the iterative scheme is adopted to implement these two subgrid models [6].There exists a fundamental requirement for subgrid models, that is, effects from subgrid should tend to vanish for laminar flows under the suitable mesh resolutions. In order to implement the investigations, Reynolds numberRe=400 is adopted to implement lid-driven cavity flows. The computational domain Ω=[0,1]×[0,1] and the top boundary velocity (u,v)=(1,0), as well as the other three boundaries are nonslip boundary conditions. The mesh scales are 1/49. Under this Reynolds number, the time-dependent Navier-Stokes equations will converge to stationary laminar solutions. The convergence tolerance is set to equal 10-8. In Figure 3, the comparisons among four methods are given based on the profiles of horizontal and vertical velocity in mid planes of x and y axes. It is obvious that the three kinds of subgrid models give the almost same results compared with the Galerkin method. In Figure 4, we show the streamline patterns by the four different computational methods. From the streamline patterns, it is observed that the center position of the bottom-left vortex by the current subgrid model coincides with that by Galerkin method very well. However, GMQ and KLR models do not show the good coincidence. Numerically, the proposed subgrid model yields 0.3% deviation from the Galerkin solution for the center position of the bottom-left vortex, but GMQ and KLR models introduce about 2%~10% errors. It is known that the lid-driven flows of Re=400 is laminar and the “ideal" subgrid models should not affect flow structures. Under the sufficient mesh resolution, the two-grid subgrid models introduce the extra-action (dissipation) to flow fields according to this simple investigation, though three subgrid models can give the almost same results about the center positions of the center vortex and the bottom-right vortex. So, the other two subgrid models tend to be “unideal."Horizontal and vertical velocity profiles.
(a)(b)Streamline patterns.
(a)
Galerkin method(b)
Current model(c)
GMQ model(d)
KLR model


### 5.3. Fluid Flows around a Cylinder by High Reynolds Numbers

In this part, we investigate the two-dimensional under-resolved fluid flow around a cylinder. In this kind of fluid flows, the flow patterns are affected by the interaction of the fluid flows with two parallel planes and the surface of the cylinder. This problem is very useful to validate the subgrid model by vortex street patterns. The success and failure of the subgrid model simulations are useful for real fluid flows and engineering applications.The geometry of fluid flows is shown in Figure5. The height H and width W of the channel are equal to 1 and 6, respectively. The origin of the cylinder is at (x,y)=(1,0.5) and the radius R is equal to 0.15. The time-dependent inlet flow velocity profiles are given by
(5.8)u1(0,y,t)=-6y(y-H)H3,u2(0,y,t)=0,
and the boundary condition of the outlet is set as
(5.9)∂u∂n=0.
The boundary conditions of the two parallel planes and the cylinder surface are set as the nonslip boundary conditions. If the diameter D of the cylinder is chosen as the characteristic length and the mean velocity Umean of inlet as the characteristic velocity, the Reynolds number is defined by Re=UmeanD/ν. In order to implement the Galerkin method (finite-element direct numerical simulation (FEDNS) under the current computer capability and limited hardware resources, the Reynolds number is set to equal 1000. By this Reynolds number, we give the comparison results among referred four different methods. The time interval Δt=0.001. The mesh scale h=1/50 for Galerkin method and h=1/16 for the three kinds of subgrid models. The mesh scale of Galerkin method is enough to resolve the small scale of current flow structures (h∝Re-1/2) [2]. From Figures 6 and 7, it is very clear that the proposed subgrid model in this paper predicts the flow structures best, compared with the other two subgrid models based on the Galerkin benchmark solutions. From these results, it is proved that the GMQ and KLR subgrid models introduce a too strong local dissipation into the flow structures, but the current subgrid model presents a suitable local dissipation behavior.Figure 5
The schematic geometry of fluid flow past a cylinder.The contour patterns of pressure fields (t=1).
(a)
Galerkin method (GM)(b)
Current model(c)
GMQ model(d)
KLR model.The contour patterns of horizontal velocity fields (t=1).
(a)
Galerkin method (GM)(b)
Current model(c)
GMQ model(d)
KLR modelIn order to assess performance of the current subgrid model, the very high Reynolds numbers (Re=105 and Re=106) are chosen to implement the computations of complex flow phenomena. In Figures 8 and 9, the snapshots of vortex shedding are given. It is obvious that under these very high Reynolds numbers, the flow structures are complex and develop into the two-dimensional turbulence. The vortex filaments are clearly visible. These results demonstrate that the proposed subgrid model is effective and flexible.The vortex shedding patterns forRe=105.
(a)
t=6(b)
t=7The vortex shedding patterns forRe=106.
(a)
t=6(b)
t=7


### 5.4. Flow around an NACA0012 Airfoil at Zero Incidence

In this part, the proposed subgrid model will be used to simulate the flow around the NACA0012 airfoil at zero incidence [5]. Two different Reynolds numbers are used to simulate fluid flows: Re=105 and Re=106. The velocity of the incoming flow at infinity is the reference velocity scale and the chord of the airfoil is the reference length scale [5]. The flow domain is [-5,6]×[-5,5] and the upstream velocity is enforced on {x=-5;-5≤y≤5}∪{-5≤x≤6;y=±5}. Natural boundary conditions are set at the downstream boundary. The mesh scale on the surface of airfoil is h=1/120. The mesh includes 45894 triangle cells. The mesh partition is given in Figure 10. The time interval Δt=10-3 is employed for implementing computations of Re=105 and Re=106.Figure 10
Mesh near the edge of airfoil.The effects of the proposed subgrid term are shown in Figures11 and 12. The proposed subgrid model suppresses the spurious oscillations plaguing the Galerkin solution [5]. The proposed subgrid model possesses the same function as the other two referred subgrid models. For the sake of brevity, we only show the comparisons between the proposed subgrid model and literature results [5]. By this, it is demonstrated that the proposed model can remove the spurious oscillations reasonably.Contour figures forRe=105 at t=12.
(a)
Velocityu(b)
Velocityv(c)
Pressurep(d)
VorticityContour figures forRe=106 at t=12.
(a)
Velocityu(b)
Velocityv(c)
Pressurep(d)
Vorticity


## 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 u=(u1,u2):
(5.4)u1(x,y,t)=-cos(nπx)sin(nπy)exp(-2n2π2νt),u2(x,y,t)=sin(nπx)cos(nπy)exp(-2n2π2νt),p(x,y,t)=-14(cos(2nπx)+cos(2nπy))exp(-2n2π2νt).
The viscosity is ν=0.001 and the corresponding Reynolds number is 1000. The mesh scales we choose are h=1/m (m=8,16,32,40,48,56,64,72). Let n=4 and α=mean(h)2, where mean(h) denotes the averaging mesh scale. The time step δt=0.001  and the final time t=2. In Table 1, the relative errors and absolute errors are given by different mesh scale h. In Figures 1 and 2, L2 and H1 convergence orders are given by two log-log plots. From these two figures, L2 and H1 convergence orders are equal to 2.31194 and 1.34385, respectively. The calculated convergence orders coincide with the theoretical analysis in Corollary 4.3. The expected convergence orders for the velocity in L2 and H1 are the second-order and the first-order, respectively. The computational convergence orders are a little higher than the expected orders. These results are mainly attributed to a specific test case which has a good regularity of the solutions.Table 1
L2 and H1 relative errors and absolute errors of numerical velocity uh.
h∥u-uh∥0/∥u∥0∥u-uh∥0|u-uh|1/|u|1|u-uh|10.1250.015750.011110.098390.617640.06250.00430.003040.060390.379140.031259.0178E-46.37147E-40.026020.163380.0255.03846E-43.55991E-40.018060.113380.020833.03229E-42.14246E-40.012810.080440.017861.95654E-41.38239E-40.009260.058140.015631.37093E-49.68631E-50.006810.042750.013891.06511E-47.52554E-50.005090.03198Figure 1
L2 Convergence rate of uh by a log-log plot.Figure 2
H1 Convergence rate of uh by a log-log plot.


## 5.2. Lid-Driven Cavity Flows

In this part, the comparisons among three kinds of subgrid models (Current model, Guermond-Marra-Quartapelle (GMQ) model [5], and Kaya-Layton-Riviere (KLR) model [6, 7]) are investigated for lid-driven cavity flows. These investigations will address the actions of subgrid models on large-scale flow structures. Now, we introduce the GMQ model and KLR model. The GMQ model is based on the two-level Lagrange finite-element setting, which is defined by [5]
(5.5)SGMQ(uhH,vhH)=cb∑K∈τk|K|∫K∇uhH·∇vhH,
where |K| is the measure of K, uhH=(1-PH)uh, and vhH=(1-PH)vh (see  [5] for the details of the functional settings and the definition of the operator PH). The KLR model is also based on two-level finite-element space, which is defined by [6, 7]
(5.6)SKLR(gH,vh)=α((∇uh,∇vh)-(gH,∇vh)),
where gH is the L2 projection of ∇uh(5.7)(gH-∇uh,lH)=0.
However, α is a user-selected stabilization parameter and typically α=𝒪(h) (see [6, 7] for the details of the functional settings). In the numerical computations, the iterative scheme is adopted to implement these two subgrid models [6].There exists a fundamental requirement for subgrid models, that is, effects from subgrid should tend to vanish for laminar flows under the suitable mesh resolutions. In order to implement the investigations, Reynolds numberRe=400 is adopted to implement lid-driven cavity flows. The computational domain Ω=[0,1]×[0,1] and the top boundary velocity (u,v)=(1,0), as well as the other three boundaries are nonslip boundary conditions. The mesh scales are 1/49. Under this Reynolds number, the time-dependent Navier-Stokes equations will converge to stationary laminar solutions. The convergence tolerance is set to equal 10-8. In Figure 3, the comparisons among four methods are given based on the profiles of horizontal and vertical velocity in mid planes of x and y axes. It is obvious that the three kinds of subgrid models give the almost same results compared with the Galerkin method. In Figure 4, we show the streamline patterns by the four different computational methods. From the streamline patterns, it is observed that the center position of the bottom-left vortex by the current subgrid model coincides with that by Galerkin method very well. However, GMQ and KLR models do not show the good coincidence. Numerically, the proposed subgrid model yields 0.3% deviation from the Galerkin solution for the center position of the bottom-left vortex, but GMQ and KLR models introduce about 2%~10% errors. It is known that the lid-driven flows of Re=400 is laminar and the “ideal" subgrid models should not affect flow structures. Under the sufficient mesh resolution, the two-grid subgrid models introduce the extra-action (dissipation) to flow fields according to this simple investigation, though three subgrid models can give the almost same results about the center positions of the center vortex and the bottom-right vortex. So, the other two subgrid models tend to be “unideal."Horizontal and vertical velocity profiles.
(a)(b)Streamline patterns.
(a)
Galerkin method(b)
Current model(c)
GMQ model(d)
KLR model


## 5.3. Fluid Flows around a Cylinder by High Reynolds Numbers

In this part, we investigate the two-dimensional under-resolved fluid flow around a cylinder. In this kind of fluid flows, the flow patterns are affected by the interaction of the fluid flows with two parallel planes and the surface of the cylinder. This problem is very useful to validate the subgrid model by vortex street patterns. The success and failure of the subgrid model simulations are useful for real fluid flows and engineering applications.The geometry of fluid flows is shown in Figure5. The height H and width W of the channel are equal to 1 and 6, respectively. The origin of the cylinder is at (x,y)=(1,0.5) and the radius R is equal to 0.15. The time-dependent inlet flow velocity profiles are given by
(5.8)u1(0,y,t)=-6y(y-H)H3,u2(0,y,t)=0,
and the boundary condition of the outlet is set as
(5.9)∂u∂n=0.
The boundary conditions of the two parallel planes and the cylinder surface are set as the nonslip boundary conditions. If the diameter D of the cylinder is chosen as the characteristic length and the mean velocity Umean of inlet as the characteristic velocity, the Reynolds number is defined by Re=UmeanD/ν. In order to implement the Galerkin method (finite-element direct numerical simulation (FEDNS) under the current computer capability and limited hardware resources, the Reynolds number is set to equal 1000. By this Reynolds number, we give the comparison results among referred four different methods. The time interval Δt=0.001. The mesh scale h=1/50 for Galerkin method and h=1/16 for the three kinds of subgrid models. The mesh scale of Galerkin method is enough to resolve the small scale of current flow structures (h∝Re-1/2) [2]. From Figures 6 and 7, it is very clear that the proposed subgrid model in this paper predicts the flow structures best, compared with the other two subgrid models based on the Galerkin benchmark solutions. From these results, it is proved that the GMQ and KLR subgrid models introduce a too strong local dissipation into the flow structures, but the current subgrid model presents a suitable local dissipation behavior.Figure 5
The schematic geometry of fluid flow past a cylinder.The contour patterns of pressure fields (t=1).
(a)
Galerkin method (GM)(b)
Current model(c)
GMQ model(d)
KLR model.The contour patterns of horizontal velocity fields (t=1).
(a)
Galerkin method (GM)(b)
Current model(c)
GMQ model(d)
KLR modelIn order to assess performance of the current subgrid model, the very high Reynolds numbers (Re=105 and Re=106) are chosen to implement the computations of complex flow phenomena. In Figures 8 and 9, the snapshots of vortex shedding are given. It is obvious that under these very high Reynolds numbers, the flow structures are complex and develop into the two-dimensional turbulence. The vortex filaments are clearly visible. These results demonstrate that the proposed subgrid model is effective and flexible.The vortex shedding patterns forRe=105.
(a)
t=6(b)
t=7The vortex shedding patterns forRe=106.
(a)
t=6(b)
t=7


## 5.4. Flow around an NACA0012 Airfoil at Zero Incidence

In this part, the proposed subgrid model will be used to simulate the flow around the NACA0012 airfoil at zero incidence [5]. Two different Reynolds numbers are used to simulate fluid flows: Re=105 and Re=106. The velocity of the incoming flow at infinity is the reference velocity scale and the chord of the airfoil is the reference length scale [5]. The flow domain is [-5,6]×[-5,5] and the upstream velocity is enforced on {x=-5;-5≤y≤5}∪{-5≤x≤6;y=±5}. Natural boundary conditions are set at the downstream boundary. The mesh scale on the surface of airfoil is h=1/120. The mesh includes 45894 triangle cells. The mesh partition is given in Figure 10. The time interval Δt=10-3 is employed for implementing computations of Re=105 and Re=106.Figure 10
Mesh near the edge of airfoil.The effects of the proposed subgrid term are shown in Figures11 and 12. The proposed subgrid model suppresses the spurious oscillations plaguing the Galerkin solution [5]. The proposed subgrid model possesses the same function as the other two referred subgrid models. For the sake of brevity, we only show the comparisons between the proposed subgrid model and literature results [5]. By this, it is demonstrated that the proposed model can remove the spurious oscillations reasonably.Contour figures forRe=105 at t=12.
(a)
Velocityu(b)
Velocityv(c)
Pressurep(d)
VorticityContour figures forRe=106 at t=12.
(a)
Velocityu(b)
Velocityv(c)
Pressurep(d)
Vorticity


## 6. Conclusion

In this paper, a subgrid model is proposed for the time-dependent Navier-Stokes equations. The corresponding error estimates are given. The numerical investigations validate that the proposed model is effective and flexible. At the same time, the proposed model only acts on the small-scale fluid flows and does not affect the large-scale flow structures. Numerical investigations also demonstrate that the proposed subgrid model is superior to recent proposed subgrid models and is effective for the simulations of very high Reynolds number fluid flows. This model is established basedP2/P1 finite-element spaces and the corresponding scale separation is achieved by L2-projection. Computationally, the implementation of the proposed model is very easy for some legacy codes of fluid flows and does not need the background coarse mesh. In future, this subgrid method will be attempted to implement simulations for 3D high Reynolds turbulence flows.



---

*Source: 494829-2009-09-06.xml*