# Finite Element Method Using a Characteristic-Based Split for Numerical Simulation of a Carbonate Fracture-Cave Reservoir

Authors: Liehui Zhang; Yuhui Zhou; Lei Zhao; Deliang Zhang
Journal: Journal of Chemistry
 (2015)

---

## Abstract

Fracture-cave carbonate reservoirs occur widely in source rocks and are prospects for exploitation worldwide. However, the presence of massive caves and multiscale fracture systems results in extremely complex fluid flow patterns. Therefore, in this paper, a discrete network model for fracture-cave reservoirs was established to study fluid flow characteristics and pressure distributions in complex flow regimes. In this study, the cave system was treated as a free-flow region, and the fluid flow in fracture systems followed the Navier-Stokes and Darcy equations, respectively. After discrete modeling, the Galerkin finite element method was used for numerical calculation of the single-phase free flow; the method maintains a high-precision result with low grid orientations during the simulation. In addition, because only one linear equation requires solving at each step, the solution is obtained quickly. Moreover, based on the proposed discrete media network model of fracture-cave reservoirs and the finite element numerical calculation method, a corresponding simulator was also developed. The finite element numerical simulation method based on the characteristic-based split (CBS) algorithm has proven to be applicable to complex flow problems in fracture-cave reservoirs.


---

## Body

## 1. Introduction

Fracture-cave reservoirs play a vital role in the petroleum industry, as these reservoirs contain large volumes of oil and are widely distributed in areas such as Saudi Arabia, Iran, Iraq, Mexico, Oman, and Syria [1]. Therefore, the development of reliable analytical solutions that consider the geometrical properties, storage capacities, and flow properties of porous media in fracture-cave reservoirs is of utmost importance. Many of the world’s giant fields contain naturally fractured and vuggy carbonate reservoirs, and the pore systems in these reservoirs are complex as a result of the characteristics of carbonate rocks, which are particularly sensitive to postdepositional diagenesis that can involve dissolution, dolomitization, and fracturing.A great deal of research and model development has been devoted to flow behavior in fractured vuggy reservoirs; the focus of such research, which has been conducted worldwide, can be categorized into efforts to develop three types of models: continuum models, equivalent continuum models, and discrete medium models. Discrete medium models describe the heterogeneity of local-scale reservoirs better than the other models do; therefore, these models describe actual reservoir characteristics more accurately. Three numerical methods are commonly used to solve flow models: the finite difference method (FDM), the finite element method (FEM), and the boundary element method (BEM). The FEM has been widely used in the field of solid mechanics and is now also widely used in the petroleum research, thus providing an important approach to numerical simulation of problems involving discrete fractured networks [2, 3].Kim and Deo [4] proposed a finite element discrete fracture model to describe multiphase flow in porous media. Their study showed that oil recovery is determined by a complex interplay of absolute matrix permeability, permeability contrasts, and flow rates; fracture capillary pressure also plays a significant role in determining water penetration into the matrix [5, 6].Liu et al. [7] discussed the discretization of coupled finite element stress equations used in commercial reservoir simulations and compared the stress predictions of these equations with the results of a commonly used commercial stress program. Sheng et al. [8] used an extended FEM to investigate an integrated and improved shale-gas transport model and characterized the main flow mechanisms and the discrete fracture network. Wei et al. [9] presented a second-order finite element algorithm that can be used to simulate the frequency-dependent laterolog response of axially symmetric, invaded, and anisotropic formations. Qin et al. [10] developed a flow calculation method based on the Stokes-Brinkman model and the discrete fracture network method, which can be used for accurate and efficient upscaling of processes in naturally fractured carbonate karst reservoirs.Gulbransen et al. [11] developed a multiscale mixed finite element (MsMFE) method that can be used to examine vuggy and naturally fractured reservoirs; the MsMFE is a first step toward developing a uniform multiscale multiphysics framework. In addition, a comparison of solutions of the MsMFE method and the fine-scale Stokes-Brinkman model for two different test cases, including both short- and long-range fractures, was also presented; the results demonstrate how fine-scale flow in fracture networks can be represented within a coarse-scale Darcy flow model using multiscale elements, with the solution computed using the Stokes-Brinkman equations. The results demonstrate that the MsMFE method can simulate flow in vuggy and naturally fractured reservoirs using a highly detailed geocellular approach.Jouini and Vega [12] estimated the elastic properties of core plug samples from a Middle East carbonate reservoir and computed tomography scanner images by X-ray. Using the FEM, they simulated the elastic properties of the reservoir rocks. Using Hertzian’s and Hookean’s contact models, the discrete element method can be used to investigate and predict the elastic properties of reservoir rocks.Jackson et al. [13] developed approaches to simulate fluid flow using the pillar-grid concept. Results showed significant improvements in the representation of multiscale geological heterogeneities and predictions of flow through those heterogeneities. In their research, which was based on more than 20 years of development of innovative numerical methods in ocean modeling, Jackson et al. refined and modified their approach to address the unique challenges associated with reservoir simulation. Chen [14] presented a new finite element in ABAQUS that incorporates the extended FEM (XFEM) in the solution of hydraulic fracture problems.Cinar et al. [15] developed a numerical model to analyze transient pressure data from naturally fractured reservoirs; their model calculated well pressures in a single well in a naturally fractured reservoir. The characteristic properties of naturally fractured reservoirs, including fracture orientation, fracture length, fracture density, and fracture distribution, were also simulated in this model. The simulation results showed that the model can predict the transient pressure response in complex naturally fractured reservoirs; in addition, it can be used to match pressure data and diagnose characteristic properties of naturally fractured reservoirs. In addition to the studies reported above, a great deal of other research related to the FEM has been carried out [16–22].In conclusion, the FEM has been widely studied and applied in petroleum engineering, especially for the analysis of flow in sandstone reservoirs and naturally fractured reservoirs, as well as in shale-gas reservoirs; however, the technique has rarely been applied to fracture-cave reservoirs. Because of the lack of a systematic theory and appropriate modeling techniques, research on theories and techniques of numerical simulation using discrete medium models to model giant fracture-cave reservoirs is still in an exploratory stage, with existing studies concentrating mostly on fractured media. The characterization of vuggy media with extremely complex internal geometries using theoretical models remains difficult. The characteristic-based split (CBS) algorithm has been used to solve problems in solid mechanics and fluid mechanics, but until now the algorithm has not been applied to a discrete media model of fracture-cave reservoirs. In this paper, we applied the FEM to the study of fracture-cave reservoirs based on a CBS discrete medium model.


## 2. Discrete Fracture-Cave Network Model Concept and Mathematical Formulation

As observed in numerous studies on the field of carbonate formation, three types of porosity are typically present in natural fracture-cave reservoirs; these are matrix porosity, fracture porosity, and cave porosity [23–26]. Fractures and caves are irregularly distributed and vary in size from microscopic to macroscopic. Figure 1 presents the conceptualized discrete fracture-cave network model. The isolated caves are connected via a discrete fracture network. Flow within the matrix and in macroscopic fracture systems follows the Darcy law and in cave systems follows a free-flow region.Figure 1
Discrete fracture-cave network model and simplified cave model.In this paper, we assume that fluid flow is isothermal, single phase, and slightly compressible and that the fluid viscosity is constant. For a laminar viscous fluid flowing through an open domain (e.g., a cave system), the flow is governed by the continuity and momentum equations, which can be expressed in vector form as [27, 28](1)
∂
ρ
∂
t
+
∇
·
ρ
u
=
0
,
ρ
∂
u
∂
t
+
ρ
u
·
∇
u
=
ρ
f
+
∇
·
σ
,where ρ is the fluid density (kg/m3), u is the fluid velocity (m/s), ∇ is a gradient operator, f is the body force per unit mass (m/s2), and σ is the fluid stress tensor (Pa). σ is given by(2)
σ
=
-
p
f
I
+
2
μ
ε
,where p
f is the pressure in the fracture system (Pa), I is a unit tensor, μ is the fluid viscosity (Pa s), and ε, the strain rate, is(3)
ε
=
1
2
∇
u
+
u
∇
.For simplicity, in this study we treat the fluid as a Newtonian fluid. The density of a slightly compressible fluid can be assumed to be nearly constant (i.e., incompressible). Therefore, the governing fluid motion equations can be reduced to the Navier-Stokes equations, which are(4)
∇
·
u
=
0
,
(5)
ρ
∂
u
∂
t
+
ρ
u
·
∇
u
=
-
∇
P
f
+
μ
∇
2
u
.The boundary conditions for (2) and (3) can be classified into two types: Dirichlet (velocity) and Neumann (traction) [29, 30], which can be expressed, respectively, as (6)
u
=
u
D
,
o
n
Γ
D
,
(7)
n
·
σ
=
t
N
,
o
n
Γ
N
,where u
D is the specific velocity at the Dirichlet boundary Γ
D, n is an outward unit normal vector at the boundary, and t
N is the specific traction.


## 3. Derivation of the Finite Element Method Based on CBS

### 3.1. CBS Procedure

The theory of the CBS algorithm is based on the split- and characteristic-Galerkin procedures, which are extended in the FEM. The CBS algorithm involves three indispensable steps: Step1 is to compute the intermediate momentum variable; Step 2 is to calculate the density or pressure field; and Step 3 involves combining the introduction of the density or pressure fields to obtain the correct momentum variables.Step 1.
The solution forΔ
U
~
∗ in matrix form is(8)
Δ
U
~
∗
=
-
M
-
1
Δ
t
C
U
~
+
K
u
~
-
f
-
Δ
t
K
u
U
~
+
f
s
n
,where the quantities identified by ~ symbols indicate nodal values, and(9)
M
=
∫
Ω
N
T
N
d
Ω
,
C
=
∫
Ω
N
T
∇
u
N
d
Ω
,
K
=
∫
Ω
B
T
μ
I
0
-
2
3
m
m
T
B
d
Ω
=
∫
Ω
μ
∂
N
T
∂
x
j
∂
N
∂
x
j
d
Ω
,
f
=
∫
Ω
N
T
ρ
g
d
Ω
+
∫
Γ
N
T
t
d
d
Γ
,where B is a matrix, N is a deformation function, m
=
[
1,1
,
0
]
T, and I
0
=
d
i
a
g
(
2,2
,
1
). After integration by parts, the expressions for K
u and f
s are (10)
K
u
=
-
1
2
∫
Ω
∇
T
u
N
T
∇
T
u
N
d
Ω
,
f
s
=
-
1
2
∫
Ω
∇
T
u
N
T
ρ
g
d
Ω
.Step 2.
The value ofΔ
p
~ can be easily determined using the expression(11)
M
P
+
Δ
t
2
θ
1
θ
2
H
Δ
P
~
=
Δ
t
G
U
~
n
+
θ
1
G
U
~
∗
-
Δ
t
θ
1
H
P
~
n
-
f
p
.The matrices are(12)
H
=
∫
Ω
∇
N
T
∇
N
d
Ω
,
M
p
=
∫
Ω
N
T
1
c
2
n
N
d
Ω
,
G
=
∫
Ω
∇
N
T
N
d
Ω
,
f
p
=
Δ
t
∫
Γ
N
T
n
T
U
~
n
+
θ
1
Δ
U
~
∗
-
Δ
t
∇
p
n
+
θ
2
d
Γ
.Step 3.
The value ofΔ
U
~ is determined from(13)
Δ
U
~
=
Δ
U
~
∗
-
M
-
1
Δ
t
G
T
p
~
n
+
θ
2
Δ
p
~
+
Δ
t
2
P
p
~
n
,where(14)
P
=
∫
Ω
∇
u
N
T
∇
N
d
Ω
.


### 3.2. Derivation Process of the Free-Flow Model

From (5), we derive(15)
∂
u
∂
t
=
-
u
∂
u
∂
x
+
v
∂
u
∂
y
-
1
ρ
∂
p
∂
x
+
μ
ρ
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
,
(16)
∂
v
∂
t
=
-
u
∂
v
∂
x
+
v
∂
v
∂
y
-
1
ρ
∂
p
∂
y
+
μ
ρ
∂
2
v
∂
x
2
+
∂
2
v
∂
y
2
.From Step 1, after the split-temporal discretization, (15) becomes (omitting the third-order derivative above)(17)
Δ
u
=
u
n
+
1
-
u
n
=
Δ
t
-
u
∂
u
∂
x
+
v
∂
u
∂
y
-
1
ρ
∂
p
∂
x
+
μ
ρ
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
+
Δ
t
2
2
u
∂
∂
x
u
∂
u
∂
x
+
v
∂
u
∂
y
+
1
ρ
∂
p
∂
x
.From Step 2, after spatial discretization, we find that(18)
Δ
u
=
N
u
Δ
u
~
,
Δ
u
∗
=
N
u
Δ
u
~
∗
,
u
=
N
u
u
~
,
v
=
N
v
v
~
,
p
=
N
p
p
~
,
N
=
N
1
,
N
2
,
N
3
,
Δ
u
~
=
Δ
u
1
,
Δ
u
2
,
Δ
u
3
T
,
Δ
u
~
∗
=
Δ
u
1
∗
,
Δ
u
2
∗
,
Δ
u
3
∗
T
,
u
~
=
u
1
,
u
2
,
u
3
T
,
v
~
=
v
1
,
v
2
,
v
3
T
,
p
=
p
1
,
p
2
,
p
3
T
,
N
u
=
N
v
=
N
p
.Then,(19)
∫
Ω
N
u
Δ
u
∗
d
Ω
=
Δ
t
-
∫
Ω
N
u
u
∂
u
∂
x
+
v
∂
u
∂
y
d
Ω
m
m
m
v
-
μ
ρ
∫
Ω
∇
·
N
u
∂
u
∂
x
+
∂
u
∂
y
d
Ω
+
Δ
t
2
2
∫
Ω
∇
N
u
u
-
u
∂
u
∂
x
+
v
∂
u
∂
y
d
Ω
+
Δ
t
μ
ρ
∫
Γ
N
u
∂
u
∂
x
+
∂
u
∂
y
d
Γ
.The left side of (19) can be expressed as(20)
∫
Ω
N
u
Δ
u
∗
d
Ω
=
∫
Ω
N
u
N
1
Δ
u
1
∗
+
N
2
Δ
u
2
∗
+
N
3
Δ
u
3
∗
d
Ω
=
∫
Ω
N
u
N
1
d
Ω
∫
Ω
N
u
N
2
d
Ω
∫
Ω
N
u
N
3
d
Ω
×
Δ
u
1
∗
Δ
u
2
∗
Δ
u
3
∗
=
∫
Ω
N
1
N
1
d
Ω
∫
Ω
N
1
N
2
d
Ω
∫
Ω
N
1
N
3
d
Ω
∫
Ω
N
2
N
1
d
Ω
∫
Ω
N
2
N
2
d
Ω
∫
Ω
N
2
N
3
d
Ω
∫
Ω
N
3
N
1
d
Ω
∫
Ω
N
3
N
2
d
Ω
∫
Ω
N
3
N
3
d
Ω
×
Δ
u
1
∗
Δ
u
2
∗
Δ
u
3
∗
.The first term on the right side of (19) can be expressed as(21)
∫
Ω
N
u
u
∂
u
∂
x
+
v
∂
u
∂
y
d
Ω
=
∫
Ω
N
u
u
∂
N
1
∂
x
u
1
+
∂
N
2
∂
x
u
2
+
∂
N
3
∂
x
u
3
+
v
∂
N
1
∂
y
u
1
+
∂
N
2
∂
y
u
2
+
∂
N
3
∂
y
u
3
=
∫
Ω
N
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
u
1
+
u
∂
N
2
∂
x
+
v
∂
N
2
∂
y
u
2
+
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
u
3
d
Ω
=
∫
Ω
N
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
d
Ω
∫
Ω
N
u
u
∂
N
2
∂
x
+
v
∂
N
2
∂
y
d
Ω
∫
Ω
N
u
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
d
Ω
u
1
u
2
u
3
=
∫
Ω
N
1
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
d
Ω
∫
Ω
N
1
u
∂
N
2
∂
x
+
v
∂
N
2
∂
y
d
Ω
∫
Ω
N
1
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
d
Ω
∫
Ω
N
2
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
d
Ω
∫
Ω
N
2
u
∂
N
2
∂
x
+
v
∂
N
2
∂
y
d
Ω
∫
Ω
N
2
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
d
Ω
∫
Ω
N
3
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
d
Ω
∫
Ω
N
3
u
∂
N
2
∂
x
+
v
∂
N
2
∂
y
d
Ω
∫
Ω
N
3
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
d
Ω
u
1
u
2
u
3
,where(22)
∫
Ω
N
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
d
Ω
=
∫
Ω
N
u
∂
N
1
∂
y
N
1
u
1
+
N
2
u
2
+
N
3
u
3
∂
N
1
∂
x
m
m
m
m
m
m
m
c
m
+
N
1
v
1
+
N
2
v
2
+
N
3
v
3
∂
N
1
∂
y
∂
N
1
∂
x
d
Ω
=
∫
Ω
N
u
N
1
∂
N
1
∂
x
d
Ω
∫
Ω
N
u
N
2
∂
N
1
∂
x
d
Ω
∫
Ω
N
u
N
3
∂
N
1
∂
x
d
Ω
∫
Ω
N
u
N
1
∂
N
1
∂
y
d
Ω
∫
Ω
N
u
N
2
∂
N
1
∂
y
d
Ω
∫
Ω
N
u
N
3
∂
N
1
∂
y
d
Ω
T
u
1
u
2
u
3
v
1
v
2
v
3
.The second term on the right side of (19) can be expressed as(23)
∫
Ω
∂
N
u
∂
x
i
∂
u
∂
x
+
∂
u
∂
y
d
Ω
=
∫
Ω
∂
N
u
∂
x
i
∂
N
1
∂
x
+
∂
N
1
∂
y
u
1
+
∂
N
2
∂
x
+
∂
N
2
∂
y
u
2
+
∂
N
3
∂
x
+
∂
N
3
∂
y
u
3
d
Ω
=
∫
Ω
∂
N
u
∂
x
i
∂
N
1
∂
x
+
∂
N
1
∂
y
∂
N
u
∂
x
i
∂
N
2
∂
x
+
∂
N
2
∂
y
∂
N
u
∂
x
i
∂
N
3
∂
x
+
∂
N
3
∂
y
u
1
u
2
u
3
d
Ω
.The third term on the right side of (19) can be expressed as(24)
∫
Ω
∂
N
u
∂
x
i
N
u
u
-
u
∂
u
∂
x
+
v
∂
u
∂
y
d
Ω
=
-
∫
Ω
∂
∂
x
i
N
u
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
u
1
+
u
∂
N
2
∂
x
+
v
∂
N
2
∂
y
u
2
+
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
u
3
d
Ω
=
-
∫
Ω
∂
∂
x
i
N
1
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
⋯
∂
∂
x
i
N
1
u
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
∂
∂
x
i
N
2
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
⋯
∂
∂
x
i
N
2
u
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
∂
∂
x
i
N
3
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
⋯
∂
∂
x
i
N
3
u
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
u
1
u
2
u
3
d
Ω
,where(25)
∂
∂
x
i
N
1
u
i
=
∂
N
1
u
∂
x
+
∂
N
1
v
∂
y
=
N
1
∂
u
∂
x
+
u
∂
N
1
∂
x
+
N
1
∂
v
∂
y
+
v
∂
N
1
∂
y
=
N
1
∂
∂
x
N
1
u
1
+
N
2
u
2
+
N
3
u
3
+
N
1
u
1
+
N
2
u
2
+
N
3
u
3
∂
N
1
∂
x
+
N
1
∂
∂
y
N
1
v
1
+
N
2
v
2
+
N
3
v
3
+
N
1
v
1
+
N
2
v
2
+
N
3
v
3
∂
N
1
∂
y
=
N
1
∂
N
1
∂
x
+
N
1
∂
N
1
∂
x
u
1
+
N
1
∂
N
2
∂
x
+
N
2
∂
N
1
∂
x
u
2
+
N
1
∂
N
3
∂
x
+
N
3
∂
N
1
∂
x
u
3
+
N
1
∂
N
1
∂
y
+
N
1
∂
N
1
∂
y
v
1
+
N
1
∂
N
2
∂
y
+
N
2
∂
N
1
∂
y
v
2
+
N
1
∂
N
3
∂
y
+
N
3
∂
N
1
∂
y
v
3
.Letting(26)
M
u
=
∫
Ω
N
u
T
N
u
d
Ω
,
C
u
=
∫
Ω
N
u
T
u
i
∂
u
i
∂
x
i
d
Ω
,
K
l
=
∫
Ω
∂
N
u
∂
x
i
∂
u
i
∂
x
i
d
Ω
,
K
u
=
∫
Ω
∂
u
N
u
T
∂
x
i
u
i
∂
u
i
∂
x
i
d
Ω
,
f
u
=
Δ
t
∫
Γ
N
u
∂
u
∂
x
+
∂
u
∂
y
d
Γ
,
f
v
=
Δ
t
∫
Γ
N
u
∂
v
∂
x
+
∂
v
∂
y
d
Γ
,then(27)
M
u
Δ
u
~
∗
=
Δ
t
-
C
u
u
~
-
μ
ρ
K
l
u
~
-
Δ
t
2
K
u
u
~
+
f
u
.Similarly, (16) can be expressed as(28)
M
u
Δ
v
~
∗
=
Δ
t
-
C
u
v
~
-
μ
ρ
K
l
v
~
-
Δ
t
2
K
v
v
~
+
f
v
.Letting(29)
Δ
U
∗
=
Δ
u
1
∗
Δ
v
1
∗
Δ
u
2
∗
Δ
v
2
∗
Δ
u
3
∗
Δ
v
3
∗
T
,we can derive the result of Δ
U
∗, and (5) becomes(30)
Δ
t
∂
u
n
+
θ
1
∂
x
i
=
Δ
t
∂
u
n
∂
x
i
+
θ
1
∂
Δ
u
∂
x
i
=
0
.According to CBS polymorphisms,(31)
Δ
u
=
u
n
+
1
-
u
n
=
Δ
u
∗
-
Δ
t
1
ρ
∂
p
n
+
θ
2
∂
x
i
+
Δ
t
2
2
u
1
ρ
∂
∂
x
i
∂
p
n
+
θ
2
∂
x
i
,and combining (31) and (30) gives(32)
Δ
t
∂
u
∂
x
i
+
θ
1
∂
Δ
u
∗
∂
x
i
-
Δ
t
θ
1
1
ρ
∂
2
p
∂
x
i
2
+
θ
2
∂
2
Δ
p
∂
x
i
2
=
0
.The weak form of (32) is written using the characteristic-Galerkin approximation as(33)
Δ
t
∫
Ω
∂
N
p
∂
x
i
u
+
θ
1
Δ
u
∗
-
Δ
t
1
ρ
∂
p
∂
x
i
+
θ
2
∂
Δ
p
∂
x
i
d
Ω
-
Δ
t
∫
Γ
N
p
u
+
θ
1
Δ
u
∗
-
Δ
t
1
ρ
∂
p
∂
x
i
+
θ
2
∂
Δ
p
∂
x
i
d
Γ
=
0
.The first term on the left side of (33) can be expressed as(34)
∫
Ω
∂
N
p
∂
x
i
u
d
Ω
=
∫
Ω
∂
N
p
∂
x
u
+
∂
N
p
∂
y
v
d
Ω
=
∫
Ω
∂
N
p
∂
x
N
1
u
1
+
N
2
u
2
+
N
3
u
3
+
∂
N
p
∂
y
N
1
v
1
+
N
2
v
2
+
N
3
v
3
d
Ω
=
∫
Ω
∂
N
1
∂
x
N
1
∂
N
1
∂
y
N
1
∂
N
1
∂
x
N
2
∂
N
1
∂
y
N
2
∂
N
1
∂
x
N
3
∂
N
1
∂
y
N
3
∂
N
2
∂
x
N
1
∂
N
2
∂
y
N
1
∂
N
2
∂
x
N
2
∂
N
2
∂
y
N
2
∂
N
2
∂
x
N
3
∂
N
2
∂
y
N
3
∂
N
3
∂
x
N
1
∂
N
3
∂
y
N
1
∂
N
3
∂
x
N
2
∂
N
3
∂
y
N
2
∂
N
3
∂
x
N
3
∂
N
3
∂
y
N
3
u
1
v
1
u
2
v
2
u
3
v
3
d
Ω
.The second term on the left side of (33) can be expressed as(35)
∫
Ω
∂
N
p
∂
x
i
Δ
u
∗
d
Ω
=
∫
Ω
∂
N
p
∂
x
Δ
u
∗
+
∂
N
p
∂
y
Δ
v
∗
d
Ω
=
∫
Ω
∂
N
p
∂
x
N
1
u
1
∗
+
N
2
u
2
∗
+
N
3
u
3
∗
+
∂
N
p
∂
y
N
1
v
1
∗
+
N
2
v
2
∗
+
N
3
v
3
∗
d
Ω
=
∫
Ω
∂
N
1
∂
x
N
1
∂
N
1
∂
y
N
1
∂
N
1
∂
x
N
2
∂
N
1
∂
y
N
2
∂
N
1
∂
x
N
3
∂
N
1
∂
y
N
3
∂
N
2
∂
x
N
1
∂
N
2
∂
y
N
1
∂
N
2
∂
x
N
2
∂
N
2
∂
y
N
2
∂
N
2
∂
x
N
3
∂
N
2
∂
y
N
3
∂
N
3
∂
x
N
1
∂
N
3
∂
y
N
1
∂
N
3
∂
x
N
2
∂
N
3
∂
y
N
2
∂
N
3
∂
x
N
3
∂
N
3
∂
y
N
3
Δ
u
1
∗
Δ
v
1
∗
Δ
u
2
∗
Δ
v
2
∗
Δ
u
3
∗
Δ
v
3
∗
d
Ω
.The third term on the left side of (33) can be expressed as(36)
∫
Ω
∂
N
p
∂
x
i
∂
p
∂
x
i
d
Ω
=
∫
Ω
∂
N
p
∂
x
∂
p
∂
x
+
∂
N
p
∂
y
∂
p
∂
y
d
Ω
=
∫
Ω
∂
N
p
∂
x
∂
∂
x
N
1
p
1
+
N
2
p
2
+
N
3
p
3
+
∂
N
p
∂
y
∂
∂
y
N
1
p
1
+
N
2
p
2
+
N
3
p
3
d
Ω
=
∫
Ω
∂
N
p
∂
x
∂
N
1
∂
x
+
∂
N
p
∂
y
∂
N
1
∂
y
∂
N
p
∂
x
∂
N
2
∂
x
+
∂
N
p
∂
y
∂
N
2
∂
y
∂
N
p
∂
x
∂
N
3
∂
x
+
∂
N
p
∂
y
∂
N
3
∂
y
p
1
p
2
p
3
d
Ω
.The fourth term on the left side of (33) can be expressed as(37)
∫
Ω
∂
N
p
∂
x
i
∂
Δ
p
∂
x
i
d
Ω
=
∫
Ω
∂
N
p
∂
x
∂
Δ
p
∂
x
+
∂
N
p
∂
y
∂
Δ
p
∂
y
d
Ω
=
∫
Ω
∂
N
p
∂
x
∂
∂
x
N
1
Δ
p
1
+
N
2
Δ
p
2
+
N
3
Δ
p
3
+
∂
N
p
∂
y
∂
∂
y
N
1
Δ
p
1
+
N
2
Δ
p
2
+
N
3
Δ
p
3
d
Ω
=
∫
Ω
∂
N
p
∂
x
∂
N
1
∂
x
+
∂
N
p
∂
y
∂
N
1
∂
y
∂
N
p
∂
x
∂
N
2
∂
x
+
∂
N
p
∂
y
∂
N
2
∂
y
∂
N
p
∂
x
∂
N
3
∂
x
+
∂
N
p
∂
y
∂
N
3
∂
y
Δ
p
1
Δ
p
2
Δ
p
3
d
Ω
,where(38)
G
=
∫
Ω
∇
N
p
T
N
u
d
Ω
,
H
=
∫
Ω
∇
N
p
T
∇
N
p
d
Ω
,
U
=
u
1
v
1
u
2
v
2
u
3
v
3
T
,
Δ
U
∗
=
Δ
u
1
∗
Δ
v
1
∗
Δ
u
2
∗
Δ
v
2
∗
Δ
u
3
∗
Δ
v
3
∗
T
,
P
=
p
1
p
2
p
3
T
,
f
p
=
-
Δ
t
∫
Γ
N
u
u
+
θ
1
Δ
u
∗
-
Δ
t
1
ρ
∂
p
∂
x
+
θ
2
∂
Δ
p
∂
x
d
Γ
.Thus, (32) can be simplified to(39)
Δ
t
2
1
ρ
θ
1
θ
2
H
Δ
P
=
Δ
t
G
U
+
θ
1
G
Δ
U
∗
-
Δ
t
1
ρ
θ
1
H
P
-
f
p
.We use a semi-implicit scheme to solve the incompressible flow problem, so0.5
≤
θ
1, θ
2
≤
1.At the completion of this stage, the value ofΔ
P is fully determined, and using the characteristic-Galerkin approximation, the weak form of (31) can be expressed as(40)
∫
Ω
N
u
Δ
u
i
d
Ω
=
∫
Ω
N
u
Δ
u
i
∗
d
Ω
-
Δ
t
1
ρ
∫
Ω
N
p
∂
p
∂
x
i
+
θ
2
∂
Δ
p
∂
x
i
d
Ω
-
Δ
t
2
2
1
ρ
∫
Ω
∇
·
N
u
i
∂
p
∂
x
i
+
θ
2
∂
Δ
p
∂
x
i
d
Ω
.The left side of (40) can be expressed as(41)
∫
Ω
N
u
Δ
u
i
d
Ω
=
∫
Ω
N
u
N
1
Δ
u
i
1
+
N
2
Δ
u
i
2
+
N
3
Δ
u
i
3
d
Ω
=
∫
Ω
N
p
N
u
N
1
N
u
N
2
N
u
N
3
Δ
u
i
1
Δ
u
i
2
Δ
u
i
3
d
Ω
=
∫
Ω
N
1
N
1
N
1
N
2
N
1
N
3
N
2
N
1
N
2
N
2
N
2
N
3
N
3
N
1
N
3
N
2
N
1
N
3
Δ
u
i
1
Δ
u
i
2
Δ
u
i
3
d
Ω
.The first term on the right side of (40) can be expressed as(42)
∫
Ω
N
u
Δ
u
i
∗
d
Ω
=
∫
Ω
N
u
N
1
Δ
u
i
1
∗
+
N
2
Δ
u
i
2
∗
+
N
3
Δ
u
i
3
∗
d
Ω
=
∫
Ω
N
p
N
u
N
1
N
u
N
2
N
u
N
3
Δ
u
i
1
∗
Δ
u
i
2
∗
Δ
u
i
3
∗
d
Ω
=
∫
Ω
N
1
N
1
N
1
N
2
N
1
N
3
N
2
N
1
N
2
N
2
N
2
N
3
N
3
N
1
N
3
N
2
N
1
N
3
Δ
u
i
1
∗
Δ
u
i
2
∗
Δ
u
i
3
∗
d
Ω
.The second term on the right side of (40) can be expressed as(43)
∫
Ω
N
p
∂
p
∂
x
i
d
Ω
=
∫
Ω
N
p
∂
∂
x
i
N
1
p
1
+
N
2
p
2
+
N
3
p
3
d
Ω
=
∫
Ω
N
p
∂
N
1
∂
x
i
N
p
∂
N
2
∂
x
i
N
p
∂
N
3
∂
x
i
p
1
p
2
p
3
d
Ω
.The third term on the right side of (40) can be expressed as(44)
∫
Ω
N
p
∂
Δ
p
∂
x
i
d
Ω
=
∫
Ω
N
p
∂
∂
x
i
N
1
Δ
p
1
+
N
2
Δ
p
2
+
N
3
Δ
p
3
d
Ω
=
∫
Ω
N
p
∂
N
1
∂
x
i
N
p
∂
N
2
∂
x
i
N
p
∂
N
3
∂
x
i
Δ
p
1
Δ
p
2
Δ
p
3
d
Ω
.The fourth term on the right side of (40) can be expressed as(45)
∫
Ω
∇
·
N
p
u
i
∂
p
∂
x
i
d
Ω
=
∫
Ω
∇
·
N
p
u
i
∂
∂
x
i
N
1
p
1
+
N
2
p
2
+
N
3
p
3
d
Ω
=
∫
Ω
∇
·
N
p
u
i
∂
N
1
∂
x
i
∇
·
N
p
u
i
∂
N
2
∂
x
i
∇
·
N
p
u
i
∂
N
3
∂
x
i
m
m
m
×
p
1
p
2
p
3
d
Ω
,where(46)
∇
·
N
p
u
i
=
∂
N
p
u
∂
x
+
∂
N
p
v
∂
y
=
N
p
∂
u
∂
x
+
u
∂
N
p
∂
x
+
N
p
∂
v
∂
y
+
v
∂
N
p
∂
y
=
N
p
∂
∂
x
N
1
u
1
+
N
2
u
2
+
N
3
u
3
+
N
1
u
1
+
N
2
u
2
+
N
3
u
3
∂
N
p
∂
x
+
N
p
∂
∂
y
N
1
v
1
+
N
2
v
2
+
N
3
v
3
+
N
1
v
1
+
N
2
v
2
+
N
3
v
3
∂
N
p
∂
y
=
N
p
∂
N
1
∂
x
+
N
1
∂
N
p
∂
x
u
1
+
N
p
∂
N
2
∂
x
+
N
2
∂
N
p
∂
x
u
2
+
N
p
∂
N
3
∂
x
+
N
3
∂
N
p
∂
x
u
3
+
N
p
∂
N
1
∂
y
+
N
1
∂
N
p
∂
y
v
1
+
N
p
∂
N
2
∂
y
+
N
2
∂
N
p
∂
y
v
2
+
N
p
∂
N
3
∂
y
+
N
3
∂
N
p
∂
y
v
3
.The fifth term on the right side of (40) can be expressed as(47)
∫
Ω
∇
·
N
p
u
i
∂
Δ
p
∂
x
i
d
Ω
=
∫
Ω
∇
·
N
p
u
i
∂
∂
x
i
N
1
Δ
p
1
+
N
2
Δ
p
2
+
N
3
Δ
p
3
d
Ω
=
∫
Ω
∇
·
N
p
u
i
∂
N
1
∂
x
i
∇
·
N
p
u
i
∂
N
2
∂
x
i
∇
·
N
p
u
i
∂
N
3
∂
x
i
m
m
m
c
×
Δ
p
1
Δ
p
2
Δ
p
3
d
Ω
,where(48)
K
p
=
∫
Ω
N
p
∂
N
∂
x
i
d
Ω
,
K
g
=
∫
Ω
∇
·
u
i
N
p
∂
N
∂
x
i
d
Ω
.Then, we derive(49)
M
u
Δ
u
i
=
M
u
Δ
u
i
∗
-
Δ
t
1
ρ
K
p
P
+
θ
2
Δ
P
+
Δ
t
2
K
g
P
+
θ
2
Δ
P
.The values of Δ
U
∗ and Δ
P can be obtained from Steps 1 and 2, and Δ
U can be calculated in Step 3 from (49). Then, U
n
+
1 can be substituted into Step 1 as an initial value for the next iterative computation.


## 3.1. CBS Procedure

The theory of the CBS algorithm is based on the split- and characteristic-Galerkin procedures, which are extended in the FEM. The CBS algorithm involves three indispensable steps: Step1 is to compute the intermediate momentum variable; Step 2 is to calculate the density or pressure field; and Step 3 involves combining the introduction of the density or pressure fields to obtain the correct momentum variables.Step 1.
The solution forΔ
U
~
∗ in matrix form is(8)
Δ
U
~
∗
=
-
M
-
1
Δ
t
C
U
~
+
K
u
~
-
f
-
Δ
t
K
u
U
~
+
f
s
n
,where the quantities identified by ~ symbols indicate nodal values, and(9)
M
=
∫
Ω
N
T
N
d
Ω
,
C
=
∫
Ω
N
T
∇
u
N
d
Ω
,
K
=
∫
Ω
B
T
μ
I
0
-
2
3
m
m
T
B
d
Ω
=
∫
Ω
μ
∂
N
T
∂
x
j
∂
N
∂
x
j
d
Ω
,
f
=
∫
Ω
N
T
ρ
g
d
Ω
+
∫
Γ
N
T
t
d
d
Γ
,where B is a matrix, N is a deformation function, m
=
[
1,1
,
0
]
T, and I
0
=
d
i
a
g
(
2,2
,
1
). After integration by parts, the expressions for K
u and f
s are (10)
K
u
=
-
1
2
∫
Ω
∇
T
u
N
T
∇
T
u
N
d
Ω
,
f
s
=
-
1
2
∫
Ω
∇
T
u
N
T
ρ
g
d
Ω
.Step 2.
The value ofΔ
p
~ can be easily determined using the expression(11)
M
P
+
Δ
t
2
θ
1
θ
2
H
Δ
P
~
=
Δ
t
G
U
~
n
+
θ
1
G
U
~
∗
-
Δ
t
θ
1
H
P
~
n
-
f
p
.The matrices are(12)
H
=
∫
Ω
∇
N
T
∇
N
d
Ω
,
M
p
=
∫
Ω
N
T
1
c
2
n
N
d
Ω
,
G
=
∫
Ω
∇
N
T
N
d
Ω
,
f
p
=
Δ
t
∫
Γ
N
T
n
T
U
~
n
+
θ
1
Δ
U
~
∗
-
Δ
t
∇
p
n
+
θ
2
d
Γ
.Step 3.
The value ofΔ
U
~ is determined from(13)
Δ
U
~
=
Δ
U
~
∗
-
M
-
1
Δ
t
G
T
p
~
n
+
θ
2
Δ
p
~
+
Δ
t
2
P
p
~
n
,where(14)
P
=
∫
Ω
∇
u
N
T
∇
N
d
Ω
.


## 3.2. Derivation Process of the Free-Flow Model

From (5), we derive(15)
∂
u
∂
t
=
-
u
∂
u
∂
x
+
v
∂
u
∂
y
-
1
ρ
∂
p
∂
x
+
μ
ρ
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
,
(16)
∂
v
∂
t
=
-
u
∂
v
∂
x
+
v
∂
v
∂
y
-
1
ρ
∂
p
∂
y
+
μ
ρ
∂
2
v
∂
x
2
+
∂
2
v
∂
y
2
.From Step 1, after the split-temporal discretization, (15) becomes (omitting the third-order derivative above)(17)
Δ
u
=
u
n
+
1
-
u
n
=
Δ
t
-
u
∂
u
∂
x
+
v
∂
u
∂
y
-
1
ρ
∂
p
∂
x
+
μ
ρ
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
+
Δ
t
2
2
u
∂
∂
x
u
∂
u
∂
x
+
v
∂
u
∂
y
+
1
ρ
∂
p
∂
x
.From Step 2, after spatial discretization, we find that(18)
Δ
u
=
N
u
Δ
u
~
,
Δ
u
∗
=
N
u
Δ
u
~
∗
,
u
=
N
u
u
~
,
v
=
N
v
v
~
,
p
=
N
p
p
~
,
N
=
N
1
,
N
2
,
N
3
,
Δ
u
~
=
Δ
u
1
,
Δ
u
2
,
Δ
u
3
T
,
Δ
u
~
∗
=
Δ
u
1
∗
,
Δ
u
2
∗
,
Δ
u
3
∗
T
,
u
~
=
u
1
,
u
2
,
u
3
T
,
v
~
=
v
1
,
v
2
,
v
3
T
,
p
=
p
1
,
p
2
,
p
3
T
,
N
u
=
N
v
=
N
p
.Then,(19)
∫
Ω
N
u
Δ
u
∗
d
Ω
=
Δ
t
-
∫
Ω
N
u
u
∂
u
∂
x
+
v
∂
u
∂
y
d
Ω
m
m
m
v
-
μ
ρ
∫
Ω
∇
·
N
u
∂
u
∂
x
+
∂
u
∂
y
d
Ω
+
Δ
t
2
2
∫
Ω
∇
N
u
u
-
u
∂
u
∂
x
+
v
∂
u
∂
y
d
Ω
+
Δ
t
μ
ρ
∫
Γ
N
u
∂
u
∂
x
+
∂
u
∂
y
d
Γ
.The left side of (19) can be expressed as(20)
∫
Ω
N
u
Δ
u
∗
d
Ω
=
∫
Ω
N
u
N
1
Δ
u
1
∗
+
N
2
Δ
u
2
∗
+
N
3
Δ
u
3
∗
d
Ω
=
∫
Ω
N
u
N
1
d
Ω
∫
Ω
N
u
N
2
d
Ω
∫
Ω
N
u
N
3
d
Ω
×
Δ
u
1
∗
Δ
u
2
∗
Δ
u
3
∗
=
∫
Ω
N
1
N
1
d
Ω
∫
Ω
N
1
N
2
d
Ω
∫
Ω
N
1
N
3
d
Ω
∫
Ω
N
2
N
1
d
Ω
∫
Ω
N
2
N
2
d
Ω
∫
Ω
N
2
N
3
d
Ω
∫
Ω
N
3
N
1
d
Ω
∫
Ω
N
3
N
2
d
Ω
∫
Ω
N
3
N
3
d
Ω
×
Δ
u
1
∗
Δ
u
2
∗
Δ
u
3
∗
.The first term on the right side of (19) can be expressed as(21)
∫
Ω
N
u
u
∂
u
∂
x
+
v
∂
u
∂
y
d
Ω
=
∫
Ω
N
u
u
∂
N
1
∂
x
u
1
+
∂
N
2
∂
x
u
2
+
∂
N
3
∂
x
u
3
+
v
∂
N
1
∂
y
u
1
+
∂
N
2
∂
y
u
2
+
∂
N
3
∂
y
u
3
=
∫
Ω
N
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
u
1
+
u
∂
N
2
∂
x
+
v
∂
N
2
∂
y
u
2
+
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
u
3
d
Ω
=
∫
Ω
N
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
d
Ω
∫
Ω
N
u
u
∂
N
2
∂
x
+
v
∂
N
2
∂
y
d
Ω
∫
Ω
N
u
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
d
Ω
u
1
u
2
u
3
=
∫
Ω
N
1
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
d
Ω
∫
Ω
N
1
u
∂
N
2
∂
x
+
v
∂
N
2
∂
y
d
Ω
∫
Ω
N
1
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
d
Ω
∫
Ω
N
2
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
d
Ω
∫
Ω
N
2
u
∂
N
2
∂
x
+
v
∂
N
2
∂
y
d
Ω
∫
Ω
N
2
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
d
Ω
∫
Ω
N
3
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
d
Ω
∫
Ω
N
3
u
∂
N
2
∂
x
+
v
∂
N
2
∂
y
d
Ω
∫
Ω
N
3
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
d
Ω
u
1
u
2
u
3
,where(22)
∫
Ω
N
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
d
Ω
=
∫
Ω
N
u
∂
N
1
∂
y
N
1
u
1
+
N
2
u
2
+
N
3
u
3
∂
N
1
∂
x
m
m
m
m
m
m
m
c
m
+
N
1
v
1
+
N
2
v
2
+
N
3
v
3
∂
N
1
∂
y
∂
N
1
∂
x
d
Ω
=
∫
Ω
N
u
N
1
∂
N
1
∂
x
d
Ω
∫
Ω
N
u
N
2
∂
N
1
∂
x
d
Ω
∫
Ω
N
u
N
3
∂
N
1
∂
x
d
Ω
∫
Ω
N
u
N
1
∂
N
1
∂
y
d
Ω
∫
Ω
N
u
N
2
∂
N
1
∂
y
d
Ω
∫
Ω
N
u
N
3
∂
N
1
∂
y
d
Ω
T
u
1
u
2
u
3
v
1
v
2
v
3
.The second term on the right side of (19) can be expressed as(23)
∫
Ω
∂
N
u
∂
x
i
∂
u
∂
x
+
∂
u
∂
y
d
Ω
=
∫
Ω
∂
N
u
∂
x
i
∂
N
1
∂
x
+
∂
N
1
∂
y
u
1
+
∂
N
2
∂
x
+
∂
N
2
∂
y
u
2
+
∂
N
3
∂
x
+
∂
N
3
∂
y
u
3
d
Ω
=
∫
Ω
∂
N
u
∂
x
i
∂
N
1
∂
x
+
∂
N
1
∂
y
∂
N
u
∂
x
i
∂
N
2
∂
x
+
∂
N
2
∂
y
∂
N
u
∂
x
i
∂
N
3
∂
x
+
∂
N
3
∂
y
u
1
u
2
u
3
d
Ω
.The third term on the right side of (19) can be expressed as(24)
∫
Ω
∂
N
u
∂
x
i
N
u
u
-
u
∂
u
∂
x
+
v
∂
u
∂
y
d
Ω
=
-
∫
Ω
∂
∂
x
i
N
u
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
u
1
+
u
∂
N
2
∂
x
+
v
∂
N
2
∂
y
u
2
+
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
u
3
d
Ω
=
-
∫
Ω
∂
∂
x
i
N
1
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
⋯
∂
∂
x
i
N
1
u
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
∂
∂
x
i
N
2
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
⋯
∂
∂
x
i
N
2
u
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
∂
∂
x
i
N
3
u
u
∂
N
1
∂
x
+
v
∂
N
1
∂
y
⋯
∂
∂
x
i
N
3
u
u
∂
N
3
∂
x
+
v
∂
N
3
∂
y
u
1
u
2
u
3
d
Ω
,where(25)
∂
∂
x
i
N
1
u
i
=
∂
N
1
u
∂
x
+
∂
N
1
v
∂
y
=
N
1
∂
u
∂
x
+
u
∂
N
1
∂
x
+
N
1
∂
v
∂
y
+
v
∂
N
1
∂
y
=
N
1
∂
∂
x
N
1
u
1
+
N
2
u
2
+
N
3
u
3
+
N
1
u
1
+
N
2
u
2
+
N
3
u
3
∂
N
1
∂
x
+
N
1
∂
∂
y
N
1
v
1
+
N
2
v
2
+
N
3
v
3
+
N
1
v
1
+
N
2
v
2
+
N
3
v
3
∂
N
1
∂
y
=
N
1
∂
N
1
∂
x
+
N
1
∂
N
1
∂
x
u
1
+
N
1
∂
N
2
∂
x
+
N
2
∂
N
1
∂
x
u
2
+
N
1
∂
N
3
∂
x
+
N
3
∂
N
1
∂
x
u
3
+
N
1
∂
N
1
∂
y
+
N
1
∂
N
1
∂
y
v
1
+
N
1
∂
N
2
∂
y
+
N
2
∂
N
1
∂
y
v
2
+
N
1
∂
N
3
∂
y
+
N
3
∂
N
1
∂
y
v
3
.Letting(26)
M
u
=
∫
Ω
N
u
T
N
u
d
Ω
,
C
u
=
∫
Ω
N
u
T
u
i
∂
u
i
∂
x
i
d
Ω
,
K
l
=
∫
Ω
∂
N
u
∂
x
i
∂
u
i
∂
x
i
d
Ω
,
K
u
=
∫
Ω
∂
u
N
u
T
∂
x
i
u
i
∂
u
i
∂
x
i
d
Ω
,
f
u
=
Δ
t
∫
Γ
N
u
∂
u
∂
x
+
∂
u
∂
y
d
Γ
,
f
v
=
Δ
t
∫
Γ
N
u
∂
v
∂
x
+
∂
v
∂
y
d
Γ
,then(27)
M
u
Δ
u
~
∗
=
Δ
t
-
C
u
u
~
-
μ
ρ
K
l
u
~
-
Δ
t
2
K
u
u
~
+
f
u
.Similarly, (16) can be expressed as(28)
M
u
Δ
v
~
∗
=
Δ
t
-
C
u
v
~
-
μ
ρ
K
l
v
~
-
Δ
t
2
K
v
v
~
+
f
v
.Letting(29)
Δ
U
∗
=
Δ
u
1
∗
Δ
v
1
∗
Δ
u
2
∗
Δ
v
2
∗
Δ
u
3
∗
Δ
v
3
∗
T
,we can derive the result of Δ
U
∗, and (5) becomes(30)
Δ
t
∂
u
n
+
θ
1
∂
x
i
=
Δ
t
∂
u
n
∂
x
i
+
θ
1
∂
Δ
u
∂
x
i
=
0
.According to CBS polymorphisms,(31)
Δ
u
=
u
n
+
1
-
u
n
=
Δ
u
∗
-
Δ
t
1
ρ
∂
p
n
+
θ
2
∂
x
i
+
Δ
t
2
2
u
1
ρ
∂
∂
x
i
∂
p
n
+
θ
2
∂
x
i
,and combining (31) and (30) gives(32)
Δ
t
∂
u
∂
x
i
+
θ
1
∂
Δ
u
∗
∂
x
i
-
Δ
t
θ
1
1
ρ
∂
2
p
∂
x
i
2
+
θ
2
∂
2
Δ
p
∂
x
i
2
=
0
.The weak form of (32) is written using the characteristic-Galerkin approximation as(33)
Δ
t
∫
Ω
∂
N
p
∂
x
i
u
+
θ
1
Δ
u
∗
-
Δ
t
1
ρ
∂
p
∂
x
i
+
θ
2
∂
Δ
p
∂
x
i
d
Ω
-
Δ
t
∫
Γ
N
p
u
+
θ
1
Δ
u
∗
-
Δ
t
1
ρ
∂
p
∂
x
i
+
θ
2
∂
Δ
p
∂
x
i
d
Γ
=
0
.The first term on the left side of (33) can be expressed as(34)
∫
Ω
∂
N
p
∂
x
i
u
d
Ω
=
∫
Ω
∂
N
p
∂
x
u
+
∂
N
p
∂
y
v
d
Ω
=
∫
Ω
∂
N
p
∂
x
N
1
u
1
+
N
2
u
2
+
N
3
u
3
+
∂
N
p
∂
y
N
1
v
1
+
N
2
v
2
+
N
3
v
3
d
Ω
=
∫
Ω
∂
N
1
∂
x
N
1
∂
N
1
∂
y
N
1
∂
N
1
∂
x
N
2
∂
N
1
∂
y
N
2
∂
N
1
∂
x
N
3
∂
N
1
∂
y
N
3
∂
N
2
∂
x
N
1
∂
N
2
∂
y
N
1
∂
N
2
∂
x
N
2
∂
N
2
∂
y
N
2
∂
N
2
∂
x
N
3
∂
N
2
∂
y
N
3
∂
N
3
∂
x
N
1
∂
N
3
∂
y
N
1
∂
N
3
∂
x
N
2
∂
N
3
∂
y
N
2
∂
N
3
∂
x
N
3
∂
N
3
∂
y
N
3
u
1
v
1
u
2
v
2
u
3
v
3
d
Ω
.The second term on the left side of (33) can be expressed as(35)
∫
Ω
∂
N
p
∂
x
i
Δ
u
∗
d
Ω
=
∫
Ω
∂
N
p
∂
x
Δ
u
∗
+
∂
N
p
∂
y
Δ
v
∗
d
Ω
=
∫
Ω
∂
N
p
∂
x
N
1
u
1
∗
+
N
2
u
2
∗
+
N
3
u
3
∗
+
∂
N
p
∂
y
N
1
v
1
∗
+
N
2
v
2
∗
+
N
3
v
3
∗
d
Ω
=
∫
Ω
∂
N
1
∂
x
N
1
∂
N
1
∂
y
N
1
∂
N
1
∂
x
N
2
∂
N
1
∂
y
N
2
∂
N
1
∂
x
N
3
∂
N
1
∂
y
N
3
∂
N
2
∂
x
N
1
∂
N
2
∂
y
N
1
∂
N
2
∂
x
N
2
∂
N
2
∂
y
N
2
∂
N
2
∂
x
N
3
∂
N
2
∂
y
N
3
∂
N
3
∂
x
N
1
∂
N
3
∂
y
N
1
∂
N
3
∂
x
N
2
∂
N
3
∂
y
N
2
∂
N
3
∂
x
N
3
∂
N
3
∂
y
N
3
Δ
u
1
∗
Δ
v
1
∗
Δ
u
2
∗
Δ
v
2
∗
Δ
u
3
∗
Δ
v
3
∗
d
Ω
.The third term on the left side of (33) can be expressed as(36)
∫
Ω
∂
N
p
∂
x
i
∂
p
∂
x
i
d
Ω
=
∫
Ω
∂
N
p
∂
x
∂
p
∂
x
+
∂
N
p
∂
y
∂
p
∂
y
d
Ω
=
∫
Ω
∂
N
p
∂
x
∂
∂
x
N
1
p
1
+
N
2
p
2
+
N
3
p
3
+
∂
N
p
∂
y
∂
∂
y
N
1
p
1
+
N
2
p
2
+
N
3
p
3
d
Ω
=
∫
Ω
∂
N
p
∂
x
∂
N
1
∂
x
+
∂
N
p
∂
y
∂
N
1
∂
y
∂
N
p
∂
x
∂
N
2
∂
x
+
∂
N
p
∂
y
∂
N
2
∂
y
∂
N
p
∂
x
∂
N
3
∂
x
+
∂
N
p
∂
y
∂
N
3
∂
y
p
1
p
2
p
3
d
Ω
.The fourth term on the left side of (33) can be expressed as(37)
∫
Ω
∂
N
p
∂
x
i
∂
Δ
p
∂
x
i
d
Ω
=
∫
Ω
∂
N
p
∂
x
∂
Δ
p
∂
x
+
∂
N
p
∂
y
∂
Δ
p
∂
y
d
Ω
=
∫
Ω
∂
N
p
∂
x
∂
∂
x
N
1
Δ
p
1
+
N
2
Δ
p
2
+
N
3
Δ
p
3
+
∂
N
p
∂
y
∂
∂
y
N
1
Δ
p
1
+
N
2
Δ
p
2
+
N
3
Δ
p
3
d
Ω
=
∫
Ω
∂
N
p
∂
x
∂
N
1
∂
x
+
∂
N
p
∂
y
∂
N
1
∂
y
∂
N
p
∂
x
∂
N
2
∂
x
+
∂
N
p
∂
y
∂
N
2
∂
y
∂
N
p
∂
x
∂
N
3
∂
x
+
∂
N
p
∂
y
∂
N
3
∂
y
Δ
p
1
Δ
p
2
Δ
p
3
d
Ω
,where(38)
G
=
∫
Ω
∇
N
p
T
N
u
d
Ω
,
H
=
∫
Ω
∇
N
p
T
∇
N
p
d
Ω
,
U
=
u
1
v
1
u
2
v
2
u
3
v
3
T
,
Δ
U
∗
=
Δ
u
1
∗
Δ
v
1
∗
Δ
u
2
∗
Δ
v
2
∗
Δ
u
3
∗
Δ
v
3
∗
T
,
P
=
p
1
p
2
p
3
T
,
f
p
=
-
Δ
t
∫
Γ
N
u
u
+
θ
1
Δ
u
∗
-
Δ
t
1
ρ
∂
p
∂
x
+
θ
2
∂
Δ
p
∂
x
d
Γ
.Thus, (32) can be simplified to(39)
Δ
t
2
1
ρ
θ
1
θ
2
H
Δ
P
=
Δ
t
G
U
+
θ
1
G
Δ
U
∗
-
Δ
t
1
ρ
θ
1
H
P
-
f
p
.We use a semi-implicit scheme to solve the incompressible flow problem, so0.5
≤
θ
1, θ
2
≤
1.At the completion of this stage, the value ofΔ
P is fully determined, and using the characteristic-Galerkin approximation, the weak form of (31) can be expressed as(40)
∫
Ω
N
u
Δ
u
i
d
Ω
=
∫
Ω
N
u
Δ
u
i
∗
d
Ω
-
Δ
t
1
ρ
∫
Ω
N
p
∂
p
∂
x
i
+
θ
2
∂
Δ
p
∂
x
i
d
Ω
-
Δ
t
2
2
1
ρ
∫
Ω
∇
·
N
u
i
∂
p
∂
x
i
+
θ
2
∂
Δ
p
∂
x
i
d
Ω
.The left side of (40) can be expressed as(41)
∫
Ω
N
u
Δ
u
i
d
Ω
=
∫
Ω
N
u
N
1
Δ
u
i
1
+
N
2
Δ
u
i
2
+
N
3
Δ
u
i
3
d
Ω
=
∫
Ω
N
p
N
u
N
1
N
u
N
2
N
u
N
3
Δ
u
i
1
Δ
u
i
2
Δ
u
i
3
d
Ω
=
∫
Ω
N
1
N
1
N
1
N
2
N
1
N
3
N
2
N
1
N
2
N
2
N
2
N
3
N
3
N
1
N
3
N
2
N
1
N
3
Δ
u
i
1
Δ
u
i
2
Δ
u
i
3
d
Ω
.The first term on the right side of (40) can be expressed as(42)
∫
Ω
N
u
Δ
u
i
∗
d
Ω
=
∫
Ω
N
u
N
1
Δ
u
i
1
∗
+
N
2
Δ
u
i
2
∗
+
N
3
Δ
u
i
3
∗
d
Ω
=
∫
Ω
N
p
N
u
N
1
N
u
N
2
N
u
N
3
Δ
u
i
1
∗
Δ
u
i
2
∗
Δ
u
i
3
∗
d
Ω
=
∫
Ω
N
1
N
1
N
1
N
2
N
1
N
3
N
2
N
1
N
2
N
2
N
2
N
3
N
3
N
1
N
3
N
2
N
1
N
3
Δ
u
i
1
∗
Δ
u
i
2
∗
Δ
u
i
3
∗
d
Ω
.The second term on the right side of (40) can be expressed as(43)
∫
Ω
N
p
∂
p
∂
x
i
d
Ω
=
∫
Ω
N
p
∂
∂
x
i
N
1
p
1
+
N
2
p
2
+
N
3
p
3
d
Ω
=
∫
Ω
N
p
∂
N
1
∂
x
i
N
p
∂
N
2
∂
x
i
N
p
∂
N
3
∂
x
i
p
1
p
2
p
3
d
Ω
.The third term on the right side of (40) can be expressed as(44)
∫
Ω
N
p
∂
Δ
p
∂
x
i
d
Ω
=
∫
Ω
N
p
∂
∂
x
i
N
1
Δ
p
1
+
N
2
Δ
p
2
+
N
3
Δ
p
3
d
Ω
=
∫
Ω
N
p
∂
N
1
∂
x
i
N
p
∂
N
2
∂
x
i
N
p
∂
N
3
∂
x
i
Δ
p
1
Δ
p
2
Δ
p
3
d
Ω
.The fourth term on the right side of (40) can be expressed as(45)
∫
Ω
∇
·
N
p
u
i
∂
p
∂
x
i
d
Ω
=
∫
Ω
∇
·
N
p
u
i
∂
∂
x
i
N
1
p
1
+
N
2
p
2
+
N
3
p
3
d
Ω
=
∫
Ω
∇
·
N
p
u
i
∂
N
1
∂
x
i
∇
·
N
p
u
i
∂
N
2
∂
x
i
∇
·
N
p
u
i
∂
N
3
∂
x
i
m
m
m
×
p
1
p
2
p
3
d
Ω
,where(46)
∇
·
N
p
u
i
=
∂
N
p
u
∂
x
+
∂
N
p
v
∂
y
=
N
p
∂
u
∂
x
+
u
∂
N
p
∂
x
+
N
p
∂
v
∂
y
+
v
∂
N
p
∂
y
=
N
p
∂
∂
x
N
1
u
1
+
N
2
u
2
+
N
3
u
3
+
N
1
u
1
+
N
2
u
2
+
N
3
u
3
∂
N
p
∂
x
+
N
p
∂
∂
y
N
1
v
1
+
N
2
v
2
+
N
3
v
3
+
N
1
v
1
+
N
2
v
2
+
N
3
v
3
∂
N
p
∂
y
=
N
p
∂
N
1
∂
x
+
N
1
∂
N
p
∂
x
u
1
+
N
p
∂
N
2
∂
x
+
N
2
∂
N
p
∂
x
u
2
+
N
p
∂
N
3
∂
x
+
N
3
∂
N
p
∂
x
u
3
+
N
p
∂
N
1
∂
y
+
N
1
∂
N
p
∂
y
v
1
+
N
p
∂
N
2
∂
y
+
N
2
∂
N
p
∂
y
v
2
+
N
p
∂
N
3
∂
y
+
N
3
∂
N
p
∂
y
v
3
.The fifth term on the right side of (40) can be expressed as(47)
∫
Ω
∇
·
N
p
u
i
∂
Δ
p
∂
x
i
d
Ω
=
∫
Ω
∇
·
N
p
u
i
∂
∂
x
i
N
1
Δ
p
1
+
N
2
Δ
p
2
+
N
3
Δ
p
3
d
Ω
=
∫
Ω
∇
·
N
p
u
i
∂
N
1
∂
x
i
∇
·
N
p
u
i
∂
N
2
∂
x
i
∇
·
N
p
u
i
∂
N
3
∂
x
i
m
m
m
c
×
Δ
p
1
Δ
p
2
Δ
p
3
d
Ω
,where(48)
K
p
=
∫
Ω
N
p
∂
N
∂
x
i
d
Ω
,
K
g
=
∫
Ω
∇
·
u
i
N
p
∂
N
∂
x
i
d
Ω
.Then, we derive(49)
M
u
Δ
u
i
=
M
u
Δ
u
i
∗
-
Δ
t
1
ρ
K
p
P
+
θ
2
Δ
P
+
Δ
t
2
K
g
P
+
θ
2
Δ
P
.The values of Δ
U
∗ and Δ
P can be obtained from Steps 1 and 2, and Δ
U can be calculated in Step 3 from (49). Then, U
n
+
1 can be substituted into Step 1 as an initial value for the next iterative computation.


## 4. Case Study and Discussion

An example calculation was performed for the initial condition of a cave with a reservoir pressure of 10 MPa, a fluid density of 800 Kg/m3, a fluid viscosity of 1 mPa s, and an axial inflow velocity of 0.01 m/s. The entrance was in the middle of the left boundary of the cave, and the exit was opposite to the entrance.The boundary conditions are nonslip, so thatu
=
v
=
0. The grid at the cave boundary is refined. The number of nodes is 700, and the number of units is 1282. The influence of mass force is ignored, and the physical parameters of the flow in the reservoir are constant with time. Figure 2 shows a two-dimensional single-phase free-flow model and Delaunay triangle subdivision. Figure 3 shows a single-grid figure.Figure 2
Two-dimensional single-phase free-flow model and Delaunay triangle subdivision.Figure 3
A single-grid figure.Based on the research results above, we used MATLAB to develop software for numerical simulation of two-dimensional single-phase free flow in a cave. The single-phase free flow of a large-scale fracture-cave reservoir can thus be simulated, with the simulation performed on any computer with MATLAB installed. The functions of the software include data input and storage, grid partitioning, simulation of single-phase flow, data output, and graphical plotting. Figures4, 5, and 6 show outputs of the numerical simulation software.Figure 4
Velocity field of the free-flow region at different times.Figure 5
Velocity vector streamline diagram of the free-flow region at different times.Figure 6
Pressure profile patterns of the free-flow region at different times.Figures4–6 represent the temporal evolution of the velocity field, velocity vector streamlines, and profile patterns of the free-flow region, respectively. The figures show that two vortices form at the two corners, beside the entrance to the cave, which is in agreement with the actual characteristics of free flow. Simulations indicate that the CBS algorithm has the advantage of flexibility of grid subdivision, and high efficiency and precision of calculations for the single-phase free-flow numerical simulation of fracture-cave carbonate reservoirs.This work provides a theoretical basis for future studies on the FEM approach for multiphase flow and Darcy flow-free flow in a discrete media network model of a fracture-cave reservoir; this is an important development in the modeling of fracture-cave flow regimes and should be of considerable use in future research in this field.


## 5. Conclusions

Based on a single-phase free-flow numerical simulation of a fracture-cave reservoir, the CBS algorithm was adopted to obtain numerical solutions accurately and quickly, in which each step only requires solution of a system of linear equations. The velocity field, velocity vector streamline diagram, and profile pattern of the free-flow region at different times were determined, and they demonstrated good conformity with the actual characteristics of free flow. The results confirm that the CBS algorithm is applicable to the numerical simulation of a discrete media model of a fracture-cave reservoir.Similarly to the traditional method, we can use statistical information and data on fractures and caves obtained from an oil field to construct a discrete fracture-cave geology geometry model. Then, we can use the approach in this study to model realistic single-phase cave free flow in fractured cave media. This paper provides an effective approach to developing a fracture-cave reservoir, including development plan design, dynamic prediction, reservoir evaluation, optimization of engineering technical problems, and decision-making.



---

*Source: 815051-2015-09-27.xml*