Software

MAFIC Software Information (continued)

The integrals which form the matrix coefficients A_nm and D_nm have been evaluated analytically using local element area (or length) coordinates as described by Pinder and Gray (1977) and Zienkiewicz (1977). This eliminates the need to perform numerical integrations. The matrices A and D are assembled locally on an element by element basis within the subroutines FELEM (for fracture elements) and MELEM (for matrix elements) and subsequently combined to form the global coefficient matrices. The matrices are calculated only at the beginning of each simulation to increase computational efficiency. Note that A and D are always symmetric and positive definite.

The right-hand-side vector contains the matrix product of the storativity matrix and the current head vector plus the prescribed nodal boundary conditions. Since these quantities generally change with each timestep they must be evaluated at the beginning of each timestep.

To conserve memory and reduce computational effort, when the matrix-block option is used equations for fracture nodes are separated from equations for matrix nodes. As shown in Figure 2-3, by placing terms expressed by Equation (2-3) which link fracture nodes and matrix nodes in the right-hand-side vector, the equations describing flow in the fractures are decoupled from those describing flow in the matrix. When solving the set of equations for fracture nodes the terms of the equations containing matrix block head values are evaluated using matrix node head values from the last iteration and placed in the right-hand-side vector. Similarly, when solving for matrix nodes, the terms containing fracture head values are placed in the right-hand-side vector. Addition of mixed nodal terms to the right-hand-side vector is performed in the subroutine MBFLUX. The iterative matrix solution procedure is repeated until the maximum head change for an iteration and the maximum unbalanced nodal flow are less than a user specified tolerance.

Convergence may be speeded following the first timestep by using a variable projection factor which estimates head values at the end of the current timestep by extrapolating head changes from the previous timestep. New head values for the current timestep are estimated once at the start of all but the first timestep. The new values are estimated using a simple linear extrapolation equation:

(2-5)

where: h_est = estimated head value at end of current timestep
h^k = head at the present timestep
h^k-1 = head at the old timestep
P = projection factor
D t_old = length of old timestep
D t_new = length of new timestep

P should always be positive and less than or equal to 1.0. When P equals zero no extrapolation is conducted. When P equals 1.0, the head gradient is extrapolated one timestep.

Modeling of Rock Matrix Flow

The diffusivity equation for three-dimensional flow in porous media can be written

(2-6)

where: S_s = Specific Storage (1/L)
h = Hydraulic Head (L)
t = Time (T)
K = Hydraulic Conductivity (L/T)
q = Source/Sink Term (1/T)
= Three-dimensional Laplace Operator

Equation (2-6) is the governing equation describing matrix flow in MAFIC. Equation (2-2) is used for the finite element solution after modifying the equation such that R is now the region volume (L³), q is the volumetric flux per unit volume (T^-1), and T and S are replaced by K and S_s.

Flow in the rock matrix may be modeled using either a fully discretized finite element solution or a generic matrix block approach. Each method offers certain advantages for particular types of problems. The relevant aspects of each are presented in the following sections.

The Fully Discretized Matrix Approach

Flow in the rock matrix may be simulated using a fully discretized rock matrix. This approach models the flow in the matrix using the Galerkin finite element solution scheme, described in Section 2.2. In order to use this approach, the rock matrix must be discretized into a finite element grid of tetrahedral elements (see Figure 2-4). Either linear elements or quadratic elements may be utilized; the program will automatically generate the additional nodes for the quadratic elements.

The fully discretized approach is suitable for small or simple problems involving only a few fractures. This method has two primary advantages. The first is that the exact geometry of the matrix is simulated. Second, this approach models flow through matrix interconnections between non-intersecting fractures. As the complexity of the fracture networks increases, however, the fully-discretized approach becomes cumbersome. One problem is the laborious task of designing finite element matrix grids for fracture systems containing many fractures of random orientation. Secondly, the computational effort required to solve a fully discretized matrix problem can quickly reach unreasonable proportions as the number of fractures increases.

The Matrix Block Approach

Considering the limitations of the fully discretized method, MAFIC provides an alternative method, utilizing generic matrix blocks, for simulating flow in the rock matrix. A similar approach has been used for conventional 'dual-porosity' simulations (Kucuk and Sawyer, 1979; Warren and Root, 1963; Streltsova-Adams, 1978). The matrix block approach defines a single generic matrix block geometry and size to represent all matrix blocks isolated by a network of surrounding fractures. Generic matrix block geometries used in the simulations include spheres and rectangular slabs. An example of how a spherical matrix block might compare to an actual fracture network is illustrated in Figure 2-5a. To proportion matrix block flow to individual fractures, a section of a matrix block with an outer surface area equal to the individual fracture areas is associated with each fracture element. For matrix block spheres, the matrix block section consists of a conical-shaped section, and for rectangular slabs the surface area of the drainage face is set equal to the fracture area. These relationships are shown in Figure 2-5b.

Flow through each matrix block is treated as a separate one-dimensional problem coupled to the average hydraulic head in the associated bounding fracture element. Two different methods are used to solve for flow within the matrix blocks, the pseudo steady-state approach and the finite element approach.

The pseudo steady-state approach (Warren and Root, 1963) assumes that flux into or out of the matrix block equals the change in storage within the block. Under this condition:

(2-7)

where: Q_mf = matrix to fracture flux (L³/T)
A_m = area of the matrix block (L²)
V_m = volume of the matrix block (L³)
K_m = hydraulic conductivity of the matrix block (L/T)
S_sm = specific storage of the matrix block (1/L)
h_m = average matrix block head (L)
h_f = average boundary fracture head (L)
t = time (T)
a = a geometric shape factor = 5/r_mb for a sphere; 3/r_mb for a slab (1/L)
r_mb = the matrix block radius or slab thickness (L)

This approximation is strictly valid only for later times when induced pressure transients at the block surface reach the center (sphere) or opposite end (slab) of the matrix block, but offers the advantage of computational simplicity.

[more]