A GENERALIZED MAC SCHEME ON CURVILINEAR DOMAINS
YIN-LIANG HUANG ∗, JIAN-GUO LIU †, AND WEI-CHENG WANG ‡
Abstract.
We propose a new finite difference scheme for Navier-Stokes equations in primitive formulation on curvilinear domains. With proper boundary treatment and interplay between covariant and contra-variant components, the spatial discretization admits exact Hodge decomposition and energy identity. As a result, the pressure can be decoupled from the momentum equation with explicit time stepping. No artificial pressure boundary condition is needed. Numerical experiments demonstrate the robustness and efficiency of the scheme.
1. Introduction. In the numerical computation of the Navier-Stokes equation ut+ (u · ∇)u + ∇p = ν∇2u+ f in Ω (1.1) ∇ · u = 0 in Ω (1.2) u= 0 on Γ (1.3)
a key difficulty is the proper implementation of boundary conditions. Since the pressure p is not described by an evolutionary equation, one has to resort to the static relation (1.2) to obtain an elliptic equation for p:
(1.4) △p = −∇ · (u · ∇u) in Ω
An extra boundary condition is needed to solve for p from (1.4). This kind of boundary condition can be derived from (1.1). For example,
(1.5) ∂p
∂n = νn · ∇
2
u on Γ
or other equivalent forms. Since the right hand side of (1.5) involves high order derivatives, direct imple-mentation of (1.5) encounters various accuracy and stability issues, especially in the finite difference setting. The main purpose for using (1.4) is to ensure incompressibility (1.2). Finding stable and accurate pressure boundary condition to realize this goal has been the focus of many research papers in the literature. An alternative pressure boundary condition
(1.6) ∂p
∂n = −ν(∇ × ∇ × u) · n on Γ
has played a key role in the development of the projection methods. This Neumann boundary condition was first proposed by in [OID] to ensure consistency with the Poisson equation in the context of projection methods. The pressure Poisson equation (1.4) with boundary condition (1.6) can be used as replacement of (1.2). Well-posedness of this reformulation and equivalency to the original Navier-Stokes was shown in [LLP]. The incompressibility constraint can be realized up to truncation error for numerical methods based on solving pressure Poisson equation (1.4) and (1.6), provided the spatial resolution is sufficient. However, when the resolution becomes marginal, this class of methods often meet some difficulties in stability. Various divergence suppression techniques were developed to overcome these difficulties. There are also some difficulty in the finite difference implementation of these pressure boundary condition in curvilinear domains as the local expansion becomes more tedious in curvilinear coordinates.
∗Department of Mathematics, National Taiwan University, Taipei, 106, Taiwan. (yin1108@gmail.com)
†Department of Mathematics and Department of Physics, Duke University, Durham, NC 27708, USA (jian-guo.liu@duke.edu)
An alternative numerical approach is to treat the pressure as a Lagrangian multiplier to enforce the discrete divergence free constraint. This requires the spatial derivative to be compatible to each other so that the discrete analogue of fundamental differential identities are preserved. In this approach, the pressure is no longer solved via an elliptic PDE and there is no boundary condition involved for the pressure. The classical Mark-And-Cell (MAC) scheme can be interpreted as a typical example of this approach [An]. The spatial compatibility is also closely related to the exact Hodge decomposition which enables the decoupling of the pressure from the momentum equation and is key to the efficiency of the scheme.
Emphasise simplicity Emphasise exact Hodge decomposition as a new approach Direct transpose has been proposed before, credit Anderson
Motivated by the merits of the classical MAC scheme, we propose a novel numerical scheme on curvilinear coordinate that preserves the desired properties. This is done by careful construction of appropriate 2nd order finite difference operators such that the differential identities such as
(1.7) curl ◦ grad ≡ 0, div ◦ curl ≡ 0
remain valid in the discrete setting. In particular, the discrete Hodge decomposition for vector fields can be performed exactly and the corresponding linear system for the pressure is symmetric and semi-definite. No pressure boundary condition is needed. Overall, the resulting scheme is robust and efficient as demonstrated by our numerical examples.
The exact Hodge decomposition also leads to a very simple and elegant error analysis for the velocity field [HLW]. Rigorous error analysis for the classical MAC scheme was first obtained in [HW] and for MAC-like schemes on Cartesian grids in [W1]. The proof in [HW, W1] is based on high order Strang’s expansion. In contrast, the argument in [HLW] explores on the special structure of the spatial discretization and makes use of both the stream function and the discrete analogue of the differential identities (1.7). As a result, an optimal O(h2) error estimate is obtained provided the exact velocity is in C4and pressure in C3. This may
be the minimal regularity requirement in finite difference setting.
The rest of the paper is organized as follows. In section 2, we review the classical MAC scheme, the boundary treatment and associated Hodge decomposition. In section 3, we describe our generalized MAC scheme in curvilinear coordinates. The velocity components are located on the same location for convenience in both programming and applications. Representation of the Navier Stokes equation in a ‘skewed’ local coordinate leads to a natural discretization that gives rise to desired crucial discrete identities. One of the key issue here is to incorporate the boundary conditions into the finite difference operators so that exact summation by parts identity holds. The key properties of MAC scheme and exact discrete Hodge decomposition mentioned above is retained here even for the non-homogeneous boundary conditions. The 3D version of our scheme, as well as possible variants are documented in the appendix. Finally, we perform a systematic numerical test and report the results in section 4.
2. Classical MAC Scheme and discrete Hodge Decomposition. In this section, we review the classical MAC scheme and fundamental discrete identities associated with it. We first recall the essential ingredients that lead to the energy estimate and the well-posedness of the Navier-Stokes equation (1.1-1.3). Take the inner product with u on both sides of (1.1) and then integrate over Ω using the following facts,
hu, ∇pi = −h∇ · u, pi, (2.1)
The vector Laplacian ∇2 is symmetric and non-positive,
(2.2)
hN(u), ui = 0, (2.3)
where N (u) = u · ∇u and hu, vi = Z
Ω
u· v dx, it is easy to obtain the basic energy estimate
(2.4) 1 2 d dtkuk 2 + νk∇uk2= 0. 2
where kuk2
= hu, ui.
It is therefore desirable if a numerical scheme can preserve discrete analogue of (2.1-2.3). In particular, the stability of such a scheme is guaranteed. A well known example satisfying (2.1, 2.2) is the classical Mark-and-Cell (MAC) scheme proposed by Harlow and Welch [MAC]. The pressure and the components of the velocity field are placed on staggered grids in such a way that 2nd order centered difference, divergence free constraint and no-slip boundary conditions all fit naturally and elegantly with the placement of the variables. The 2D case is illustrated in Figure 2.1.
Fig. 2.1.
The time continuous version of the MAC scheme is given by
(2.5) ∂tu+ Nh(u) + ∇hp = ν∇2hu on ‘ → ’ and ‘ ↑ ’, ∇h· u = 0 on ‘ • ’, u= 0 on Γ.
The success of the MAC scheme may be attributed to the fact that the discrete analogue of (2.1,2.2) and (1.7) remain valid. These identities not only ensures the stability of the scheme, they also lead to exact Hodge decomposition of discrete vector fields which decouples the pressure equation from the momentum equation with explicit time stepping. The resulting scheme is robust and efficient. However, the use of staggered and uniform grids have limited the applicability and popularity of the MAC scheme. In addition, issues of high order time discretization and cell Reynolds number constraint have raised controversy and were not fully understood until the 90’s [EL].
In this paper, we propose a generalized MAC scheme for Navier-Stokes equation on curvilinear domains. The velocity field and pressure are placed on cell centers and grid points respectively. With a set of ‘skewed’ coordinates, the Navier-Stokes equation can be discretized naturally in such a way that the discrete analogue of (2.1-2.3) remain valid. As a consequence, the resulting scheme admits a discrete energy estimate. In addition, the scheme also preserves the vector identities (1.7) in discrete settings. This is key to discrete Hodge decomposition and plays an essential role in the efficiency, as well as rigorous error analysis of our scheme. Spatial discretizations with this property play important roles in many other fields such as computational elasticity and computational electromagnetic.
One of the main issues in accomplishing the discrete analogue of (2.1) and (2.2) is the boundary treatment of the spatial discretization. This is done by retaining a fraction of the original difference operators at the boundary and was introduced by Anderson [An] as ‘reduced’ operator in the Cartesian case. In the curvilinear case, an analogue of the reduced Laplacian in the skewed variables was proposed for elliptic interface problems in [Wa, HWa].
Another key to the success of the MAC scheme is the decoupling of the pressure from the momentum equation. Here we illustrate this fact with first order forward Euler discretization. Proper high order time-discretization and its connection with the cell Reynolds number can be found in [EL]. The fully discrete MAC scheme with explicit treatment for the viscous term is given by:
(2.6) un+1− un ∆t + Nh(u n ) + ∇hpn+1= ν∇2hun on ‘ → ’ and ‘ ↑ ’, ∇h· un+1= 0 on ‘ • ’, un+1= 0 on Γ.
Note that the divergence free constraint and the no-slip boundary condition all fit naturally with standard 2nd order differencing due to the staggered placement of the variables (u, v) and p.
Given (un, pn) at time tn with ∇
h· un = 0, the procedure of updating (un+1,pn+1) in (2.6) can be
decomposed into 2 steps: Step 1. Solve for u∗
from (2.7) u ∗ − un ∆t + Nh(u n ) = ν∇2hun
in the interior nodes of ‘→’ and ‘↑’. The tangential component of un on the ghost points is needed
in both the evaluation of ν(∇2
hun) · n and the realization of no-slip boundary condition un× n = 0.
Step 2. Find un+1with ∇
h· un+1= 0 and un+1· n|Γ= 0 from
(2.8) u
n+1− u∗
∆t + ∇hp
n+1= 0.
Step 1 is plain evaluation. The solvability of the step 2 is equivalent to the existence of discrete Hodge decomposition for u∗:
(2.9) u∗= un+1+ ∆t∇
hpn+1, ∇h· un+1= 0.
In practice, step 2 can be solved efficiently and accurately by first solving the pressure Poisson equation
(2.10) △hp = ∇h· u
∗
∆t .
One can think of (2.10, 2.12) as discretization of their continuum counterpart. From standard PDE theory, an extra pressure boundary condition needs to be supplied. This is usually done by an Neumann boundary condition
(2.11) Dnp =
u∗
· n ∆t on Γ which follows from (2.8).
After pn+1is solved, one then substitute back in (2.8) to get
(2.12) un+1= u∗− ∆t∇hpn+1.
The value of p on the ghost points is needed both in the pressure Poisson equation (2.10-2.11) and the updating (2.12). An alternative interpretation for (2.12) is to treat p pressure as a discrete Lagrangian multiplier purely to enforce the discrete divergence free constraint. The following observation explains why this is possible in practice.
Proposition 1. For classical MAC scheme (2.6), the result of un+1 from solving (2.10, 2.11) and
(2.12) is independent of the value of u∗
· n on Γ.
As a direct consequence of Proposition 1, the discrete pressure on cell centers within the computational domain can be solved as a Lagrangian multiplier without referring to the ghost values outside the computa-tional domain. Motivated by Proposition 1, we will construct proper finite difference operators in curvilinear coordinates that preserves discrete analogue of (2.1-2.3) and (1.7). In addition, the discrete pressure can be solved purely as Lagrangian multipliers without introducing artificial pressure boundary conditions.
3. Generalized MAC Scheme on Curvilinear Domains. Our scheme is based on discretization of the Navier-Stokes equation in rotational form:
(3.1)
ut+ ω × u + ∇p = −ν∇ × ω + f
ω = ∇ × u in Ω
∇ · u = 0
u = 0 on Γ
On a curvilinear domain, we place all three components of the velocity on cell centers x(ξ1 i+1 2, ξ 2 j+1 2, ξ 3 k+1 2),
while the pressure and the vorticity components are placed on the grid points x(ξ1
i, ξj2, ξk3). Here (ξ1, ξ2, ξ3)
is the coordinate in the computational domain while x is the position vector in physical domain Ω.
In addition to generalization to curvilinear coordinates, our scheme differs from the classical MAC scheme in the placement of the velocity components. One advantage of placing all three components of velocity at the same place is that (2.3) can be naturally realized in the discrete setting. Together with other vector identities, this ensures the stability of our scheme. In addition, the resulting discrete Laplacian for p is self-adjoint and non-positive definite as long as the cells are non-singular, regardless of the regularity of the grids. See the discussion in section 3.4. For the fully staggered case, where different components of u are placed on different positions as in the classical MAC case, positivity of the pressure equation on curvilinear domains may require extra assumption on smoothness of the grids ([BD], p147-150).
For simplicity of presentation, we start with the 2D case:
(3.2) ut+ ωu⊥+ ∇p = ν∇⊥ω + f ω = ∇⊥ · u in Ω ∇ · u = 0 u = 0 on Γ
The discretization of (3.2) and boundary treatment in curvilinear coordinate will be explained in detail in following subsections.
3.1. Differential operators in curvilinear coordinate. In the 2D case, x = (x, y) is the position vector in the physical domain and (ξ1, ξ2) the coordinate in the computational domain with mesh size
∆ξ1= h
1 and ∆ξ2= h2. We further introduce a new set of coordinates in the skewed direction by
(3.3) ‰ξ1,h2ξ1+ h1ξ2 ph2 1+ h22 , ‰ξ2, −h2ξ1+ h1ξ2 ph2 1+ h22 . as illustrated in Figure 3.1.
Once the a local coordinate is chosen, the intrinsic differential operators can be determined following standard procedure. Through out this paper, we use ‘‰’ to emphasis that the corresponding quantities are computed in the skewed variables ‰ξα. Denote by
(3.4) ‰ e1= ∂x ∂‰ξ1, ‰e2= ∂x ∂‰ξ2, ‰ e1= ∇‰ξ1, e‰2= ∇‰ξ2,
‰ξ2 = constant ‰ξ2= constant ξ2= constant ξ2 = constant ‰ξ1 = constant ‰ξ1= constant ξ1= constant ξ1= constant
Fig. 3.1. The computational domain (left) and the physical domain (right).
the metric tensors with respect to the skewed coordinate (‰ξ1, ‰ξ2) is then given by
(3.5) ‰gµν = ‰eµ· ‰eν, ‰gµν = ‰eµ· ‰eν, µ, ν = 1, 2.
and
(3.6) ‰g , det(‰gµν).
The following identities follow immediately from definition:
(3.7) p‰g = det ∂x ∂‰ξ (3.8) 2 X γ=1 ‰gµγ‰gγν= δνµ.
We use (‰u1, ‰u2) and (‰u1, ‰u2) to denote the covariant and contra-variant components of a vector field u in
the ‰ξαcoordinate:
u= ‰u1‰e1+ ‰u2‰e2= ‰u1e‰1+ ‰u2‰e2.
The transformation between co-variant and contra-variant components is given by
(3.9) u‰µ= 2 X γ=1 ‰gµγu‰γ, ‰uν = 2 X γ=1 ‰gγνu‰γ, µ, ν = 1, 2. 6
We summarize relevant formulae for (3.2) as follows: ∇p =∂p ∂‰ξ1‰ e1+ ∂p ∂‰ξ2‰ e2, ∇⊥ ω =√1 ‰g − ∂ω ∂‰ξ2‰ e1+ ∂ω ∂‰ξ1‰ e2 ! , ∇ · u =√1 ‰g ∂ ∂‰ξ1 (p‰g‰u1) + ∂ ∂‰ξ2 (p‰g‰u2) ! , ∇2p =∇ · ∇p = √1 ‰g 2 X µ,ν=1 ∂ ∂‰ξµ p ‰g‰gµν ∂p ∂‰ξν ! , u⊥ =p‰g− ‰u2‰g11+ ‰u1‰g21e‰1+ p ‰g− ‰u2‰g12+ ‰u1‰g22‰e2, ω =∇⊥· u = √1 ‰g ∂‰u2 ∂‰ξ1 − ∂‰u1 ∂‰ξ2 ! .
Whenever necessary, the co-variant and contra-variant components can be converted to each other using (3.9). For example, ∇p = ∂p ∂‰ξ1‰g 11 + ∂p ∂‰ξ2‰g 21 ! ‰ e1+ ∂p ∂‰ξ1‰g 12 + ∂p ∂‰ξ2‰g 22 ! ‰ e2.
3.2. Spatial Discretization. The choice of the skewed coordinate in motivated by the jump condition capturing scheme developed for elliptic interface problems [Wa]. It was shown that, using the skewed coordinate as independent variables, the interface jump conditions can be naturally incorporated into the finite difference operator. The resulting scheme is symmetric, definite and second order accurate even when the diffusion coefficient has a jump continuity across the material interface.
We now give detailed description on the spatial discretization of our scheme, starting with the metric tensor. In case the coordinate mapping (ξ1, ξ2) 7→ (x, y) is explicitly given, one can compute the metric
tensors from (3.4-3.8) and use it in the discretizations (3.29-3.35) below. Alternatively, one can compute the numerical metric tensors ‰ghαβfrom straightforward centered difference:
‰gh11 i+1 2,j+ 1 2 , xi+1,j+1‰ − xi,j h · xi+1,j+1− xi,j ‰ h , ‰gh22 i+1 2,j+ 1 2 , xi,j+1− xi+1,j ‰ h · xi,j+1− xi+1,j ‰ h , ‰gh12 i+1 2,j+ 1 2 =‰gh21 i+1 2,j+ 1 2 , xi+1,j+1− xi,j ‰ h · xi,j+1− xi+1,j ‰ h . where (3.10) ‰h, 2h1h2 ph2 1+ h22 = ∆‰ξ1= ∆‰ξ2
is the mesh size in the skewed directions ‰ξ1 and ‰ξ2. Note that the indices i, j refer to the ξ variables, not the ‰ξ ones.
The numerical contra-variant components (3.8) are defined through analogue of (3.8): (‰gαβh ), (‰g
h µν)
that is 2 X λ=1 ‰gµλh ‰g h λν= δµν.
The numerical Jacobian on cell centers and grids are given by (‰gh)i+1 2,j+ 1 2 , det ‰gh11 ‰gh12 ‰gh21 ‰g h 22 i+1 2,j+ 1 2 and p ‰ghi,j , 1 4 p ‰ghi+1 2,j+ 1 2 + p ‰ghi+1 2,j− 1 2 + p ‰ghi−1 2,j+ 1 2 + p ‰ghi−1 2,j− 1 2 .
Next, we introduce the discrete grad, div, curl and Laplacian. With the semi-staggered placement of the variables u, p and ω as in Figure 3.4, it is straight forward to discretize (3.2) using centered difference. The crucial issue here is the boundary treatment of these operators which plays an essential role in the stability of our scheme. We define the 1D reduced finite difference operator by
(3.11) (D′ b)i , b1 2 −0 h/2 i = 0 bi+ 1 2 −bi− 1 2 h 1 ≤ i ≤ M − 1 0−b N − 12 h/2 i = M
A direct consequence of the reduced differencing is the exact summation by parts identity as observed in [An]: Proposition 2. (3.12) h M X i=1 bi−1 2(Da)i− 1 2 = −h M X i=0 ′ ai(D′b)i
where the primed sum denotes half weight on the boundary: (3.13) M X i=0 ′ ai= 1 2a0+ M −1 X i=1 ai+ 1 2aN. We now introduce the following notations
Ωc , { (ξi−1 1 2, ξ 2 j−1 2) |1 ≤ i ≤ M, 1 ≤ j ≤ N } (3.14) ˚Ωg , { (ξ1 i, ξj2) |1 ≤ i ≤ M − 1, 1 ≤ j ≤ N − 1 } (3.15) ¯ Ωg , { (ξi1, ξ 2 j) |0 ≤ i ≤ M, 0 ≤ j ≤ N } (3.16) Γg , ¯Ωg\ ˚Ωg (3.17) ¯ Ωge , { (ξ 1 i, ξ 2 j) ∈ ¯Ωg| i + j is even } (3.18) ¯ Ωgo , { (ξ 1 i, ξ 2 j) ∈ ¯Ωg| i + j is odd } (3.19) and denote by L2( ¯Ω
g, R) the collection of real valued functions on ¯Ωg and L2(Ωc, R2) the vector fields on Ωc:
L2( ¯Ω g, R) ,{a : ¯Ωg → R} (3.20) L2(Ω c, R) ,{b : Ωc→ R} (3.21) L2(Ω c, R2),{ u : Ωc→ R2} (3.22) 8
and L2( ¯Ω g, R)/R2, n a ∈ L2( ¯Ω g, R) X ¯ Ωge ′ (p ‰gha)i,j = 0 = X ¯ Ωgo ′ (p ‰gha)i,j o (3.23) L2 c( ¯Ωg, R) ,{a ∈ L2( ¯Ωg, R) a = constant on Γg} (3.24)
In view of Proposition 2, it is natural to define the (reduced) difference operator in the skewed variables as Definition 1. 1. For a ∈ L2( ¯Ω g, R), (3.25) ( ‰D1a)i−1 2,j− 1 2 , ai,j− ai−1,j−1 ‰ h , ( ‰D2a)i− 1 2,j− 1 2 , ai−1,j− ai,j−1 ‰ h . 2. For b ∈ L2(Ω c, R), (3.26) ( ‰D′1b)i,j, bi+1 2,j+ 1 2 − bi− 1 2,j− 1 2 ‰ h on ˚Ωg; bi+1 2, 1 2 ‰ h/2 0 < i < M, j = 0; b1 2, 1 2 ‰ h/4 (i, j) = (0, 0); 0 (i, j) = (M, 0). (3.27) ( ‰D′2b)i,j, bi−1 2,j+ 1 2 − bi+ 1 2,j− 1 2 ‰ h on ˚Ωg; bi−1 2, 1 2 ‰ h/2 0 < i < M, j = 0; 0 (i, j) = (0, 0); bM −1 2,12 ‰ h/4 (i, j) = (M, 0).
Here in Definition 1 and the rest of the paper, we use the prime to denote that the reduced difference is applied on boundary when the standard finite difference requires grid points outside the computational domain.
The (reduced) finite difference in the ‰ξ variables extends naturally to the discrete grad, div, curl and Laplacian operators as follows:
Definition 2. For a ∈ L2( ¯Ω g, R), we define (3.28) ∇‰h: L2( ¯Ωg, R) 7→ L2(Ωc, R2), ∇‰ha, ( ‰D1a) ‰e1+ ( ‰D2a) ‰e2 and (3.29) ∇‰⊥h : L 2 ( ¯Ωg, R) 7→ L2(Ωc, R2), ∇‰ ⊥ ha, − ‰ D2a p‰gh ‰ e1+ ‰ D1a p‰gh ‰ e2.
Definition 3. Let u = ‰u1e‰1+ ‰u2‰e2∈ L2(Ω
c, R2). We define (3.30) ∇‰′h· : L2(Ω c, R2) 7→ L2( ¯Ωg, R), ∇‰ ′ h· u = 1 p‰gh ‰ D′1( p ‰gh‰u 1 ) + ‰D′2( p ‰ghu‰ 2 )
Fig. 3.2. Schematic illustration of ‰∇′h· u. Fig. 3.3. Schematic illustration of ‰∇ ⊥′ h · u. and (3.31) ∇‰⊥′h · : L2(Ω c, R2) 7→ L2( ¯Ωg, R), ∇‰ ⊥′ h · u = 1 p‰gh ( ‰D′1u‰2− ‰D ′ 2u‰1). That is (3.32) ( ‰∇′h· u)i,j= 1 (p‰gh)i,j (p‰ghu‰ 1 )i+1 2,j+ 1 2 − (p‰ghu‰ 1 )i−1 2,j− 1 2 ‰ h +(p‰ghu‰ 2 )i−1 2,j+ 1 2 − (p‰ghu‰ 2 )i+1 2,j− 1 2 ‰ h ! on ˚Ωg; 2 (p‰gh)0,j (p‰ghu‰ 1 )1 2,j+12 − (p‰gh‰u 2 )1 2,j−12 ‰ h , i = 0, 0 < j < N ; 4 (p‰gh)0,0 (p‰ghu‰ 1 )1 2, 1 2 ‰ h , (i, j) = (0, 0), and (3.33) ( ‰∇⊥′h · u)i,j= 1 (p‰gh)i,j (‰u2)i+1 2,j+ 1 2 − (‰u2)i− 1 2,j− 1 2 ‰ h −(‰u1)i− 1 2,j+ 1 2 − (‰u1)i+ 1 2,j− 1 2 ‰ h ! on ˚Ωg; 2 (p‰gh)0,j (‰u2)1 2,j+ 1 2 + (‰u1) 1 2,j− 1 2 ‰ h , i = 0, 0 < j < N ; 4 (p‰gh)0,0 (‰u2)1 2, 1 2 ‰ h , (i, j) = (0, 0).
See also Figure 3.2 and Figure 3.3.
Finally, the discrete Laplacian is defined in a similar way: Definition 4. For a ∈ L2( ¯Ω g, R), define (3.34) △‰′h: L2( ¯Ωg, R) 7→ L2( ¯Ωg, R), △‰ ′ ha = 1 p‰gh 2 X µ,ν=1 ‰ D′µ( p ‰gh‰g µν h D‰νa). 10
That is, (3.35) ( ‰△′ha)i,j = 1 p‰ghi,j (‰qh 11‰ D1a + ‰qh 12‰ D2a)i+1 2,j+ 1 2 ‰ h + (‰qh 21‰ D1a + ‰qh 22‰ D2a)i−1 2,j+ 1 2 ‰ h −(‰qh 11‰ D1a + ‰qh 12‰ D2a)i−1 2,j− 1 2 ‰ h − (‰qh 21‰ D1a + ‰qh 22‰ D2a)i+1 2,j− 1 2 ‰ h ! on ˚Ωg; 2 p‰gh0,j (‰qh 11‰ D1a + ‰qh 12‰ D2a)1 2,j+ 1 2 ‰ h − (‰qh 21‰ D1a + ‰qh 22‰ D2a)1 2,j− 1 2 ‰ h ! , 0 < j < N ; 4 p‰gh0,0 (‰qh 11‰ D1a + ‰qh 12‰ D2a)1 2, 1 2 ‰ h , (i, j) = (0, 0), where ‰qαβh =p‰gh‰g αβ h .
We are now ready to state the key Lemma associated with the reduced difference operators. Define the discrete inner products
h u , v iΩc= h1h2 M −1 X i=0 N −1 X j=0 (u · v)√gh i+1 2,j+ 1 2 = h1h2 M −1 X i=0 N −1 X j=0 (‰u1‰v1+ ‰u2‰v2)√gh i+1 2,j+ 1 2 = h1h2 M −1 X i=0 N −1 X j=0 ‰ u1‰v1+ ‰u2‰v2 √ gh i+1 2,j+12 , u, v ∈ L2(Ωc, R2). (3.36) (3.37) h a , b iΩ¯g = h1h2 M X i=0 ′ N X j=0 ′ a b√gh i,j, a, b ∈ L 2( ¯Ω g, R)
and the corresponding norms
kuk2Ωc= h u , u iΩc, kak 2 ¯ Ωg = h a , a iΩ¯g, where √gh, 2h1h2 h2 1+ h 2 2
p‰ghis the numerical Jacobian with respect to the default coordinate (ξ1, ξ2). Applying
Proposition 2 in the skewed directions ‰ξ1and ‰ξ2, it is easy to derive the following discrete identities:
Lemma 3.1. Let u ∈ L2(Ω c, R2) and a ∈ L2( ¯Ωg, R). We have 1. (3.38) h u , ‰∇ha iΩc = −h ‰∇ ′ h· u , a iΩ¯g 2. (3.39) h u , ‰∇⊥ha iΩc = −h ‰∇ ⊥′ h · u , a iΩ¯g 3. (3.40) ∇‰′h· ‰∇ha = ‰∇⊥′h · ‰∇ ⊥ ha = ‰△ ′ ha on ¯Ωg; 4. If a ∈ L2( ¯Ω g, R), then (3.41) ∇‰′h· ‰∇⊥ha = ‰∇⊥′h · ‰∇ha = 0 on ˚Ωg. In addition, if a ∈ L2 c( ¯Ωg, R), then (3.42) ∇‰′h· ‰∇⊥ha = ‰∇⊥′h · ‰∇ha = 0 on ¯Ωg.
Fig. 3.4. Positions of velocity field (տր), vorticity and pressure (•) for the generalized MAC scheme (3.43).
3.3. Generalized MAC Scheme. On a curvilinear domain, we place both components of the velocity on cell centers x(ξ1
i+1 2, ξ
2 j+1
2), while the pressure and the vorticity are both placed on the grid points x(ξ
1 i, ξ2j)
as shown in Figure 3.4.
Our scheme can be summarized as follows: Solve for u ∈ C1([0, T ]; L2(Ω
c, R2)) and p ∈ C0([0, T ]; L2( ¯Ωg, R)) such that
(3.43) ut+ ¯ωu⊥+ ‰∇hp = ν ‰∇ ⊥ hω + f on Ωc ω = ∇‰⊥′h · u on ¯Ωg ‰ ∇′h· u = 0 on ¯Ωg where ¯ωi−1 2,j− 1 2 = 1
4(ωi,j+ ωi−1,j+ ωi,j−1+ ωi−1,j−1).
As a direct consequence of (3.38-3.39), we have the following discrete energy estimate for (3.43): Lemma 3.2. Let (u, ω) be a solution to (3.43) with f = 0, then
1 2 d dtkuk 2 Ωc+ νkωk 2 ¯ Ωg = 0.
It is worth noting that the no-slip condition u|Γ = 0 has been implicitly incorporated into the operators
‰
∇⊥′h · and ‰∇ ′
h· in (3.43). In case of non-homogeneous data, u = ubon Γ, the corresponding reduced operators
should be modified accordingly. Let F ∈ L2(Ω
c, R) and f ∈ L2(Γ, R). With slight abuse of notations, we define the extended reduced
operators, still denoted by ‰D′1, ‰D ′ 2, as follows: (3.44) D‰′1, D‰ ′ 2 : L 2 (Ωc, R) × L2(Γ, R) 7→ L2( ¯Ωg, R), 12
‰ D′1(F ⊕ f)i,j, Fi+1 2,j+ 1 2 − Fi− 1 2,j− 1 2 ‰ h on ˚Ωg; Fi+1 2,12 − fi,0 ‰ h/2 0 < i < M, j = 0; F1 2, 1 2 − (2f0,0+ f0,1+ f1,0)/4 ‰ h/4 (i, j) = (0, 0). (3.45) ‰ D′2(F ⊕ f)i,j, Fi−1 2,j+ 1 2 − Fi+ 1 2,j− 1 2 ‰ h on ˚Ωg; Fi−1 2, 1 2 − fi,0 ‰ h/2 0 < i < M, j = 0; f0,1− f1,0 ‰ h (i, j) = (0, 0). (3.46)
One can define ‰∇′h· and ∇‰ ⊥′
h · : L2(Ωc, R2) × (L2(Γ, R))27→ L2( ¯Ωg, R) in a natural way.
When the boundary data is identically zero, (3.45-3.46) reduce to the original reduced operators,
(3.47) D‰′α(F ⊕ 0) = ‰D ′ αF. (3.48) ∇‰′h· u = ‰∇ ′ h· (u ⊕ 0), ∇‰ ⊥′ h · u = ‰∇ ⊥′ h · (u ⊕ 0).
The corresponding scheme for inhomogeneous boundary velocity is given by
(3.49) ut+ ¯ωu⊥+ ‰∇hp = ν ‰∇ ⊥ hω + f on Ωc ω = ∇‰⊥′h · (u ⊕ ub) on ¯Ωg ‰ ∇′h· (u ⊕ ub) = 0 on ¯Ωg
subject to the following compatibility condition for the boundary velocity ub:
(3.50) h 1Ω¯ ge, ‰∇ ′ h· (0 ⊕ ub) iΩ¯ g = 0 = h 1Ω¯go, ‰∇ ′ h· (0 ⊕ ub) iΩ¯ g.
To see this, we first note that the kernel of ‰∇hconsists of linear combinations of indicator functions of ¯Ωge
and ¯Ωgo:
(3.51) ker( ‰∇h) = span{1Ω¯
ge, 1Ω¯go}.
From (3.47, 3.38, 3.51), it follows that for all c1, c2∈ R,
h c11Ω¯ ge+ c21Ω¯go, ‰∇ ′ h· (u ⊕ ub) iΩ¯ g =h c11Ω¯ge + c21Ω¯go, ‰∇ ′ h· (u ⊕ 0) + ‰∇ ′ h· (0 ⊕ ub) iΩ¯ g =h c11Ω¯ge + c21Ω¯go, ‰∇ ′ h· (0 ⊕ ub) iΩ¯g (3.52)
In view of (3.52) and the 3rd equation of (3.49), the compatibility condition (3.50) follows. It is easy to see that both conditions in (3.50) can be interpreted as discrete analogues of
(3.53)
Z
Γ
3.4. Time stepping and Hodge Decomposition. As in section 2, we illustrate the time stepping by the forward Euler method:
(3.54) un+1− un ∆t + (¯ωu ⊥)n+ ‰∇ hpn+1 = ν ‰∇ ⊥ hωn+ fn on Ωc ωn = ∇‰⊥′ h · (un⊕ unb) on ¯Ωg ‰ ∇′h· (un+1⊕ un+1b ) = 0 on ¯Ωg
In high Reynolds number calculations, a high order Runge-Kutta method such as RK4 is needed for stability consideration [EL].
Given (un, pn) with ‰∇′h· (un⊕ unb) = 0, (3.54) is solved via the following steps:
Step 1. Evaluate ωn up to the boundary
(3.55) ωn= ‰∇⊥′h · (un⊕ unb) on ¯Ωg.
Step 2. Evaluate u∗ on cell centers
(3.56) u ∗ − un ∆t + (¯ωu ⊥ )n= ν ‰∇⊥hωn+ f n on Ωc.
Step 3. Solve for (un+1, pn+1) such that
(3.57) un+1− u∗ ∆t + ‰∇hp n+1= 0 on Ω c, ‰ ∇′h· (un+1⊕ un+1b ) = 0 on ¯Ωg,
This is the (inhomogeneous) Hodge decomposition for u∗
and can be performed as follows: Step 3-1. Solve for pn+1 up to the boundary from
(3.58) △‰′hpn+1= 1 ∆t∇‰ ′ h· (u ∗ ⊕ un+1b ) on ¯Ωg.
Step 3-2. Update un+1from
(3.59) u
n+1− u∗
∆t + ‰∇hp
n+1= 0 on Ω
c.
It follows from (3.40, 3.48, 3.58, 3.59) that ‰ ∇′h· (un+1⊕ un+1b ) = ∇‰ ′ h· (u∗⊕ un+1b ) + ‰∇ ′ h· (un+1− u∗) ⊕ 0 = ∇‰′h· (u ∗ ⊕ un+1b ) + ‰∇ ′ h· (un+1− u ∗ ) = 0.
Step 1 can be viewed as a vorticity boundary condition that incorporates the tangential component of ub. On the other hand, step 3 depends on ub only through its normal component. In step 3-1, the pressure
is solved as a Lagrangian multiplier without introducing the ghost value and artificial pressure boundary conditions. From (3.38) and (3.40), it is easy to see that
(3.60) h q , ‰△′hp iΩ¯g = −h ‰∇hq , ‰∇hp iΩc = −h1h2 M X i=1 N X j=1 √g h X µ,ν=1,2 ‰gµνh ( ‰Dµp)( ‰Dνq)i−1 2,j−12. 14
Ω
Fig. 3.5. A domain patched by multiple non-overlapping coordinate charts.
In other words, ‰△′h is self-adjoint with respect to the inner product (3.37) and non-positive definite as long
as each of the 2 × 2 matrices {‰gµνh }2µ,ν=1
i−1
2,j−12 is positive definite on the cell centers. This only requires
the cells to be non-singular. No further grid regularity is needed for stability concerns. In addition, (3.61) ker( ‰△′h) = ker( ‰∇h) = span{1Ω¯
ge, 1Ω¯go}
thus (3.58) is solvable if and only if (3.62) h c11Ω¯ge + c21Ω¯go, ‰∇
′ h· (u
∗
⊕ unb) iΩ¯g = 0 for all c1, c2∈ R.
In view of (3.52), it follows that the solvability condition for (3.58) is exactly the compatibility condition (3.50). In addition, it is worth noting that (3.58) can be decoupled and solved on ¯Ωge and ¯Ωgo separately.
Remark 1. In case of complex geometry, it may be necessary to patch the domain with several non-overlapping coordinate charts as illustrated in Figure 3.6. In this case, the metric tensor has a jump dis-continuity across coordinate boundaries. We define the discrete divergence operator at a grid point P on the coordinate interface by (3.63) ( ‰∇h· u)P = p‰gh 2 ∇‰ ′ h· u P(+)+ p‰gh 2 ∇‰ ′ h· u P(−) (p‰gh)P(+)+ (p‰gh)P(−)
At a multiple coordinate junction Q, such as the one given in Figure 3.6, the formula becomes (3.64) ( ‰∇h· u)Q = X Q∈ith chart (p‰gh)Q(i) −1 X Q∈ith chart p‰gh 4 ∇‰ ′ h· u Q(i)
A similar formula applies to discrete curl and discrete Laplacian. In this way, it is easy to see that the summation by parts identities in Lemma 3.1 remains valid. This same treatment for discrete Laplacian on material interface has been proposed for elliptic interface problems in [Wa, HW].
Remark 2. A natural approach for realizing the no-slip boundary condition associated with the semi-staggered grid is the reflection boundary condition
(3.65) u−1 2,j+ 1 2 = −u 1 2,j+ 1 2
In general, the reflection boundary condition (3.65) gives rise to a different boundary treatment from the reduced differencing. For example, if Ωh= Ω = (0, 1) × (0, 1) with h1= h2= h, then the reduced curl gives
(3.66) ∇‰⊥′h · u0,j = 2(‰u2)1 2,j+ 1 2 + (‰u1) 1 2,j− 1 2 ‰ h .
Ω− Ω + Γ Ω1 Ω2 Ω3 P Q
Fig. 3.6. Schematic illustration of the discrete divergence at a point P on the coordinate interface Γ (left) and a triple coordinate junction Q (right).
while the full curl operator, with the reflection boundary condition (3.65) gives
(3.67) ∇‰⊥h · u0,j = (‰u2)1 2,j+ 1 2 + (‰u1) 1 2,j− 1 2 ‰ h + (‰u1)1 2,j+ 1 2 + (‰u2) 1 2,j− 1 2 ‰ h .
4. Numerical Result. In this section, we report several numerical test results including standard accuracy check and benchmark test problems with classical RK4 as time stepping.
Example 1: Convergence Test
We start with standard accuracy check for the spatial discretization. We take the exact solution to be (4.1) pe= cos(t) cos(x) sin(y), ψe(x, y) = cos(t) xy(1 − x)(k(x) − y)2
, ue= ∇⊥
ψe
on a curvilinear domain Ω =(x, y) | 0 ≤ x ≤ 1, 0 ≤ y ≤ k(x) with k(x) = 1 + 0.2 sin(2πx), and generate the corresponding forcing term f with ν = 10−3.
We take a simple coordinate mapping
(4.2) x(ξ1, ξ2) = ξ1, y(ξ1, ξ2) = ξ2 1 + 0.2 sin(2πξ1)
and RK4 for time discretization. The result at T = 3.2 is given in Table 4.1 with clean second order accuracy as expected.
Table 4.1
Absolute error and rate of convergence of computed velocity, vorticity and pressure in L∞and L2norms for Example 1.
M × N kue− uk
∞ Order kue− ukL2 Order kωe− ωk∞ Order kpe− pk∞ Order
100 × 100 3.163E − 5 8.012E − 6 2.094E − 3 2.890E − 5
200 × 200 8.767E − 6 1.85 2.011E − 6 1.99 5.265E − 4 1.99 7.224E − 6 2.00 400 × 400 2.297E − 6 1.93 5.033E − 7 2.00 1.317E − 4 2.00 1.806E − 6 2.00 800 × 800 5.876E − 7 1.97 1.258E − 7 2.00 3.294E − 5 2.00 4.515E − 7 2.00
Example 2: Benchmark Test: Lid-Driven Cavity Flow
We continue with the benchmark problem of lid-driven cavity flow [Bu, GGS, BP]. The results for Re = 1, 000 with 256 × 256 cells and 10, 000 with 360 × 360 cells are presented in Figure 4.1 and Figure 4.2
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. 4.1. Vorticity contour plot of the lid-driven cavity flow (Example 2) at Re=1,000.
respectively. Both are simulated till T ≈ 1.5√Re: T1000= 48.0, T10000= 255.0. The Re = 1, 000 case agrees
well with the benchmark result in [BP]. The singularities around the corners are quite well resolved. Example 3: Lid-Driven Cavity Flow with Bottom Topography
In this example, we consider the cavity with a bottom topography Ω = (x, y) | − 1 ≤ x ≤ 0, −1 + 0.2 sin(2πx) ≤ y ≤ 0 at Re = 2, 000. The result at T = 80.0 with 250 × 250 cells is shown in Figure 4.3.
Acknowledgment. This work is partially supported by National Science Foundation, Center of Math-ematical Sciences at NTHU, National Science Council and the National Center for Theoretical Sciences in Taiwan.
REFERENCES
[An] C. R. Anderson, Derivation and solution of the discrete pressure equations for the incompressible Navier-Stokes equa-tions, preprint LBL-26353, Lawrence Berkeley Laboratory, Berkeley, CA, 1988.
[BCG] J. B. Bell, P. Colella and H. Glaz, A second-order projection method for the incompressible Navier-Stokes equations, J. Comp. Phys. 85 (1989) 257–283.
[BD] F. Bertagnolio and O. Daube Solution of the div-curl problem in generalized curvilinear coordinates J. Comp. Phys. 138(1997) 121–152.
[BK] R. S. Bernard and H. Kapitza, How to discretize the pressure gradient for curvilinear MAC grids, J. Comput. Phys., 99(1992) 288 – 298.
[BP] O. Botella and R. Peyret, Benchmark spectral results on the lid-driven cavity flow, Comput. & Fluids, 27 (1998) 421–433.
[Bu] O. R. Burggraf, Analytical and numerical studies of the structure of steady separated flows, Journal of Fluid Mechanics, 24(1966) 113-151.
Driven Cavity Flow, contour of numerical ω, 360x360, Re=10000, time=255 3 2 −1 −2 3 0 0 −0.5 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
Fig. 4.2. Vorticity contour plot of the lid-driven cavity flow (Example 2) at Re=10,000.
[Ch] A. J. Chorin, Numerical solution of the Navier-Stokes equations, Math. Comp., 22 (1968) 745–762
[EL] W. E and J.-G. Liu, Vorticity boundary condition and related issues for finite difference schemes, J. Comput. Phys. 124 (1996) 368–382.
[FPT] M. Fortin, R. Peyret and R. Temam, R´esolution num´erique des ´equations de Navier-Stokes pour un fluide incompressible. J. M´ecanique 10 (1971), 357–390.
[GGS] U. Ghia, K. N. Ghia and C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comp. Phys. 48, (1982) 387–411.
[HW] Thomas Y. Hou and Brian T. R. Wetton, Second-Order Convergence of a Projection Scheme for the Incompressible Navier-Stokes Equations with Boundaries SIAM. J. Num. Anal., 30 (1993) 609 – 629.
[HLW] Y.-L. Huang, J.-G. Liu and W.-C. Wang, Error Analysis of the Generalized MAC Scheme, Preprint (aXiv XXX). [KM] J. Kim and P. Moin, Application of a fractional-step method to incompressible Navier-Stokes equations. J. Comp. Phys.
59(1985) 308–323.
[MAC] F.H. Harlow and J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids, 8 (1965), 2182–2189.
[JL] H. Johnston and J.-G. Liu, Accurate, stable and efficient Navier-Stokes solvers based on explicit treatment of the pressure term, J. Comp. Phys. 199 (1) (2004) 221–259.
[LLP] J.-G. Liu, J. Liu and R. Pego, Stability and convergence of efficient Navier-Stokes solvers via a commutator estimate, Comm. Pure Appl. Math. 60 (2007) 1443–1487.
[LHZW] J. Li, M. Hesse, J. Ziegler and A. W. Woods, An arbitrary Lagrangian Eulerian method for moving-boundary problems and its application to jumping over water, J. Comput. Phys., 208 (2005) 289 – 314.
[LW1] J.-G. Liu and W.C. Wang, Energy and helicity preserving schemes for hydro- and magnetohydro-dynamics flows with symmetry, J. Comput. Phys., 200 (2004) 8 – 33.
[LW2] J.-G. Liu and W.C. Wang, An energy preserving MAC-Yee scheme for the incompressible MHD equation, J. Comput. Phys. 174 (2001), 12–37.
[Ne] J.-C. N´ed´elec, A new family of mixed finite elements in R3, Numer. Math. 50, (1986) 57–81.
[OID] S. A. Orszag, M. Israeli and M. Deville, Boundary conditions for incompressible flows. J. Sci. Comput. 1, (1986) 75–111. 18
[RT] P.-A. Raviart and J. M. Thomas, A mixed finite element method for 2nd order elliptic problems, in Mathematical Aspects of Finite Element Methods, Lecture Notes in Mathematics, Springer, Berlin, Vol. 606 (1975) 292–315. [TC] D. P. Trebotich and P. Colella, A projection method for incompressible viscous flow on moving quadrilateral grids, J.
Comput. Phys., 166 (2001) 191 – 217.
[N1] R.A. Nicolaides, Analysis and Convergence of the MAC Scheme I. The Linear Problem SIAM Journal on Numerical Analysis, vol 29 (1992) 1579-1591.
[N2] R.A. Nicolaides and X. Wu, Analysis and convergence of the MAC scheme II, Navier-Stokes equations, Math. Comp. 65 (1996), 29-44.
[HW1] H Huang, B.R. Wetton, Discrete compatibility in finite difference methods for viscous incompressible fluid flow Journal of Computational Physics, 1996
[HWa] Y.-L. Huang, W.-C. Wang, A Monotone Finite Difference Scheme for Elliptic Interface Problems on Arbitrary Domains, preprint.
[Wa] W.-C. Wang, A Jump Condition Capturing Finite Difference Scheme for Elliptic Interface Problems, SIAM J. SCI. COMPUT. Vol. 25, No. 5 (2004) 1479-1496.
[W1] B.R. Wetton. Analysis of the spatial error for a class of finite difference methods for viscous incompressible flow, SIAM J. NUMER. ANAL. Vol. 34, No. 2 (1997) 723-755,
Appendix A. A Variant of the Generalized MAC Scheme. An alternative placement of the unknown variables
(A.1) u∈ L2(˚Ωg, R2), p, ω ∈ L2(Ωc, R)
leads to a variant of the scheme (3.43). In this setting, we place both components of the velocity on grid points, pressure and vorticity on cell centers. Note that the boundary velocity ub∈ L2(Γg, R2) is prescribed
data and not listed in (A.1) as part of active variables.
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
Associated with this new placements are the following function spaces
(A.2) L2(˚Ωg, R2), {u : ˚Ωg → R2}, L2(Γg, R2), {ub: Γg → R2}, L2(Ωc, R) , {a : Ωc→ R}.
and finite difference operators ‰ ∇◦h:L 2 (Ωc, R) 7→ L2(˚Ωg, R), ∇‰ ◦ ha, ( ‰D1a) ‰e1+ ( ‰D2a) ‰e2 on ˚Ωg; (A.3) ‰ ∇◦⊥h :L 2(Ω c, R) 7→ L2(˚Ωg, R), ∇‰ ◦⊥ h a, − ‰D2a p‰gh ‰ e1+ ‰ D1a p‰gh ‰ e2 on ˚Ωg. (A.4)
Note that the operators ‰∇◦hand ‰∇ ◦⊥
h are exactly the standard finite difference operators as (3.28) and (3.29),
except (A.3) and (A.4) are defined for vector fields on Ωc. We use the superscript ‘◦’ to emphasize that ‰∇ ◦ h
and ‰∇◦⊥h are ‘interior’ gradient operators whose range are vector fields defined only on the interior grids.
For inhomogeneous boundary velocity, the operators ‰∇◦′h· and ‰∇ ◦⊥′
h · involve both the active variable
u∈ L2(˚Ω
g, R2) and the data ub∈ L2(Γg, R2). Denote by U = u ⊕ ub, that is
(A.5) U‰α= ( ‰ uα on ˚Ωg, ‰ uαb on Γg, ‰ Uα= ( ‰ uα on ˚Ωg, (‰ub)α on Γg, , α = 1, 2, we define ‰ ∇◦′h· : L 2 (˚Ωg, R2) ⊕ L2(Γg, R2) 7→ L2(Ωc, R), ∇‰ ◦′ h · (u ⊕ ub) = 1 p‰gh ‰ D1( p ‰ghU‰ 1 ) + ‰D2( p ‰ghU‰ 2 ) (A.6) ‰ ∇◦⊥′h · : L2(˚Ωg, R2) ⊕ L2(Γg, R2) 7→ L2(Ωc, R), ∇‰ ◦⊥′ h · (u ⊕ ub) = 1 p‰gh ( ‰D1U‰2− ‰D2U‰1) (A.7) As before, for u ∈ L2(˚Ω g, R2), (A.8) ∇‰◦′h · u , ‰∇ ◦′ h · (u ⊕ 0), ∇‰ ◦⊥′ h · u , ‰∇ ◦⊥′ h · (u ⊕ 0). Finally, (A.9) △‰◦′h : L2(Ωc, R) 7→ L2(Ωc, R), △‰ ◦′ ha = ‰∇ ◦′ h · ‰∇ ◦ ha = ‰∇ ◦′ h · ( ‰∇ ◦ ha) ⊕ 0.
Fig. A.1. Positions of velocity field (տր), vorticity and pressure (•) for the scheme (A.17). 20
Define the discrete inner products h u , v i˚Ωg = h1h2 M −1 X i=1 N −1 X j=1 (u · v)√ghi,j, u, v ∈ L2(˚Ωg, R2), (A.10) h a , b iΩc = h1h2 M X i=1 N X j=1 a b√ghi−1 2,j−12, a, b ∈ L 2(Ω c, R). (A.11)
The counterpart of Lemma 3.1 is also valid: Lemma A.1. If u ∈ L2(˚Ω
g, R2) and b ∈ L2(Ωc, R). Then we have
1. (A.12) h u , ‰∇◦ha i˚Ω g = −h ‰∇ ◦′ h · u , a iΩc 2. (A.13) h u , ‰∇◦⊥h a i˚Ωg = −h ‰∇ ◦⊥′ h · u , a iΩc 3. (A.14) ∇‰◦′h · ( ‰∇ ◦ ha) = ‰∇ ◦⊥′ h · ( ‰∇ ◦⊥ h a) = ‰△ ◦′ ha on Ωc; 4. If a ∈ L2(Ω c, R), then (A.15) ∇‰◦′h · ( ‰∇ ◦⊥ h a) = ‰∇ ◦⊥′ h · ( ‰∇ ◦ ha) = 0 on ˚Ωc, where ˚Ωc =(ξi−1 1 2, ξ 2 j−1 2) ∈ Ωc : 2 ≤ i ≤ M − 1, 2 ≤ j ≤ N − 1 . If in addition, a is constant on Ωc\ ˚Ωc, then (A.16) ∇‰◦′h · ( ‰∇◦⊥h a) = ‰∇◦⊥′h · ( ‰∇ ◦ ha) = 0 on Ωc.
The alternative scheme is thus given by: Solve for u ∈ C1 [0, T ]; L2(˚Ω
g, R2) and p ∈ C0 [0, T ]; L2(Ωc, R) such that
(A.17) ut+ ¯ωu⊥+ ‰∇ ◦ hp = ν ‰∇ ◦⊥ h ω + f on ˚Ωg ω = ∇‰◦⊥′h · (u ⊕ ub) on Ωc ‰ ∇◦′h · (u ⊕ ub) = 0 on Ωc
The compatibility condition for the boundary velocity ub is
(A.18) h 1Ωce, ‰∇ ◦′ h · (0 ⊕ ub) iΩc= 0 = h 1Ωco, ‰∇ ◦′ h · (0 ⊕ ub) iΩc, where (A.19) Ωce, (ξ 1 i−1 2, ξ 2 j−1 2) ∈ Ωc, i + j is even , Ωco , (ξ 1 i−1 2, ξ 2 j−1 2) ∈ Ωc, i + j is odd . Furthermore,
(A.20) ker( ‰△◦′h) = ker( ‰∇ ◦
Appendix B. The 3D Scheme.
Generalization of our scheme with the semi-staggered variables placement as either (B.1) u∈ L2(Ωc, R3), p ∈ L2( ¯Ωg, R), ω∈ L2( ¯Ωg, R3).
or
(B.2) u∈ L2(˚Ωg, R3), p ∈ L2(Ωc, R), ω∈ L2(Ωc, R3).
are both straight forward once a skewed coordinate is chosen. We only elaborate with the version (B.1). Take, for example, (
ξ1,ξ2,ξ3) as shown in left-top corner of Figure B.1:
(B.3) ξ1, h3ξ2+ h2ξ3 ph2 2+ h23 , h1= ∆ ξ1= 2h2h3 ph2 2+ h23 ξ2, h1ξ 3+ h 3ξ1 ph2 3+ h21 , h2= ∆ ξ2= 2h3h1 ph2 3+ h21 ξ3, h2ξ1+ h1ξ2 ph2 1+ h22 , h3= ∆ ξ3= 2h1h2 ph2 1+ h22
The discrete metric tensors {
gh αβ}3α,β=1, { gαβh }3 α,β=1 and q gh with respect to ( ξ1,ξ2,ξ3) as well as the
reduced difference operators
D′ αand
Dαcan be defined in a similar way as in the 2D case.
We therefore arrive at the first version of our scheme for 3D Navier-Stokes equation (3.1): Scheme 3D-skew: (B.4) ut+ ¯ω× u + ∇hp = ν ∇h× ω + f on Ωc ω = ∇′ h× u on ¯Ωg ∇′ h· u = 0 on ¯Ωg where ∇h:L2( ¯Ωg, R) 7→ L2(Ωc, R3), ∇hp, ( D1p) e1+ ( D2p) e2+ ( D3p) e3, (B.5) ∇′ h· :L2(Ωc, R3) 7→ L2( ¯Ωg, R), ∇′ h· u , 1 q gh 3 X α=1 D′ α( q gh uα), (B.6) ∇′h× :L 2 (Ωc, R3) 7→ L2( ¯Ωg, R3), ∇′h× u , 1 q gh e1 e2 e3 D′ 1 D′ 2 D′ 3 u1 u2 u3 , (B.7) ∇h× :L2( ¯Ωg, R3) 7→ L2(Ωc, R3), ∇h× ω , 1 q gh e1 e2 e3 D1 D2 D3 ω1 ω2 ω3 , (B.8) and (B.9) △′ h: L2( ¯Ωg, R) 7→ L2( ¯Ωg, R), △′ hp = 1 q gh 3 X µ,ν=1 D′ µ q gh gµνh Dνp . 22
Fig. B.1.
Alternatively, one could have chosen {
ξα}3 α=1, { ξα}3 α=1, or { ξα}3
α=1 as coordinates and discretized
accordingly (Figure B.1).
The resulting schemes differ by each other slightly. It should be noted that the averaged version of the four skew operators, for example,
(B.10) ∇¯hp, 1 4( ∇hp + ∇hp + ∇hp + ∇hp)
is isotropic, or more precisely, invariant under the coordinate transformation (ξ1, ξ2, ξ3
) 7→ (ξ2, ξ3, ξ1). We
therefore propose the following isotropic version of the 3D scheme and elaborate with more details below. Scheme 3D-isotropic:
We first give an alternative expression for the isotropic gradient operator. To this end, we define the following operators with respect to the default coordinate {ξ1, ξ2, ξ3
}: ¯ D1,D1A2A3: L2( ¯Ωg, R) → L2(Ωc, R), (B.11) ¯ D2,D2A3A1: L2( ¯Ωg, R) → L2(Ωc, R), (B.12) ¯ D3,D3A1A2: L2( ¯Ωg, R) → L2(Ωc, R), (B.13) ¯ D′1,D ′ 1A ′ 2A ′ 3: L2(Ωc, R) → L2( ¯Ωg, R), (B.14) ¯ D′2,D ′ 2A ′ 3A ′ 1: L 2 (Ωc, R) → L2( ¯Ωg, R), (B.15) ¯ D′3,D ′ 3A ′ 1A ′ 2: L 2(Ω c, R) → L2( ¯Ωg, R). (B.16)
where D is the standard short-stencil finite difference operator, D′
is the reduced finite difference operator as defined in (3.11), and the averaging operators are defined by
(B.17) (Ag)i−1 2 = 1 2(gi+ gi−1), 1 ≤ i ≤ N, (A ′ f )i= f1 2 i = 0 1 2(fi+1 2 + fi− 1 2) 1 ≤ i ≤ N − 1 fN −1 2 i = N
With lengthy but elementary calculations, it can be shown that the isotropic gradient defined in (B.10) can be recast in terms of the operators in the default coordinate
(B.18) ∇¯hp = ( ¯D1p)e1+ ( ¯D2p)e2+ ( ¯D3p)e3∈ L2(Ωc, R3).
where eα= ∇ξα, e α=
∂x
∂ξα. We therefore recast our isotropic 3D scheme as
(B.19) ut+ ¯ω× u + ¯∇hp = ν ¯∇h× ω + f on Ωc ω = ∇¯′ h× u on ¯Ωg ¯ ∇′ h· u = 0 on ¯Ωg with ¯ω= (¯ω1, ¯ω2, ¯ω3) and ¯ ω× u ∈ L2(Ω c, R3) defined by (B.20) ω¯× u ,√gh e1 e2 e3 ¯ ω1 ω¯2 ω¯3 u1 u2 u3 , ω¯α, A1A2A3ωα∈ L2(Ωc, R), and ¯ ∇′ h· : L2(Ωc, R3) 7→ L2( ¯Ωg, R), ∇¯′h· u = 1 √g h X3 α=1 ¯ D′α(√ghuα) , (B.21) ¯ ∇′ h× : L 2 (Ωc, R3) 7→ L2( ¯Ωg, R3), ∇¯′h× u = 1 √g h e1 e2 e3 ¯ D′1 D¯ ′ 2 D¯ ′ 3 u1 u2 u3 , (B.22) ¯ ∇h× : L2( ¯Ωg, R3) 7→ L2(Ωc, R3), ∇¯h× ω = 1 √g h e1 e2 e3 ¯ D1 D¯2 D¯3 ω1 ω2 ω3 . (B.23)
It can be shown that analogue of (B.10) remains valid for the discrete operators ¯∇′ h·, ¯∇
′
h× and ¯∇h× as well.
In (B.20-B.23), the co-variant components and contra-variant components of a vector field v in L2(Ω c, R3)
or L2( ¯Ω
g, R3) are transformed to each other through
(B.24) vµ= 3 X ν=1 ghµνvν and vµ= 3 X ν=1 ghµνvν
whenever necessary. Here ghµν and gh
µν are (numerical) metric tensors with respect to the default coordinate
(ξ1, ξ2, ξ3) defined on cell centers and satisfy 3
X
γ=1
gµγh gh γν = δνµ.
The corresponding isotropic Laplacian is given by (B.25) △¯′ h: L 2 ( ¯Ωg, R) 7→ L2( ¯Ωg, R), △¯′hp = ¯∇ ′ h· ¯∇hp = √1 gh 3 X µ,ν=1 ¯ D′µ(√ghghµνD¯νp).
The 3D analogue of Lemma 3.1 remains valid for the isotropic difference operators Lemma B.1. For u ∈ L2(Ωc, R3), v ∈ L2( ¯Ωg, R3) and p ∈ L2( ¯Ωg, R), we have
1. (B.26) h u , ¯∇hp iΩc = −h ¯∇ ′ h· u , p iΩ¯g 2. (B.27) h u , ¯∇h× ω iΩc = h ¯∇ ′ h× u , ω iΩ¯g 3. (B.28) △¯′ hp = ¯∇ ′ h· ¯∇hp 4. If ψ ∈ L2( ¯Ω g, R3) with ψ × n = 0 on Γg, then (B.29) ∇¯′ h· ¯∇h× ψ ≡ 0
In addition, from (B.26) and (B.28), it is easy to see that the null space of Laplacian ¯△′
h is spanned by
grid functions which are constant along the diagonals within each cell. That is
(B.30) ker( ¯△′ h) = span{1ee, 1eo, 1oe, 1oo} where 1eei,j,k,(1 + (−1) i+j)(1 + (−1)j+k) 4 (B.31) 1eoi,j,k,(1 + (−1) i+j)(1 − (−1)j+k) 4 (B.32) 1oei,j,k,(1 − (−1) i+j)(1 + (−1)j+k) 4 (B.33) 1ooi,j,k,(1 − (−1) i+j)(1 − (−1)j+k) 4 . (B.34)
The corresponding versions for the alternative grid placement (B.2) or inhomogeneous boundary velocity and discrete compatibility conditions can be derived as in the 2D case. We omit the details.
