FULLY CONSISTENT AND ACCURATE BOUNDARY TREATMENT FOR THE GENERALIZED MAC SCHEME ON CURVILINEAR
DOMAINS (V06)
Abstract. In this paper, we introduce a new finite difference scheme to deal with the unsteady Navier-Stokes equations in velocity-pressure formulation. The scheme is based on the curvilinear coordinate and the MAC scheme. The main issue to be concerned is how to impose an appropriate numerical pressure boundary condition for the pressure Poisson equation to avoid the numerical boundary layer. We also provide stability and convergence analysis for our scheme that achieves second order of accuracy. Finally, we will demonstrate the numerical results when applying our scheme to the Lid-Driven Cavity Flow.
1. Introduction. In the study of the Navier-Stokes equations for incompressible fluid dynamics, the pressure has long been a main source of trouble for understanding and computing solutions to these equations. Pressure playes a role like a Lagrange multiplier to enforce the incompressibility constraint. This leads to computational difficulties, typically related to the lack of an evolution equation for updating the pressure dynamically, and a lack of useful boundary conditions for determining the pressure by solving boundary-value problems.
The MAC scheme of Harlow and Welch in 1965 is a direct discretization of the Navier-Stokes equations using second order finite differences implimented on a stag-gered grid. It results in a discrete Poisson equation for the pressure by applying the incompressibility constraint to the discretized momentum equation. However, use of a staggered grid generally limits the application oh the MAC scheme to simple geometries.
In the late 60s Chorin and Temam introduced the projection method. Projection methods are time splitting schemes in which an intermediate velocity is first com-puted and then projected onto the space of incompressible vector fields by solving a Poisson equation for pressure. Unfortunately, the time splitting introduces numerical boundary layers in the pressure and intermediate velocity field.
In the finite element approach, the fast Stokes solver schemes for fully coupled systems have proven highly successful. To ensure stability, the inf-sup compatibility condition (or Babuska-Brezzi condition) between the pressure and velocity finite el-ement spaces has to be satisfied. This generally complicates the implel-ementation of fast solver for the resulting linear systems.
2. Incompressible Navier-Stokes Equation. Let Ω be a bounded domain in R2 or R3 with boundary ∂Ω = Γ. The incompressible Navier-Stokes Equation with
primitive variables takes the form
ut+ (u · ∇)u + ∇p = ν△u momentum
∇ · u = 0 incompressibility u|Γ = 0 no-slip B.C.
where u is the velocity, p the pressure, and ν the kinematic viscosity. Define
ω= ∇ × u (vorticity) and use the vector identity (Bernoulli’s Law)
(u · ∇)u = ω × u + ∇(|u|2/2),
the incompressible NSE are then given by ut+ ω × u + ∇˜p = −ν∇ × ω in Ω ∇ · u = 0 in Ω u|Γ = 0 on Γ
where ˜p = p + |u|2/2 is the total pressure. But p will stand for ˜p in the rest of this
paper.
2.1. Pressure Poisson Equation. Find p in the Navier-Stokes Equation
ut+ ω × u + ∇p = ν∆u in Ω
u = 0 on Γ
so that ∇ · u = 0 in Ω.
Since p is not described by an evolutionary equation, it needs to be solved via static relation utilizing the incompressibility constraint
∇2p = −∇ · (ω × u) in Ω
From standard PDE theory, one needs to supply a boundary condition for p, e.g. ∂p
∂n|Γ = ν∇
2
u· n or ∂p
∂n|Γ= −ν(∇ × ∇ × u) · n
The main question is how to impose an appropriate boundary condition for p in discrete setting?
3. Classical MAC Scheme and Hodge decomposition. We recall the Mark and Cell (MAC) scheme (Welch and Harlow, 1965) on staggered grids in this section. We consider the case of 2D NSE. The components of u = (u, v) and p on staggered grids are placed as
3.1. The Spatial Discretization of MAC. MAC uses the following centered difference operators: 1. Discrete gradient Dxpi,j+1 2 = pi+1 2,j+1 2− pi− 1 2,j+ 1 2 h Dypi+1 2,j = pi+1 2,j+1 2− pi+ 1 2,j− 1 2 h 2. Discrete Divergence ∇h· u = Dxu + Dyv 3. Discrete Laplacian ∇2hp = Dx2p + D2yp = ∇h· ∇hp 4. Discrete Curl ∇h× u = Dxv − Dyu
3.2. Time discretization of MAC. Given (un,pn) with ∇ h· un = 0, find (un+1,pn+1) so that ∇ h· un+1= 0 and un+1− un ∆t + (u n· ∇ h)un+ ∇hpn+1 = ν∇2hun in Ω un+1| Γ = 0 on Γ
We decompose the MAC scheme into 2 steps: 1. Solve for u∗ from
u∗− un ∆t + (u n· ∇ h)un= ν∇2hun. 2. Find un+1 with ∇ h· un+1= 0 and un+1|Γ = 0 from un+1− u∗ ∆t + ∇hp n+1= 0.
The above is done by first solving the pressure Poisson equation (PPE) ∇2hp = ∇h· u
∗
∆t , then substitute back to get
un+1 = u∗− ∆t∇ hpn+1.
This is the Hodge decomposition of u∗:
u∗= un+1+ ∆t∇ hpn+1.
Moreover, the divergence free part and the gradient part are perpendicular to one another.
3.3. Consistent and stable boundary treatment. In NSE, p is treated as a Lagrangian multiplier purely to enforce the ∇h· u = 0 constraint.
Proposition 1. For standard MAC on the staggered grid, the result of un+1= u∗− ∆t∇
hpn+1 in Ω
is independent of the choice/value of u∗· n on Γ in solving the Neumann
boundary-value problem ∇2hp = ∇h· u ∗ ∆t in Ω, n· ∇p = u∗· n ∆t on Γ.
This explains why n · ∇p = 0, formally inconsistent, still gives 2nd order accuracy for MAC scheme on the staggered grid.
The advantages of staggered grids are
• the discretization operator defined naturally on the staggered grids so that the usual vector identities still valid in the discretization. These vector identities are key important in NSE, Maxwell system, etc.
• exact Hodge decomposition • divergence free can be achieved
ξ = 1 constant ξ = 2 constant η =2 constant 1 η =constant ξ = 1 constant ξ = 2 constant η =2 constant 1 η =constant
Fig. 4.1: The computational coordinate
• avoid the pressure boundary condition
If non-staggered grid is used, then one may need consistent pressure boundary condi-tion, otherwise there will be numerical boundary layer. It is not trivial.
The main feature of the MAC is that the discretization operator defined naturally on the staggered grid so that the usual vector identities related to ∇h, ∇h·, ∇h× still
valid in the discretization. These vector identities are key important in many of physical system such as NSE, Maxwell etc.
4. Generalized MAC on Curvilinear Domains. In this paper, we mainly aim to solve the following 2D NSE on a curvilinear domain.
(4.1) ut+ ωu⊥+ ∇p = ν∇⊥ω + f in Ω, ∇ · u = 0 in Ω, u = 0 on ∂Ω.
Here u is the velocity field, f is force, p is pressure, and the vorticity field ω is the z-component of ∇ × u.
4.1. curvilinear coordinate. First we introduce the curvilinear coordinate sys-tem. Every point x in the physical domain with Euclidean coordinate (x, y) has a curvilinear coordinate (ξ1, ξ2) refer to the computational domain, i.e., ξ1 = ξ1(x, y)
and ξ2= ξ2(x, y).
After a coordinate system has been determined, we place an uniform mesh on the computational domain with mesh size ∆ξ1 = h
1 and ∆ξ2= h2. Then a skewed
coordinate system will be introduced by using two straight lines L1 and L2
L1: ξ2= h2
h1
ξ1 and L2: ξ2= −h2
h1
ξ1
to be the new coordinate axis in the computation plan, as illustrated in the left of Figure 4.1. Let
ϕ1(ξ1, ξ2) = h1ξ2− h2ξ1 and ϕ2(ξ1, ξ2) = h1ξ2+ h2ξ1.
Then, the skewed variables (η1, η2) are defined as η1:= ∇(ξ1,ξ2)ϕ2 k∇(ξ1,ξ2)ϕ2k · (ξ1, ξ2) = (h2, h1) p h2 1+ h22 · (ξ1, ξ2) (4.2)
and η2:= ∇(ξ 1 ,ξ2 )ϕ1 k∇(ξ1,ξ2)ϕ1k · (ξ1, ξ2) =p(−h2, h1) h2 1+ h22 · (ξ1, ξ2). (4.3)
Hence, the step size along the η1and η2direction is
∆η1:= η1(h
1, h2) − η1(0, 0) and ∆η2:= η2(−h1, h2) − η2(0, 0).
And we can define
h := ∆η1= ∆η2=p2h1h2 h2
1+ h22
.
The metric tensor with respect to the new skewed variables (η1, η2) is then given
by (4.4) bgµν:= ∂x ∂ηµ · ∂x ∂ην and (4.5) bg := det(bgµν).
Then bgµγ is defined by the inverse of bg
γν, that is
(4.6) bgµγbgγν= δµν.
An obvious conclusion is that
(4.7) pbg = det(∂x
∂η) and
(4.8) bgµν = ∇ηµ· ∇ην.
5. Approximation of Metric Tensors. Let (bgij)h denote the second order
approximation of the metric tensor corresponded to the skewed curvilinear coordinate (η1, η2). Then (bg
i,j)h in cells can be easily computed by the difference quotients of
the position vectors xij’s. That is
((bg11)h)i+1 2,j+ 1 2 := xi+1,j+1− xi,j h · xi+1,j+1− xi,j h , ((bg22)h)i+1 2,j+ 1 2 := xi,j+1− xi+1,j h · xi,j+1− xi+1,j h , ((bg12)h)i+1 2,j+1 2 = ((bg21)h)i+ 1 2,j+ 1 2 := xi+1,j+1− xi,j h · xi,j+1− xi+1,j h .
It should be noted that these quantities are all evaluated in cells. At interior grids, we define ((bgαβ)h)i,j := 1 4 n ((bgαβ)h)i+1 2,j+ 1 2 + ((bgαβ)h)i+ 1 2,j− 1 2 + ((bgαβ)h)i−1 2,j+ 1 2 + ((bgαβ)h)i− 1 2,j− 1 2 o . 5
At boundary grids, say ξ0,j, we define ((bgαβ)h)0,j := 1 2 ((bgαβ)h)1 2,j+ 1 2 + ((bgαβ)h) 1 2,j− 1 2 . At corner grids, say ξ0,0, we define
((bgαβ)h)0,0 := ((bgαβ)h)1 2,
1 2.
By Equation (4.5), the numerical analogous formula is b
gh:= det ((bgh)µν) .
According to Equation (4.6), the numerical quantities bgµνh ’s can be obtained from
solving the linear system
b
ghµλ(bgλν)h:= δµν.
6. Vector Expression and Differential Operator in Curvilinear Coordi-nate. From now on, all componentwise expression are all referred to the curvilinear (η1, η2)-coordinate. Let (u1, u2) stand for the contra-variant component of a vector
field u, that is
u= u1∂x ∂η1 + u
2∂x
∂η2.
And let (u1, u2) be the co-variant component of a vector field u, that is
u= u1∇η1+ u2∇η2.
The transformation between co-variant component and contra-variant component is u1= u1bg11+ u2gb21 and u2= u1bg12+ u2bg22,
or
u1= u1bg11+ u2bg21 and u2= u1bg12+ u2bg22.
The contra-variant expression of u⊥ is
u⊥ =pgb− u2bg11+ u1bg21 ∂x ∂η1 + p b g− u2gb12+ u1bg22 ∂x ∂η2.
The gradient of the pressure field p can be written in its co-variant component as ∇p = ∂p
∂η1∇η 1+ ∂p
∂η2∇η 2,
or its contra-variant component as ∇p = ∂p ∂η1bg 11+ ∂p ∂η2bg 21 ∂x ∂η1 + ∂p ∂η1bg 12+ ∂p ∂η2bg 22 ∂x ∂η2. And ∇⊥ω = 1 p b g −∂ω ∂η2 ∂x ∂η1 + ∂ω ∂η1 ∂x ∂η2 ,
∇ · u = p1 b g ∂ ∂η1( p b gu1) + ∂ ∂η2( p b gu2) . Therefore, the Laplacian of p is
∇2p = ∇ · ∇p = p1 b g 2 X µ,ν=1 ∂ ∂ηµ( p b gbgµν∂η∂pν).
Finally, the vorticity field ω is the z-component of ∇ × u, hence, ω = ∇⊥· u = 1 p b g ∂u 2 ∂η1 − ∂u1 ∂η2 .
7. Space Discretization. We define the skewed difference operators in the fol-lowing way:
Definition 1. If a is a field defined at grids, then ( bD1a)i−1 2,j− 1 2 = ai,j− ai−1,j−1 h and ( bD2a)i−12,j− 1 2 = ai−1,j− ai,j−1 h .
If b is a field defined in cells, then ( bD1b)i,j= bi+1 2,j+ 1 2 − bi− 1 2,j− 1 2 h and ( bD2b)i,j = bi−1 2,j+ 1 2 − bi+ 1 2,j− 1 2 h
at interior grids. But, at boundary or corner grids, they only involve those terms in the interior cells. For example,
( bD1b)0,j = b1 2,j+1 2 h and ( bD2b)0,j = −b1 2,j− 1 2 h at ξ0,j, ( bD1b)0,0 = b1 2,1 2 h and ( bD2b)0,0 = 0 at ξ0,0.
7.1. Discrete Operator: ∇h, ∇h·, ∇2h, ∇⊥h and ∇⊥h· . Let G be the space of
all scalar fields defined at grids and C the space of all vector fields defined in cells. We introduce the following discrete operators:
∇h and ∇⊥h : G → C,
∇h· and ∇⊥h· : C → G,
and
∇2
h: G → G.
7.1.1. Discrete Gradient. Definition 2. If a is a scalar field defined at grids, then the discrete gradient of a , in its co-variant componentwise expression, is
(∇ha)i−1 2,j− 1 2 = ( bD1a, bD2a)i− 1 2,j− 1 2
at each cell center. And its contra-variant componentwise expression is (∇ha)i−1 2,j−1 2 = ( bD1a bg 11 h + bD2a bgh21, bD1a bg12h + bD2a bg22h )i−1 2,j− 1 2.
But the discrete gradient-perp of a, in its contra-variant componentwise expression, is (∇⊥ ha)i−1 2,j− 1 2 = 1 p b gh ( − bD2a , bD1a )i−1 2,j− 1 2 7
7.1.2. Discrete Divergence. Definition 3. If b is a vector field defined in cells and (b1, b2) is its contra-variant expression, then the discrete divergence of b is
a scalar fields defined at grids in the following way: 1. At interior grids: (∇h· b)i,j= 1 p b gh b D1( p b ghb1) + bD2( p b ghb2) i,j,
2. At boundary grids, say ξ0,j:
(∇h·b)0,j = 2 (pbgh)0,j (pbghb1)1 2,j+ 1 2 − ( p b ghb2)1 2,j− 1 2 h ! , (0 < j < Ny),
3. At corner grids, say ξ0,0:
(∇h· b)0,0 = 4 (pbgh)0,0 (pbghb1)1 2, 1 2 h . Fig. 7.1:
7.1.3. Discrete Divergence-perp. Definition 4. Let b be a vector field that vanishes on the boundary, and (b1, b2) its co-variant component. Then the discrete
divergence-perp of b is defined in the following way: 1. At interior grids: (∇⊥ h · b)i,j= 1 p b gh ( bD1b2− bD2b1)i,j
2. At boundary grids, say ξ0,j:
(∇⊥h · b)0,j = 2 (pgbh)0,j (b2)1 2,j+ 1 2 + (b1) 1 2,j− 1 2 h ! , (0 < j < Ny),
3. At corner grids, say ξ0,0:
(∇⊥ h · b)0,0 = 4 (pbgh)0,0 (b2)1 2, 1 2 h .
Fig. 7.2:
7.1.4. Discrete Laplacian. Definition 5. If a is a scalar fields defined at grids, then the numerical Laplacian of a is defined at grids in the following way:
1. At interior grids: (∇2ha)i,j= 1 p b ghi,j 2 X µ,ν=1 b Dµ( p b ghbghµνDbνa)i,j
2. At boundary grids, say ξ0,j: for 0 < j < Ny,
(∇2ha)0,j =p2 b gh0,j ( (pbghbg11h Db1a + p b ghbgh21Db2a)12,j+ 1 2 h − ( p b ghbg12h Db1a + p b ghbgh22Db2a)1 2,j− 1 2 h ) , 3. At corner grids, say ξ0,0:
(∇2ha)0,0= 4 p b gh0,0 p b ghbgh11Db1a + p b ghbg21h Db2a 1 2, 1 2 h .
Lemma 7.1. If a is a scalar fields defined at grids, and vanishes (or is constant ) at boundary, then
∇h· ∇⊥ha = 0
at each grid (include boundary and corner).
Lemma 7.2. Let b be a scalar field defined at grids. At interior and boundary grids, we have
∇h· ∇hb = ∇2hb.
Lemma 7.3. Let b be a scalar field defined at grids. At interior and boundary grids, we have
∇⊥
h · ∇⊥hb = ∇2hb.
The proofs is in Appendix.
The advantages of these Difference Operators are
• No-slip boundary condition is incorporated into operators.
• No boundary condition is needed for the pressure Poisson Equation.
• These difference operators are all consistent in the sense that local truncation error is O(h2) in interior and O(h) on boundary.
7.2. IMAC Scheme. We place the contra-variant components of the velocity field u, (u1, u2), in cell centers. The vorticity field ω and pressure field p are placed
at grids. Explicitly speaking, u1 i−1 2,j− 1 2 and u2 i−1 2,j− 1 2 are at ξi−1 2,j− 1
2, but ωi,j and
pi,j are at ξi,j.
Because ω is the z-component of ∇ × u, we have the following definition. Definition 6. We define the numerical vorticity field eω at grids via the co-variant components (eu1, eu2) of numerical velocity euby
e
ω = (∇h)⊥· eu.
Consequently, our numerical scheme is to discretize the NSE at the cell center ξi−1 2,j− 1 2 in the form (7.1) e ut+ eω eu⊥+ ∇hp = ν∇e ⊥hω + fe in cells, e
ω = (∇h)⊥· eu at grids (including boundary),
∇h· eu = 0 at grids (including boundary),
e
u = 0 at boundary grids.
where ep is the numerical pressure field and e ωi−1 2,j− 1 2 = 1
4(eωi,j+ eωi−1,j+ eωi,j−1+ eωi−1,j−1).
7.3. Nonlinear Stability. Let G be the space of all grid scalar fields and C the space of all cell vector fields.
W.L.O.G., we assume the computation domain is { ξ1 ≥ 0 , ξ2≥ 0 }, and hence
the indices are
0 ≤ i < ∞, and 0 ≤ j < ∞. We also assume f is zero.
Define the discrete inner product of two vector fields in the space C as
h a , b ih,C= h1h2 ∞ X i=0 ∞ X j=0 n a1b1+ a2b2 pgbh o i+1 2,j+1 2 = h1h2 ∞ X i=0 ∞ X j=0 n a1b 1+ a2b2 pgbh o i+1 2,j+1 2 . (7.2)
The induced discrete norm in the space C is
kak2h,C= h a , a ih,C.
Define the discrete inner product of two scalar fields in the space G as
(7.3) h a , b ih,G= h1h2 ∞ X i=1 ∞ X j=1 a bpbgh i,j+ 1 2 ∞ X i=1 a bpbgh i,0 + 1 2 ∞ X j=1 a bpbgh 0,j+ 1 4 a bpbgh 0,0 .
The induced discrete norm in the space G is
kak2h,G= h a , a ih,G.
Lemma 7.4. If a is a vector filed defined in cells and b, c are scalar fields defined at grids. Then we have
1. h a , ∇⊥
hb ih,C= −h ∇⊥h · a , b ih,G
2. h a , ∇hb ih,C = −h ∇h· a , b ih,G
3. h ∇hb , ∇hc ih,C = −h ∇2hb , c ih,G
The proof is in Appendix.
Lemma 7.5. For the numerical solution eu and eω to (7.1) with f = 0, we have the stability 1 2 d dtkeuk 2 h,C+ νkeωk2h,G= 0.
Proof. By Lemma 7.4, we have
h eu, ∇⊥hω ie h,C = −h ∇h⊥· eu, eω ih,G= −keωk2h,G
Taking the discrete inner product with eu, we have h eu, eut+ ¯eωeu ⊥ + ∇hp ie h,C= h eu, ν∇⊥hω ie h,C. By Lemma 7.4, we have h eu, ∇hp ie h,C = −h ∇h· eu, ep ih,G= 0 Noting that h eu, ¯eωeu⊥ih,C = 0,
the stability is implied.
8. Time splitting. For the time splitting, we explain the idea by using forward Euler in this section. But in the practice, Fourth order Runge-Kutta time stepping is used. We consider the following time splitting scheme: Given u0, find eun+1by
1. Direct evaluate u∗ in cells by
u∗− eun ∆t + (eωeu
⊥
)n= ν∇⊥ hωen
2. solving epn+1 with mean zero from
∇2hpen+1= 1 ∆t∇h· u
∗ at all grids (including boundary).
REMARK: NO need for pressure boundary condition!! 3. Find eun+1 from e un+1− u∗ ∆t + ∇hpe n+1= 0 such that ∇h· eun+1= 0. 4. Update eωn+1= ∇⊥ h · eu n+1 at grids. 11
9. Error Analysis. Define k · kCk( ¯Ω)= k · kCk.
Theorem 1. Let Ω = [0, 1] × (0, 2π) ⊂ R2 with periodic boundary condition in y. Assume u ∈ C1([0, T ], C4( ¯Ω)), ω and p ∈ C0([0, T ], C3( ¯Ω)) are the exact solution
of Equation (4.1) with f = 0. Let eu, eω and ep be the numerical solution of Equation (7.1). Then we have sup 0≤t≤Tke u− uk2 h,C+ ν Z T 0 ke ω − ωk2 h,G ≤ C1h4, (9.1) Z T 0 k∇h(ep − p)k2h,C ≤ C2h2. (9.2) where C1= C1(T, ν, kukC2, k∂tuk C2, kωk C3, kpk C3). Proof.
Step 1. Strange’s expansion and truncation error analysis
Let ψ be the stream function, that is, ∇⊥ψ = u. ψ can be solved by the Poisson
equation
∇2ψ = ∇⊥· u with homogeneous Neumann boundary condition.
We recall that, for a function ϕ ∈ C3( ¯Ω), we have
∇2 hϕ|x=0,1= ±2 h1 ∂ϕ ∂n+ ∇ 2ϕ + R(ϕ, 3)h,
where R(ϕ, k) is a/some linear combination of the derivatives of ϕ up to kth order. We now construct an auxiliary function φ ∈ C1([0, T ], C3( ¯Ω)) that satisfies
1. φ vanishes at the boundary of Ω. 2. ∂φ ∂n|x=0,1= ∓h1 2 R(ψ, 3). Define (9.3) ( U = ∇⊥ h(ψ + h2φ), in cells, W = ∇⊥
h · U = ∇2h(ψ + h2φ), at grids up to the boundary of Ω.
Then W = ω + hR(ψ, 3) + ±2 h1 ∂φ ∂n + {R(ψ, 4) + R(φ, 2)}h2 at x = 0, 1, ω + {R(ψ, 4) + R(φ, 2)}h2 at interior grids.
By the construction of φ, we have 1. |U − u| ≤ C3(kψkC3, kφkC1)h2.
2. |W − ω| ≤ C4(kψkC4, kφk C2)h2.
3. ∇h· eu= 0 and ∇h· U = 0.
By the local truncation error analysis, we have
∂tU+ ωU⊥+ ∇hp = ν∇⊥hω + E in cells,
where
Since
ωU⊥− ωu⊥ = ω(U − u)⊥+ (ω − ω)u⊥= O(h2),
we arrived
|E| ≤ C5h2,
where C5= C5(kωkC3, kpkC3, kukC2, k∂tukC2).
Step 2. Error equation and energy estimate Define e = eu− U , we have ∇h· e = 0, and
∂te+ (eωeu ⊥
− ωU⊥) + ∇h(ep − p) = ν∇⊥h(eω − ω) − E.
By Lemma 7.4 and the fact ∇h· e = 0,
h e , ∂teih,C+ h e , eωeu⊥− ωU⊥ih,C = νh e , ∇⊥h(eω − ω) ih,C− h e , E ih,C.
The viscous term, by Lemma 7.4, we have h e , ∇⊥
h(eω − ω) ih,C = h eu− U , ∇⊥h(eω − ω) ih,C
= − h eω − W , eω − ω ih,G
= − h eω − ω , eω − ω ih,G− h ω − W , eω − ω ih,G
The non-linear term is
h e , eωeu⊥− ωU⊥ih,C= h e , eω(eu− U )⊥+ (eω − ω)U⊥ih,C
= h e , (eω − ω)U⊥i h,C.
Here we use the fact that e = eu− U is perpendicular to eω(eu− U )⊥. Hence
h e , eωeu⊥ − ωU⊥i h,C ≤1 νkek 2 h,C(max |U ⊥|2) +ν 4keω − ωk 2 h,C ≤C6 ν kek 2 h,C+ ν 4keω − ωk 2 h,G. Here we use keω − ωk2 h,C≤ keω − ωk2h,G and C6= max |U ⊥ |2. Then, 1 2∂tkek 2 h,C+ νkeω − ωk2h,G = − νh ω − W , eω − ω ih,G− h e , (eω − ω)U⊥ih,C− h e , E ih,C ≤ν 4keω − ωk 2 h,G+ νkω − W k 2 h,G+ C6 ν kek 2 h,C+ ν 4keω − ωk 2 h,G+ kek 2 h,C+ 1 4kEk 2 h,C ≤ (1 + C6 ν )kek 2 h,C+ ν 2keω − ωk 2 h,G+ (νC42+ C2 5 4 )h 4. Therefore, we obtain 1 2∂tkek 2 h,C+ ν 2keω − ωk 2 h,G≤ (1 + C6 ν )kek 2 h,C+ (νC42+ C2 5 4 )h 4.
By Gronwall’s inequality, we have sup 0≤t≤T kek2h,C+ ν Z T 0 ke ω − ωk2h,G≤ C1h4. 13
Step 3. Error estimate of pressure
Take inner product with ∇hǫ on the both sides. After using Lemma 7.4 and
∇h· (U − eu) = 0, we obtain h ∇hǫ , ωU⊥− eωeu ⊥ ih,C+ k∇hǫk2h,C= νh ∇hǫ , ∇⊥h(ω − eω) ih,C+ h ∇hǫ , E ih,C. Then (9.4) k∇hǫk2h,C≤ 3 4k∇hǫk 2 h,C+ ν2k∇ ⊥ h(ω − eω)k2h,C+ kEk2h,C+ kωU ⊥ − eωeu⊥k2h,C. Note that k∇⊥h(ω − eω)k2h,C≤ 16 h2kω − eωk 2 h,G and kωU⊥ − eωeu⊥k2 h,C≤ kωk2C0kU − euk2h,C+ (keuk2L∞)k¯ω − ¯eωk 2 h,C. Therefore 1 4 Z T 0 k∇hǫk2h,C≤ 16ν2 h2 Z T 0 kω − e ωk2h,G+ Z T 0 kEk2h,C+ sup [0,T ] kωk2C0 ! Z T 0 kU − e uk2h,C + sup [0,T ]ke uk2L∞ ! Z T 0 kω − e ωk2h,G (9.5) keuk2L∞ ≤ 1 h2keuk 2 h,C ≤ 2 h2(keu− U k 2 h,C+ kU k 2 h,C) ≤ C7 h2,
where C7 depends on C1, C6 and |Ω|.
(9.6)
Z T 0
k∇hǫk2h,C ≤ C2h2
where C2 is a constant that depends on C1.
10. Numerical Result.
10.1. Convergence Rate Test. In this subsection, we report the accuracy check results. We use a stream functionsψ(x, y) to determine the velocity field u(x, y) and the vorticity ω by
(10.1) u= ∇⊥ψ and ω = ∇2ψ
as exact solutions of (4.1). The pressure field is given by
(10.2) p = cos(t) cos(x) sin(y).
The force term f (time dependent) is chosen to ensure that u, ω and p are solutions of (4.1).
The physical domain is the region bounded by the following four boundary seg-ments B1 : x = 0, B2 : y = 0, B3 : x = 1 and B4 : y = 1 + α sin(2πx). When the
parameter α 6= 0, B4 is a curvilinear boundary segment.
The stream function: ψ(x, y) = cos(t) (x − x2)2{y + α sin(2πx)y − y2}2. Re = 10000; Final time = 3.2; α = 0.2;
−1 −0.9 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 −1 −0.8 −0.6 −0.4 −0.2 0
The grid distribution on the domain with a curvilinear bottom
Fig. 10.1: Grid distribution on a curvilinear domain
Table 10.1: Absolute errors and order of accuracy of the numerical velocity eu for Case 3.
Nx× Ny ku − euk∞ Order ku − eukL1 Order
100 × 100 2.336E − 5 5.656E − 6
200 × 200 7.269E − 6 1.684 1.425E − 6 1.989 400 × 400 2.118E − 6 1.779 3.569E − 7 1.998
Table 10.2: Absolute errors and order of accuracy of the numerical vorticity eω and numerical pressure field ep for Case 3.
Nx× Ny kω − eωk∞ Order kp − epk∞ Order
100 × 100 5.543E − 3 2.920E − 5
200 × 200 1.349E − 3 2.04 7.346E − 6 1.99 400 × 400 3.325E − 4 2.02 1.836E − 6 2.00
10.2. Lid-Driven Cavity Flow.
10.2.1. On the Square Domain – Benchmark Problem. Mesh size 256 × 256, ∆t = 3.8E − 3.
Reynolds number = 1000
Final time of numerical simulation T = 1.5/Re ≃ 48.0 Real time cost = 46.0 min.
The vorticity contour agrees with the benchmark result of [1].
0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 x 10−5
Error of pressure, Re=10000, T=3.2
Fig. 10.2: Error of pressure in Example 1, Re=10000
10.2.2. On the Square Domain – High Reynolds Number. 1. Reynolds number = 8000
Final Time T = 135
The domain = [−1, 0] × [−1, 0].
10.2.3. On the Domain with a Curvilinear Bottom. The physical domain in this numerical simulation is the region bounded by the four boundary segments x = 0, y = 0, x = −1 and y = −1 + 0.2 sin(2πx).
Final Time T = 80.0. Reynolds number = 2000.
Lid−Driven Cavity Flow, numerical vorticity contour, 256x256, Re=1000, time=48.00 −3 −4 −5 −2 −1 −0.5 0 1 2 3 1 0 0 2 3 3
Fig. 10.3: Numerical vorticity contour of Lid-Driven Cavity Flow at Re=1000
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Driven Cavity Flow, contour of numerical ω, 256x256, Re=8000, time=135
−1 −2.23 2.8 2 0 2.8 0 0 2.8
Driven Cavity Flow, contour of numerical ω, 256x256, Re=2000, time=80.00 −2 −1 −0.5 0 0.5 1 2 3 3 0 0 0 0.5 1 2 3 −3 −4 −5 −0.5 −1 −0.9 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 −1 −0.8 −0.6 −0.4 −0.2 0
Fig. 10.5: Numerical vorticity contour of Lid-Driven Cavity Flow on the domain with a curvilinear bottom at Re=2000
Appendix A. Some Proofs.
Proof. [Proof of Lemma 7.2] At interior grids, the result is straightforward by the definition. We exam the result at boundary points, say ξ0,j. Note that the
contra-variant componentwise expression of ∇hb is
(∇hb)i−1 2,j− 1 2 = ( bD1b bg 11 h + bD2b bgh21, bD1b bgh12+ bD2b bgh22)i−1 2,j− 1 2.
Then, according to the definition of the discrete divergence at ξ0,j, we have
(∇h· ∇hb)0,j =p2 b gh0,j p b ghbg11h Db1b + p b ghbg21h Db2b 1 2,j+ 1 2 −pbghbgh12Db1b + p b ghbgh22Db2b 1 2,j− 1 2 h = (∇2hb)0,j.
According to the definition of the discrete divergence at ξ0,0, we have
(∇h· ∇hb)0,0= 4 p b gh0,0 p b ghbgh11Db1b + p b ghbgh21Db2b 1 2, 1 2 h = (∇ 2 hb)0,0.
Proof. [Proof of Lemma 7.3] At interior grids, ∇⊥ hb = − bD2b p b gh ∂x ∂η1 + b D1b p b gh ∂x ∂η2 = − bpD2b b gh (bg11)h+ b D1b p b gh (bg21)h ! ∇η1+ − bpD2b b gh (bg12)h+ b D1b p b gh (bg22)h ! ∇η2 =−pbghgb22h Db2b − p b ghgb12h Db1b ∇η1+pbghbgh12Db2b + p b ghbgh11Db1b ∇η2. Therefore, ∇⊥ h · ∇ ⊥ hb = 1 p b gh n b D1 p b ghbg12h Db2b + p b ghgb11h Db1b + bD2 p b ghbg22h Db2b + p b ghbg12h Db1b o = ∇2hb.
At boundary points, say ξ0,j,
(∇⊥ h · ∇ ⊥ hb)0,j =p2 b gh0,j p b ghbg12h Db2b + p b ghbg11h Db1b 1 2,j+ 1 2 −pbghbgh22Db2b + p b ghbgh12Db1b 1 2,j− 1 2 h = (∇2hb)0,j.
In the case of corner grids,
(∇⊥ h · ∇ ⊥ hb)0,0= 4 p b gh0,0 p b ghbg12h Db2b + p b ghbg11h Db1b 1 2, 1 2 h = (∇ 2 hb)0,0.
Proof. [Proof of Lemma 7.4] 1. h a , ∇⊥ hb ih,C = h1h2 ∞ X i=0 ∞ X j=0 −a1Db2b + a2Db1b i+1 2,j+ 1 2 = h1h2 ∞ X i=0 ∞ X j=0 −(a1)i+1 2,j+ 1 2 bi,j+1− bi+1,j h + (a2)i+12,j+ 1 2 bi+1,j+1− bi,j h = h1h2 ∞ X i=0 ∞ X j=1 (−a1)i+1 2,j− 1 2 bi,j h + ∞ X i=1 ∞ X j=0 (a1)i−1 2,j+1 2 bi,j h + ∞ X i=1 ∞ X j=1 (a2)i−1 2,j−1 2 bi,j h − ∞ X i=0 ∞ X j=0 (a2)i+1 2,j+ 1 2 bi,j h = − h1h2 ∞ X i=1 ∞ X j=1 ( bD1a2− bD2a1)bi,j+ ∞ X i=1 (a2)i+1 2, 1 2 − (a1)i− 1 2, 1 2 h bi,0 + ∞ X j=1 (a2)1 2,j+ 1 2 + (a1) 1 2,j−1 2 h b0,j+ (a2)1 2, 1 2 h b0,0 = − h1h2 ∞ X i=1 ∞ X j=1 (∇⊥h · a p b ghb)i,j+ ∞ X i=1 (∇⊥ h · a p b ghb)i,0 2 + ∞ X j=1 (∇⊥ h · a p b ghb)0,j 2 + (∇⊥ h · a p b ghb)0,0 4 = − h ∇⊥ h · a , b ih,G
2. The proof is similar to the previous one. 3. By Lemma 7.2, we have
h ∇hb , ∇hc ih,C= −h ∇h· ∇hb , c ih,G= −h ∇2hb , c ih,G
REFERENCES
[1] O. Botella and R. Peyret, Benchmark spectral results on the lid-driven cavity flow, Comput-ers and Fluids, 27 (1998), pp. 421–433.
