Software

MAFIC Software Information (continued)

The finite element approach models one-dimensional transient flow within the matrix block as governed by Equation (2-6). Although this approach uses the same Galerkin finite element method as presented previously, advantage is taken of the resulting tri-diagonal matrix. A Thomas algorithm (Pinder and Gray, 1977) is used to solve for the matrix heads.

In order to retain this computational advantage, linear one-dimensional elements are always used to discretize the matrix blocks regardless of whether linear or quadratic elements are used to discretize the fracture network. Note that this approach is essentially identical to the pseudo steady-state approach when a single element is used to discretize the matrix block. For this reason the pseudo steady-state approach is always used when a single element is used to discretize the matrix block and the finite element approach is used when more than one element is prescribed.

An iterative solution is used to obtain both fracture and matrix heads in a manner analogous to that used for the fully discretized matrix approach. Fracture nodal heads are solved after calculating the flux from the matrix blocks (using matrix head values from the last iteration) and placing it in the right-hand-side vector. For the pseudo steady-state approach this flux is given by the left hand side of Equation (2-7). For the finite element approach this flux is determined by Equation (2-3), where n is an adjacent fracture node and m is the matrix node closest to the matrix block boundary. Matrix nodal heads are subsequently solved using the new fracture head estimates. For the pseudo steady-state approach, the average matrix block head is obtained by solving for m in Equation (2-7). Solution convergence is obtained in the same manner as previously described for the fully discretized matrix approach.

The matrix block equations are formulated in the subroutine MBLOCK. Solution of matrix heads and calculation of matrix to fracture flux is performed in the subroutine MBFLUX. Examples of the discretization used for these two methods are shown in Figure 2-6.

As mentioned above, the matrix block method requires defining a generic matrix block characteristic of the simulated region. Several options are available for simulating various fracture geometries. Matrix blocks may be defined as slabs or as spheres. Slabs are suitable for predominantly parallel fractures and spheres are suitable for networks composed of orthogonal fracture sets or randomly oriented fractures.

In addition to specifying the geometry of the matrix blocks, the user must also specify the radius or thickness of the generic matrix block. The size of the generic matrix block should be specified such that: 1) the total surface drainage area of the matrix blocks equals the available fracture surface area, and 2) the total volume of the region equals the total volume of the matrix blocks. When using spheres, for example, the radius should be chosen such that:

(2-8)

and

(2-9)

are both satisfied, where:

n = number of matrix blocks
rmb = radius of matrix block
Afr = total fracture area contained in simulation region
Vr = volume of simulation region

In Equation (2-8), total surface area is twice the total fracture area because each fracture will usually have flow from both sides. Solving these equations for rmb results in the following expression:

(2-10)

(2-11)

Although the matrix block radius could be calculated automatically by the computer program it has been left as a user specified parameter to facilitate accommodation of special circumstances. For example, any fractures which lie on the boundary of the simulation region will only have flow from one side and the area equation will have to be adjusted accordingly to account for the reduced total fracture area.

Using a generic matrix block means that the matrix block associated with every fracture has the same geometry, length, and flow parameters. This approach, while simplistic, offers several computational advantages over the fully discretized approach:

  • Eliminates need to create complex matrix grid
  • Reduces matrix problem to one-dimensional finite elements (tri-diagonal matrix solution), or analytical solution
  • Allows a finer discretization near the fracture/matrix boundary without large increases in computation effort
  • Reduces storage and computation requirements since all matrix blocks have an identical coefficient matrix.

The matrix block method is conceptually applicable for situations in which flow occurs predominantly in the fracture network. Consequently, this approach is ideally suited for fractured crystalline rock, shale, low permeability sandstone, and surficial basalt. These materials commonly contain fractures with conductivities much larger than the matrix conductivity and fracture storativities much smaller than the matrix storativity. Therefore, although most of the water is held in the matrix, the high fluid velocities in the fractures account for the majority of flow. The matrix block method does not account for flow between the matrix blocks; thus, it will not effectively simulate any situation where matrix flow dominates fracture flow. The dominance of the matrix flow increases as the hydraulic conductivity of the matrix approaches the conductivity of the fractures, and/or as fractures become less interconnected. If matrix flow is dominant, however, the fractures may be ignored. MAFIC may be used to model matrix flow alone (i.e., as a standard porous media simulator) if fully discretized matrix elements are used and no fracture elements are specified.

Matrix Storage

The matrix equation solver only stores terms in the matrix which are non-zero. Terms representing fixed-head nodes are not stored in the main matrix.

Nodal Groups

Nodal groups are a user-specified option which may be used to identify nodes with common characteristics. For example, a nodal group designation may be used to identify nodes having identical boundary conditions. In addition, a nodal group specification is required to apply time-varying flux or head conditions to a node or group of nodes. Since MAFIC calculates the total flux into each nodal group, a nodal group designation allows the user to observe the flux into any collection of nodes, including a group of boundary or well nodes. Furthermore, if summary output is specified (see Section 3.4) only results for nodes in nodal groups are written to the output file; thus, the user may create brief output files which provide information only on certain nodes.

Incorporation of Flow Boundary Conditions

Three different types of nodal boundary conditions are used in MAFIC: 1) specified head, 2) specified flux, and 3) specified nodal group flux.

Specified head boundary conditions may be either constant or time varying. The specified head condition at a node is implemented by eliminating the corresponding matrix row and column from the coefficient matrices. This is accomplished by eliminating the corresponding nodal equation (matrix row) and transferring the terms containing the fixed head nodes (matrix column) to the right hand side vector. The specified flux boundary condition is implemented by adding the flux term to the right-hand-side vector in Equation (2-4).

The group flux boundary condition allows for the specification of the total flux from all nodes within a nodal group. This allows for the simulation of a commonly occurring well boundary condition whereby the total flux into or out of the well is known but the production of individual nodes within the well is not (see Figure 2-7). If it is assumed that all nodes in the nodal group have the same hydraulic head value (i.e., head loss along the well is negligible) then the individual equations for each node in the nodal group can be combined into one expanded equation (i.e., for matrix solution purposes all nodes in the group are treated as a single node). In MAFIC this expanded equation replaces the equation for the first node in the nodal group and the equations for all other nodes in the nodal group are eliminated. The specified flux boundary condition for the combined nodal equation then becomes the total flux from the nodal group.

Time-varying prescribed nodal boundary conditions may be specified for nodes belonging to nodal groups. Time-varying boundary conditions are specified for each nodal group as chronologically tabulated data pairs of elapsed time and boundary value. For each timestep, the corresponding boundary value is determined from linear interpolation of the tabulated values using the elapsed time at the end of the timestep, as shown in Figure 2-8. Interpolation is performed by the subroutine INTERP. Beyond the last tabulated time, the boundary condition is taken as constant and equal to the last tabulated volume for all subsequent timesteps up until the end of the simulation.

[more]