Page 95 - Fundamentals of Computational Geoscience Numerical Methods and Algorithms

82 4 Algorithm for Simulating Fluid-Rock Interaction Problems

e φ 0 ∂C Ci φ 0 ∂C Ci
R = u x0 − + u yo − (i = m + 1, m + 2,..., n),
i u x u y
φ ∂x φ ∂y
(4.20)
e
D = φ 0 D i (i = m + 1, m + 2,..., n), (4.21)
0i
where φ 0 is the initial porosity of the porous rock mass; u x0 and u y0 are the reference
e
velocities, which are the initial Darcy velocities, in the x andy directions; R is the
i
equivalent source/sink term due to the variation in both the velocity of pore-fluid
and the porosity of the rock mass.
It is clear that the proposed term splitting algorithm consists of the following
two main steps. First, by means of some rigorous mathematical manipulations,
Eq. (4.9) with variable coefficients in front of the first and second derivatives has
been changed into Eq. (4.19), which has constant coefficients in front of the first
and second derivatives. From the mathematical point of view, this means that we
have changed a partial differential equation with variable coefficients (Eq. (4.9))
into another one with constant coefficients in front of the first and second deriva-
tives. Second, Eq. (4.19) is directly used to obtain a numerical solution for Equation
(4.9) in the finite element analysis, because both Eqs. (4.19) and (4.9) are mathemat-
ically equivalent. However, from the computational point of view, Eq. (4.19) can be
solved much easier than Eq. (4.9) in the finite element analysis. The reason for this
is that in order to solve Eq. (4.9), we need to deal with a problem with variable mesh
Peclet number and Courant number. But in order to solve Eq. (4.19), we need only
to deal with a problem with constant mesh Peclet number and Courant number. On
the other hand, the numerical solvers presently available are more stable and robust
when they are used to solve Eq. (4.19), instead of Eq. (4.9).
Note that the equivalent source/sink term presented here has a very clear physical
meaning. For a representative elementary volume/area, a change either in the poros-
ity or in the pore-fluid velocity is equivalent to the addition of a source/sink term
into the representative elementary volume/area. Since the constants involved in both
the advection and dispersion terms are basically independent of time, it is possible to
use the initial mesh and time step, which are determined using the initial conditions
at the beginning of computation, throughout the whole process of the transient anal-
ysis of fluid-rock interaction problems. As a result, the finite element method with
the proposed term splitting algorithm can be efficiently used to solve both the first
and second type of transient transport equations (i.e. Eqs. (4.8) and (4.19)) for fluid-
rock interaction problems in pore-fluid saturated hydrothermal/sedimentary basins.

4.3 Application Examples of the Term Splitting Algorithm

In order to illustrate the usefulness and applicability of the newly proposed concepts
and numerical algorithms, we have built them into a finite element code so that fluid-
rock interaction problems in pore-fluid saturated hydrothermal/sedimentary basins
can be solved effectively and efficiently. As shown in Fig. 4.1, the application exam-