Page 94 - Fundamentals of Computational Geoscience Numerical Methods and Algorithms

4.2 Development of the Term Splitting Algorithm 81

to the porosity of the aquifer. In this regard, we have to deal with the problem of
variable Peclet and Courant numbers in the transient analysis of fluid-rock inter-
action problems. To solve this problem effectively and efficiently, we develop the
term splitting algorithm in this section. The basic idea behind this algorithm is that
through some mathematical manipulations, we invert the problem of variable Peclet
and Courant numbers into that of constant Peclet and Courant numbers so that a
fixed finite element mesh and fixed time step of integration can be safely used in the
numerical modelling of fluid-rock interaction problems.
If the pore-fluid saturated porous medium has a non-zero initial porosity, then
we can use this initial porosity as a reference porosity and change Eq. (4.9) into the
following form:

φ ∂C Ci ∂C Ci ∂C Ci e 2 φ
φ 0 +u x +u y = D ∇ C Ci + φ 0 R i (i = m+1, m+2,..., n).
i
φ 0 ∂t ∂x ∂y φ 0
(4.16)
Dividing both the left and right hand sides of Eq. (4.16) by φ/φ 0 yields the fol-
lowing equation:

∂C Ci φ 0 ∂C Ci φ 0 ∂C Ci φ 0 e 2
φ 0 + u x + u y = D ∇ C Ci +φ 0 R i (i = m+1, m+2,..., n).
i
∂t φ ∂x φ ∂y φ
(4.17)
Using the initial Darcy velocities, u x0 and u y0 , as reference velocities in the x and
y directions respectively, Eq. (4.17) can be rewritten as follows:

∂C Ci φ 0 ∂C Ci φ 0 ∂C Ci
φ 0 + u x0 + u x − u x0 + u y0 + u y − u y0
∂t φ ∂x φ ∂y
φ 0 e 2
= D ∇ C Ci + φ 0 R i
i
φ
(i = m + 1, m + 2,..., n). (4.18)
It is noted that in Eq. (4.18), the coefficient in front of ∂C Ci/∂x has been split
into a constant coefficient, u x0 , and a variable term, (u x φ 0 )/φ − u x0 . Similarly, the
coefficient in front of ∂C Ci/∂y has also been split into a constant coefficient, u y0 ,
and a variable term, (u y φ 0 )/φ − u y0 . This term splitting process, although it is car-
ried out mathematically, forms a very important step in the proposed term splitting
algorithm.
If the first derivative terms with constant coefficients are kept in the left hand side
of Eq. (4.18) and the rest terms are moved to the right hand side, Eq. (4.18) can be
rewritten as follows:

∂C Ci ∂C Ci ∂C Ci e 2 e
φ 0 +u x0 +u y0 = D ∇ C Ci +φ 0 R i + R i (i = m +1, m +2,..., n),
0i
∂t ∂x ∂y
(4.19)