Boundary Element Formulation

Numerical modelling is an attempt to mathematically simulate the way the rock mass responds to mining. This is done by accommodating the effects of:

•Loading conditions - Prior to mining the rock mass is loaded by overburden and tectonic forces.

•Structural support elements - There may also be loads do to mechanical ground support, backfill, heating, etc.

•Geometry - When you make excavations, these pre-existing loads are redistributed around the excavations, concentrating in abutments and pillars. Stiff dykes and soft ore zones will also influence the response.

•Elasticity - At locations where the stresses do not exceed the strength, the rock deforms in a more or less elastic manner. This simply means that the deformations are relatively small in magnitude and are mostly recoverable.

•Non-linearity - At locations where the stresses are concentrated to the point where they exceed the strength, the rock will yield to these loads and deform. This means that the deformations can be relatively large in magnitude and are mostly non-recoverable. This can include the effects of fault slip, cracking and generalized 3D non-linear yielding.

How does numerical modelling work?

Numerical modelling achieves this simulation by using certain physical constraints on how the rock mass can respond:

Equilibrium – Applied forces must always balance one another at all locations in the model. If you cut out a small cube of material and examine the stresses acting on this cube these stresses must be in equilibrium.

∂σxx /∂x + ∂τxy /∂y + ∂τxz /∂z + X = 0

∂τxy /∂x + ∂σyy/∂y + ∂τyz /∂z + Y = 0

∂τxz /∂x + ∂τyz /∂y + ∂σzz /∂z + Z = 0

Note that terms X, Y and Z, represent body forces that can be used to apply any sort of external field loading. Such field loading can be the result of heating, fluid pressures material non-linearity etc. In addition, these can be determined from many forms of in situ monitoring including for example fluid pressures (e.g. well drawdown, dams, hydrofracturing), heating (e.g. natural heating, nuclear waste storage), deformations (e.g. monitored with extensometers), and seismic activity (e.g. definition of fault slip or material non-linearity from seismic data). It is these terms that will be used for the integration of numerical modelling with seismic monitoring in Map3Di.

Continuity – In the rock mass continuum, the mass of material must be maintained. You cannot have material disappearing or being created.

∂²εxx /∂y² + ∂²εyy/∂x² = 2 ∂²εxy/∂x/∂y

∂²εyy /∂z² + ∂²εzz/∂y² = 2 ∂²εyz/∂y/∂z

∂²εzz /∂x²  + ∂²εxx/∂z² = 2 ∂²εzx/∂z/∂x

∂²εxx /∂y/∂z = ∂/∂x ( -∂εyz /∂x + ∂εzx /∂y + ∂εxy /∂z )

∂²εyy/∂z/∂x  = ∂/∂y ( -∂εzx /∂y + ∂εxy/∂z + ∂εyz /∂x )

∂²εzz /∂x/∂y = ∂/∂z ( -∂εxy/∂z + ∂εyz /∂x + ∂εzx /∂y )

Elasticity - At locations where the stresses do not exceed the strength, the rock deforms in a linear elastic manner: stresses varying in direct proportion to the elastic strains.

σxx  = σ°xx + [ (1-υ) εxx + υ ( εyy + εzz  ) ] E/[(1+υ)(1-2υ)]

σyy = σ°yy + [ (1-υ) εyy + υ ( εzz + εxx  ) ] E/[(1+υ)(1-2υ)]

σzz  = σ°zz  + [ (1-υ) εzz  + υ ( εxx + εyy ) ] E/[(1+υ)(1-2υ)]

τxy  = τ°xy + εxy E/(1+υ)

τyz  = τ°yz + εyz E/(1+υ)

τzx  = τ°zx  + εzx E/(1+υ)

Note that this is where the pre-mining or far-field stresses (σ°) are incorporated into the physics.

Non-linearity - At locations where the stresses are concentrated to the point where they exceed the strength, the rock will yield to these loads and deform. Deformations are all allowed to proceed until the stresses relax down to the strength. This may be accompanied by some dilation. This can include the effects of fault slip, cracking and generalized 3D non-linear yielding.

Note that the equations of equilibrium and continuity are expressed as differential equations. What we need to do to solve these is to integrate them over the rock mass volume such that the appropriate boundary conditions are satisfied. There are many ways of accomplishing this. Boundary Element Methods integrate the equations analytically, and then use a numerical approximation to satisfy the boundary conditions. Finite Element Methods and Finite Difference Methods use a numerical integration scheme to integrate the differential equations.

No matter which method is used, in the end they all use some sort of numerical approximation and you end up with a large set of equations that describes how various parts of the rock mass interact upon one another. In BEM the resulting equations appear as

σsurface = σfar_field + M P

where M is the set of simultaneous equations and P represents the loads or deformations that need to be applied to cause the stresses at the excavation surfaces to be zero. This constitutes a mathematical description of how the rock mass responds. All numerical models use some variation on this approach.

These equations need to be solved simultaneously such the boundary conditions are satisfied. By solving these equations throughout the rock mass, along with the loading conditions and geometry, you conduct a "stress analysis".

Map3D is based on a very efficient Indirect Boundary Element Method (Banerjee and Butterfield, 1981), and incorporates simultaneous use of both fictitious force and displacement discontinuity elements. Special proprietary boundary elements are incorporated for the thermal and non-linear analysis versions.

This Boundary Element formulation offers many advantages over other stress analysis techniques. Direct Boundary Element formulations require approximately twice the computational effort to assemble and solve the boundary element matrix, compared to the indirect method used in Map3D.

Matrix lumping techniques are used to reduce the matrix size, and as a result, very large problems (more than 333333 elements - 1000000 degrees of freedom) can be accommodated on PC platforms. Without lumping, a 64000 element problem would take 147 GB of disk space. With lumping, this can be reduced to less than 1 GB. This permits users to specify existing mining geometry in detail, and add new mining as required. This greatly reduces the effort required to set up and run analyzes and permits the whole mine to be considered when necessary.

With the Boundary Element formulation, one starts with an infinite homogeneous elastic medium (rock mass). The process of model building consists of making excavations and superimposing non-homogeneous zones (dykes, ore zones or yielding zones) and any faults or joints upon which slip may occur. Since one starts with an infinite medium, far field boundaries are automatically accommodated.

This is unlike domain formulations such as the Finite Element or Finite Difference methods, where one starts with empty space. For these latter methods, model building consists of assembling the entire rock mass and all its components. Elements must be assembled out to some far field boundary many diameters away from the excavations.

While Finite Element and Finite Difference formulations are very well developed for non-linear problems such as plasticity, transient heat and fluid flow, and dynamic simulations, there are also very difficult to use and require very long run times. The difficulty in use comes from the large amount of information required to discretize the host rock mass. While this can be minimized by use of mesh generators and powerful front-end graphics, the cost and effort to learn the interface outweighs the simplicity with which the Map3D Boundary Element analyzes can be conducted.

These fundamental differences make the boundary element method much more suitable and economic (in terms of analysis time) for rock mass problems.

Many of the supposed limitations of the boundary element method have been overcome in Map3D. The program can accommodate multi-step mining sequences and multiple material zones with different material properties and stress states. These zones are permitted to behave non-linearly. Multiple intersecting fault planes can slip and open according to user specified shear strength. Steady state thermal/fluid stress analysis can be simulated.

The built-in solid modelling technology makes model construction very straightforward by building complex intersections internally. The discretization routines built into Map3D relieve the user of the burden of attempting to optimize the use of elements by automatically concentrating them only where results are requested. The internal lumping procedures reduce both matrix size and computational effort to a nearly linear dependence on problem size.