ORIGINAL PAPER
A simple projection method for the coupled Navier-Stokes and Darcy flows
Ming-Chih Lai1· Ming-Cheng Shiue1· Kian Chuan Ong2
Received: 26 September 2017 / Accepted: 25 September 2018 / Published online: 1 October 2018
© Springer Nature Switzerland AG 2018
Abstract
The flow through both an incompressible fluid region and a porous medium occurs in a wide range of applications encompassing industrial processes and geological phenomena. Recently, the fluid transport phenomena across the interface between the distinct regions have received increasing attention both from the mathematical and the numerical points of view.
The primary objective of the present study lies in the development of a simple and efficient method for the computations of coupled Navier-Stokes and Darcy flows with complex interface conditions. A finite difference projection method is developed within a staggered grid framework to solve the coupled system in a segregated manner using primitive variables.
Numerical simulations are carried out to demonstrate the order of convergence and its capability. The proposed method renders the versatility in solving the coupled system, and it is readily extendible to multi-physics fluid flows and turbulent flows for a broad range of applications.
Keywords Coupled Navier-Stokes and Darcy flows· Beavers-Joseph-Saffman interface conditions · Projection method · Finite difference· Staggered grid
1 Introduction
The fluid transport phenomena across the interface between an incompressible fluid region and a porous medium are often encountered in various science and engineering applications. A plethora of physicochemical phenomena, such as surface reaction, may occur at this interface.
Those coupled problems have received increasing attention both from the mathematical and the numerical points of view. The mathematical model consists of the Navier- Stokes or Stokes equations in the fluid region and the Darcy law in the porous medium. Despite the fact that
Kian Chuan Ong kcong@ncts.ntu.edu.tw Ming-Chih Lai
mclai@math.nctu.edu.tw Ming-Cheng Shiue mshiue@math.nctu.edu.tw
1 Department of Applied Mathematics, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu 300, Taiwan
2 National Center for Theoretical Sciences, No. 1, Sec. 4, Roosevelt Road, Astronomy-Mathematics Building, National Taiwan University, Taipei 10617, Taiwan
the numerical schemes of each sub-problem are well- known and well-developed, the numerical analysis of the coupled problem is still interesting to investigate. Therefore, it is prudent that attention is focused on developing a highly efficient numerical method for solving the system.
The numerical algorithm is required to compute different sets of conservation laws on different regions, on top of the coupling with complex interface conditions. Over the last decade, various approaches have been developed in the literature. For instance, based on whether these different sets of variables in different regions are solved simultaneously, finite element methods can be classified into two categories: coupled finite element methods and decoupled finite element methods. Coupled finite element methods include the coupled mixed finite element methods [1–9], coupled discontinuous Galerkin (DG) methods [10–
14], coupled hybrid methods [14], and coupled multi-grid methods [15]. Decoupled finite element methods include the domain decomposition methods [16–21], Lagrange multiplier methods [22–24], two-grid methods [25, 26], and partitioned time stepping methods [27, 28]. Besides finite element methods, other approaches such as boundary integral methods [29,30], a spectral method [31], a pseudo- spectral least squares method [32], and a finite volume
method [33] have already been studied and recently, a finite difference method based on an immersed interface method [34] was considered (as for more sophisticated variations of Navier-Stokes/Darcy or Stokes/Darcy systems that consider two-phase flow, poroelasticity, dual porosity media and so on, we refer the interested readers to see [35–42] for more details). Comprehensive reviews of the numerical algorithms in the literature revealed that the vast majority are based on finite element method instead of other discretization practice. Furthermore, the transient solutions [43] are rarely investigated as opposed to the steady-state solutions. Contrary to Stokes and Darcy system, there are also significantly less effort invested in the development on Navier-Stokes and Darcy system [12,26, 44–47].
The primary objective of the present study lies in the development of a simple and efficient method for the computations of coupled Navier-Stokes and Darcy flows with complex interface conditions. The proposed method is based on the projection method [48] in conjunction with appropriate interface boundary conditions. The fundamental discretization approach is the finite difference method based on the staggered grid framework. The proposed projection method solves the velocities, Darcy velocities, pressure, and scaled pressure in a segregated procedure.
The mathematical analysis of the numerical algorithm for coupling Stokes and Darcy flows based on the staggered grid framework is studied in [49]. In [49], the authors considered the steady-state solution. In the present article, nonstationary Navier-Stokes and Darcy flows based on the first-order projection method in the temporal variable and finite difference approximation in spatial variables are performed. Numerical simulations are conducted for the model problem with a designed analytical solution and applications in order to investigate the numerical accuracy and functionality of the method. The aforementioned methodology is flexible due to its segregated algorithm, and therefore it can be extended for a broad range of applications with relative ease. Numerical experiments suggest that the numerical scheme for coupling Navier-Stokes and Darcy flows has the same theoretical results as the case of only Navier-Stokes flow.
2 Governing equations
Consider a model problem consisting of an incompressible flow in the fluid region Fwith boundary ∂Fand a Darcy flow in the porous medium D with boundary ∂D. The bounded domains F and D ⊂ R2 are separated by an interface with the unit normal vector n pointing out of the fluid domain F. The schematic description of the model is depicted in Fig.1. In the fluid region F, the incompressible
Navier-Stokes equations for the velocity u= (u, v) and the pressure p can be read as
∂u
∂t + (u · ∇) u + ∇p = νu + F, in F, (1)
∇ · u = 0, in F, (uis given on ∂F) (2) where ν is the kinematics viscosity and F = (f, g) stands for the external forcing term. On the other hand, the porous medium can be defined by a model where the fluid and solid occupy the whole region on the macroscopic scale. The porous medium is a homogeneous continuum domain where a representative volume element is larger than the average pore size but much smaller than the length scale of the system. For modeling this saturated flow in homogeneous porous media domain D, the Darcy law can be used, i.e.,
uD= −k∇φ, in D, (3)
∇ · uD= S, in D, (uDis given on ∂D) (4) where uD = (uD, vD) is the Darcy velocity, k is the hydraulic conductivity coefficient, φ denotes the scaled pressure, and S is the external source/sink term.
To close the system (1)–(4), the governing equations must be coupled across the interface by suitable condi- tions [22] in the following. Firstly, the mass conservation across must be hold by
u· n = uD· n on . (5)
The second interface condition is the balance of normal force across as
2νn· D · n = p − φ on , (6)
where D = 12
∇u + ∇Tu
is the deformation tensor. The pressure is allowed to be discontinuous across the interface.
Lastly, since the fluid model is viscous, a condition on the tangential fluid velocity on must be given. In general, the simplest assumption of no-slippage along the interface
is invalid due to the large deviation from experimental measurements [22]. Therefore the Beavers-Joseph-Saffman (BJS) interface condition [50,51] is applied which states that the tangential component of the fluid velocity is proportional to the shear stress from the free fluid and the proportionality constant depends linearly on the square root of the permeability, i.e.,
u· τ = −2
˜k
α n· D · τ on , (7)
where τ is the unit tangent vector along the interface ,
˜k = νk and α are the positive constants. Here, the form
˜k/α has the physical meaning of friction coefficient.
Fig. 1 Schematic representation of the coupled Navier-Stokes and Darcy system
3 Finite difference discretization
To solve the coupled Navier-Stokes and Darcy flows with complex interface conditions, an underlying finite difference numerical scheme based on staggered grid is first described in this section. For the realistic applications with curved interface, the original curve can be approximated by a union of horizontal and vertical line segment interfaces.
For simplicity, we assume to be a rectangular domain [0, Lx] ×
−Ly, Ly
, where F = [0, Lx] × 0, Ly
, and D = [0, Lx] ×
−Ly,0
in which Lx, Ly are positive constants. Therefore, the fluid-porous interface is a horizontal straight line at y= 0. Consequently, the interface outward normal vector is simply written as n = (0, −1)T and the interface conditions on are simplified to
v= vD= v (8)
2ν∂v
∂y = p − φ (9)
u=
˜k α
∂u
∂y +∂v
∂x
. (10)
Similarly, in case of a vertical straight interface, the outward normal vector is n= (1, 0)T, then the interface conditions on are given by
u= uD= u (11)
Fig. 2 Finite difference staggered grid
2ν∂u
∂x = p − φ (12)
v=
˜k α
∂v
∂x +∂u
∂y
. (13)
For the clarity of presentation, the case with just a horizontal straight interface will be considered in the remainder of this section.
Let xi = (i − 1/2)x, yj = (j − 1/2) y where the mesh sizes x and y are equal to Lx/M and Ly/N, respectively, with M and N as positive integers.
The schematic representation of the finite difference discretization within the staggered grid framework is depicted in Fig.2. As illustrated in Fig.2, all the primitive variables are defined at different locations, with pressure p and scaled pressure φ defined at the cell centers, while the velocity components u, v, uD, and vD defined at the center of the cell faces. It should be noted that due to the staggered grid arrangement, the interface separating the fluid flow region and the porous medium is placed to be coinciding with the vertical velocity component v = v = vD. Therefore, to satisfy the interface conditions, an additional equation, i.e., Eq.9, is required to solve for v.
Before discussing the spatial discretization for the governing equations and the interface conditions, an
Fig. 3 Decoupled multi-step temporal integration scheme
extended projection method in light of the original projection method [48] is developed to calculate the primitive variables u, p, uD, φ at time level n+ 1 on the staggered grid. The projection method begins from a temporal discretization of the Navier-Stokes equations and Darcy law. Using the Euler implicit time stepping results in
un+1− un
t + un· ∇
un−νun+1+∇pn+1= Fn+1, (14)
∇ · un+1= 0, (15)
unD+1= −k∇φn+1, (16)
∇ · unD+1= Sn+1, (17)
2ν∂vn+1
∂y | =
pn+1− φn+1
|. (18)
The system of discretized equations is split into two sub- steps, i.e.,
u∗− un
t + un· ∇
un− νu∗= Fn+1, (19)
u∗=
˜k α
∂u∗
∂y +∂vn
∂x
on . (20)
and⎧
⎪⎪
⎪⎪
⎪⎪
⎨
⎪⎪
⎪⎪
⎪⎪
⎩
un+1−u∗
t + ∇pn+1= 0, (21a)
∇ · un+1= 0, (21b)
unD+1= −k∇φn+1, (21c)
∇ · unD+1= Sn+1, (21d)
2ν∂v
n+1
∂y | =
pn+1− φn+1
|. (21e)
By taking the divergence of Eq. 21a and 21c, and invoking Eq.21band21d, system of Eq.21can be further reduced to
⎧⎪
⎨
⎪⎩
pn+1= t1∇ · u∗, (22a)
−kφn+1= Sn+1, (22b)
2ν∂v
n+1
∂y |=
pn+1− φn+1
|. (22c)
Equation22ato22care to be solved by coupling each other with appropriate boundary conditions. Subsequently, un+1 and unD+1 are calculated through the following correction sub-step
un+1= u∗− t∇pn+1, (23a)
unD+1= −k∇φn+1. (23b)
The time-marching algorithm consists of three sub-steps for each time step: (i) the first sub-step accounts for viscous effects and computes the provisional velocity in the incompressible flow region. (ii) The second sub-step accounts for the incompressibility, and coupling of p, φ, v via interface conditions. (iii) Finally, the provisional velocity is corrected, and the velocity in the Darcy law region is updated. The present version of the projection method for Navier-Stokes equations shares the same spirit of SIMPLE algorithm [52,53] but differs in some aspects where the latter algorithm is popularly used in some commercial CFD softwares. In the aforementioned sub- step (i), SIMPLE algorithm treats the nonlinear advection term semi-implicitly as (un · ∇)u∗ and an explicit term
∇pn is added in the momentum equation. Thus, in step (2), the SIMPLE algorithm solves the incremental pressure Poisson equation rather than solving the pressure Poisson equation used in our present scheme. Since we have to solve the fluid pressure p and Darcy’s scaled pressure φ together, SIMPLE algorithm for Navier- Stokes solver is not ultilized in the present study. For completeness of the numerical formulation, the spatial discretizations in both flow domains and their coupling at the interface are articulated in the remainder of this section.
Let ui+1/2,j, vi,j+1/2, and pi,j be denoted the dis- crete approximations of the flow velocity u
xi+1/2, yj
,
Table 1 Grid refinement analysis of the fluid horizontal velocity u and Darcy horizontal velocity uDat T = 1
M× 2N uexact− u∞ Rate uD,exact− uD∞ Rate
32× 64 8.484×−03 – 5.285×−03 –
64× 128 4.963×−03 0.77 2.673×−03 0.983
128× 256 2.694×−03 0.88 1.298×−03 1.04
256× 512 1.398×−03 0.95 6.148×−04 1.08
Time step is set as t= h/2 with h = 1/M
Table 2 Grid refinement analysis of the fluid vertical velocity v and Darcy vertical velocity vDat T = 1
M× 2N vexact− v∞ Rate vD,exact− vD∞ Rate
32× 64 8.764× 10−3 – 7.457× 10−3 –
64× 128 4.604× 10−3 0.93 4.200× 10−3 0.83
128× 256 2.317× 10−3 0.99 2.221× 10−3 0.92
256× 512 1.163× 10−3 0.99 1.134× 10−3 0.97
Time step is set as t= h/2 with h = 1/M
v
xi, yj+1/2
, and the pressure p (xi, yi), respectively. In the similar context, let (uD)i+1/2,j, (vD)i,j+1/2, and φi,j
denote the discrete approximations of the flow velocity uD
xi+1/2, yj
, vD
xi, yj+1/2
, and the scaled pressure φ (xi, yi), respectively. First, using the central differencing scheme, the momentum equation, i.e., Eq. 19, for provi- sional velocities u∗and v∗are discretized as
u∗i+1/2,j− uni+1/2,j
t +uni+1/2,j
uni+3/2,j− uni−1/2,j 2x
+vni+1/2,j
uni+1/2,j+1− uni+1/2,j−1 2y
−ν
u∗i+3/2,j− 2u∗i+1/2,j+ u∗i−1/2,j (x)2
+ u∗i+1/2,j+1− 2u∗i+1/2,j+ u∗i+1/2,j−1 (y)2
= fi+1/2,jn+1 , (24)
for (i, j )∈ Su
=
(i, j )∈ Z2|1 ≤ i ≤ M − 1, 1 ≤ j ≤ N , and
v∗i,j+1/2− vi,jn +1/2
t +vi,jn +1/2
vni,j+3/2− vi,jn −1/2
2y
+uni,j+1/2
vni+1,j+1/2− vi−1,j+1/2n 2x
−ν
vi+1,j+1/2∗ − 2vi,j+1/2∗ + v∗i−1,j+1/2
(x)2
+ vi,j∗ +3/2− 2v∗i,j+1/2+ vi,j∗ −1/2
(y)2
= gni,j+1+1/2, (25)
for (i, j )∈ Sv
=
(i, j )∈ Z2|1 ≤ i ≤ M, 1 ≤ j ≤ N − 1 , respectively. The non-linear advection terms are dis- cetized explicitly while the diffusion terms are treated implicitly to assist stability. In this configuration, the
coupling interface is located along j = 1/2 and is coincided with the vertical velocity component (v)i,1/2. Therefore, additional treatments are required to impose on Eqs. 24 and25 in order to include the coupling interface conditions.
Considering Eq.24, the u∗component at j = 1 adjacent to the interface is given by
u∗i+1/2,1− uni+1/2,1
t +uni+1/2,1
uni+3/2,1− uni−1/2,1
2x
+vin+1/2,1
uni+1/2,2− uni+1/2,0
2y
−ν
u∗i+3/2,1− 2u∗i+1/2,1+ u∗i−1/2,1 (x)2
+ u∗i+1/2,2− 2u∗i+1/2,1+ u∗i+1/2,0 (y)2
= fi+1/2,1n+1 ,
i = 1, · · · , M − 1, (26)
where the ghost nodes u∗i+1/2,0and uni+1/2,0are required to be evaluated. The Beavers-Joseph-Saffman condition (10) is approximated and discretized at j = 1/2 as
u∗i+1/2,1+ u∗i+1/2,0
2 =
˜k α
u∗i+1/2,1− u∗i+1/2,0
y
+
˜k α
(v)ni+1,1/2− (v)ni,1/2
x
. (27)
Table 3 Grid refinement analysis of the fluid vertical velocity at the interface vat T = 1
M× 2N v,exact− v∞ Rate
32× 64 9.372× 10−03 –
64× 128 4.830× 10−03−03 0.96
128× 256 2.402× 10−03−03 1.01
256× 512 1.189× 10−03−03 1.02
Time step is set as t= h/2 with h = 1/M
Table 4 Grid refinement analysis of the pressure p and scaled pressure φ at T = 1
M× 2N pexact− p∞ Rate φexact− φ∞ Rate
32× 64 3.149× 10−1 - 1.590× 10−3 -
64× 128 2.228× 10−1 0.50 7.901× 10−4 1.01
128× 256 1.562× 10−1 0.51 3.778× 10−4 1.06
256× 512 1.093× 10−1 0.52 1.779× 10−4 1.09
M× 2N pexact− p2 Rate φexact− φ2 Rate
32× 64 9.434× 10−2 – 3.746× 10−4 –
64× 128 5.322× 10−2 0.83 1.641× 10−4 1.19
128× 256 2.935× 10−2 0.86 7.165× 10−5 1.20
256× 512 1.601− ×10−2 0.87 3.150× 10−5 1.19
Time step is set as t= h/2 with h = 1/M
This yields the ghost node u∗i+1/2,0
u∗i+1/2,0 =
2
˜k/αy − 1 2√
k/αy+ 1
u∗i+1/2,1
+ 2
˜k/αx 2
˜k/αy + 1
(v)ni+1,1/2− (v)ni,1/2 . (28) Similarly, the ghost node uni+1/2,0 can be approximated in the identical manner.
Now, considering Eq. 25 on the other hand, the v∗ component at j = 3/2 adjacent to the interface is given by v∗i,3/2− vi,3/2n
t +uni,3/2
vni+1,3/2− vin−1,3/2
2x
+vni,3/2
vni,5/2− (v)ni,1/2 2y
−ν
vi∗+1,3/2− 2v∗i,3/2+ vi∗−1,3/2
(x)2 +vi,5/2∗ − 2vi,3/2∗ + (v)∗i,1/2
(y)2
= gi,3/2n+1,
i= 1, · · · , M. (29)
The (v)ni,1/2prevailing at previous time step n is employed for the explicit discretization of advection term. However, the (v)∗i,1/2 term in the implicit diffusion term can be
obtained from Eq.9and approximated by first-order one- sided difference as
2ν
v∗i,3/2− (v)∗i,1/2
y
= pi,1n −φi,0n , i= 1, · · · , M. (30)
Next, for the coupled system of Eq. 22, the spatial discretization yields
pin+1,j+1 − 2pni,j+1+ pin−1,j+1
(x)2 +pni,j+1+1− 2pni,j+1+ pi,jn+1−1
(y)2
= 1
t
u∗i+1/2,j − u∗i−1/2,j
x +v∗i,j+1/2− vi,j∗ −1/2
y
, (31) for (i, j )∈ S
=
(i, j )∈ Z2|1 ≤ i ≤ M, 1 ≤ j ≤ N ,
φin+1,j+1−2φni,j+1+φni−1,j+1 (x)2
+φi,j+1n+1−2φ(y)i,jn+12+φn+1i,j−1 = −Si,jn+1k , (32) for (i, j )∈ D
=
(i, j )∈ Z2|1 ≤ i ≤ M, −N + 1 ≤ j ≤ 0 ,
Table 5 Grid refinement analysis of the scaled pressure gradient φxand φyat T = 1
M× 2N φx,exact− φx2 Rate φy,exact− φy2 Rate
32× 64 1.253× 10−3 – 1.154× 10−3 –
64× 128 5.748× 10−4 1.12 5.466× 10−4 1.08
128× 256 2.618× 10−4 1.13 2.535× 10−4 1.11
256× 512 1.200× 10−4 1.13 1.174× 10−4 1.11
Time step is set as t= h/2 with h = 1/M
Table 6 Grid refinement analysis of the fluid horizontal velocity u and Darcy horizontal velocity uDat T = 1
M× 2N uexact− u∞ Rate uD,exact− uD∞ Rate
32× 64 4.145× 10−4 – 1.584× 10−4 –
64× 128 1.080× 10−4 1.94 3.446× 10−5 2.20
128× 256 2.754× 10−5 1.97 7.911× 10−6 2.12
256× 512 6.951× 10−6 1.99 2.133× 10−6 1.89
Time step is set as t= h2/2 with h= 1/M
2ν
vi,3/2n+1−(v)ni,1/2+1
y
=pni,1+1− φni,0+1, i=1, · · · , M.
(33) Note that at the current stage, vi,3/2n+1 is unknown. In this case, the correction sub-step Eq.23must be invoked, i.e.,
vni,3/2+1 = vi,3/2∗ − t
pni,2+1− pi,1n+1
y
,
i= 1, · · · , M. (34)
Hence, Eq.33is written as 2νvi,3/2∗ − (v)ni,1/2+1
y = pni,1+1− φi,0n+1
+2νt (y)2
pi,2n+1− pni,1+1
,
i = 1, · · · , M. (35)
To impose the coupling interface conditions and bound- ary conditions, additional treatments are required for Eqs.31and32. Note that along the interface j = 1/2 in i direction, where all the variables of pn+1are located at the cell center along j = 1, Eq.31is coupled with (v)ni,1/2+1 as
pi+1,1n+1 − 2pni,1+1+ pni−1,1+1
(x)2 +pi,2n+1− 2pni,1+1+ pi,0n+1
(y)2
= 1
t
u∗i+1/2,1− u∗i−1/2,1
x + vi,3/2∗ − (v)ni,1/2+1
y
. (36)
Here, the ghost value pni,0+1 can be obtained easily by approximating zero Neumann boundary condition on the interface . On the other hand, φn+1 located at the cell center along j = 0 adjacent to the interface is given by
φni+1,0+1 − 2φni,0+1+ φin−1,0+1
(x)2 +φi,1n+1− 2φi,0n+1+ φni,+1−1
(y)2 = −Si,0n+1 k ,
(37)
where the ghost value φi,1n+1 is evaluated by direct substitution of the gradient at the interface as
φni,1+1− φi,0n+1
y = −1
k(v)ni,1/2+1 , (38)
which is also coupled with (v)ni,1/2+1.
Finally, in the correction sub-step, un+1 and vn+1 in incompressible fluid region are computed as
uni+1/2,j+1 = u∗i+1/2,j− t
pni+1,j+1 − pni,j+1
x
, (39)
for (i, j )∈ Su,
vni,j+1+1/2= v∗i,j+1/2− t
pni,j+1+1− pni,j+1
y
, (40)
for (i, j ) ∈ Sv. The correction sub-step ensures that un+1 and vn+1in the incompressible fluid flow region satisfy the divergence-free kinematic constraint. Concurrently, unD+1 and vnD+1in the porous medium region are also updated, i.e.
(uD)ni+1/2,j+1 = −k
φin+1,j+1 − φni,j+1
x
, (41)
Table 7 Grid refinement analysis of the fluid vertical velocity v and Darcy vertical velocity vDat T = 1
M× 2N vexact− v∞ Rate vD,exact− vD∞ Rate
32× 64 3.287× 10−4 – 2.571× 10−4 –
64× 128 8.396× 10−5 1.97 6.966× 10−5 1.88
128× 256 2.114× 10−5 1.99 1.788× 10−5 1.96
256× 512 5.275× 10−6 2.00 4.398× 10−6 2.02
Time step is set as t= h2/2 with h= 1/M
Table 8 Grid refinement analysis of the fluid vertical velocity at the interface vat T = 1
M× 2N v,exact− v∞ Rate
32× 64 3.699× 10−4 –
64× 128 9.458× 10−5 1.97
128× 256 2.378× 10−5 1.99
256× 512 5.900× 10−6 2.01
Time step is set as t= h2/2 with h= 1/M
for (i, j )∈ Du
=
(i, j )∈ Z2|1 ≤ i ≤ M − 1, −N + 1 ≤ j ≤ 0 ,
(vD)ni,j+1+1/2= −k
φi,jn+1+1− φni,j+1
y
, (42)
for (i, j )∈ Dv
=
(i, j )∈ Z2|1 ≤ i ≤ M, −N + 2 ≤ j ≤ 0 .
The aforementioned discretization framework still holds in case of straight vertical interface n = (1, 0)T by switching the “role” between the horizontal u and vertical v velocity components adjacent to the interface. The vertical interface separating the fluid flow region and the porous medium is placed to be coinciding with the horizontal velocity component u = u = uD. Therefore, Eq. 12 is required to solve for u. The discretization for u follows exactly Eq.33. The Beavers-Joseph-Saffman condition (13) is discretized on and imposed as an implicit boundary condition for v∗equations similar to Eq.26.
In the majority of real-world situations, the fluid velocity magnitude in the incompressible fluid region is marginally higher than the fluid velocity magnitude in the porous medium. Since the coupled system evolves on different time scales, the multi-step temporal integration scheme [33] can
be adopted. By introducing two temporal grids, the coupled system of Eqs.31and33is decoupled from Eq.32, i.e., 2νv∗i,3/2− (v)ni,1/2+1
y = pni,1+1− φmi,0
+2νt (y)2
pni,2+1− pni,1+1
, (43) where the time level n represents the fine time step and time level m represents the coarse time step. In this case, the incompressible fluid flow solutions are computed with a fine time step t and the porous medium solutions are computed with a coarse time step mt where m≥ 1 is the ratio between the coarse and fine time steps. Figure3shows the two temporal grids schematically. The coupled system of Eqs. 31to33is only solved once at every coarse time level mt. Note that the ratio m is a problem-dependent parameter which depends on the local characteristic of the flows. It is possible that the ratio m can be set dynamically and adaptively throughout the simulation.
4 Numerical results
4.1 Accuracy check of the present method
To validate and to check the numerical convergence of the present numerical method, a grid refinement analysis is conducted by constructing an analytic solution. Consider the Navier-Stokes flow in 0≤ x ≤ 1, 1 ≤ y ≤ 2 and the Darcy flow in 0≤ x ≤ 1, 0 ≤ y ≤ 1, the exact solution is given by u= e−t
(y− 1)2+ x(y − 1) + 3x − 1
, (44)
v= e−t
x(x− 1) − 0.5(y − 1)2− 3y + 1
, (45)
p= e−t(2x+ y − 1) , (46)
uD= e−t[(x− 1)(y − 1) + x(y − 1) − 2] , (47)
Table 9 Grid refinement analysis of the pressure p and scaled pressure φ at T = 1
M× 2N pexact− p∞ Rate φexact− φ∞ Rate
32× 64 7.446× 10−2 – 5.155× 10−5 –
64× 128 3.693× 10−2 1.01 1.217× 10−5 2.08
128× 256 1.835× 10−2 1.01 2.832× 10−6 2.11
256× 512 9.144× 10−3 1.01 7.792× 10−7 1.86
M× 2N pexact− p2 Rate φexact− φ2 Rate
32× 64 3.229× 10−2 – 2.028× 10−5 –
64× 128 1.553× 10−2 1.06 4.892× 10−6 2.05
128× 256 7.613× 10−3 1.03 1.187× 10−6 2.04
256× 512 3.768× 10−3 1.01 2.914× 10−7 2.03
Time step is set as t= h2/2 with h= 1/M
Table 10 Grid refinement analysis of the scaled pressure gradient φxand φyat T = 1
M× 2N φx,exact− φx2 Rate φy,exact− φy2 Rate
32× 64 3.321× 10−5 – 3.025× 10−5 –
64× 128 6.781× 10−6 2.29 6.426× 10−6 2.24
128× 256 1.212× 10−6 2.48 1.168× 10−6 2.46
256× 512 2.746× 10−7 2.14 2.672× 10−7 2.13
Time step is set as t= h2/2 with h= 1/M
vD= e−t
x(x− 1) − (y − 1)2− 2
, (48)
φ= e−t
x(1− x)(y − 1) +(y− 1)3
3 + 2x + 2y + 4
, (49) with parameters ν = 1, k = 1, and α = 1. In addition, the exact solution also satisfies the interface conditions Eqs.8 to10. It should be noted that along the interface y = 1, the aforementioned solution has discontinuous pressure and non-zero velocity. The convergence tests are carried out for all the primary variables using four levels of grid-resolution starting from M× 2N = 32 × 64 with uniform mesh size
x = y = h, time step t = h/2 and t = h2/2.
The computations are performed with time-marching up to time T = 1.0. For the case with t = h/2, the errors and convergent rates in L∞-norm for u and uD, v and vD, v, p and φ are summarized in Tables1,2,3, and4, respectively.
In addition, the errors and convergent rates in L2-norm for p, φ, φx, and φyare summarized in Tables4and5.
In a same manner, the corresponding errors and convergent rates in L∞-norm for the case with t = h2/2 are summarized in Tables6,7,8, and 9. In addition, the corresponding errors and convergent rates in L2-norm are summarized in Tables9and10. For the case t proportional to h, the fluid velocity u and v, Darcy velocities uDand vD, vertical velocity at the interface v, and scaled pressure φ are first-order accurate in L∞-norm, whereas the pressure p is only half-order accurate in L∞-norm. Pressure p, scaled
Fig. 4 Initial and boundary conditions which represent the percolation of waters of hydrological basins through rocks and sand
pressure φ, and its gradients φx, φy are first-order accurate in L2-norm. This first-order accuracy for the velocity is due to the choice of first-order temporal discretization of the projection method. However, for the case of t proportional to h2, the fluid velocities u and v, Darcy velocities uD
and vD, vertical velocity at the interface v, as well as scaled pressure φ are second-order accurate in L∞-norm, whereas the pressure p is first-order accurate in L∞-norm and L2-norm. Scaled pressure φ and its gradients φx, φyare second-order accurate in L2-norm.
4.2 Percolation of fluid through porous medium A physical example is used to demonstrate the capability of the developed numerical algorithm. Consider the domains
F= [0, 8]×[1, 2] and D= [0, 8]×[0, 1], the initial and boundary conditions are depicted in Fig.4. This canonical case may represent the percolation of fluid through porous medium [33]. The inlet boundary condition at the left boundary of the incompressible fluid flow region ∂F = {0} × [1, 2] is given by
u= U
1− 4 (y − 1.5)2
, v= 0 (50)
where U = 10 in this case. The outflow condition is prescribed at the right boundary and no-slip boundary is prescribed at the top boundary of the incompressible fluid flow region ∂F = {8} × [1, 2]. For the remaining boundaries of the porous medium domain, the no-flux boundary condition ∇φ · n = 0 is imposed. The initial velocities in F and D are zero. The computation is performed with time-marching, starting from the initial condition to the steady-state using mesh size of x= y = h = 1/60 and time step of t = h2/2. The steady-state solution is justified by assessing the convergence of the primitive variables.
Figure5depicts the comparison of the flow field using different hydraulic conductivity coefficients k = 1.0 (Fig.5a) and k = 0.01 (Fig.5b). In both cases, the fluid from the incompressible flow region F flows into the porous medium region D until the mass flux across the interface has reached the steady-state. Subsequently, the
Fig. 5 Comparison of velocity vectors with relative magnitude between the case a k= 1.0 and the case b k = 0.01
boundary layer along the top boundary due to the no-slip condition, and along the interface in the incompressible flow region Fdue to the BJS interface condition, continue to develop toward steady-state. The primary discrepancy between the two cases is the magnitude of the mass flux across the interface . Figure6illustrates the comparison of normalized velocity v/Uprofiles along the interface using different k values starting from k = 0.0001 to k = 1.0. Due to the conservation of mass, the sum of the mass flow across the interface must be zero. Table11justifies that the projection method enforces the discrete divergence- free conditions, and the mass balance error
vdh across the interface is up to the machine accuracy. Since k
Fig. 6 Normalized velocity v/Uprofiles along the interface using different k values starting from k= 0.0001 to k = 1.0
is proportional to permeability, the magnitude of normalized velocity flows from the incompressible fluid region into the porous medium region increases as the value of k increases.
To illustrate this, Fig. 7 depicts the mass flux across the interface normalized by the inlet mass flux. It is noted that the normalized mass flux increases significantly with permeability proportional to k.
The accuracy of the multi-step temporal scheme is eval- uated herein using the preceding case study with identical boundary conditions and initial conditions. Table12sum- marizes the comparison of normalized mass flux across the interface using the coarse to fine time step ratio m = 1, 2, 3, 5 with k = 1.0, 0.001, 0.0001. The predicted nor- malized mass fluxes for m = 2, 3, 5 are consistent with the fully coupled case m = 1. By decoupling the system of Eqs.31to33and only solve once at every coarse time level mt, the computational cost can be reduced without degrading the accuracy.
Table 11 Mass balance error
vdh across the interface
k Mass balance error on
k= 1.0 8.750× 10−13
k= 0.1 1.145× 10−13
k= 0.01 4.297× 10−14
k= 0.001 6.886× 10−15
k= 0.0001 8.058× 10−16
