of the stresses does become important, and it is essential
to work with total stresses during the equation-solving
process.

Following the notation used by Sinha (2), the funda-
mental boundary-element equation relating elemental
stress components (T) to elemental displacement discon-
tinuities (D) is

t- > [w+ |p
mey

a4)

where [TC] and [TC"] are m x n matrices of displacement
coefficients for the infinite body and half-space solutions,
respectively, and T and D are n dimensional vectors. In
principle, equation 14 is identical to equation 5 given in
the prior simplified explanation of BEM workings.

Induced-Stress Boundary Conditions

For an induced-stress formulation, the boundary con-
ditions for the mined-out areas are

(T] + [P] = [0] (as)
and for the elastic areas are
[T] + [K}7[D] = [0] (16)

where P is the far field primitive stress and K is the ele-
ment stiffness. For linear elastic elements, the stiffnesses

.K are constants and do not depend on the elemental
displacements D or the induced-elemental stresses T in
any respect. Combining equations 14, 15, and 16 results in
a linear system of equations of the form

Ax=b. (7)
Direct solution of this system is not feasible since the
unknown elemental displacements D must satisfy the ad-
ditional restriction

[D] = [t] (as)
where t is a constant vector of elemental thicknesses. This
restriction necessitates use of an indirect method like the
Gauss-Seidel scheme provided by Dahlquist and Bjork (33)
for n = 1,2,... N:

-1
-Da

m=1

28? — aga? -

me=n+l1

where n is the current unknown index, N is the maximum
number of unknowns, k is the iteration number, m is a
summation index, and w is the overrelaxation factor. w is
defined by the user and typically ranges from 1.0 to 1.5.
A new variable for summation of element stress influences,
T, is defined as

n-1
m-- Y [ Fn + Tq [DEP
m=1

- [Fem + Fa | De

N
- 2

m=n+1

[ Fm * Term | D®. 2)

For the mined-out areas, combining equations 14 and 15
and then forming a recursive relationship after 19 results
in

-P,
| 1)

D&D 2 pH 4 | a
Hy Hy =
TCan * TCin

Similarly, in the elastic areas, combining equations 14 and
16 results in

-kK, ©

pD&) -p® +0 lL (22)

TCan * TCan * Kn

Basically, Sinha’s original MULSIM, and its enhancement
MULSIM/BM, implement and solve equations 20, 21, and
22 for a linear elastic induced-stress approach. Figure 14,
modified from Sinha (2), provides a flowchart for this
critical area of the program - namely, the equation solver.
These three equations are noted in this flowchart.

Total-Stress Boundary Conditions

With nonlinear materials, however, the elemental
stiffness K depends on the absolute magnitude of the
elemental displacements and stresses. (In fact, K becomes
a function of D.) It is therefore no longer sufficient to

Se)
DL amakm + bp

OD LO Ly

ann

(19)


approaches zero (i.e., the number of mining steps used to
create an excavation becomes a continuum), the kinetic
energy release component becomes zero, while the strain
energy release component approaches the total energy
release.

The step-size effect tests seek to observe the theoretical
energy release relationships with numerical experiments
using MULSIM/NL. Figure 36 shows the four meshes
used for the experiments. In each of these tests, an
800-m-wide opening is created with an increasing number
of effective mining steps; since, in each case, the number
of actual mining steps is always just two. The size of this
second step decreases such that its size is the same size as
the last step of an actual 4-, 8-, then 16-step model.

Figure 37 plots the kinetic and strain ERR normal-
ized to the total energy release rate versus the number of
steps used to excavate the opening. The figure shows ana-
lytic derivations by Salamon (40) for the excavation of a

‘Step 1
rStep 2
Pz ‘483
i sac
A 2-step
rt model
i
i
‘Step 4
[ rStep 2
56
apes 4-step
sg model
ster, 1
Step 2
7
7] 8-step KEY
o A 98
Ea Ei model Fine- mesh
i Elastic
‘ 2; area
4
Step KN open area
rStep 2
—~_!
46 50-m blocks
st 30 by 30
ge 4 |6-ste coarse mesh
3 4 mode! 3,750-MPa
mei modulus
a I-m thickness

Figure 36.—Four MULSIM/NL grids for step-size effect tests.

spherical cavity. When that cavity is created in a single
step, half the total energy release rate comes from the ki-
netic component and the other half derives from the strain
energy component. Increasing the number of steps used
to create the spherical cavity causes the strain ERR to
dominate the total ERR, while the kinetic ERR tends to
zero. Again, the total ERR remains constant.

Figure 37 also shows the numerical calculations with
MULSIM/NL for the excavation of an 800-m-wide tabular
opening. In this example, the kinetic ERR dominates the
total ERR. For practical purposes, the strain ERR is neg-
ligible. However, the same underlying trend still exists
such that increasing the number of mining steps used to
create an excavation causes the kinetic ERR component to
decrease relative to the strain ERR component, while the
total ERR remains constant.

NONLINEAR EFFECTS ON ENERGY
CALCULATIONS

The ERR derivations showed that total ERR is com-
prised of a strain ERR component and a kinetic ERR
component, which, in turn, is composed of a linear and a
nonlinear part. The following studies examine these three
components to discern what factors affect their relative
magnitude, hence, relative importance in the total ERR.
These parametric studies examine, systematically, the fol-
lowing important factors using several different small
MULSIM/NL models:

(1) The size of the open area relative to the total
backfill area in the model.

(2) The amount of open area in the whole model
which must either increase, decrease, or remain constant
between mining steps.

(3) The average stiffness of the backfill material.

(4) The degree of nonlinearity of the backfill material.

10

ase ~~
We key
« 3 — == MULSIM/NL solution
£ for longwail panel
= — Analytic solution for
x 6 Um/Wr spherical cavity
3 Wr Total energy release
2 | rate
2 W/ Wr Wk Kinetic energy release
4 rate
e. | Uy. Strain energy release
rote
Uy/ Wr ~
ol 4 saath as
V 2 4 e 16

NUMBER OF STEPS USED TO
EXCAVATE OPENING

Figure 37Normalized energy release rate components ver-
‘sus number of excavation steps for analytic case and numerical
case.


widths. For these models, stresses agree to within about
5% of analytic. Displacements can also show reasonable
agreement with analytic (within 5%) depending on how
well the two-dimensional solution approximates the model.

The important point of this exercise is that the choice
of clement size can have a profound effect on stress dis-
placement calculations. Figures 29 and 30 show that near
an opening boundary, the clement size can affect displace-
ment calculations by over 30%, stress calculations by al-
most 50%. The choice of element size also impacts the
rigid boundary effect. In fact, the two are inseparable.
Smaller element sizes which tend to improve the accuracy
of stress and displacement calculations also tend to show
deleterious rigid boundary effects because of the smaller
model size. In summary, because of these somewhat fickle
element size and rigid boundary effects on stresses and
displacements computed with MULSIM/NL, caution is re-
quired in all comparisons especially when different ele-
ment sizes and different model geometries are involved.
Again, as shown in figures 29 and 30, large differences can
arise that may be artificial.

35

30F

STRESS (c), MPa

—— Analytical solution

°
950

1,000

__DISTANCE, m _

Figure 30—Computed and analytic stresse:
opening with different block and element sizes.

1,050

1,100

1,150

edge of 100-m

25-m opening 100-m opening 400-m opening
35 ey 35 ; r 35 7 ,
+ -—--12.5-m block,
eg * 2'5"m elements 50
--50-m block, 4
= 28 10m elenients 25
rs 20 —— Analytic solution 4 20
Q 1s N\ 4 15
Yio \ H 10
5 i 4
QD 5 \ H 5
° roan L °
: O15 —— 16
= 14
2 12
5 1.0
= 8
8 6
2 4
ao 2
a _/ RA} og . i
1000 1050 1100 150 800 1,000 1200 1400
DISTANCE, m

Figure 31.—Comparison of computed and analytic stresses and displacements Illustrating element size effects.


displacement (convergence) changes as opposed to total or
absolute stresses and displacements.

The energy calculations constitute the last unique
feature of MULSIM/NL. Salamon (6) provides the theo-
retical basis for the comprehensive energy release rate
(ERR) computations performed by MULSIM/NL. The
total energy release is comprised of three basic terms:
(1) a strain energy release from the mined-out material,
(2) a linear kinetic energy release due to the change in
gravitational potential over the mined-out area, and (3) a
nonlinear kinetic energy release due to nonlinear materials
in the total backfill and/or gob area. All terms are in-
cluded in the implementation; however, the strain energy
release term is modified slightly to account for nonlin-
ear stress-strain behavior in the unmined materials,
MULSIM/NL also computes various useful strain energy
quantities for each unmined element throughout the BEM.
model.

Comparisons of numerical ERR calculations from
MULSIM/NL with known analytic solutions show excel-
lent agreement. The ratio between open-area-to-backfill
or gob area influences the relative magnitude and im-
portance of the three energy release components. With
large open-area-to-gob-area ratios, the kinetic energy re-
lease component (change in gravitational potential) dom-
inates the total energy release. With the small ratios that
are more typical of most practical BEM models of mines,
then all three components of the energy release (strain,
linear kinetic, and nonlinear kinetic) have importance.

These extensive enhancements necessitated complete
restructuring and reprogramming to create MULSIM/NL.
Earlier versions of MULSIM used complicated indexing
and one-dimensional arrays throughout the program. The
new MULSIM/NL takes advantage of today’s larger mem-
ory computers and uses numerous multidimensional arrays
and vectors that better reflect the underlying mathemat-
ical foundations of MULSIM/NL. The reprogramming
resulted in a much more understandable, reliable, and
modifiable program that should serve users well in the
foreseeable future.

Development of MULSIM/NL stems from a need to
alleviate safety hazards in certain U.S. coal mines plagued
with coal mine bumps. With its six nonlinear in-seam
material models, multiple mining step capabilities, and
comprehensive energy release and strain energy calcu-
lations, MULSIM/NL provides a means to calculate stress,
displacement, and energy changes for various mining con-
figurations in bump-prone or rock burst conditions. In

addition, MULSIM/NL has many practical uses beyond
coal mine bump research. Stresses and displacements cal-
culated by MULSIM/NL can help an engineer design pil-
lar sizes and understand multiple-seam mining interactions.
This program is one tool out of many that engineers can
use in the fight to alleviate coal mine bumps and rock
bursts. Through judicious use of numerical models such
as MULSIM/NL, engineers can choose mining methods
that decrease the risk of these violent failures.

RECOMMENDATIONS

Many avenues exist to improve MULSIM/NL. The
three possibilities suggested here may be especially useful
in the near term,

First, MULSIM/NL could benefit from further restruc-
turing and reprogramming to take advantage of modern
vectorizing compilers and future parallel computers. One
way to render MULSIM/NL highly vectorizable is to elim-
inate use of the coarse mesh,

The second suggestion is to calculate a condition num-
ber for the system of equations solved by MULSIM/NL.
Operating experience with MULSIM/NL shows that
certain models lead to an ill-conditioned system of equa-
tions with slow solution convergence and an unreliable
final solution. It is very easy to create such misbehaved
models. A condition number for the system of equations
provides a simple, quantitative measure of the solution
reliability.

The third suggestion is to include topographic influence
in the boundary-element calculations. One approach to
these calculations would define an upper displacement dis-
continuity at the topographic surface. Each element of
this surface would require elevation data. The solution
procedure would then require calculation of semi-infinite
body influence coefficients between the topographic sur-
face and the main seam plane. Unfortunately, the com-
puted topographic influence coefficients would require a
lot of memory, and the computational effort may prove
very costly.

Another, less rigorous, approach would define in situ
stresses at all elements in the seam plane. The existing
boundary-element calculations can include these element-
specific, in situ stresses very easily. The problem then
becomes one of estimating the seam plane in situ stresses
that include the topographic influence. Horizontal seam
planes only require an estimate of the vertical stress com-
ponent, which considerably simplifies the problem.


energy release per unit area mined. (As will be discussed

later, the previous ERR relation corresponds to the domi-
nant term in an ERR expression derived and discussed by
Salamon (6).)

In addition to ERR, Itasca also provided calculations
for the energy release due to element modulus softening
given by

Was = 3 (41 + 2)) (P= Dia), @

where Was = energy release due to modulus soften-
ing

and o; = element stress in current mining step.

With the multiple mining step feature, the user can ap-
proximate a yielding or softening behavior in the seam
materials by reducing the linear elastic modulus from a
high value in the (i-1) step to a low value in the i step.
Such a modulus reduction causes the energy release given
by equation 2 and shown in figure 4 (bottom). The theo-
retical basis of this procedure seems problematic. Fur-
thermore, Ways is also step size and path dependent. All
theoretical problems aside, the energy release due to mod-
ulus softening does provide a useful approximation to a
nonviolent component of the energy dissipation occurring
due to mining.

NONLINEAR MATERIAL PROPERTIES IN MULSIM/NL

GENERAL MECHANICS OF BOUNDARY-
ELEMENT METHOD

This section provides a simplistic discussion on the ba-
sics of how BEM work, beginning with certain mathemati-
cal concepts and analytic solutions central to the BEM.
These fundamental analytic solutions lead to the formation
of a large system of algebraic equations. Boundary con-
ditions applied to this system of equations then fully define
the particular problem at hand, Finally, the equations are
solved for stresses and displacements at each element.

Numerical methods to solve complex problems in elas-
ticity fall into two broad categories: differential methods,
such as finite-element and finite-difference methods, and
integral methods, such as BEM. With differential meth-
ods, the entire problem domain (i.c., the rock mass) is
divided into elements or node points, and certain boundary
conditions are specified on the edges of the problem do-
main, Stresses and displacements are then calculated at
each element or node point in the subdivided rock mass.
In contrast, integral methods only require subdivision of
the problem boundary into elements; therefore, the in-
tegral methods require much less computational effort
than the equivalent differential methods. Stress and
displacement calculations occur along the boundary and at
requested points in the interior. In many practical rock
mechanics problems, boundary stresses and displacements
are the only ones of interest; therefore, integral methods
are well suited.

The BEM applies quite well to the stress analysis of
thin tabular deposits such as coal seams, deep gold reefs,
and narrow veins. Sinha (2), and later Wardle (14), dis-
cuss that the key abstraction underlying BEM is visualizing

the tabular deposit as a crack or discontinuity in an other-
wise homogeneous, isotropic, linear elastic rock mass. Fig-
ure 5 shows this crack or seam plane in the surrounding
rock mass or host media. The top and bottom surfaces of
the crack plane form the problem boundary (i.c., the roof
and floor in coal mining or the hanging wall and footwall
in metal mining). The rock mass is the internal region for
this problem, whereas the coal seam or vein material is an
external region for the problem.

The next critical step in the BEM is dividing the crack
or seam plane into several square boundary elements.
Figure 6 shows the seam plane subdivided into elements.
The actual mine plan within this seam plane may have a
very complex geometry. The individual clements are then

Earth's surface
A

Depth

Figure 5.—Boundary-element-method seam plane in surround-
ing rock mass.


APPENDIX E.—NONLINEAR STRESS-STRAIN RELATION FOR GOB

INTRODUCTION

This appendix develops a nonlinear stress-strain re-
lationship for gob material, which exhibits strain-hardening
behavior, and discusses implementation of that model
in MULSIM/NL. Gob formation occurs via a two step
mechanism: first, the overburden rock behind a longwall
face fractures, collapses, and swells to fill the newly
created void, and second, the gob consolidates and follows
a nonlinear, strain-hardening, constitutive relationship as
vertical stress on the gob increases. Derivation of the
model begins by assuming that the tangent modulus of the
gob increases linearly with stress from some initial value
to a final value. Based on this assumption, a nonlinear
stress-strain relation for gob results is:

oe a

where  E,_ = initial modulus,

Ep, = final modulus,
o, = final stress,
D = seam closure,

t = seam thickness,

and 2

W

gob height factor.

MECHANISM OF GOB FORMATION

Several geologic studies have examined rock move-
ments in the overburden behind an advancing longwall
face (42-44). Figure E-1 shows the three zones generally
found as one progresses from the mined-out seam level up
through the geologic column. First, a totally broken zone
exists for a height ranging from two to six times the seam
thickness. Rock within this zone exhibits considerable in-
elastic behavior, since it has first fractured, then moved
downward, and then rotated randomly during its downward
displacement. Above the totally broken zone is a fractured
zone, whose top may lie 20 to 60 times above seam thick-
ness. In this zone, considerable fracturing occurs with the
downward displacement and mechanical properties such as
modulus certainly decrease; however, gross rotation of the
fragments does not occur and individual rock strata remain
recognizable and correlatable, Finally, above the broken
zone lies an elastic zone where no new fracturing has oc-
curred and mechanical properties remain the same as the
original rock mass.

Subsidence data also provide information on the delin-
eation and behavior of these three zones. At risk of gross
oversimplification, surface subsidence generally amounts to
about 50% of seam thickness. Displacements at the base
of the elastic zone shown in figure E-1 may equal 50% of
the seam thickness as well. At the base of the fractured
zone, displacements are approximately 60% of the seam
thickness, Downward movement at the base of the totally
broken zone is 100% of the seam thickness.

Mechanistically, the overburden rock first swells to fill
the mined-out void, then consolidates to induce surface
subsidence. During this swelling and consolidation, the
totally broken zone sees a net vertical strain on the order
of (1.0t - 0.6t)/(4t), which equals 0.100 or 10% strain, The
fracture zone experiences a strain of (0.60t - 0.50t)/(40t)
which equals 0.0025 or 0.25% strain. The elastic zone sees
zero vertical strain. Since the strains in the broken zone
are much larger than the strains in the fractured zone,
then the mechanical behavior of the broken zone probably
dominates overall gob behavior.

In this research, the fractured zone is effectively ne-
glected by including it mechanistically with the clastic
zone. Gob formation occurs via the two-step mechanism
hypothesized in figure E-2. In the first step, the overbur-
den rock behind the longwall face fractures, collapses, and
swells to fill the newly created void. During this swelling
phase, the broken material is essentially stress free and

Surface

—>

2-6T

amas

—
©

Figure E-1.—Three mechanical zones found above mined-out
coal seam. Tis seam thickness. Lowestis totally fractured zone;
middle is fractured zone; and top is elastic zone.


APPENDIX C.—SUPPLEMENTAL DATA ON ELEMENT SIZE EFFECT TESTS

35 , 0.16 : : r 1
A A KEY
30 KEY, 4p en block . 4
- =m blo
IT Som Bloc — Analytic solution
— Analytic solution dab ee - 4
25 _
: ' tok —~ 4
20f i
! it ost 4
18 ifs \
yi ost \ 4
ol \ I \
L j oaf i
5 i
| i af
° it a es ae eo
985 995 1,005 1,015 1,025 1,035. ° 4 4
35 985 99500540258
T T T o4
B KEY B
30 12-m block 4
‘25-m block
f=" 50-m block H
25 f— Anolytical solution 4 ak 4
é q { é
= ry
2 5
g o2b 4
8 Fd
& 8
5 3
” a
5
Can |
— Analytica’
solution
Lo : Nee
iso %50 ooo 10504100, 150
16 T T T
c
14 4
J 2 4
10
{I 8
| er
i 4 KEY
! sone 50-m block
i 4 ~1-!= 100-m block
| 2| — Analytic
[ ‘solution
i ° 1
1200 14001600 800 000 ——«200~=S=«,AOO—=~S=«*GO
DISTANCE, m _ _ _ DISTANCE,m
Figure C-1.—Computed and analytic stresses at edge of open- Figure C-2—Computed and analytic displacements across

ing with different block and element sizes. A, 25 m;B,100m;C, opening with different block and element sizes. A, 25 m; B,
400 m. 100 m; C, 400 m.


This expression can be changed to

1 1 1
3Js,Tiaess = -3 foTouds + 5 fo Tues

1 1
= ARu,ds - = RyAud:
2 bo mes A vc

+

1 1
5 Js, T4 - 5 fXidudv. (16)

Finally, substituting A-5, A-6, A-9, and A-12 into A-1 and
simplifying with A-16 results in the following expression
for the total energy release:

1 1
Ware [El ds - ahi]

1 1
+ 5 fo,7ia ves + Zaye 4B uss} (A-17)

45

The first term in brackets is called the strain energy re-
lease from the mined material (Uy), and the second term
is called the kinetic energy release (W,). In most practical
circumstances, the body force term X; in the strain energy
release is negligible compared with the surface traction
term T,. The energy subroutine in MULSIM/NL does
neglect this body force term. In further notation, the first
term of the kinetic energy release is referred to as the
linear part of Wx and basically represents a change in
gravitational potential. The second term in Wy, referred
to as the nonlinear part of W,, arises from the nonlinear
stress-displacement behavior of the backfill and/or gob
material. In certain practical circumstances, it, too, is
negligible compared with the linear part of Wx (ie, the
change in gravitational potential), The energy subroutine
in MULSIM/NL does retain this nonlinear term. In addi-
tion, the relative value of the first and third terms (i.e., the
strain energy release Uy and the kinetic energy release
W,) is step-size dependent. In many practical cases, the
linear kinetic energy release may dominate the strain
energy release from the mined material. Again, that is not
always the case, hence MULSIM/NL retains the strain
energy release term for completeness.


In both cases, locating this intercept on the prescribed
stress-strain curve provides an updated modulus (or
stiffness) value for use in the next iteration of the equation
solver. Upon returning to the solver with this secant
modulus, the program continues the current iteration using
equation 27 to calculate updated displacements.

Material Property Subroutine Structure

MULSIM/NL uses a 26 by 10 array called EPROP to
store parameters used by the six nonlinear stress-strain
material models. The program permits up to 26 different
materials labeled A through Z which in turn can follow
any one of the six material models shown in figure 12.
EPROP allows each material model to have up to 10 de-
scriptive parameters. The value of the first parameter
directs the program to the proper stress-strain model for
each material defined. These are the assignments for the
six models currently permitted: 1 - linear elastic coal, 2 -
strain-softening coal, 3 - elastic-plastic coal, 4 - bilinear
gob, 5 - strain-hardening gob, and 6 - linear elastic gob.
Table 1 shows the required parameters and their location
in EPROP for the six material models. Quite naturally,

Table 1.—Material model paramete!

when defining materials for a given model and reading the
input file, MULSIM/NL expects the material descriptors
in a compatible format.

The six material property subroutines have very similar
structures. Associated with each element or block is a ma-
terial code, A through Z, which corresponds to positions
1 to 26 in EPROP. When MULSIM/NL requires updated
material properties, it reads the first parameter in EPROP
for that material which, in turn, directs MULSIM/NL to
the proper material property subroutine. Each subroutine
can then read the remaining parameters in EPROP for
that material and interpret those parameters accordingly.
As discussed earlier, each subroutine returns updated se-
cant moduli for the normal direction of the element and
the two tangential directions.

Finally, each nonlinear material property subroutine
incorporates a convergence check. The old total stress
values from the initial stage of the iteration step must
equal (or nearly equal) the new total stress value com-
puted by the material model subroutine. If the old and
new total stresses agree within an acceptable user-defined
limit, then the solution has converged. Otherwise, the
iteration steps continue up to a user-defined maximum,

and material property array EPROP structure

Model 1 2 3 4 5 6

1... Linear elastic Young's Shear modulus NAp........+ NAP. os... ee NAp.
coal-(1). modulus (G).

2... Strain soft- Peak stress. Peak strain... Residual stress Residual strain Poisson's ratio
ening coal ().

3. | rr Plastic modulus Poisson's ratio NAp.

).

4. Offset Offset strain Hardening mod-..do ....... Gob height
hardening _ stress. ulus. factor (n).
gob-(4).

5... Strain. Initial mod- Final modulus Final stress (c). Gob height Poisson's ratio
hardening ulus (E). (EJ. * factor (n). W.
gob-(5).

6... Linear elastic Young's Shear modulus Gob height fac. | NAp........ NAp.
gob-(6). modulus (6). tor (n).

Nap Not applicable.

NOTE.—For these 6 material models, there are no parameters 7 through 10.


Constant open area tests Decreasing Increasing
Small area Middle area Large area area test area test
A y
Step 1 oe i Be,
ZN
Step 2 me 8 LR, toe 23
Pp Zz in

Step 3 fa a a as 23

ZA Elastic area

KEY
Fine-mesh area

(SX Backfill (gob) area

GB Open area

Figure 38.—Five different three-step MULSIM/NL models used to study constant, decreasing, and increasing area effects on total

energy release and its components.

Figure 38 shows the five, three-step MULSIM/NL grids
used for these studies. The first three meshes have a
constant open area that moves in steps across the grid.
These three meshes examine final ratios between open
area and backfill (gob) area of 0.05 (small), 0.10 (middle)
and 0.33 (large). The last three meshes vary the amount
of open area between steps two and three. In the large,
constant area test, of course, that open area amount
remains constant. In one mesh, that area decreases, and
in the other, it increases.

These models compare two basic materials for the
backfill. The first is simple linear elastic, with a 70-MPa
modulus. The second is strain hardening with a 7-MPa
initial modulus and a 560-MPa final modulus at 7 MPa
stress. Figure 39 shows the stress-strain curves for the
linear and nonlinear backfill materials used in these mod-
els. The figure shows typical values only; however, during
these studies, numerous other linear and nonlinear mate-
rial parameters received consideration.

Figures 40 and 41 show the computed total ERR and
its three components (linear kinetic, nonlinear kinetic, and
strain energy release rate) for the five models shown in
figure 38. The calculations in figure 40 use the linear
backfill properties whereas those in figure 41 use the
nonlinear strain-hardening model.

STRESS (co), MPa

Figure 39.—Linear and nonlin
fill material used in model stud!

o

+

n L 1
0.02 0.04 0.06 008 O10
STRAIN (€), m/m

stress-strain curves for back-


CONTENTS—Continued

Summary and recommendations
Summary
Recommendations

References . .

Appendix A.—Rederivation of Salamon’s energy release equation oe

Appendix B.—Supplemental data on rigid boundary effect tests .

Appendix C.—Supplemental data on element size effect tests . .

Appendix D.—Supplemental data on nonlinear effect tests .

Appendix E.—Nonlinear stress-strain relation for gob . .

ILLUSTRATIONS

1, Multiple-seam grids for program MULSIM (2) . . .
2. Programming structure of MULSIM (2) . .

3. Modified program structure of MULSIM wit mul ipl le , mining steps and linear energy calculations
4. Stress-displacement curves showing strain energy release

5. Boundary-element-method seam plane in surrounding rock mass . .
6.

1.

8.

Boundary-element-method seam plane subdivided by elements ......
Actual mine plan in-seam plane and gridded modeling approximation .
Pressurized crack problem in two dimensions .......... eens

9. Displacement-discontinuity problem in two dimensions

10. Source element and field element within seam plane .

11. Three boundary conditions applied to elemental displacement discontinuity

12. Six different stress-strain models available in MULSIM/NL

13, Assumed linear relationship between tangent modulus and stress over range 0 to oy «eee eee ee eee

14. Flowchart of iteration scheme using induced-stress approach to solve for unknown stresses and

displacements 2.0.0.0... 00 e cece cece e eee eee eee

15. Flowchart of iteration scheme using total-stress approach to solve for ui unknown stresses and

displacements 6.0.6.6. 6c eect e eee tener ene e ee en eee bebe e nee enene

16. Two-stage iteration process for nonlinear material models in MULSIM/NL

17. Basic notation for energy calculations as mining progresses from state I to state II . :

18. Strain energy calculations for various linear and nonlinear material behaviors showing how parameter a

accounts for degree of nonlinearity

19. Step-size dependence of energy release components for radial expansion of circular tunnel .

20. Basie strain energy relationships for linear and nonlincar material models in MULSIM/NL that apply to.

unmined seam material

21, Three special cases of mining progress from state I to state II showing how elements can undergo six

different status changes between state I and state II

22. Simplified flowchart for subroutine ENERGY showing six possible logic paths .

23. Elliptic crack of width 2c subject to stress P perpendicular to crack plane .

24, MULSIM/NL models for rigid boundary effect studies of infinitely long planar crack of width 2c

25. Computed and analytic displacements across 100-m opening with different coarse-mesh sizes . .

26. Computed and analytic stresses at edge of 100-m opening with different coarse-mesh sizes ........

27. Comparison of computed and analytic stresses and displacements illustrating rigid boundary effects .

28. MULSIM/NL models for clement size effect studies

29. Computed and analytic displacements across 100-m opening with different block and element sizes

15
16
19

21

Ra

27
27
27

29
29


25-m opening
6.25-m blocks, 1.25-m elements

100-m opening
25-m blocks, 5-m elements

400-m opening
100-m blocks, 20-m elements

35 —

KEY
30F ---Middie-size coarse
L mesh (20 by 20cm)

— Analytic solution
20F H 4

STRESS (co), MPa

T T 35 T T

Sr) ‘
4 20h

tJ

!
———|

DISPLACEMENT (D), m

1.6 r T
1.4

° 4 °
985 995 1,005 1,015 1,025 1,035 950 1,000

_ DISTANCE, m

@O WV Sb

1050 1,100 1,150 300 «1,000 1,200 1,400 1,600

Figure 27.—Comparison of computed and analytic stresses and displacements illustrating rigld boundary effects.

of a 100-m-wide opening. Appendix B shows similar fig-
ures for the 25- and 400-m openings. Again, numerical
computations agree with analytic solutions, and the com-
putations show similar rigid boundary effects.

Finally, figure 27 compares numerical calculations with
a middle-size coarse mesh to the analytic solutions for the
25-, 100-, and 400-m openings. With the displacement cal-
culations, the magnitude of the peak displacement and the
shape of the displacement profile scales with increased
opening width. However, with the stress calculations, the
magnitude of the peak stress increases slightly and only the
shape of the stress profile scales with opening width, In
all cases, agreement with the analytic solution appears
adequate.

Several interacting factors affect the peak stress mag-
nitudes shown in figure 27. Intuitively, the peak stress and
the stress profile should increase with opening width, and
indeed consideration of the analytic solution given by
equation 46 shows that it does. Block and element size

can also have a significant effect on peak stress computa-
tion, Smaller blocks and elements approximate actual
stress profiles more accurately, and, therefore, models with
smaller blocks and elements should have higher peak
stresses at opening edges. However, figure 27 exhibits just
the opposite trend, with peak stresses increasing slightly as
block and/or clement and opening size increase. The
seeming discrepancy is a rigid boundary effect. The mod-
els made with smaller block and/or element and opening
widths are also smaller in overall dimensions; hence, the
rigid boundary is physically closer to the region of interest.
Again, presence of the rigid boundary tends to decrease
the magnitude of peak stresses and displacements in
elements too close to that boundary.

There are no strict guidelines that one can apply to
eliminate rigid boundary effects from BEM computations.
The magnitude of these effects depends on the specific
problem geometry under consideration. The simple prob-
lem shown here, coupled with experience, suggests that the


Rigid eleaent

Begin equation solver
in Tterati

Unaines
clonent

Mined-out
‘element

[compute total stresses] [Coapute total stresses

Reset convergence one sot

indicator
Wo = 0

i i

‘Coapute changes in Compute changes in
isplocenent ri i

a=

Begin block Toop 2

Tan + Tan t K

Gest
eanvargence
he

< Yer

Ostput
Tota! stron
Disptacener

1

[eer]

i

‘nd noraal ste

from blocks outside
ohborhood

[Compute displacenents
One = Dots + MD

[Compute displacenente
Dnow = Doig + 8D

Ty TLE + Te]0

neighborhood
Te = EIT + Teo

i

Coapute induced
stresses

Tahth

fend equation
solver

Key

Displocenent
Change in disp
Displacement conver ger
Element stiff
Priaitive or
Stress com

criterion

body disploceaent influence coefficient
Hal f-space displaceaent influence coefficient

Compute priaitive
str

Figure 14.Flowchart of iteration scheme using Induced-stress approach to solve for unknown stresses and displacements.

consider just the induced stresses in the elastic areas;
hence, total stresses are required. For a total-stress form-
ulation, the boundary conditions for the mined-out areas
are

(T] + [P] = [0] (23)
and for the elastic areas are
[T] + [P] + [K}"[D] = [0]. (24)

As noted earlier, for nonlinear materials, K is a function
of D. Therefore, combining equations 14, 23, and 24, now
results in a nonlinear system of equations of the form:

bo)
where the right-hand side is now some nonconstant func-

tion of x. In addition, the restriction given by equation 18,
limiting D to less than or equal to the element thickness,

Alx = (25)

still applies. Fortunately, however, the Gauss-Seidel itera-
tion procedure given by equation 19 also applies to non-
linear systems after appropriate modifications.

For the mined-out areas, combining equation 14 with 23
and then forming a recursive relation after 19 results in

p&) - p® to

TCon * TCon

which is, of course, identical to equation 21, Similarly, in
the elastic areas, combining equations 14 and 24 results in

3)
-T, -K,D,” -P,
pe 2p +4 a =Pa | 7

TCan + TCan * Ky

As stated earlier, with the nonlinear material models for
the elements, K is no longer constant, but depends on D.


Substituting into equation 34, the energy balance becomes

W = Uc- Uy + Wop + Wr. (38)
Solving for the energy release, this step yields
Wr = Uy + (W - Uc - Wgp). 39)

The energy release this step is comprised of the strain
energy released from the mined-out rock, Uy, and the
kinetic energy release, Wy, where

Wx = W- Uc - Wsp- (40)
Considering the stresses in the rock mass T, body forces
in the rock mass X, stresses in the backfill (gob) and sup-
ports R, and the associated displacements u, Salamon (6)
derived the following general expression for the energy
release at this step based on equation 39:

1 1
Wp = [5 fa, Teds - 3 hme" |

1 1
+ [3 Ja,t#ess + Zf,G-aara was Jan

where Ty = rock mass stress in state I,
uy = displacement in state I,
Au = change in displacement from state I to
state II,
AR = change in gob (backfill) stress from
state I to state II,
X,_ = body force in state I,
Sy = surface area mined this step,
Vy = volume mined this step,
Son = surface area of gob (backfill) in state II,
and a@ = nonlinearity factor for the gob or back-

fill where 0 < a,

Specifically, for linear materials, « = 1, and for strain-
hardening materials, « < 1.

‘An abbreviated derivation of this relation is pre-
sented in appendix A. As with the more generic energy

conservation relation given by equation 39, the first brack-
eted term is the strain energy release from the mined-out
rock, or Uy, and the second bracketed term is called the
kinetic energy release, or Wy. Two basic parts makeup
the kinetic part of the energy release namely, a linear part
and a nonlinear part. The linear part (i.c., the first term)
represents a change in gravitational potential from the
newly mined-out area. The nonlinear term accounts for an
energy release in the total gob and/or backfill area due to
the nonlinear stress-strain behavior of this material.
Figure 18 illustrates the significance of the nonlinearity
parameter a. Introduction of this parameter accounts for
material nonlinearity when calculating strain energy at the
point (R, u). For linear gob and/or backfill materials, a
equals 1 and the nonlinear part of the kinetic energy re-
lease is zero. For materials with & less than 1, a certain
positive kinetic energy release occurs in the gob and/or
backfill. For materials with a between 1 and 2, a negative
kinetic energy release occurs in the gob and/or backfill;
however, such materials do not seem to occur in nature.
Again, Salamon (6) derived equation 41 on the basis of
a small mining step that excavates a certain volume of rock
(or coal) Vy, over a mining area Sy. If the final excava-
tion geometry such as that shown in figure 17 state II, is
achieved in a single large mining step, then the concepts
advanced by Cook (4) apply and an energy release equal
to one-half the change in gravitational potential occurs as
per equation 33. Considering equation 41, Salamon (6)
proves analytically that for large mining steps, the first
term (the strain energy release from the mined-out rock)
is negligible compared with the second term (the kinetic
energy release). Therefore, for large mining steps with
linear elastic backfill and/or gob, equation 41 reduces to

Wr * 5 Ja, Tad (42)

GOB STRESS (R)

DISPLACEMENT (u)
a=1 O<a<t

Strain energy = 2 a Ru

Figure 18.—Strain energy calculations for various linear and
nontinear material behaviors showing how parameter « accounts
for degree of nonlinearity.


STRESS (oc)

STRAIN (e)

Figure 16.—Two-stage iteration process for nonlinear material
models in MULSIM/NL. Top, Strain-softening and elastic-plastic
materials; bottom, strain-hardening and bilinear hardening ma-
terials. Alphanumeric labels indicate two-stage iteration progress
toward final solution for stress and displacement.

FORTRAN IMPLEMENTATION OF MATERIAL
NONLINEARITIES

Equation-Solver Logic

Figure 15 shows the flowchart-logic diagram for the
MULSIM/NL equation solver that incorporates mate-
rial nonlinearities. Much of the FORTRAN code appears
in the MULSIM/NL subroutine called SOLVER. This
solver uses three sets of nested DO loops for this process.

The first set (outermost) is the main iteration loop for the
Gauss-Seidel scheme. The iterative calculations in the
loop repeat themselves until satisfactory convergence
occurs, or the user-specified maximum number of itera-
tions is exceeded. The second set of DO loops increments
through each block in the mesh, while the optional third
set of DO loops (innermost) increments through all 25 ele-
ments in a fine-mesh block. The second and third sets of
DO loops evaluate equation 20 which provides the induced
and total stress components for the block or element.

After evaluation of the induced and total stress com-
ponents for an element or block, the program branches
three ways according to the applied boundary conditions.
In the special case of rigid blocks or clements (i.e., infinite
stiffness), the displacement discontinuity components are
identically zero as are the induced stresses. The total
stresses remain identically equal to the primitive stresses.

For the mined-out areas, the program first calculates
displacement discontinuity components for the clement ac-
cording to equation 26 and then applies a series of con-
vergence tests. By definition, total stresses in the open
elements are zero, and the induced stresses are equal and
opposite to the primitive stresses. Acceptable convergence
occurs when these total stresses are approximately zero.
An exception can occur in the normal direction of the ele-
ment if the calculated normal displacement or closure ex-
ceeds the thickness of the clement. In this case, the nor-
mal displacement is set equal to the thickness, and the
total stresses are allowed to have positive or compressive
values. With complete normal closure, the stresses begin
to increase, and ultimately, the total normal stress com-
ponent may approach the primitive normal stress compo-
nent. If an element fails either a displacement or a stress
convergence criterion, then the equation solver automat-
ically proceeds with another iteration up to a uscr-set
limit.

For the elastic areas, the program first calculates
displacement-discontinuity components according to equa-
tion 27. It then retrieves the most recent total stress and
displacement values for the element or block and enters
the appropriate material property subroutine to obtain a
new, updated modulus (or stiffness) estimate. Obtaining
updated values for total stress is the first stage of the itera-
tion process. Once inside the proper material subroutine
(from a current choice of six material models), the second
stage of the iteration process begins. As discussed earlier,
and shown in figure 16, the strain-softening and elastic-
plastic models start from the current displacement to find
the intercept on the stress-strain curve of the material,
whereas the strain-hardening and bilinear hardening mod-
els start from the current total stress to find this intercept.



model edge must be about 10 block widths away from a
region of interest. Alternatively, the coarse mesh should
extend for 10 blocks beyond the fine mesh to minimize
rigid boundary effects everywhere within the fine mesh.

ELEMENT SIZE EFFECTS

As alluded to in the prior section, element size impacts
the numerical accuracy of BEM programs, particularly
near openings where stresses tend to become large at a
very fast rate. This section will examine element size ef-
fects systematically and show the effect element size can
have on the accuracy of stress and displacement calcu-
lations with MULSIM/NL. Figure 28 shows the model

/— 20 blocks —+}
T 25-m opening with
£ 12.5-m blocks
r=)
2 100-m opening with
3 50-m blocks
nN

— 20 blocks >|
T 25-m opening with
» 6.25-m blocks
2
3 FertttH |!OO-m opening with
ray 25-m blocks
8 400-m opening with
{ 100-m blocks

1-—20 blocks —>1 :
Tt 100-m opening with
o 12.5-m blocks
SHH 400-m opening with
2 50-m blocks
a

KEY

Fine-mesh blocks
Open blocks

Figure 28.—MULSIM/NL models for element size
ies.

geometries used in these studies. As in the rigid boundary
effects studies, this section considers an infinitely long
planar opening of width 2c. The coarse-mesh size is 20 by
20 blocks with a central fine-mesh 4 blocks by 12 (ie.,
20 elements by 60). The models examine three different
opening widths (25, 100, and 400 m) with block (clement)
sizes of 6.25 (1.25), 12.5 (2.5), 25 (5), 50 (10), and 100
(20) m.

Figure 29 shows displacement profiles across a 100-m
opening modeled with 12.5-, 25-, and 50-m blocks (2.5-, 5-,
and 10-m clements). Normally, accuracy improves with
decreasing block size; however, once again, rigid boundary
effects enter the problem and produce the opposite effect.
Close examination of the 100-m-opening model geometries
in figure 28 shows why rigid boundary effects will impact
the 12.5-m-block model more than the 50-m-block model.
The latter case approximates a two-dimensional problem
more closely with the analytic displacement solution given
by equation 47, whereas the former requires a true three-
dimensional solution,

The stress profiles away from the opening edge do show
the expected trend, that is smaller block sizes capture high
peak edge stresses more accurately. Indeed, figure 30
shows that the peak edge stress increases dramatically as
block size decreases. (Appendix C shows additional stress
and displacement comparisons for 25- and 400-m- opening
widths. Similar observations are noted for those numerical
experiments.)

Figure 31 compares stress and displacement calculations
made with different element sizes for different opening

0.4

ow
T

KEY
12-m block
25-m block
50-m block
Analytical
solution

DISPLACEMENT (D), m
iv

T

ok===— L Se os
950 1,000 1,050 1,100
DISTANCE, m

1,150

Figure 29.—Computed and analytic displacements across
100-m opening with different block and element sizes.


Linear Elastic for Coal (Model 1)

This basic model requires little discussion since it is the
basis for MULSIM (2) and MULSIM/BM (3). The re-
quired parameters for the model are E (Young’s modulus)
and G (shear modulus) which may or may not be related
by

__E
S* samy ®

In the BEM calculations, E relates the normal stress to
the normal displacement (closure) across an element,
whereas G relates the associated shear stresses to the
corresponding shear displacements (rides). As discussed
earlier, the three boundary conditions applied to each
element behave independently of one another.

With the linear elastic model, MULSIM/NL works
internally with a stiffness K given by

_E
x-= 0)

where E is the modulus and t is the seam thickness.

As will be discussed later, the linear elastic model
works equally well in either induced-stress or total-stress
approaches to the BEM.

Linear Elastic for Gob (Model 6)

This model is the counterpart to model 1 and is in-
tended for the gob material: left in wake of the mining.
Required parameters are again E and G plus a "gob height
factor," n, The factor n is the ratio between the height of
broken, rotated gob fragments and the unmined seam
height. The factor n typically ranges from 2 to 6 and av-
erages about 4. Beyond a range of 2 to 6 seam thick-
nesses, the rock mass surrounding mined-out gob areas
remains linear elastic. With linear elastic gob elements,
MULSIM/NL uses a stiffness given by

Ez. (10)
nt

In principle, the factor n accounts for the larger effective
seam thickness present in extracted areas that are now
filled with gob material. Another way to consider n is as
a modulus reduction factor as done by Beckett and Madrid
(3). For further details on the meaning of n, the gob
height factor, see appendix E of this report.

Strain-Softening for Coal (Model 2)

This model, after Crouch and Fairhurst (28), approxi-
mates the complete stress-strain curve observed during

laboratory strength tests on coal conducted under true
displacement control. In principle, it also describes the
yielding behavior of moderately sized pillars or the perim-
eters of large pillars. Field observations by Wang (34) and
Iannacchione (35) strongly support the use of a strain-
softening model to simulate the behavior of full-size pillars
in many mines.

Required input parameters for this stress-strain model
are a peak stress and peak strain plus a residual stress and
residual strain. In addition, the model requires a Poisson’s
ratio, v. For the normal components of stress and dis-
placement, the strain-softening model uses the peak and
residual stress-strain points. However, for the two shear
components, MULSIM/NL scales the specified peak and
residual stresses by a factor of 1/[2(1+v)]. In the linear
elastic portion of the strain-softening model, MULSIM/
NL relates normal stress and displacement with an elastic
modulus computed as E = o,/¢,, where o, is peak stress
and ¢, is peak strain. The shear stresses and displace-
ments satisfy a similar relation in the initial linear por-
tion where shear modulus is computed as G = (2,/E,)/
[2@+v)]. Recall that E and G are related via equation 8.

One restriction on the strain-softening model is that
residual strain must exceed peak strain and peak stress
must exceed residual stress. For strains greater than
residual, stress remains constant at the residual level. As
a final point, implementation of the model requires a total-
stress approach to the BEM since the model is defined in
terms of absolute or total, stresses and strains.

Elastic-Plastic for Coal (Model 3)

This stress-strain model, closely akin to the strain-
softening model, approximates a "pseudo-ductile" behavior
believed to occur in pillar cores (36). Required input for
the model is again a peak stress and peak strain, the slope
of the postyield portion (in degrees), and a Poisson’s ratio,
v. As with the previous model, the factor 1/[2(1+v)]}
scales the amplitude of the normal stress-strain curve to
obtain the shear stress-strain relations. E and G satisfy
equation 8 in the initial linear portion. Since the material
model is specified in terms of total stress and strain, as
opposed to induced stress and strain, implementation again
requires a total-stress formulation of the BEM.

Bilinear Hardening for Gob (Model 4)

This stress-strain model, analogous to the elastic-plastic
model in certain respects, permits a certain amount of de-
formation to occur prior to introducing significant element
stiffness. Input requirements for the model are stress and
strain at the hinge point, modulus in the hardening region
past the hinge point, and Poisson’s ratio, v. Again, the


may provide some indication of the burst-bump potential
of specific regions in the BEM model.

Last, MULSIM/NL features multiple mining steps to
simulate various stages of mine development, This feature
enables the user to examine stress and displacement
changes as the mine development advances. Such changes
are generally much more readily comparable to field
measurement programs that tend to measure stress and
displacement (convergence) changes as opposed to total or
absolute stresses and displacements. The energy com-
putations also need the multiple mining step feature since
they require stress and displacement changes.

These extensive enhancements necessitated complete
restructuring and reprogramming to create MULSIM/NL.
Earlier versions of MULSIM used complicated indexing
and one-dimensional arrays throughout the program. The
new MULSIM/NL takes advantage of today’s larger mem-
ory computers and uses numerous multidimensional arrays
and vectors that better reflect the underlying mathemat-
ical foundations of MULSIM/NL. The reprogramming

resulted in a much more understandable, reliable, and
modifiable program that should serve users well in the
foreseeable future.

SCOPE

The remainder of this manual discusses theoretical and
programming details of these enhancements that form the
MULSIM/NL program. The enhancements are largely
based on Sinha’s (2) original program MULSIM; hence,
many details on vast portions of the BEM program are
best addressed in that Ph.D. dissertation.

Discussion first focuses on certain BEM developments
leading up to this Bureau version. Next, sections are
devoted to the material nonlinearities and the energy cal-
culations. Subsequent sections document checks on the
nonlinear BEM and energy computations. Finally, rec-
ommendations are suggested for future development of
MULSIM/NL.

ORIGINS AND PREVIOUS VERSIONS OF MULSIM/NL

HISTORICAL REVIEW

BEM applied to solid mechanics and geomechanics
emerged into a mature engineering discipline during the
1970’s and 1980’s. Numerous text books now exist on the
subject (10-12), as well as a steady stream of specialty
conferences (13). Very basically, BEM evolved as a nu-
merical method for solving the integral equations of
elasticity theory, at first in two dimensions, and later in
three. Two fundamentally different formulations of BEM.
exist: namely, a direct method and an indirect method.
According to Wardle (14), no one formulation is clearly
superior for all problems, rather each method may have
mathematical and numerical advantages depending on the
problem geometry and boundary conditions. Research and
development work have resulted in various schemes to
improve accuracy through better element formulations (15)
and in methods to include anisotropy in the media (16).
Other significant developments have occurred in elasto-
plasticity of the continuum as reviewed by Telles (17) and
in coupling with other methods, such as finite elements
(18-19).

Much of the initial practical development work in BEM
stemmed from mining applications and the study of the
tabular orebodies, like coal seams and gold-bearing reefs.
The key notion for analysis of planar excavations was
to treat the excavation plane as a discontinuity in dis-
placement within an infinite clastic medium. Berry (20)

and Berry and Sales (27-22) used this displacement-
discontinuity approach to derive solutions for subsidence
profiles over tabular excavations. Salamon (23-25) form-
ulated the "face element method" for arbitrary tabular ex-
cavation geometries in which the displacement disconti-
nuity at one clement was the sum of the contributions
from every other clement. Lacking a computer, he used
a graphical method to evaluate the summations and com-
pute closures across complicated excavation shapes. Even-
tually, Salamon’s face element approach served as the ba-
sis for the first computer program using an indirect (the
displacement discontinuity) formulation of the BEM (26).

The carly BEM formulations only considered the clo-
sure distribution across the displacement-discontinuity
plane that is subject to a normal stress field. Later de-
velopments by Starfield and Crouch (27) permitted an in-
clined excavation plane relative to the primitive stress
directions. Further work by Crouch and Fairhurst (28) in-
cluded certain nonlinear material models as boundary con-
ditions across the plane of the displacement discontinuity,

In a sense, the initial indirect BEM formulations only
considered a one-component displacement discontinuity
(that is, one normal to the excavation plane). Research
to develop more general displacement discontinuities that
included shear components began in two dimensions.
Crouch (29) developed a two-dimensional BEM analysis of
a tabular scam near the Earth’s surface by making use of
elementary displacement-discontinuity solutions within an



MULSIM/NL PROGRAM CHECKS AND BEHAVIOR

This section provides checks on the accuracy of the
stress and displacement calculations conducted by
MULSIM/NL and illustrates certain behavioral aspects of
the model. The first check compares the linear stress-
displacement calculations of MULSIM/NL with known
analytic solutions. Simultaneously, these calculations also
examine rigid boundary effects. The presence of com-
pletely rigid material outside the problem grid can have a
severe impact on the quality of numerical results. Another
series of MULSIM/NL models examines element size ef-
fects on the numerical reliability of the stress-displacement
calculations. Finally, the last study examines the new non-
linear material models in MULSIM/NL. While it is not
possible to check the accuracy of the nonlinear models
against an exact solution, certain qualitative behavioral
checks are feasible. The intent of all these checks is to
provide some essential quality control on this complicated
numerical model.

SOME SIMPLE ANALYTIC SOLUTIONS

Jaeger and Cook (41) provide certain analytic solutions
essential for checking BEM programs. The solutions are
for an elliptic crack of width 2c subject to a stress P
perpendicular to this crack plane. Figure 23 shows this
two-dimensional crack problem. In more practical terms,
this problem corresponds to an infinitely long slot of width
2c in an infinite medium, i.e., a longwall panel of width 2c
and infinite length. The vertical stress distribution along
the edge of this opening is given by

oy = P coth(e), (46)

where € = cosh" (x/c).

Y

——- x
[+ 2¢ |
to ft tf ff ft -

Figure 23.—Elliptic crack of width 2c subject to stress P per-
pendicular to crack plane.

This equation applies along the x axis for x/e > 1. Nor-
mal convergence across this opening is given by

2

Ve iGy) P (2x2), (Co)
which applies across the region -1 < x/e < +1. These re-
lationships provide important checks on the accuracy of
MULSIM/NL. The need for this checking arises because
of the radical changes to the equation solver in MULSIM/
NL brought about by the addition of nonlinear material
models.

RIGID BOUNDARY EFFECTS

Figure 24 shows the first simple problem geometry
studied with MULSIM/NL which is an infinitely long pla-
nar opening of width 2c. The objective of these initial
studies is to compare the accuracy of MULSIM/NL with
known linear elastic solutions and to illustrate rigid
boundary effects on the computed stresses and displace-
ments. Analytic solutions presented in the previous sec-
tion give computed stresses and displacements along the
symmetry line of the meshes shown in figure 24.

Rigid boundary effects arise in BEM because the in-
seam material outside the defined mesh is infinitely stiff or
rigid. The presence of rigid material just outside the de-
fined mesh can affect stress and displacement calculations
in elements near the mesh boundary. Sinha (2) considered
this rigid boundary when designing the original linear
MULSIM program. In MULSIM, a coarse-mesh sur-
rounds a central fine mesh which is the area where ac-
curate stress and displacement calculations are desired.
The coarse mesh improves the accuracy of numerical
calculations within the fine mesh by pushing the rigid
boundary far away from it. The following studies will
illustrate these rigid boundary effects on the BEM calcula-
tions and suggest some guidelines for minimizing their
effects in practical modeling.

The numerical model used has a central fine-mesh area
4 blocks by 12 (ic., 20 fine-mesh elements by 60). As
shown in figure 24, three different coarse-mesh sizes are
considered: a small coarse mesh 20 by 10; a middle mesh
20 by 20, and a large mesh 20 by 30. In addition, three
different block sizes are used in the model; namely, 6.25,
25, and 100 m that have corresponding element sizes of
1.25, 5.0, and 20 m. Varying the block and/or element
size also changes the opening width 2c from 25 to 100 m


As shown in appendix A, which summarizes Salamon’s
work (6), the strain energy release for the mined-out
material (equation A-6) assumes linear stress-strain be-
havior. The modified strain energy release calculations in
subroutine ENERGY account for certain nonlinear stress-
strain behaviors, Figure 20 shows the linear and nonlinear
material models available in MULSIM/NL for the un-
mined seam material (coal or rock) and the assumed
strain energy relationships. On loading to a point (T,u),
each material requires a total strain energy input equal to
the total area under the given stress-strain curve up to that
point. Upon unloading, the modified strain energy release
calculations assume that the unloading modulus equals the
initial loading modulus. Thus, for nonlinear materials, the
strain energy release from the mined-out material must
equal the "recoverable strain energy" depicted in figure 20.
The difference between the total and recoverable strain
energy is the “dissipated strain energy" which is energy
presumed to dissipate nonviolently through cracking and
pseudoplastic flow of the material. As can be seen in fig-
ure 20, the dissipated strain energy is zero for the linear
material model, and the total strain energy equals the re-
coverable strain energy; hence, the strain energy release.

In summary, the strain energy release component of
equation 45 is modified as shown in figure 20 to account
for nonlinear (dissipative) behavior in the material mined-
out this step over the area Sy. In MULSIM/NL energy
calculations, the recoverable strain energy calculated as
shown in figure 20 replaces the strain energy release term
of equation 45.

Finally, subroutine energy writes the final output files
for MULSIM/NL. For each block and/or element, these
files contain block and/or element centroidal coordinates,
the material property code, three stress components, three
displacement components, and the following six energy
quantities for the block and/or element: total strain en-
ergy, recoverable strain energy, dissipated strain energy,

KEY
[SDissipated strain energy
EEE) Recoverable strain energy

(Tu),

(7, u)

3

STRESS (o)-—

DISPLACEMENT (u)

Linear elastic Strain softening Elastic plastic

Figure 20.—Basic strain energy relationships for linear and
nonlinear material models in MULSIM/NL that apply to unmined
‘seam material (either coal or rock). (T, u) is stress-displacement
‘state of element.

strain energy release, kinetic energy release, and total
energy release. The three energy release quantities for
the block and/or element may be nonzero over the area
mined this step (S,,) and the total backfill and/or gob area
(Sou); whereas, the strain energy quantities may be non-
zero only in areas outside of Sy, and Soy.

Global Energy Calculations

As subroutine ENERGY proceeds through the local el-
emental energy calculations, it also performs summations
of the elemental energy releases that, in effect, approxi-
mate the integrations prescribed in equation 45. The
strain energy release is the summation over the mined-out
area this step (S,,) of the elemental recoverable strain
energy. The prior section discusses this minor modifica-
tion from strict adherence to equation 45. The linear
component of the kinetic energy release is the summation
over S,, of the elemental change in gravitational potential
energy. The nonlinear component of the kinetic energy
release is the summation over the total backfill and/or gob
area (Sou) of energy dissipation due to material nonlinear-
ity. The total kinetic energy release is simply the sum of
the linear and nonlinear components, and finally, the total
energy release is the sum of the strain energy and kinetic
energy releases. Upon exiting ENERGY, the aforemen-
tioned energy release summations are written to a print
file. For each mining step, this file contains the different
energy release components for the coarse-mesh blocks, the
fine-mesh elements, and their total. Finally, these energy
release quantities are normalized by the area mined this
step (area Sy) to obtain the corresponding energy release
rates.

Three Special Cases

Earlier discussions alluded to certain special cases that
can arise which complicate the programming of equa-
tion 45. Figure 21 shows the three possible cases as min-
ing progresses from state I to state II. In case A, the new-
ly mined-out area S,, exactly balances the newly backfilled
area AS, and the open area in state I equals the open
area in state II. In case B, the backfilling lags the mining
and the open area increases, whereas in case C, the back-
filling exceeds the mining and the open area decreases.
Also, as shown in figure 21, individual elements can un-
dergo several types of changes as mining progresses in a
step from state I to state II.

In the simplest case A, four types of change exist.
Those elements in region 1 are unmined in state I and re-
main unmined in state II. Total strain energy, recoverable
strain energy, and dissipated strain energy are computed


UNIT OF MEASURE ABBREVIATIONS USED IN THIS REPORT
cm centimeter MJ/m? —_megajoule per square meter
ft foot m/m meter per meter
Tb pound MPa megapascal
m meter psi pound per square inch

Disclaimer of Liability

The U.S. Bureau of Mines expressly declares that there are no warranties expressed or implied that apply to the software
contained (or described) herein. By acceptance and use of said software, which is conveyed to the user without
consideration by the Bureau of Mines, the user hereof expressly waives any and all claims for damage and/or suits in-

cluding special, consequential, or other similar damages arising out of or in any way connected with the use of the soft-
ware contained (or described) herein.



38.

39.
. Total energy release from mining step 2 to 3 and its components for five different test models with

41,

A-1.
A-2.

A3.
B-1,
B-2.
C1.
C2.
D-1,
D-2.
E-1.

E-2.
E-3,

. Computed and analytic stresses at edge of 100-m opening with different block and clement sizes
31.
32.
. Comparison of computed stresses and displacements using linear and nonlinear material properties
. MULSIM/NL grid for widening-crack problem
35.
. Four MULSIM/NL grids for step-size effect tests
31.

ILLUSTRATIONS—Continued

Comparison of computed and analytic stresses and displacements illustrating element size effects .... .
MULSIM/NL models for nonlinear effect studies .

Comparison of analytic and numerical energy release rates for widening-crack problem

Normalized cnergy release rate components versus number of excavation steps for analytic case and
mumerical case 0... eee net n eee n eben eee tbenetetaees

Five different three-step MULSIM/NL models used to study constant, decreasing, and increasing area
effects on total energy release and its components ........ 20.0. 0ee cece cece cues

Linear and nonlinear stress-strain curves for backfill material used in model studies

linear backfill material... 6.0.6... eee cece cece cece cece n ence een eee
Total energy release from mining step 2 to 3 and its components for five different test models wit

nonlinear backfill material
Basic problem geometry for energy release derivations
Strain energy calculations for various linear and nonlinear material behaviors showing how parameter a

accounts for degree of nonlinearity 2.1.6.6... cece eee eee eet e teen teen nen eee
Typical strain-hardening stress-displacement behavior for backfill or gob with basic notation
Computed and analytic stresses at edge of opening with different coarse-mesh sizes ......
Computed and analytic displacements across opening with different coarse-mesh sizes ....
Computed and analytic stresses at edge of opening with different block and element sizes ...........
Computed and analytic displacements across opening with different block and clement sizes .........
Computed stresses at edge of opening using various lincar and nonlinear material models ...........
Computed displacements across opening using various linear and nonlinear material models
Three mechanical zones found above mined-out coal seam
Swelling and consolidation mechanism of gob formation
Assumed linear relationship between tangent modulus and stress 22.2... eee eee eee eee eee

TABLES

Material model parameters and material property array EPROP structure .
Strain-softening properties used in nonlinear models



stress-strain relation for the normal direction of the ele-
ment is scaled by 1/2[(1+v)} to obtain the stress-strain
relation for the tangential directions of the element.

This gob model also requires the gob height factor n.
This factor increases the effective seam thickness for those
elements employing the various material models for gob.
An alternative way to view n is as a modulus reduction
factor to account for effective seam thickness in the gob
areas that exceeds seam thickness in the coal areas.

Strain-Hardening for Gob (Model 5)

This model allows the gob material to increase in stiff-
ness as it consolidates under increasing load, Appendix E
provides the complete derivation of this model. The mod-
el assumes that the gob formation occurs via a two-step
mechanism: first, the overburden rock behind the active
face fractures, collapses, and swells to fill the newly
created void, and second, the gob consolidates and follows
a nonlinear, strain-hardening constitutive relationship ver-
tical stress on the gob increases. As shown in figure 13,
derivation of the model begins by assuming that the tan-
gent modulus of the gob increases linearly from an initial
modulus value E; to a final modulus value Ef over the
stress range 0 to oy, where oy is the final stress. Beyond
a stress level of oy, the tangent modulus of the gob may
continue to increase, but for practical purposes, stress
levels beyond oy are unlikely except in rare multilevel
mining cases. The relation describing the linear increase
in tangent modulus is

Ep-E,
E= ea k +E.
Oy
The differential relation between stress (2) and strain (¢)
is

(11)

do = Ede. (12)

Combining these relations and integrating produces the
following stress-strain relation for gob:

_E;[ no, Ep-E; |[p]_
-3feslb Es Rh]

where E, = initial modulus,

Ep

final modulus,

oy = virgin vertical stress,

‘TANGENT MODULUS (E)
m
£

a

STRESS (co)

Figure 13.—Assumed linear relationship between tangent
modulus (E) and stress (2) over range 0 to cy E, Is initial
modulus at zero vertical stress, and Epis final modulus at vertical
stress oy

D = seam closure,
t = seam thickness,

and n= gob height factor.

The parameter (Ey - E,)/(n o,) controls the degree of
nonlinearity in the model. Input parameters are E,, Ep,
o,,n, and v. As with the prior models, scaling the above
stress-strain relation for the normal direction of the
element by a factor 1/[2(1+v)] provides an appropriate
stress-strain relation for the tangential directions of the
clement. Again, a total-stress formulation of the BEM is
required.

TOTAL-STRESS FORMULATION FOR MULSIM/NL

In making the change from strictly linear elastic mate-
tial properties, it becomes necessary to shift the basic
boundary-element formulation from an induced-stress ap-
proach to a total-stress approach. The original MULSIM
and MULSIM/BM programs conduct computations strictly
with the induced stresses. After the equations are solved,
the induced stresses are added to the primitive stresses to
obtain the total stresses for output. With linear elastic
materials, the elastic modulus serves as the proportionality
constant between the induced stress and displacement - a
proportionality unaffected by the magnitude of total stress.
With nonlinear materials, however, the absolute magnitude


INTRODUCTION

PURPOSE AND OBJECTIVE

The program described in this U.S. Bureau of Mines
report is part of the MULSIM/NL package, which fea-
tures the actual BEM program described herein, as well
asa preprocessor program called MULPRE/NL and
a plotting postprocessor program called MULPLT/NL.
The preprocessor helps the user generate the requisite
input file for the main program MULSIM/NL, whereas
the postprocessor assists the user in studying the calcu-
lated stresses and displacements in a very rapid graphical
manner.

Documentation and instructions for the MULSIM/NL
package are divided into two related reports with each
providing essential elements toward an understanding of
the whole. The objective of this report, "MULSIM/NL -
Theoretical and Programmer’s Manual," is to provide
certain mathematical and programming details to those
engineers and programmers who need to fully understand
the FORTRAN program or desire to alter and enhanc
A related report, "MULSIM/NL - Application and Practi-
tioner’s Manual," provides users with detailed operating in-
structions for MULSIM/NL along with several practical
examples to better illustrate the capabilities of the
program,

Development of MULSIM/NL stems from a need to
reduce mine accidents due to bumps in U.S. coal mines,
A bump (very similar to a rock burst in many respects) is
the violent failure of highly stressed in-place coal. They
can range in size from small face bumps, which violently
eject several pounds of coal across an entry, to massive
longwall panel bumps, which violently eject hundreds of
tons of coal and damage or destroy millions of dollars
worth of longwall mining equipment. Naturally, coal mine
bumps pose a severe hazard to the miners. The smaller
bump events may severely injure a miner unlucky enough
to be in its path, whereas the larger events may result in
multiple fatalities (2).2 Even if a bump does not result in
costly lost-time accidents or fatalities, recovery of produc-
tion may take days to complete while damaged mining
equipment is repaired. Furthermore, sterilization of
large coal reserves may occur if the bump risk becomes
too great and sections of the mine are abandoned.
MULSIM/NL is an important tool for bump control be-
cause it provides a means to calculate stress, displace-
ment, and energy changes for various mining configura-
tions in bump-prone conditions. Through judicious use of

Italic numbers in parentheses refer to items in the list of references
preceding the appendixes at the end of this report.

numerical models such as MULSIM/NL, engincers can
choose alternate mine designs that may decrease the bump
tisk,

In addition, MULSIM/NL has many practical uses be-
yond coal mine bump research. Stresses and displace-
ments calculated by MULSIM/NL can help an engineer
design pillar sizes and understand multiple-seam mining
interactions,

MULSIM/NL ENHANCEMENTS

Three significant enhancements to earlicr versions of
MULSIM resulted in MULSIM/NL, namely (1) nonlinear
material models, (2) multiple mining steps, and (3) com-
prehensive energy release and strain energy computations.
These enhancements stem from the needs of current Bu-
reau research on coal mine bumps. Field measurements
made on this and numerous other coal pillar research
projects clearly indicate that coal follows highly nonlinear
stress-strain paths of various forms. In addition, the gob
or backfill material left in the wake of coal extraction
follows yet other nonlinear stress-strain paths as it con-
solidates and closure occurs. Prior versions of MULSIM
(2-3) only permitted linear stress-strain relations for in-
seam materials such as coal or gob. MULSIM/NL now
has six material models from which to choose including
(1) linear elastic for coal, (2) strain softening, (3) elastic
plastic, (4) bilinear hardening, (5) strain hardening, and
(6) linear elastic for gob.

MULSIM/NL also calculates detailed energy changes
due to mining. The basis for these calculations is the
energy release rate (ERR) concept advanced by Cook (4),
later by Walsh (5), and most recently by Salamon (6) and
Brady and Brown (7). Early South African research on
rock bursts in deep gold mines found that the ERR, i.e.,
the energy dissipation per unit area mined, correlated well
to the incidence and severity of devastating rock bursts
(8-9). Based on the South African experience with ERR
and rock bursts, the hypothesis is put forth that ERR also
correlates with the incidence and severity of coal bumps.
MULSIM/NL implements the theoretical developments of
Salamon (6) into its ERR computations. The implementa-
tion is also complete in that it correctly considers the
nonlinear component of the energy release arising from
nonlinear backfill or gob material. In addition, MULSIM/
NL also computes various useful strain energy quantities
for the unmined, nonlinear material in the BEM mesh.
These energy quantities include the total strain energy
delivered to an clement, the recoverable elastic strain
energy available from an element, and the dissipated strain
energy from an element. Examination of these quantities


NONLINEAR ENERGY CALCULATIONS IN MULSIM/NL

This section begins by briefly tracing the origins of
energy release rate (ERR) concepts applied to rock burst
control in the mines. Salamon (6) re-examined the theo-
retical basis of the ERR concept. The analytic results of
his work provide the basis for the nonlinear ERR cal-
culations in MULSIM/NL. The discussions summarize
Salamon’s (6) derivations and provide key details of the
FORTRAN implementation in MULSIM/NL.

EVOLUTION OF ENERGY RELEASE
RATE CONCEPTS

Cook (37) first noted the enormous energy changes that
result from underground mining and postulated "excess po-
tential energy causes the damage noticed as rock bursts."
Subsequent papers by Cook (4, 38) provided complete der-
ivations and discussions of what became the ERR concept.

When an excavation is created underground, an enor-
mous gravitational potential energy, Wgp, source becomes
available and is given approximately by

Wop = % Vex (28)

where g, is the in situ vertical stress and V,, is the total
excavation volume. Immediately after excavation creation,
the excavation walls converge to some extent and a portion
of the available gravitational potential energy given by
equation 28 is transformed into work. This change in
gravitational potential energy, Wg, is
We = 2 Ves (29)
where V, is the excavation closure volume. Part of the
gravitational potential energy change given by equation 29
transforms into an increase in the strain energy storage,
Wg, of the rock mass surrounding the excavation. Cook
(4) showed that this strain energy storage necessarily has
the following bounds:
O< We<% Wg. (30)
Thus, creation and subsequent partial closure of an ex-
cavation results in a certain positive strain energy storage
increase in the rock mass that never exceeds one-half the
gravitational potential energy change.
Energy conservation requires that
Wg = Ws + Wr, (1)
where We is the energy release due to mining. Note,
that the energy release, Wa, is introduced to preserve

conservation. Considering equations 30 and 31, the follow-
ing additional inequality must also apply:
Ws < Wr < Wo. (32)
In the special case of a linear elastic rock mass, then equa-
tions 30 and 32 become:
Wp = Ws = % Wo. 33)
Stated semantically, upon creation of an underground ex-
cavation, the increase in strain energy storage in the
surrounding rock mass must equal the energy release iden-
tically, which in turn, equals exactly one-half the change in
gravitational potential energy. Cook (4) concluded

“the energy release must either be in the form of
nonviolent dissipation in the course of crushing of
rock or supports, or in the form of violent events. In
the former case, the released energy is transformed
mainly into heat through friction, and in the latter
case, into kinetic energy. A portion of the violent
events is manifested as rock bursts.”

Thus, the energy release may provide a quantitative
measure of the damage potential from rock bursts.

Associated with the energy release is the ERR defined
as the energy release per unit area mined. (Note that the
use of rate in ERR is unfortunate since there is no
relation to time involved.) The energy release is fixed by
the depth and size of excavation, and engineers have
relatively little control over its magnitude. However, the
ERR depends on the layout and extraction sequence for
an excavation, and engineers can exert considerable in-
fluence over its magnitude. Researchers in South Africa
have demonstrated a clear correlation between the number
and incidence of damaging rock bursts and the ERR (8-9,
39). Thus, the ERR provides an indicator to assess rock-
burst potential. By analogy, a similar correlation may exist
between ERR and the risk of coal mine bumps. When
evaluating mine layouts and extraction sequences with
ERR, the objective is to minimize the peak ERR. This
objective is achieved by keeping the ERR relatively con-
stant over an extraction sequence. In general, the most
uniform ERR results when extraction begins at the most
highly stressed areas and proceeds toward the low stress
areas,

Subsequent researchers continued to develop un-
derstanding of the energy changes due to mining. Walsh
(5) presented general theoretical relations describing
mining-induced energy changes. These included terms to


BUREAU OF MINES
INFORMATION CIRCULAR/1992

MULSIM/NL Theoretical and
Programmer’s Manual

By R. Karl Zipf, Jr.

UNITED STATES DEPARTMENT OF THE INTERIOR



39. Brady, B. HL G., and E. T. Brown. Energy Changes and Stability
in Underground Mining: Design Applications of Boundary Element
Methods. Trans. Inst. Min, and Metall (London), Section A, v. 90, Apr.
1981, pp. A61-A68.

40. Salamon, M. D. G. Rockburst Hazard and the Fight for Its
Alleviation in South African Gold Mines. Paper in Rockbursts: Predic-
tion and Control. Inst. Min, and Metall. (London), 1983, pp. 11-36.

41. Jaeger, J.C, and N. G. W. Cook. Fundamentals of Rock Mech-
anics. Chapman and Hall (London), 1979, 593 pp.

42. Dahl, H. D., and H. S. Von Schonfeldt. Rock Mechanics Elements
of Coal Mine Design. Paper in Proceedings 17th U.S. Symposium on
Rock Mechanics, Univ. UT, Salt Lake City, UT, 1976, pp. 4A1-1- 41-9.

43. Owili-Eger, A. S.C. Dynamic Fractured Flow Simulation Model.
Min. Eng., AIME, 1989, pp. 110-114.

4, Singh, M. M., and F. S. Kendorski. Strata Disturbance Prediction
for Mining Beneath Surface Water and Waste Impoundments. Paper in
Proceedings 1st Conference on Ground Control in Mining, WV Univ,
Morgantown, WV, 1981, pp. 76-89.



MULSIM/BM

During the period 1984-85, the Bureau upgraded and
restructured Sinha’s original program to create MULSIM/
BM (3). All the original capabilities of the program
remain in the enhancement since MULSIM/BM still cal-
culates stresses and displacements in single or multiple
parallel seams (up to four) that are inclined at any angle
to the Earth’s surface and oriented at any angle to the
principal stress field. The program neglects the effect of
the Earth’s surface and assumes that the seam planes are
at great depth (i.c., depth greater than opening width).
MULSIM/BM includes these major new features:

(1) Additional Linear In-Seam Material Properties -
The original program permits only one linear property set
for the in-seam material, whereas the new version permits
up to 26 different linear elastic materials.

(2) Inserted Material or Gob Model - Since the orig-
inal program permits only one material, no method exists
to have a different lower modulus material remain in the
wake of mining. In addition to allowing numerous linear
elastic materials, the new version features an inserted
material or gob model to represent the layer of broken
tock, such as gob left behind in coal mining, or other
artificial supports such as backfill, packwalls, or cribs.
Certain in-seam materials can be identified as inserted
materials, in which case, the model uses a linear stress-
strain relationship with zero stress at zero seam closure
@).

(3) Increased Coarse- and Fine-Mesh Size - The orig-
inal MULSIM allows up to four seams each with a 12 by
12 array of coarse-mesh blocks, each of which can be
subdivided into a 5 by 5 array of fine-mesh elements.
MULSIM/BM expands the permissible coarse-mesh array
size up to 40 by 40. The allowable number of paral-
lel seams remains at four. The number of coarse-mesh
blocks divisible into fine-mesh elements also increased
slightly from a 12 by 12 coarse-mesh block array to a 20 by
20 array, thereby permitting a 100 by 100 array of fine-
mesh elements.

(4) Coarse-Mesh Block Stress and Displacement Cal-
culations - In the original version, MULSIM only prints
stress and displacement calculations within the array of
fine-mesh elements. Usually, only the fine-mesh cal-
culations are of interest; however, for certain applications,
having the coarse-mesh stress and displacement calcula-
tions available is desirable. MULSIM/BM now displays
those calculations in the print file. Data files retain block
and element locations and stress and displacement calcu-
lations separately for the coarse-mesh blocks and fine-
mesh clements.

(5) Addition of Extraction Ratios in Coarse Mesh -
The new version, MULSIM/BM, allows the user to specify
the extraction ratio of each coarse-mesh block. This fea-
ture enables the user to include effects of prior mining in
the coarse-mesh area without having to specify all the de-
tails of the mine layout. The program then decreases the
stiffness of a partially extracted block as a function of the
extraction ratio.

MULSIM/BM has seen considerable use in various Bu-
reau projects. Beckett and Madrid (3) describe its utility
in general coal mine modeling. Kripakov (32) also pre-
sents examples of its utility in coal mine modeling and
discusses a coupling procedure between MULSIM/BM
and a finite-clement analysis. Recently, Kripakov and
Rockwell (32) used MULSIM/BM to study various coal
mine gate road pillar configurations under bump-prone
conditions. Among other conclusions, they recommended
further developments of Bureau BEM programs to include
nonlinear in-seam material properties and various energy
calculations. These developments have been included in
the Bureau’s latest BEM development, MULSIM/NL.

ITASCA ADDITIONS

During 1988, the Itasca Consulting Group, Minneapolis,
MN, provided certain upgrades to MULSIM/BM while
under contract to the Bureau. These modifications, name-
ly multiple mining steps and basic linear energy release
computations satisfied portions of the Bureau’s goals for
a BEM-based numerical model with (1) nonlinear mate-
rial properties, (2) multiple mining steps, and (3) com-
prehensive nonlinear energy release and strain energy
computations.

Implementing multiple mining steps is an essential pre-
requisite for implementing energy computations. By itself
though, the multiple mining step feature provides a useful
tool for the analyst and/or engineer. With it, the user can
simulate a continuous mining process by several discrete
steps. For example, a mining-cut sequence might split a
barrier pillar with 20 10-ft-deep continuous miner cuts.
Each cut might be one step in a detailed BEM analysis,
To simulate progression of a longwall face, each mining
step in the BEM analysis might remove an entire row of
elements along the longwall face. Then in a sequence of
about 10 mining steps, the BEM analysis might simulate
the mine-by cycle of a tailgate pillar. With the multiple
mining step option, the analyst and/or engineer can ex-
amine progressive changes in stress and displacement.
Load transfers can be observed numerically as they occur
during mine development.


isotropic infinite half-space. Further two-dimensional gen-
eralizations by Crouch (30) enabled analysis of nontabu-
lar and multiple seams arbitrarily close to the Earth’s
surface. Finally, Sinha (2) developed general three-
dimensional solutions for stresses and displacements both
normal and parallel to the displacement-discontinuity
plane. He presented programs for (1) tabular seams near
the Earth’s surface, (2) multiple seams at great depth,
and (3) a nontabular seam at depth. It is the second pro-
gram that formed the essence of MULSIM/BM (3) and
MULSIM/NL.

SINHA’S MULSIM

The original MULSIM was intended for stress and dis-
placement analysis of multiple seams at great depth below
the Earth’s surface. It could accommodate up to four
parallel seams where cach seam is modeled by a 12 by 12
array of coarse mesh blocks. A 5 by 5 grid of fine-mesh
elements can subdivide each block if requested. Figure 1,
taken from Sinha, shows the general modeling-grid
scheme. The seams are parallel to each other and share
a common Z axis. The local X and Y axes in each seam
also remain parallel (i.¢., there is no rotation about the
Z axis), and there is no origin offset.

The surrounding rock mass is a linear elastic medium.
The in-seam material is also linear elastic; however,
certain anelastic behavior is also considered. For example,
one of the seams can be considered a bedding or fault
plane parallel to the other seams. MULSIM models any
gouge material in a seam plane with a Mohr-Coulomb
failure criterion that allows the element to yield and
deform in response to normal and shear stresses.

The general structure of the original program MULSIM
is shown in figure 2. In principle, program flow is straight-
forward and simple, although the details are complex. In-
put consists of a file defining the problem and its ge-
ometry. After calculating primitive stresses and influence
coefficients, the program enters an iterative equation-
solving loop to calculate the unknown displacement com-
ponents. When it has found all the elemental displace-
ments, the program then calculates induced and total
stresses and outputs the results. All output is simply
printed as tables in a rather lengthy line printer file.

In his dissertation, Sinha (2) derives all the three-
dimensional displacement-discontinuity influence coeffi-
cients from potential theory, and he presents the basic
programs for analyzing common mining situations in seam
deposits. He also conducts rigorous tests of the program
to ensure its accuracy.

Element “mee

Ye
id eam |
iy
Bila)
Ee,
Hy

Enter

?

Read ond print Read seam geometry
input variables and material properties

i i

Obtain ride and closure
correction factors

i i

Obtain in situ stresses Solve for unknown rides
‘and gradients ‘and closures

i i

Obtain influence Obtain element stresses,
coefficients - stress ‘and displacements,
‘nd displacement and output

i i

Obtain requested off-seam
Store influence stresses and displacements

i i

Initialize random access
Tass storage file for
rides and closures

Define constants

Exit

Figure 2Programming structure of MULSIM (2).


Kinetic

eneray release rate, linear component

Kinetic energy release rate, nonlinear component

1.6 pt : : : r r r r : : 1 T
iat 4 4 8 4
Lab 4, + 4
of J ft |
ef 14. |
6 4 4
4 J
Ww
Book 4 4
= wz fara
a, 0 J J
we
25 -2l L L 1 1 1 L 4 1 1 .
a2
go Strain energy release rate Total energy release rate
= Ber 1 + + + 1 . 1 1 T
&~ 14t ¢ J Lo |
Pa
Sek 4 of J
1.0 +f 4
+]
ot | |
4b 4 +f 4
ab | 1
of am | | L |
-plu 1 fi ‘ fi fl 1 fl f 1 f A
Constant open area Decrease Increase Constant open area Decrease Increase
Small Middle Lorge open area open area Small Middle Large openarea openarea

Figure 41.—Total energy release from mining step 2 to 3 and its components for five different test models with nonlinear backfill

material. Note: graph A + graph B + graph C = graph D.

ERR and the strain ERR should be retained in total ERR
calculations.

(4) Alternatively, there is a possibility that the above
trends which are credited to the open-area-to-backfill-area
ratio may be due to a larger overall backfill area. This
possibility requires further study.

(5) Additional MULSIM/NL tests not supported by
tables and figures suggested that increasing the modulus of
the backfill (final modulus for the nonlinear tests with the
strain-hardening backfill model), causes the linear kinetic
ERR and the total ERR to decrease. The relative magni-
tude of the ERR decrease depends on the open-area-to-
backfill-area ratio. The small ratio models had larger
relative decreases in ERR. In these example problems,

varying the modulus in the linear backfill from 50 to
150 MPa caused a 25% decrease in ERR for the small
area ratio model and a 10% decrease from the large area
ratio model. Varying the final modulus in the nonlinear
strain-hardening backfill from 140 to 560 MPa caused 15%
and 10% ERR decreases for the small and large area ra-
tios. Present measurements of the stress-strain behavior
of gob limit the reliability of gob modulus estimates to an
order of magnitude only. Such poor reliability in gob
modulus estimates adversely impacts ERR estimates.
Thus, research studies on the stress-strain behavior of gob
material may substantially improve the quality of ERR
calculations.


case, elements in region 6 exist which have gone from un-
mined coal to backfill (gob) in one step without remaining
open during an intermediate step. Such elements belong
to both the set S,, and Sgy. For these elements, the strain
energy release is computed as the recoverable strain ener-
gy; the linear component of the kinetic energy release is
computed as % T, Au, and the nonlinear component of the
kinetic energy release is computed % AR Au. The total
energy release for these special elements is the sum of
these three parts, The strain energy release, total kinet-
ic energy release, and total energy release are recorded

DO ALL ELEMENTS

Type 1 element
compute strain

Type yes | energies

for the element, and all energy release summations are
updated.

Subroutine Energy Logic Diagram

The three special cases just discussed give rise to six
element classes or types which require separate considera-
tion during the energy calculations. The logic diagram
shown in figure 22 reflects these six logic paths and pro-
vides the essence of subroutine ENERGY.

this. step

unnined, ~ Total

~ Recoverable
= Dissipated

Compute non=tineor
kinetic energy
releose

update sunnations

Type
this step
backfill

Type 2 element
compute strain
energies

= Total

~ Recoverable

~ Dissipates
Conpute energy
releases

= Strain

= Lineor kinetic
= Total

Update suamations

had natericl
Jost step

Type 6 elenent
compute strain
energies

~ Total

= Recoveroble

Yes | ~ Dissipated
Conpute energy
releases

~ Strain

~ Linear kinetic
~ Toto!

Update suasat ions

Etenent

Type
last step
pockfil

Type 3 elenent

Type 5 element

Write elenent
energies
END D0

Figure 22.—Simplified flowchart for subroutine ENERGY showing six possible logic paths.


APPENDIX A.—REDERIVATION OF SALAMON’S ENERGY RELEASE EQUATION

This appendix gives an abbreviated rederivation of Sala-
mon’s (6)! energy release equation presented as equa-
tion 41 in the text. The rederivation here uses a notation
similar to Salamon’s that is improved and, hopefully, easier
to follow. The analysis begins with equation 39, which is
the energy release relationship derived from basic energy
conservation considerations.

Wr = Uy + (W-Uc-Wep), (A-1)

where Wp = total energy release for this step,
Uy. = strain energy release from the mined-
out rock,
W_ = work done on the system by boundary
and body forces,
Uc = change in strain energy in the rock
mass,
and Wgp = change in strain energy in the backfill

(gob) and support systems.

The term in parentheses is called the kinetic energy re-
lease Wx.

Figure A-1 (same as figure 17) illustrates the basic
problem and notation considered in this derivation, Min-
ing proceeds in a small step from state I to state II during
which a volume of rock Vy, with surface area S,, is extract-
ed. The open area and volume change from So, and Vo,
to Son and Voy. The backfill (gob) area and volume
change from So, and Vg; to Sy and Voy. The volume of
the rock mass and its outer boundary area remain un-
changed at V and S, respectively. A summary of basic
volume and surface area definitions follows.

V_ = rock mass volume

Vo, = volume open in state I

Von = volume open in state IL
Vg = volume gob in state I

Von = volume gob in state IL
Vu = volume mined this step

Ntalic numbers in parentheses refer to items in the list of references
preceding the appendixes.

4Vo = You - Yor

AVG = Von- Var

$ = rock mass boundary

Soq = surface open in state I

Son = surface open in state II

Sg = surface backfilled in state I
Sou = surface backfilled in state I
Sy = surface mined this step
ASo = Son - Sor

48g = Son - Sor

In this notation, subscripts have the following meanings:
I - state I
TI - state II
G_ - gob (backfill) area or volume
© - open area or volume

M_ - mined-out area or volume this step

none - rock mass area or volume.
State 1 fa Key
Vo pence V_ Rock mass volume
1 Vw_i-Sm Vor Volume open in state I
Vor_}-Sor} Von Volume open in state IE
LL. Ver Volume of backfill in state I
Ver [Set Ver Volume of backfill in state IT
Vy Volume mined this step
S-_ Rock mass boundary area
Stote £8 Sor Surface area open in state I
Vv Som Surface area open in state IL
Son, Seq Surface area of backfill in state I
Sen Surface area of backfill in state IE
* Sen! Sm Surface area mined this step
Vex

Figure A-1.—Basic problem geometry for energy release derl-
vations.


1 STATE IL
2 Som
CASE A
ASo=Sm 3 Ser
Sox* Son \N
Set 4
1
{
Ss, KEY
CASE B 2 [| _ 0% | Sy Mined-out area this step
ASg< Sy Sor $ Sor Open area ~ state I
So1<Somt NS Som Open area - state IZ
Se 4 Sen Sor Backfill area - state 1
Sgn Backfill area - state IZ
1 ASg Sem~ Set
1 Possible changes in element status
Ss state I to state IL
2 OIL .
z 4 Unmined —* Unmined
CASE C 2  Unmined —* Open
BSg>Sm Sor 3 Sor 3 Open — Backfilled
Sot > Som SN 4 Backfilled —* Backfilled
Sar 4 5 Open  — Open
6 Unmined —» Backfilled
1

Figure 21.—Three special cases of mining progress from state | to state Il showing how elements
can undergo six different status changes between state | and state Il.

and recorded for type 1 elements with the methods shown
in figure 20. Elements in region 2 are mined out in going
from state I to state II. As discussed earlier, strain energy
release for these elements is computed as the recoverable
strain energy shown in figure 20. This method for nonlin-
ear materials replaces the first term in equation 45 where
the strain energy release from a linear element is com-
puted as % T, u,. In addition, elements in region 2 also
contribute a linear component of the kinetic energy release
given as % T, Au. Together, the strain energy release and
linear kinetic energy release comprise the total energy
telease for the element. These three energy releases are
recorded for the element, and the energy release summa-
tions for this model step are updated. Elements in
regions 3 and 4 begin as open or backfilled, respectively,
in state I and change to, or remain, backfilled in state II.
For these elements, a nonlinear component of the kinetic

energy release is calculated as ¥ (1-0) AR Au. (The only
difference between regions 3 and 4 is the state I stress,
which is zero in region 3 and nonzero in 4.) In these ele-
ments, the nonlinear kinetic energy release makes up the
whole of the total energy release. These two energy re-
leases are recorded for the element, and the energy re-
lease summations for the model are updated.

Case B in figure 21 illustrates the first exception. As
before, case B has elements in regions 1 through 4; how-
ever, when the backfilling lags the mining, a region (5)
exists with elements that were mined-out and open in
state I and remained so in state II. Such elements contrib-
ute zero energy to the energy release and contain no strain
energy. Therefore, such elements require no computations
that is easily accommodated in the program!

Finally, case C in figure 21 shows the other exception,
which arises when backfilling exceeds the mining, In this


MULSIM/NL THEORETICAL AND PROGRAMMER’S MANUAL

By R. Karl Zipf, Jr."

ABSTRACT

MULSIM/NL (multiple seams, nonlinear) is a three-dimensional boundary-element method (BEM)
program developed at the U.S. Bureau of Mines for stress and displacement analysis of coal mines and
thin metalliferous veins. It can analyze one to four parallel seams that can have any orientation with
respect to the Earth’s surface. Three main new features distinguish MULSIM/NL from its predeces-
sors: (1) nonlinear material models, (2) multiple mining steps, and (3) energy release and strain energy
computations. MULSIM/NL has six material models for the in-seam material: (1) linear elastic for
coal, (2) strain softening, (3) elastic plastic, (4) bilinear hardening, (5) strain hardening, and (6) linear
elastic for gob. Detailed checks show that numerical stress and displacement calculations compare well
to known analytic solutions for simple problems. The total energy release calculated by MULSIM/NL
is comprised of three basic terms: (1) strain energy release from the mined-out material, (2) linear
kinetic energy release due to change in gravitational potential over the mined-out area, and (3) nonlinear
kinetic energy release due to nonlinear materials in the total backfill and gob area. ‘The strain energy
release term is modified to account for nonlincar stress-strain behavior in the unmined materials.
Comparisons of these numerical energy release rate calculations with known analytic solutions show
excellent agreement.

‘Structural engineer, Denver Research Center, U.S. Bureau of Mines, Denver, CO.


The solution procedure must continuously update the val-
ues of K according to the latest estimates of D until stable,
final values of D are achieved.

The new nonlinear equation solver in MULSIM/NL
invokes equations 26 and 27 and continuously updates K
as a function of D. Figure 15 presents the new flowchart
and provides key details of the new solver.

Illustration of Convergence Schemes

Each iteration step for the nonlinear system proceeds
in two stages. The first stage executes the summation
given by equation 20 and provides total stress at that block
or element. Based on the new total-stress estimate and
the displacement estimate from the prior iteration, the
second stage of the iteration process computes a new es-
timate of the stiffness K at that block or element following
the specified material model and its material properties.
Using this updated stiffness value K, MULSIM/NL then
computes a new displacement estimate based on cither
equation 26 or 27. Repeating this two-stage iteration proc-
ess for total stress then stiffness and displacement results

Begin equation solver
Begin acin iteration
1009

Rigid efenent

Vnained
leaent

in rather rapid convergence to stable final values for the
unknown stresses and displacements at each block and
element.

Figure 16 shows this two-stage iteration process for two
kinds of material models. In both cases, total stress and
displacement (or strain) lie at the point 1A for the current
iteration. For the strain-softening or elastic-plastic ma-
terial models, MULSIM/NL uses the current displacement
estimate to determine the intercept labeled 1B on the
stress-strain curve, It then determines a new stiffness
estimate for use in the next iteration from that intercept.
In contrast, with the strain-hardening or bilinear harden-
ing material models, MULSIM/NL uses the current total
Stress estimate to determine the intercept labeled 1B
on the stress-strain curve. In both cases, the iterative
equation-solving process repeats itself until the calculated
stresses and displacements (strains) converge to the pre-
scribed stress-strain model within acceptable limits. The
number of iterations required can range from 10 to 200
depending on the size and complexity of the model. A
typical problem may need 20 to 50 iterations.

Mined-out Yer

elonent

Get actericl type
ond enter subroutine
with TS ond D

Reset. convergence
indicators
110 ='0 UN's 0

Compute changes in
Gisptacenent

Output
Total stresses ~ 1S
Displacenents ~ 0

[Compute displacenente
Drew = Dota + Gd

fron blocks outeid
‘neighborhood
Ty = FTE + TD

sonvergence criterion
atr

eaten P

Coupute total stresses
S=P4+T

Te = £10 + TD IND = IND +1 Te body displacement influence coefficient
TE? rteapace sapiveent In lores setiitet
T Tote sess
Capit Tete thick Conn theta
‘tres e_brrtetontionfstr
ree
Compa rive

Figure 15.Flowchart of iteration scheme using total-stress approach to solve for unknown stresses and displacements.


50-m opening 100-m opening 200-m opening
40 —+ 407-— — 40 T —
KEY
o 35[--Linearelements only | 35F 35
é — Both strain- softenin L
= 307 "ond strain-hardeniag 3° ! 7
oe ast elements 25
® ' 20
8 ! 15
#
i | lok
5
°
1.0
E 9
= 25 7 5 J
a 3
20 4 at 7°
5 oy 4 8
Ww L
@ is 3 5
ir 4
Q 10 2 3
az ost 2
g° | Jl
So 6 - Cy) °.
970 990 1010 1,030 050 1070 950 990 1,030,070 I,10 4150 900 1,000 1,100 1,200 1,300

DISTANCE, m

Figure 33.—Comparison of computed stresses and displacements using linear and nonlinear material properties.

ANALYTIC CHECKS ON ENERGY
RELEASE CALCULATIONS

The same problem shown in figure 23 also provides
analytic checks on the accuracy and reliability of the
energy release rate implementation in MULSIM/NL.
Analytic energy release rate expressions exist for a thin,
two-dimensional crack of width 2c subject to a farfield
stress P perpendicular to the crack plane. For the two-
dimensional crack plane shown in figure 23, equation 47
gives the convergence across the opening as a function of
position, modulus, and farfield stress. As the opening
width (W=2c) increases, a critical width W., is attained at
which the convergence just equals the opening thickness.
Solving equation 47 for this critical width gives:

HE

HE 48
2(1-v)P &

Wer =

where H

crack or seam thickness,

E = Young’s modulus,

v= Poisson’s ratio,
and P = farfield stress perpendicular to plane.

Brady and Brown (7) provide the required energy release
rate expressions for the two-dimensional crack shown in
figure 23 which approximates thin tabular excavations such
as longwall panels. For a panel width less than critical
(W < W,,), the energy release rate is

(-
2E

py 2.
Wet Pw. (49)

For a panel width greater than critical (W > W,), the
energy release rate approaches the following value
asymptotically:

Wp = HP. (50)
Equation 50 shows that, well beyond a certain critical

width, this mining geometry can approach steady-state
conditions and achieve a uniform energy release rate.


Figure 34 shows the coarse-mesh grid used to simu-
late extension of a two-dimensional crack plane with
MULSIM/NL. In this model, the crack width extends
from 50 to 500 m in 10 equal 50-m steps. Figure 34 shows
just a few of these steps for clarity. Since the crack length
always exceeds the crack width by a factor of 3 or more,
the three-dimensional numerical model approximates the
two-dimensional analytic solution for the energy release
rate. Both the analytic solutions and the numerical model
use the following input parameters:

Young’s modulus - 3,750 MPa
Poisson’s ratio - 0.25

Farfield stress - 10 MPa
Thickness - 1 m.

Figure 35 compares the ERR calculated numerical-
ly with MULSIM/NL to those determined analytically.
Equation 48 gives a critical opening width of 200 m. For
widths less than critical, the ERR calculated with equation
49 increases linearly to about 8 MJ/m?. For widths great-
er than 200 m, the ERR approaches an asymptote of
10 MJ/m*. As shown in figure 35, the 10-step numerical
ERR agree reasonably well with the analytic values. The
numerical values fall short of the analytic by about 10%.
These discrepancies are attributed to end effects. The
two-dimensional analytic solutions apply to an opening of
width W and infinite length, whereas the three-dimensional

#30 blocks ———>|

KEY

. 9 Fine-mesh
area

(7 Elastic

GA area

RN Open area

/k—30 blocks

50-m blocks
~ MP
I lOth step modulus °
I 7th step \-m thickness
—4th step
r— Ist step

Figure 34.—MULSIM/NL grid for widening-crack problem.

numerical models have a finite length. Calculated dis-
placement profiles at the ends of the opening are some-
what lower than those in the middle; hence, ERR are also
somewhat low. The agreement between the numerical and
analytic results provides good evidence that the energy
calculations in MULSIM/NL function properly. In addi-
tion, this test also serves to confirm the accuracy of the
stress and displacement calculations within MULSIM/ NL,
since these calculations must be correct in order for the
subsequent energy calculations to be correct.

STEP-SIZE EFFECT TESTS

As discussed in a prior section, the total energy release
is composed of a kinetic energy release component and a
strain energy release component. Salamon (6, 40) showed
theoretically that the relative proportion between these two
components depends on the mining step size. When using
smaller mining areas per step to achieve a given state (ie.,
using more mining steps to reach a given excavation geom-
etry), the relative proportion between kinetic and strain
energy release changes, although the total energy release
remains constant. In the limit, as the mined areas per step

OPENING WIDTH (W), m
100 200 300 400
T | TT

ERR asymptote
for W>Wor

500

iy
°

6

. Numerical ERR with
Analytic for MULSIM/NL

WS Wer

TOTAL ENERGY RELEASE RATE (ERR), Md/m2
a

ab 4
L —Critical opening
2 width (Wer)
1 1 1 1
fe) 2 4 6 8 10

MINING STEP NUMBER _

Figure 35.—Comparison of analytic and numerical energy
release rates for widening-crack problem.



STRESS (¢), MPa

Figure D-1.—Computed stresses at edge of opening using
various linear and nonlinear material models. A, 50 m; B, 100 m;

APPENDIX D.—SUPPLEMENTAL DATA ON NONLINEAR EFFECT TESTS

40 T T T T
4 KEY
35 Linear only
Softening only
30

Hardening only
Both

1010-4030 1,050

4070

1150

C, 200 m.

1,100
DISTANCE, m

4,200

1300

030

25

T T T T

KEY
=== Linear only
—-— Softening only
Hardening only
Yt!

1,070

wu >

DISPLACEMENT (D), m.
iv

1030-4070 Ilo

1,150

a @ oN ® © 0 oO

°
900

T T T

1,000 1,100 1,200
DISTANCE, m

1,300

Figure D-2.—Computed displacements across opening using
various linear and nonlinear material models. A, 50 m; B, 100 m;

C, 200 m.


ie., the total energy release equals half the change in
gravitational potential again, in concurrence with Cook (4).

However, mining does not occur in giant steps, rather
it occurs in small steps that approach zero size at the limit.
In this case, Salamon (6) proves that the kinetic energy re-
lease tends to zero, and, therefore, the total energy release
is nothing more than the strain energy release from the
mined-out rock! This result is plainly visible if one con-
siders a weightless rock mass (X, = 0) and a linear elastic
gob and/or backfill (2 = 1) in which case equation 41
reduces to

1 1
Wpe 3Js,74s + a (43)

At the limit, as the mining step size approaches zero, Au
tends to zero, and equation 43 can be approximated by

1
Wr X als, Tyujds. (4)

Figure 19, modified from Salamon (40), illustrates the
step-size dependence of the energy release components,
namely strain energy release U,, and kinetic energy release
W that are normalized by the total energy release. For
this example, Figure 19 shows that the kinetic energy term
decays rapidly as the mining step sizes become smaller.
As the mining step size becomes small, the kinetic energy
release tends to zero, and the total energy release becomes
equal to the strain energy release from the mined-out ma-
terial. Thus, the kinetic energy release component of the

RATIO OF STARTING TUNNEL RADIUS
TO FINAL TUNNEL RADIUS

0 05 0750875 09375096875 0984375

|ORMALIZED

Un/ We

ok.
ried

<es

Fou

3

eet:

nEa

g2q J
aye

B23

so <

>J20

Bae ab N 4
eee N

§

Sea

aB2

&

ZZO oF

gar 2

ox

bE®

a

2 °

x i) 2

NUMBER OF MINING STEPS TO EXCAVATE
CIRCULAR TUNNEL,

Figure 19.—Step-size dependence of energy release compo-
nents for radial expansion of circular tunnel.

total energy release cannot serve as an energy source for
devastating rock bursts or coal bumps since that term dis-
appears with small step sizes. Note, however, that the
total energy release, Wa which is the sum of U,, and Wi,
remains constant and independent of mining step size! In
the limiting case of small step sizes, the energy release is
nothing more than the strain energy release from the rock
excavated during the current step. That strain energy is
dissipated completely as heat during the excavation process
and is, therefore, totally unavailable to drive rock bursts or
coal mine bumps as was originally postulated by Cook (4).

These insights into total energy release and the associ-
ated ERR do not diminish its utility as a practical indica-
tor of rock burst or coal mine bump potential, The old in-
terpretation of energy release made the supposition that
high total energy releases corresponded to high kinetic en-
ergy releases which, in turn, might drive rock bursts. As
was shown by Salamon (6), this is not true since in the
limit, the total energy release approaches the strain energy
release from the mined-out material. The correlation be-
tween energy release and rock bursts arises because rock
bursts are most frequently associated with highly stressed
rock. Mining out highly stressed rock must result in a
high strain energy release which, in turn, must equal a
high total energy release. Thus, the total energy release
merely serves as an indicator of the stress level in the
rocks mined-out this step. Higher energy releases cor-
relate to higher rock stress states which in turn correlate
to greater rock burst hazards.

PROGRAMMING ENERGY RELEASE
CALCULATIONS IN MULSIM/NL

The step-wise energy release equation (41 or A-17)
derived by Salamon (6) provides the basis for energy cal-
culations in MULSIM/NL. This section discusses the pro-
gramming of this relationship into a practical boundary-
element program. The discussion begins with an overview
of the intermediate stress-displacement files written by
MULSIM/NL. These files serve as the essential input to
subroutine ENERGY which performs all MULSIM/NL
energy calculations. These energy calculations proceed on
two different levelslocal (i.c., for each element) and
global (ie., for the entire BEM). Subroutine ENERGY
then generates final output files from MULSIM/NL that
contain calculated stresses, displacements, and the local
energy quantities on an element-by-element basis. Last,
ENERGY reports global energy changes for the entire
model in accordance with equation 41 (i.e. total energy
release, strain energy release, and kinetic energy release
for each mining step).

Implementation of equation 41 is not as straightforward
as it seems. Three different cases can arise that com-
plicate the programming logic. These cases are discussed


true strains. Derivation of the connection between true
and effective moduli and strain follows.

Consider a zone with thickness t subject to a compres-
sive displacement of 0.1t. Suppose the NT is the true
thickness of the zone and it is subject to a vertical stress
g,. The true strain is

ae” Ton e7
True stress is simply
oy = Eyer = Ey [ia] (E-8)
Ta
where By = true modulus.
The effective strain is
@ = Se ONLS. (E-9)
Effective stress is
oy = Ene = Ep Fal (E-10)
10
where Ey = effective modulus.
Equating stresses yields
mbsf] =
or
Ey = E,/n.

(E-12)

Therefore, the true modulus values require modification by
factor n which is equal to the broken zone height divided
by the seam thickness.

In terms of true elastic moduli, the nonlinear stress-
strain relation for gob is

= initial modulus,

iv]
cs
0

final modulus,

o, = virgin vertical stress,

i=}
u

seam closure,
t = seam thickness,

and n

gob height factor.
IMPLEMENTATION INTO MULSIM/NL

Subroutine MODHARD implements this nonlinear
stress-strain model for gob in MULSIN/NL. An array
called EPROP stores material property data for up to 26
different materials. EPROP has dimensions 26 by 10 so
each material can have up to 10 descriptors in its material
model. The first of the 10 material descriptors in EPROP
acts as a flag identifying the material model to use, which
is five for strain-hardening. For the strain-hardening case,
the EPROP descriptors are as follows:

EPROP(*,1) = 5 ~ material type strain-hardening

EPROP(*,2) =~ initial modulus
EPROP(*,3) = - final modulus at stress
EPROP(*,4) =~ virgin vertical stress
EPROP(*,5) = - gob height factor

EPROP(*,6) = — Poisson’s ratio
Descriptors 7 through 10 remain unused.

Subroutine MODHARD returns values for the Young’s
modulus perpendicular to the seam plane and the shear
moduli parallel to the seam plane based on current esti-
mates of the stresses. The user specifies material prop-
erties normal to the seam only. The associated shear
moduli are calculated as

E
S- oa (E-14)

While the material property specifications for the nor-
mal and shear directions for an element are linked, the
boundary conditions for these directions are decoupled and
act independently. Coupling these boundary conditions
may have certain intuitive appeal; however, it will over-
specify the system (i.e., more equations would exist than
unknowns). Such a procedure may lead to numerical diffi-
culties of unknown proportions.


CONTENTS

Page

Purpose and objective
MULSIM/NL enhancements . .
SCOPE oe eee ee cece cceee eee eee eee e es

Origins and previous versions of MULSIM/NL .
Historical review
Sinha’s MULSIM .
MULSIM/BM .
Itasca additions

Nonlinear material properties in MULSIM/NL
General mechanics of boundary-element method .

Basic equations of boundary-clement method .
Application of boundary conditions .
Some notes on BEM boundary conditions
Six material models for MULSIM/NL
Linear clastic for coal .. .
Linear elastic for gob .. .
Strain-softening for coal .
Elastic-plastic for coal . .
Bilinear hardening for gob .
Strain-hardening for gob . .
Total-stress formulation for MULSIM/NL see
Induced-stress boundary conditions .
Total-stress boundary conditions . .
Illustration of convergence schemes

FORTRAN implementation of material nonlinearities 16
Equation-solver logic 16
Material property subroutine structure .... 17

Nonlinear energy calculations in MULSIM/NL .. 18

Evolution of energy release rate concepts . 18

Recent energy release rate concepts 19

Programming energy release calculations in MULSIM/NL 2
Intermediate stress-displacement files . 22
Local energy calculations . . 2
Global energy calculations . 23
Three special cases . 2B
Subroutine energy logic diagram : 25

MULSIM/NL program checks and behavior . 26

Some simple analytic solutions 26

Rigid boundary effects . 26

Element size effects 29

Nonlinear effects . . 31

Subroutine ENERGY checks and behavior : . 31

Analytic checks on energy release calculations . 32

Step-size effect tests 33

Nonlinear effects on energy calculations . 34



REFERENCES

1, Campoli, A. A., C. A. Kertis, and C. A. Goode. Coal Mine
Bumps: Five Case Studies in the Eastern United States. BuMines
IC 9149, 1987, 34 pp.

2. Sinha, K. P. Displacement Discontinuity Technique for Analyzing
Stresses and Displacements Due to Mining in Seam Deposits. Ph.D.
‘Thesis, Univ. MN, Ann Arbor, MI, 1979, 311 pp.

3. Beckett, L. A,, and R. S. Madrid. MULSIM/BM - A Structural
Analysis Computer Program for Mine Design. BuMines IC 9168, 1988,
302 pp.

4. Cook, N. G. W., E. Hoek, J. P. G. Pretorius, W. D. Ortlepp, and
M.D. G. Salamon. Rock Mechanics Applied to the Study of Rock
Bursts. J. S. Afr. Inst. Min, Metall., v. 66, No. 10, 1966, pp. 435-528.

5. Walsh, J. B. Energy Changes Due to Mining. Int. J. Rock Mech.
Min. Sci., v. 14, 1977, pp. 25-33.

6. Salamon, M. D. G. Energy Considerations in Rock Mechanics:
Fundamental Results. J. S. Afr. Inst. Min. Metall., v. 84, No. 8, 1984,
pp. 233-246.

7. Brady, B. H. G., and E. T. Brown. Rock Mechanics for Under-
ground Mining. George Allen and Unwin (London), 1985, 527 pp.

8. Hodgson, K., and N.C. Joughin. The Relationship Between En-
ergy Release Rate, Damage, and Seismicity in Deep Mines. Paper in
Proceedings of 8th Symposium on Rock Mechanics, Univ. MN, ed. by
C. Fairhurst. AIME, 1967, pp. 194-203.

9. Cook, N.G. W. Rockburst and Rockfalls, Chamber of Mines of
S. Afr., Johannesburg, SA, Publ. No. 216, 1978.

10. Crouch, S. L., and A. M. Starfield. Boundary Element Methods
in Solid Mechanics. George Allen and Unwin (London), 1983, 322 pp.

11. Banerjee, P. K., and R. Butterfield, Boundary Element Methods
in Engineering Science. McGraw-Hill (New York), 1982, 452 pp.

12. Brebbia, C. A., and S. Walker. Boundary Element Techniques in
Engineering. Newnes-Butterworth (London), 1980, 323 pp.

13. Brebbia, C. A., and J. J. Connor. BEM-11. Proceedings of 11th
International Conference on Boundary Element Methods in Engineering,
Springer-Verlag (Berlin FRG), 3 volumes, 1989, 1,217 pp.

14, Wardle, L. J. Boundary Element Methods for Stress Analysis of
Tabular Excavations. Ph.D. Thesis, Univ. Queensland, Australia, 1986,
11 pp.

15. Lachat, J. C., and J. O. Watson, Effective Numerical Treatment
of Boundary Integral Equations: A Formulation for Three-Dimensional
Elastostatics. Int, J. Numer. Meth. Eng,, v. 10, 1976, pp. 991-1005.

16. Vogel, S. M., and F. J. Rizzo. An Integral Equation Formulation
of Three-Dimensional Anisotropic Boundary Value Problems. J. Elast.,
v. 3, 1973, pp. 203-216.

17, Telles, J. C.F. The Boundary Element Method Applied-to In-
elastic Problems. Lecture Notes in Engineering, v. 1, Springer-Verlag
(Berlin), 1983, 243 pp.

18. Zienkiewicz, O. C. Finite Element Methods in Engineering Sci-
ence. 3rd ed., McGraw-Hill (London), 1977, 787 pp.

19. Brady, B. H. G., and A. Wassyng. A Coupled Finite Element-
Boundary Element Method of Analysis. Int. J. Rock Mech. Min. Sci.
and Geomech. Abstr., v. 18, 1981, 475-485 pp.

20. Berry, D. S. An Elastic Treatment of Ground Movement Due to
Mining - Part I. Isotropic Ground. J. Mech. Phys. Solids, v. 8, 1960,
pp. 280-292.

21. Berry, D. S., and T. W. Sales. An Elastic Treatment of Ground
Movement Due to Mining - Part II. Transversely Isotropic Ground.
J. Mech. Phys. Solids, v. 9, 1961, pp. 52-62.

22. Berry, D. S., and 'T. W. Sales. An Elastic Treatment of Ground
Movement Due to Mining - Part III. Three-Dimensional Problem,
Transversely Isotropic Ground. J. Mech. Phys. Solids, v. 10, 1962,
pp. 73-83.

23, Salamon, M. D. G. Elastic Analysis of Displacements and Stresses
Induced by the Mining of Seam or Reef Deposits. J. S. Afr. Inst. Min.
Metall, v. 64, Nov. 1963, pp. 128-149.

Elastic Analysis of Displacements and Stresses Induced
ig of Seam or Reef Deposits - Part II. Practical Methods
of Determining Displacement, Strain and Stress Components from a
Given Mining Geometry. J. S. Afr. Inst. Min. Metall, v. 64, Jan. 1964,
pp. 197-218.

25. Elastic Analysis of Displacements and Stresses Induced by
the Mining of Seam or Reef Deposits - Part IV. Inclined Reef. J. S.
Afr. Inst. Min, Metall., v. 65, Dec. 1964, pp. 319-338.

26. Plewman, R. P., F. H. Deist, and W. D. Ortlepp. The Develop-
‘ment and Application of a Digital Computer Method for the Solution of
Strata Control Problems. J. S. Afr. Inst. Min. Metall., v. 70, Sept. 1969,
pp. 33-44.

27. Starfield, A. M., and S.L. Crouch. Blastic Analysis of Single Seam
Extraction. Paper in Proceedings of 14th U.S. Symposium on Rock
Mechanics, Penn. State Univ. ASCE (New York), 1973, pp. 421-439.

28. Crouch, S. L., and C. Fairhurst, The Mechanics of Coal Mine
Bumps and the Interaction Between Coal Pillars, Mine Roof, and Floor
(contract H0101778, Univ. MN). BuMines OFR 53-73, 1973, 88 pp.;
NTIS PB 222 898.

29. Crouch, S. L. Two-Dimensional Analysis of a Near Surface,
Single-Seam Extraction. Int. J. Rock Mech, Min. Sci. and Geomech.
Abstr., v. 10, 1973, pp. 85-96.

30. ___. Analysis of Stresses and Displacements Around Under-
ground Excavations: An Application of the Displacement Discontinuity
Method. Univ. MN Geomech. Rep. Dept. of Civ. Miner. Eng., Univ.
MN, Minneapolis, MN, 1976, 268 pp.

31. Kripakov, N. P., L. A. Beckett, D. A. Donato, and J. S, Durr.
Computer-Assisted Mine Design Procedures for Longwall Mining. Bu-
Mines RI 9172, 1988, 38 pp.

32. Kripakov, N. P., and L. E. Rockwell. Application of Numerical
Modeling Techniques for Gate Road Design in Bump Prone Mines.
Paper in Proceedings of Longwall USA, Pittsburgh, PA, Industrial
Presentations, Inc., Aurora, CO, 1989, pp. 87-99.

33. Dahlquist, G., and A. Bjork. Numerical Methods. Prentice-Hall,
1974, S71 pp.

34, Wang, F. D., W. A. Skelly, and J, Wolgamott. In Situ Coal Pillar
Strength Study (contract H0242022, CO Sch. Mines), BuMines OFR.
107-19, 1976, 243 pp.

35. Tannacchione, A. Behavior of a Coal Pillar Prone to Burst in the
Southern Appalachian Basin of the United States. Paper in Proceedings
of Second International Symposium of Rockbursts and Seismicity in
Mines. Univ. MN, June 1988, pp. 427-439.

36. Heasley, K. A., and K. Barron. A Case Study of Gate Pillar Re-
sponse to Longwall Mining in Bump-Prone Strata. Paper in Proceedings
of Longwall USA, Pittsburgh, PA, 1988, Industrial Presentation, Inc.,
Aurora, CO, 1988, pp. 91-105.

37. Cook, N. G. W. The Basic Mechanics of Rockbursts. J. S. Afr.
Inst. Min. Metall, v. 64, Oct. 1963, pp. 71-81.

38. ___. The Design of Underground Excavations. Paper in Pro-
ceedings of 8th Symposium on Rock Mechanics, Univ. MN, ed. by
C. Fairhurst. AIME, 1967, pp. 167-193.



(6) Changing the modulus of the backfill and/or gob,
whether linear or nonlinear, had no perceptible effect on
the strain ERR values. In all cases considered, the strain
ERR component of the total ERR remained negligible.

The changing open area tests (i.e., the last three models
in figure 38) led to many observations similar to those of
the constant area tests. However, the strain ERR and the
nonlinear kinetic ERR are negligible in the total ERR due
to the large open area to backfill and/or gob area ratio in
these models. Increasing the modulus of the backfill (gob)
always tends to decrease the total ERR and its dominant

component, the linear kinetic ERR. Also, changes in
backfill modulus or its degree of nonlinearity has an
almost negligible effect on the calculated strain ERR.

In addition to the prior observations, the changing area
tests suggest that decreasing the open area between mining
steps (causing the backfilling or gob formation activities to
accelerate) has little effect on the total ERR. Conversely,
increasing the open area between steps (allowing the back-
filling or gob formation activities to lag) causes the total
ERR to increase. In this example, that increase is on the
order of 10%. Such an increase may suggest greater burst
and/or bump potential.

SUMMARY AND RECOMMENDATIONS

SUMMARY

MULSIM/NL is a new Bureau boundary-clement-
method program for stress analysis of coal mines and thin
metalliferous veins. This Information Circular provides
certain theoretical documentation on that program. Ac-
cordingly, its intended audience is those engineers and
programmers who need to understand the FORTRAN
program fully, or who need to alter and enhance it,

MULSIM/NL is a three-dimensional boundary-
element-method (BEM) program capable of performing
stress analyses of planar structures in coal seams or thin
veins. It can analyze one to four parallel seams that can
have any orientation with respect to the Earth's surface.
Topographic or free surface effects are neglected; there-
fore, the seam depth should exceed the modeling block
width by a factor of 3 or more. MULSIM/NL utilizes a
coarse-mesh grid of blocks to approximate the larger
model bounds and a fine-mesh grid of elements embedded
within the coarse mesh to provide greater numerical detail
within a region of interest. The allowable coarse-mesh
block array size is 40 by 40. Up to a 20 by 20 array of
these blocks may be subdivided into a 100 by 100 array of
fine-mesh elements.

Three main new features distinguish MULSIM/NL
from its predecessors MULSIM (2) and MULSIM/BM
(), namely: (1) nonlinear material models, (2) multiple
mining steps, and (3) comprehensive energy release and
strain energy computations, MULSIM/NL now has six
material models for the in-seam material from which to
choose: (1) linear elastic for coal, (2) strain softening,
(3) elastic plastic, (4) bilinear hardening, (5) strain hard-
ening, and (6) linear elastic for gob. The first three mod-
els are intended for the in-place coal or vein material

prior to mining, whereas the latter three models are in-
tended for the gob or backfill material left in the wake of
full extraction mining. The strain-hardening model for
gob and/or backfill material derived and implemented as
part of this research may be a unique addition to the
variety of material models available in mining-oriented
BEM programs.

Unlike its predecessors, MULSIM/NL uses a total
stress approach when solving the basic BEM system of
equations. Such radical changes in the equation-solving
procedure, plus the addition of the six nonlinear material
models, necessitated detailed checks on the MULSIM/NL
stress and displacement calculations. Numerical stress and
displacement calculations compare reasonably well to
known analytic solutions for simple problems. MULSIM/
NL models exhibit rigid boundary and element size effects.
Until studies with a particular model dictate otherwise,
having at least 10 coarse-mesh blocks beyond the fine-
mesh boundary will usually keep the rigid boundary out-
side the model far enough away and lead to good quality
calculations in the fine mesh. Smaller element sizes tend
to improve the accuracy of stress and displacement cal-
culations; however, they can also lead to deleterious rigid
boundary effects because of the smaller model size. Exact
checks on the nonlinear material models are not feasible;
however, stresses and displacements calculated with these
models behave qualitatively as expected.

The second distinguishing feature of MULSIM/NL is
the multiple mining step capability which enables simu-
lation of the various stages of mine development. This
feature allows the user to calculate stress and displace-
ment changes as the mine development advances. These
changes are much more directly comparable with field
measurement programs that typically measure stress and


assigned material properties (or boundary conditions) to
approximate this geometry. Figure 7 shows an actual mine
plan and the modeling approximation where individual ele-
ments are considered either mined out or unmined with
assigned material properties.

The central foundations to BEM rest on certain analytic
solutions from basic elasticity theory. In two dimensions,
one of these solutions is for a crack of width 2a subjected
to uniform pressure components P, and P,. Figure 8
shows this pressurized crack problem. By integrating ana-
lytic solutions to the Kelvin problem, expressions are
obined for the stress and displacement components of the
pressurized crack in terms of the crack pressure compo-
nents P, and P, and the spatial coordinates x and z.

Another even more important solution is that of a
"displacement discontinuity” in an infinite solid as shown
in figure 9. Instead of applying constant pressure com-
ponents to the crack walls, constant displacement com-
ponents are applied resulting in a so-called "displace-
ment discontinuity.” Again, it is possible to find analytic

Figure 6.—Boundary-element-method seam plane subdivided
by elements.

expressions for the stress and displacement components in
terms of the displacement discontinuity components D,
and D, and the spatial coordinates x and z.

Basic Equations of Boundary-
Element Method

Figure 9 shows a two-dimensional displacement discon-
tinuity; however, it is also possible to extend the concept
and solutions to three dimensions. The solution forms,
following Sinha’s (2) notation, are

1, = bD,+eDy,

4

aD, +eDy,

=aD,, @)

Figure 8.—Pressurized crack problem in two dimensions. 2a
Is crack width, and P, and P, are pressure components.

=N

Dz

U

Ul

oOo

Actual mine plan Modeling approximation

Figure 7.—Actual mine plan in-seam plane and gridded model-
Ing approximation.

20

Figure 9.—Displacement-discontinulty problem In two dimen-
sions. 2a is crack width, and D, and D, are displacement com-
ponents,


does not transmit any significant vertical stresses to the
floor rock. In the second step, as the longwall face moves
farther away, the gob changes gradually from a swelling
phase to a consolidation phase as it begins to experience
higher vertical stresses. Ultimately, the vertical stress in
the mined-out region may approach the virgin vertical
stress and further vertical deformation or consolidation
ceases,

DERIVATION OF NONLINEAR STRESS-STRAIN
RELATION FOR GOB

Two major assumptions underlie the gob model pre-
sented here. First, all the vertical strain that results in
surface subsidence occurs within a broken zone whose
height varies from two to six times the seam thickness.
Justification for this assumption was presented carlier.
The second assumption states that the elastic modulus of
the gob material in the broken zone increases linearly
from some initial value E, to some final value Ep. The
initial modulus is that of the swelled gob material at zero
vertical pressure, whereas the final modulus is that of the
consolidated gob material at the virgin vertical stress o,.
The assumed linear relationship between tangent modulus
and vertical stress is shown in figure E-3.

The equation describing this linear relationship is

Ep -E,
E= F/ +E.
oy

(2)

The differential relation between stress and strain is
de = Lac. (B3)
E

Substituting E-2 into E-3 and integrating between the
limits 0 to € and 0 to o yields

PO? ail, ea
lz} es
Rearranging and exponentiating gives

le} ©

Ep -E,
ala a
FO Ey ey

From figure E-2, the strain in the broken zone is calcu-
lated as

e=

2B (E-6)
t
where D is the seam closure and t is the seam thickness.

THE GOB HEIGHT FACTOR “n"

The above formulation utilizes effective modulus values
and effective strain as opposed to true modulus values and

| vind} pen im "Swelling — Consolidating |

Figure E-2.—Swelling and consolidation mechanism of gob
formation. T is seam thickness; NT is totally fractured zone
thickness where N is gob height factor; and S is surface sub-
sidence assumed as 0.5 T.

ry
° -
3 “
5 a
8 ERp-------5
= 1
_ I
z 1
ai
8 !
4 1
z
= |
Ey !
1
L
ov

STRESS, (o)

Figure E-3.—Assumed linear relationship between tangent
modulus (E) and stress (2) over range 0 to oy E; is initial
modulus at zero vertical stress, and Eyis final modulus at vertical
stress oy


For the mined-out or open elements, the elemental
stresses are zero, and the boundary conditions are

0)

Thus, the boundary conditions given by equations 6 and
7 provide the additional 3N relations, which together with
the boundary-element equations 5, enable computation of
the remaining unknown elemental stresses and displace-
ment discontinuities. Numerical solution of this system of
linear algebraic equations is normally accomplished with
a Gauss-Seidel iteration procedure. See Dahlquist and
Bjork (33) for details of this method.

Some Notes on BEM Boundary Conditions

The three boundary conditions on an element must act
independently of one another. Thus, the normal and shear
components of stress and displacement applied to an ele-
ment are uncoupled. Coupling the normal and shear com-
ponents with some complex constitutive relation may result
in an overdetermined system of equations with unknown
consequences on the solution, Also note (fig. 11) that hor-
izontal components of stress do not exist in the in-seam
elements, Therefore, the engincer and/or analyst using
BEM must find artificial means to account for the effects
of horizontal stress or confinement on the behavior of an
element. One common procedure involves the use of a
family of stress-strain curves ranging from highly strain
softening through elastic plastic to represent materials
under increasing confinement. The highly strain-softening

Figure 11.—Three boundary conditions applied to elemental
ement discontinulty. D,, D,, and D, are elemental dis-
placement discontinuity components.

material is used near an opening where the confinement
is low, and the elastic-plastic material is used in the
interior of a model where the confinement is high.

SIX MATERIAL MODELS FOR MULSIM/NL

As alluded to in the previous discussions on the formu-
lation of BEM nonlinear in-seam material properties are
incorporated into the model via the boundary conditions
on the system of boundary-clement equations. This sec-
tion discusses theoretical and programming details of the
nonlinear material models added to MULSIM/NL., Initial
discussions focus on the derivation and formulation of the
six material models now available in MULSIM/NL. Addi-
tional mathematical details provide the essence of total-
stress versus induced-stress BEM formulations. Imple-
menting the material nonlinearities required a shift from
the induced-stress approach used in MULSIM to the total-
stress approach used in MULSIM/NL. Subsequent discus-
sions then focus on the logic, algorithms, and convergence
schemes necessary to actually program the models. Final-
ly, the section closes with certain programming details of
the material property subroutines and the total-stress
BEM formulation.

Figure 12 shows the stress-strain models for MULSIM/
NL that include (1) linear elastic coal, (2) strain-softening
coal, (3) elastic-plastic coal, (4) bilinear hardening gob, (5)
strain-hardening gob, and (6) linear elastic gob. The first
three models are for the unmined in-seam coal material,
while the latter are for the broken gob material left in the
wake of mining. As will be discussed later, the distinction
between coal material models and gob material models
becomes crucial in the subsequent energy calculations.

Linear elostic coo!

Strain-softening coal _Elastic plastic coal

7 2 3

Linear elastic gob __Strain-hardening gob Bilinear hardening gob

6 $ 4

STRESS (o)

STRAIN (¢)

Figure 12.-Six different stress-strain models available in
MULSIM/NL.


Library of Congress Cataloging in Publication Data:

Zipf, R. Karl.
MULSIM/NL theoretical and programmer’s manual / by R. Karl Zipf, Jr.

p. cm. — (Information circular; 9321)
Includes bibliographical references (p. 40).
Supt. of Docs. no.: I 28.27:9321.-
1, MULSIM/NL. I. Title. II. Series: Information circular (United States. Bu-

reau of Mines); 9321.

TN295.U4 [TN288] 622 s—dc20 [622’.28'0285369] 91-47152 CIP



NONLINEAR EFFECTS

Figure 32 shows the model geometry used to illustrate
the qualitative behavior of the nonlinear material models
in MULSIM/NL. Again, the problem considered is an in-
finitely long, planar opening of width 2c. The model uses
a 30- by 40-block coarse-mesh with a central fine mesh
measuring 4 by 20 blocks (20 by 100 elements). Opening
width increases from 50 to 100 m and then to 200 m. Two
different material property sets are compared. The first
set is completely linear. Young’s modulus for the rock
mass and the seam is set to 3,450 MPa. The second set
includes the linear elastic properties of the first and also
introduces strain-hardening and strain-softening elements.
For the hardening elements placed in the mined-out area,
initial modulus is 28 MPa and final modulus is 140 MPa at
a stress of 7 MPa. Table 2 gives the four types of soft-
ening elements used by the second set. The first type is
used in the first row of elements adjacent to the opening;
the second type in the second row and so forth.

Table 2.—Strain-softening properties used In nonlinear models

Peak Residual
Type Stress, Strain Stress, Strain
MPa MPa
1. 6.90 0.0020 275 0.0083
2 13.80 0040 6.90 0080
3 24.15 0070 20.70 0120
4. 34.50 .0100 34.50 1000

Finally, figure 33 compares calculated stresses and dis-
placements using the strictly linear material property set
with the nonlinear set. Using strain-softening elements
along the opening edge decreases the peak edge stress
dramatically. The magnitude of the peak stress tends to
increase with opening width as expected. The increase is

b+————-40 blocks ————4}

50-, 100-, and 200-m
openings with 12.5-m
blocks

30 blocks ——-

KEY
Fine-mesh blocks

Open blocks

Figure 32—MULSIM/NL models for nonlinear effect studies.

ificant (about 40%) for the linear model and somewhat
less (about 25%) for the nonlinear model. Presence of the
hardening elements causes the mined-out area to increase
its stress level; however in this case, the stress levels are
still much less than the in situ vertical stress of 7 MPa.
Hardening elements also cause the vertical displace-
ment to decrease (up to 25%) in the mined-out area as
expected.

In summary, while exact checks are not readily possible
on the nonlinear capabilities of MULSIM/NL, the model
does behave as expected when these nonlinearities are
utilized. Stress profiles in the strain-softening areas re-
flect the strain-softening material input parameters. Cal-
culated stresses and displacements for an element always
lie on or near the prescribed nonlinear stress-strain curve
of that element. These validation checks, while far from
exhaustive, lend strong credibility to the correctness of
MULSIM/NL, both in concept and in programming.

SUBROUTINE ENERGY CHECKS AND BEHAVIOR

This section examines the energy subroutine in
MULSIM/NL with three different test sets. The purpose
of these tests is to validate the energy calculations and to
examine the many factors affecting them, The first tests
enable a rigorous comparison between MULSIM/NL en-
ergy calculations and the known analytic solution to a very
simple problem. The second tests study the influence of

mining step size on the relative proportion of energy
components comprising the total energy release. Finally,
the last tests examine quantitatively the effect of certain
key geometric and material property values on the total
energy release and its components. These final tests are
a simple parametric study of the energy release.


account for energy dissipation due to cracking in the rock
mass, Brady and Brown (39) provided practical mine de-
sign applications of BEM incorporating energy change cal-
culations. Finally, energy release concepts have been
included in rock mechanics text books (7).

RECENT ENERGY RELEASE RATE CONCEPTS

Recent theoretical work by Salamon (6) demonstrated
that certain ERR concepts require serious revision. In the
1960's, Cook (4, 37-38) recognized that when an excavation
is created, half of the change in gravitational potential
energy (i.c., half the work done by the external body
forces) is stored as elastic strain energy in the rock mass.
Therefore, as seen in the prior derivations, to preserve
energy conservation, the other half of the gravitational
potential energy change is necessarily released.

These early researchers hypothesized that the released
energy was cither dissipated nonviolently as fracturing and
other inelastic behavior in the rock mass, or it was trans-
formed into kinetic energy (i.e., seismic energy). Some of
these kinetic or seismic events might, in turn, manifest
themselves as damaging rock bursts. The early derivations
of energy release assumed, erroneously, that opening crea-
tion occurs suddenly or in one large step.

Salamon’s recent theoretical findings dispel the notion
that the energy release serves as the source for kinetic,
seismic, and/or rock burst energy—it is not. Other sources
for seismic and rock burst energies are required. This ma-
jor revision in thinking concerning energy release, how-
ever, does not necessarily diminish its practical utility as an
indicator of rock burst and/or coal mine bump potential.

Salamon (6, 40) reconsidered the basic energy balance
for mining and, most important, considered the step-wise
nature of mining. Retaining Salamon’s notation, the fol-
lowing derivation captures the essence of his work. As
shown in figure 17, mining progresses from state I to state
II in small steps. The small mined-out volume of rock this
mining step is V,,, and the surface area of the new mined-
out area is Sy. Work, W, is done on the system, as shown
by external forces acting on the system surface, and by
internal body forces acting within the system. Basic energy
conservation requires that this energy supply W is im-
mediately transformed into three energy sinks: strain
energy in the rock mass, strain energy in the opening sup-
ports and backfill, and energy dissipation. Stated sym-
bolically, energy conservation in going from state I to state
lis

W = AU + Wep + Wp, (34)

State 1 Ss

KEY

V__ Rock mass volume

Vor Volume open in state I
Volume open in state IE

Volume of backfill in state I

Vem Volume of backfill in state I

Vm Volume mined this step

S$ _ Rock mass boundary area

State I £s. Sor Surface area open in state I
v Son Surface area open in state IL
[Vor | seal Scr Surface area of backfill in state 1
Sen Surface area of backfill in state IL
~~~~BSoql Sw Surface area mined this step
Ver

Figure 17.—Basic notation for
progresses from state | to state Il.

nergy calculations as mining

where W = total work done on system,
AU = change in strain energy in rock mass
volume V,
Wep = change in backfill and support strain
energy,
and Wa = total energy release going from state I
to state II.

The total strain energy in the rock mass in state I, Uy, is

U,=U+ Uy, 35)
where U_ = total strain energy in rock mass volume
V outside the mined-out volume Vyq
and Uy = total strain energy in rock mass volume

Vm mined-out this step.
The total strain energy in the rock mass in state II, Uy is
Uy = U + Ug, (36)

where induced strain energy in rock mass

volume V for this step.

The change in strain energy for this step is thus
AU = Uy- Uy =U + Up

-(U + Uy) = Uc- Uy. (en)


Information Circular 9321

MULSIM/NL Theoretical and
Programmer’s Manual

By R. Karl Zipf, Jr.

UNITED STATES DEPARTMENT OF THE INTERIOR
Manuel Lujan, Jr., Secretary

BUREAU OF MINES
T S Ary, Director


The following notations are used for displacements,
stresses, surface tractions, and body forces:

u_ = displacement

T = stress

T = surface traction over rock mass areas
R= surface traction over backfilled areas
X = body forces.

In going from state I to state Il, these displacements,
stresses, surface tractions and body forces must obey the
following relations:

Ug =u, + Au

Ty = 7) + Ar

Ty = T, + AT (A-2)
Ry = Ry + aR
Xy = X.

Equilibrium of the rock mass in state I satisfies

fxav + fotias + ft + fis = 0, (A3)

whereas in state II it satisfies

fXnav + fetuds + a Pus =0. (A-4)

The work W done on the system in going from state I
to state II is found by considering the tractions on, the
bounding surface $ and the body forces within the rock
mass V and backfill (gob) Vg as displacements undergo a
change Au. The work done on the system is

_ ca
we fene uds ah TAuds
- fsduav - fuXaeuse

1
“3 JavgeGaudy.

(A5)

The negative signs arise from the sign convention dis-
cussed by Salamon (6).

As indicated in equation A-1, the work done on the sys-
tem depicted in figure A-1 results in strain energy changes
in the three components of that system: the rock mass
mined-out during this step, the surrounding rock mass, and
the backfill and/or gob. The strain energy changes in each
of these areas are considered next.

In going from state I to state II, the strain energy
change (or release) over the mined-out volume V,, bound-
ed by surface Sy is found directly as

1
Um = 5 [Jatin - Sure } (A6)

In state I, the total strain energy in the rock mass vol-
ume V bounded by surfaces $, Sy, and So, is

1
Uy = -4 [I Xyuydv + fetinds

+f Tyyds + fs Ryu; ds } (A-7)
M or
Again, the negative sign arises due to the sign convention
used. The strain energy in state I is one-half the product
of the equilibrium forces in state 1 (equation A-3)
multiplied by the displacements in state I u,.

In state Il, the total strain energy in the rock mass
volume V bounded by surfaces $ and Soy is

1
Uva = “5 [fue

+ ftnds + Soq,8utnds } (48)

Substituting equations A-2 into A-8, expanding the result,
and subtracting A-7 results in the following expression for
the change in strain energy in the rock mass volume V:

Uc = = [feta + Jotiduds + J,erues |

3 [formes +f Rituds + ftRauss]
1 1

1
3 [. gh Reds + fh gPibuds + f. siRauss|

1 1
-5 fXtuav + 5 fot (A-9)


and then to 400 m. All models are subject to a vertical
stress of 7 MPa and use a Young’s modulus of 3,450 MPa
and a Poisson’s ratio of 0.25.

Figure 25 shows computed and analytic displacements
(ie, convergence) across an opening width 2c of 100 m.
Going from the small coarse-mesh (20 by 10 blocks) to the
middle-size coarse-mesh (20 by 20) gives a marked im-
provement to the accuracy of the numerical computations.
Going from the middle mesh size to the large size (20 by
30 blocks) does not result in significant improvement. As
shown in figure 25, numerical displacement calculations
agree with the analytic to within 5%. The discrepancies
are attributed to normal discretization and truncation
errors and do not suggest any programming error.

Figure 26 shows computed and analytic stresses away
from the edge of the 100-m opening. For all different

L+—20 blocks —4} Small coarse mesh
TEE
10
blocks
Bae
+ —20 blocks > Middle coarse mesh
20
blocks
/k—20 blocks—>| Large coarse mesh
FF Crack width 2C
30
es blocks
ail
z KEY
Fine -mesh
+ Baild blocks
Open blocks

Figure 24.—MULSIM/NL models for rigid boundary effect
studies of Infinitely long planar crack of width 2c.

coarse-mesh sizes, the agreement with analytic appears
acceptable. Coarse-mesh size does have a marked effect
on the peak edge stress which increases from 22.5 to
25 MPa (or 10%) as coarse mesh size increases.

Figures 25 and 26 illustrate numerical accuracy and
rigid boundary effects for a simple MULSIM/NL analysis

0.40 — T
25-m blocks
.35f 5m elements

30F

25

DISPLACEMENT (D), m
8

Small, 10 by 20 em |L
Middle, 20 by 20 em |
~ Large, 30 by 20 cm
—— Analytic solution

|
er

O°
950 1,000

1,050
DISTANCE, m

1,100 1,150,

Figure 25.—Computed and analytic displacements across
100-m opening with different coarse-mesh sizes.

35 SEE EEE ER
KEY
30 =-— Small, 10 by 20 em J
25+ P Analytic . 4
Fd ‘solution
=
= 20F il 4
2 ij jl
2 i iF
i Isr ! | z|
= i I
5
lof | 25-m blocks ! |
= | S-metements |
sb i |
|
° deo is
950 1,000 1,050 100 1150

DISTANCE, m

Figure 26.—Computed and analytic stresses at edge of 100-m
opening with different coarse-mesh sizes.


Kinetic

energy release rate, linear component

Kinetic energy release rate, nonlinear component

1.6 T T T T T : T ‘ T T
lap 4 4 a 4
2p 4
oF 7 +f 1
BF i+
6b 1
4h 1
w
2 4 4
F: 7
ea ° J [Le ——— — 4
@S -2l fl ! 1 1 fl 1 1 1 A 1
Ws
dl Strain energy release rate Total energy release rate
© 1.6.4 T T T T T T ¥ T T T
wi
2 at ¢ 7 + 2 7
2
G  o.e2b 1 oF 1
Lor 4 + 4
+° |
.6F 4 4
4h 4 4
2F 4 1
Ob ———_—— a | L & 4
-2lu 1 fl \ 1 1 1 1 fl 1 1 1
Constant open area Decrease —_ Increase Constant open area Decrease Increase
Small Middle Lorge  openarea open area Small Middle Large ‘open area

open area

Figure 40.—Total energy release from mining step 2 to 3 and its components for five different test models with linear backfill

material. Note: graph A + graph B + graph C = graph D.

The constant open area numerical experiments (i.e., the
first three models in figure 38) suggest the following
observations:

(1) The linear kinetic ERR and the strain ERR ap-
pear to vary approximately as the square of the open area
to backfill area ratio for both the linear and nonlinear
backfill models.

(2) The linear kinetic ERR dominates the strain ERR
by an order of magnitude or more. However, neglecting
the strain ERR is not recommended, since it can become

significant with low open-area-to-backfill and/or gob-area
ratios, Such low ratios are perhaps more typical of prac-
tical MULSIM/NL modeling.

(3) The nonlinear kinetic ERR may vary inversely as
the area ratio. Large backfill or gob areas may contrib-
ute a nonlinear kinetic ERR component that actually ex-
ceeds the linear kinetic ERR. Again, small open-area-to-
backfill-area ratios (ie., large backfill and/or gob areas
relative to the open area) probably characterize practical
mining situations better than the large area ratios. Thus,
until further studies dictate otherwise, the nonlinear kinetic


in relation to a flowchart diagramming the energy calcu-
lations in MULSIM/NL.

Equation 41 is repeated below as 45 for more direct
reference.

1
Wr = [5 f, trod - sh 10 |

1 1
+ [3 Janos + 5 Jog: (ta) 4Rauds |~

where T; = rock mass stress in state I,
u, = displacement in state I,
Au = change in displacement from state I to
state II,
AR = change in gob (backfill) stress from

state I to state HI,
X= body force in state I,
Sy = surface area mined this step,
Vy = volume mined this step,
Sgu = surface area of gob (backfill) in state II,
and @ = anonlinearity factor for the gob and/or
backfill (for linear materials, @ = 1,

and for strain-hardening materials,
a< 1).

The first two terms comprise the strain energy release
from the mined-out material. The second term, reflecting
the body forces, is always small relative to the first and is
neglected in all subsequent calculations. As discussed
earlier, the third and fourth terms are called the linear and
nonlinear components of the kinetic energy release. To-
gether, the strain energy release and the kinetic energy
release form the total energy release.

Intermediate Stress-Displacement Files

To evaluate equation 45, stresses and displacements are
required from MULSIM/NL for each coarse-mesh block
and fine-mesh element, and for every mining step in
the sequence. As developed in the previous section,
MULSIM/NL calculates each prescribed mining step using
a variety of nonlinear stress-strain models for the in-
place coal and/or vein material and the gob and/or

backfill material left in the wake of mining. Upon com-
pletion of each mining step, MULSIM/NL writes the
stress-displacement results to intermediate files. Each
seam requires a coarse-mesh and a fine-mesh data file;
thus, two stress-displacement data files are written for
each seam. These files contain the block and/or element
centroidal coordinates, the material property code, three
components of calculated stresses, and three displacement
components. These files provide a convenient way to ac-
cess stresses and displacements for the current mining step
and the previous mining step.

Local Energy Calculations

Subroutine ENERGY enters its main loops and con-
siders cach block and element in the model. For each
block and/or clement, six local energy quantities arc
computed as follows:

(1) Strain energy release.
(2) Kinetic energy release.
Linear part (gravitational).
Nonlinear part (backfill and/or gob).
(3) Total energy release.
(4) Total strain energy.
(5) Recoverable strain energy.
(6) _Dissipated strain energy.

The first three quantities are found by applying equa-
tion 45 on an elemental basis. For those elements within
the mined-out area this step, S,,, ENERGY calculates %
T, u, (the strain energy release from the mined-out mate-
rial for a linear element) and % T, Au (the linear part of
the kinetic energy release for that element). For those
elements within the area Sgy, i.e., the current backfill
and/or gob area, ENERGY evaluates % (1-2) AR Au
(the nonlinear part of the kinetic energy release). Next,
the total energy release for that element is computed.
Normally, this quantity is, either the sum of %T, u, plus
T, Au if the element is within S,,, or just % (1-2) AR Au
if the element is within $,,. Sometimes the total energy
release for that clement is the sum of all three of these
quantities for the special case of an element within both
Sy and Sou.

The last three strain energy quantities are calculated
over each unmined (nongob and/or backfill) element in
the model. These energy quantities, especially the recover-
able strain energy, may provide a very useful indicator of
potentially dangerous highly stressed areas, Moreover, the
procedure for calculating recoverable strain energy pro-
vides a basis for modifying equation 45 to account for
strain energy release from nonlinear materials in the
mined-out area Sy.


Figure 3 shows the necessary program structure modifi-
cations to implement multiple mining steps. After initial
data are read and influence coefficients are computed and
stored, the program enters a FORTRAN "Do loop,” where
it reads the seam geometry and seam materials and pro-
ceeds with stress, displacement, and finally energy com-
putations. Then, it re-reads the seam geometry and seam
materials and proceeds with calculations for the next min-
ing step. This procedure is repeated for all mining steps.

As an aside, the multiple mining step option has con-
siderable computational advantages over running a series
of separate single-step problems. MULSIM uses a simple
Gauss-Seidel iteration procedure to solve the system of
equations for unknown elemental displacements. For the
first iteration of the first mining step, the unknown dis-
placements are customarily set to zero. Once convergence
is achieved and the elemental displacements are known for
the current mining step, then these displacements serve as
the initial guess for the solution at the next step. If the
mining steps are small and the change in mining geometry
is not too radical, then using the displacements of the
prior step as the initial guess ensures rapid convergence to
the solution of the next step.

‘As shown in figure 3, calculation of linear energy
changes follows immediately the stress and displacement

calculations, The energy release for an element mined-out
this step is

a)» (e)

1
Wr = 57-1 (D; -

where Wp = total energy release going from state 1
to state 2,
9;.; = element stress during prior mining step,
element displacement in current mining
step,
and —_-D,_ = element displacement in prior mining

step.

Graphically, this energy quantity is shown in figure 4 (top).
Total energy release is simply the sum of all the elemental
energy releases. ERR is just the total energy release di-
vided by the total elemental area mined this step, ie., the

4

STRESS (c)
9

Strain energy
[released

Enter @
Read and print Read seam geometry
input variables cand material properties
Obtain ride and closure
Define constonts correction factors,
Obtain in situ stresses Solve for unknown rides
‘and gradients and closures
q i

‘Obtain influence
coefficients - stress.

Obtain element stresses.
‘and displacements

rides and closures

‘and displacement and output
Obtain requested off-seam
Store influence
‘coefficients stresses and cuapigeements
Initialize random access
imoys storage ie for Compute linear eneray

6

Repeat for multiple
mining steps

Exit

Figure 3.—Modified program structure of MULSIM with mul-
tiple mining steps and linear energy calculations.

Di-1

y

Dj

DISPLACEMENT (D)

ci-|

yi

Strain energy
/ T released

STRESS (co)

ca

L-- "Dit

Li

—_

DISPLACEMENT (D)

Figure 4.-Stress-displacement curves showing strain energy

softening.

@. Top, From element extraction; bottom, from element


The backfill and/or gob material may exhibit certain
nonlinear stress-displacement behavior. Accordingly, the
parameter a is introduced to account for this nonlinearity
in the strain energy calculations. Figure A-2 (same as
figure 19) illustrates the meaning of this parameter with
linear elastic and strain-hardening materials. Most typ-
ically, backfill and/or gob materials exhibit a strait
hardening or linear elastic behavior; therefore, a has val-
ues in the % to 1 range.

Figure A-3 shows a typical stress-displacement relation-
ship expected for backfill and/or gob. Based on this fig-
ure, the total strain energy in state I of the backfill and/or
gob volume V<; bounded by surface Sq; is

1
Uo = 5 [formes - Ju,Xor (A-10)

In state I, the total strain energy in the backfill and/or
gob volume Vou bounded by surface Sou is

1
Ugn = = Ryuyds - X dy |. (A-11)
ou = 5 ( Tha Svcs Gun "| (11)

Substituting equations A-2 into A-11, expanding the result,
and subtracting A-10 results in the following expression for
the change in strain energy in the gob volume Voy:

Ws = f. aR,Auds ~f XgAudv
Sor Vor

1 1
+3f. aaRAuds - Sav, Xortuev. (A-12)

2JSon

Equations A-5, A-6, A-9, and A-12 provide expressions for
W, Uj Uc and Wg, respectively, necessary to evaluate the
change in energy release given by equation A-1. However,
to simplify the expression that results after substituting
A-5, A-6, A-9, and A-12 into A-1 requires another rela-
tionship derived from the reciprocity theorem. This im-
portant theorem from elasticity theory states that, the
energy associated with the first equilibrium force system
acting through the displacements of the second system is
equal to the energy associated with the second equilibrium
force system acting through the displacements of the first
system. Stated mathematically, the reciprocity theorem is

Fy uy = Fy up, (A-13)

or alternatively,

F, Au = OF y. (A-14)

GOB STRESS (R)

ESS
DISPLACEMENT (u)
O<a<t

Strain energy = V2 a Ru

Figure A-2.—Strain energy calculations for various linear and
nonlinear material behaviors showing how parameter a accounts
for degree of nonlinearity.

o
Sop

oe “DO T
Sz !
49 or |
ze i
z
Se Rp
ge du
oOo

ur un

DISPLACEMENT (u) ———>

Figure A-3.—Typical strain-hardening stress-displacement be-
havior for backfill or gob with basic notation.

With the equilibrium force system given by equation A-3,
the reciprocity theorem gives

fXeuav + foTiAuds
+ Ja,To eas + Jo,,Fau4s
= faXurav + Joo Tras

+ Jao Tards + Jogo Benes. (A-15)


APPENDIX B.—SUPPLEMENTAL DATA ON RIGID BOUNDARY EFFECT TESTS

35 : : 0.10, r T , ,
A 625-m blocks ‘ A
key O9f 125-m elements
30f
— Small, 10 by 20 em
Middle, 20 by 20 cm o8r 4
L + Large, 30 by 20 6m
25 — Analytic solution o7r 4
j o6F 4
20
08
15P ! 4 7 04F 4
!\
1 o3- 4
‘0 625-m blocks | = Sd
125-m elements | 02 4
sr |
! or Analytic soltion J
° i ° f 1 L
985 995 1,005 015,025 1,038. 985 995 1005 1015 1025 1035
35 : : . 040 —
B 25-m blocks 8
30F 35f B-m elements
30 4
25 ; €
2 =
S | 8 2
a 20 t Pa
6 Hy i
2 i i % 20
8 j a
is i i 8
& i i a 15
® H H & y
10 | 25-m blocks ! 3 iol j |
= | “S-melements — j i
5 ! ! |
i i
° rn A i =
950 1,000 4050-4100 _—1180 4,150
35 r
c c
30
‘ J
i J
i
i i |
! !
| 100-m blocks i |
i 20-m elements j
5 i ! 4
\ |
6 i 1 I .
00 1000 —+4,200.~——«,400 1600 200 17000 1,200 1400 ,600

_ DISTANCE, m

____ DISTANCE, m

Figure B-1.Computed and analytic stresses at edge of Figure B-2.Computed and analytic displacements across
‘opening with different coarse-mesh sizes. A, 25m;B,100m;C, opening with different coarse-mesh sizes. A, 25 m; B, 100 m; C,
400 m. 400 m,


where ry = shear stress in local x direction of ele-

ment,

D, = displacement discontinuity in x direc-

tion,
Dy = displacement discontinuity in y direc-
tion,
Ty = shear stress in local y direction of ele-
ment,

¢,, = normal stress in local z direction of
element,

D, = displacement discontinuity in z direc-
tion,

and a, b,c, d, and ¢ are influence functions depending on
the material properties (E,v) of the host medium and the
spatial coordinates (x, y, z). Equation 3 relates an applied
displacement discontinuity at the origin to the induced
stresses at a point (x, y, z) due to that applied displace-
ment discontinuity. While equation 3 applies to a point
displacement discontinuity, it also applies approximately
over an element as shown in figure 10. Here, constant dis-
placement discontinuities D,, D,, and D, over the source
element (k, 1) induce stresses 7,, 7,, and o, at the field
clement (i, j). Equation 4 reflects the change from point
stresses and displacement discontinuities to elemental
stresses and displacement discontinuities.

ij ij i,j
tr, = bD, + cD
‘wl © kl”
ij ij i,j
r= dD, + 6D,
yg kl”
ij ij
o,= aD, @)
kl

Again, equation 4 gives induced stress components at ele-
ment (i, j) due to displacement discontinuities at (k, 1).
Therefore, the total induced stress at element (i, j) is the
sum of the induced stresses from all the displacement
discontinuities at all the elements (k, 1). Taking this
summation on equation 4 results in

i,j ij ij
a> ° D+ »,|

* KI ty kl

6)

If N represents the number of elements, then equation 5
represents a system of 3N equations with 6N unknowns,
3N stress components 7,, 7), and ¢,, and 3N displacement
discontinuities D,, D,, and D,.

Application of Boundary Conditions

Obviously, a unique solution for the unknowns requires
3N additional independent equations which are generated
from the boundary conditions. Figure 11 shows an ele-
mental displacement discontinuity and the nature of the
applied boundary conditions. For unmined clastic ele-
ments (i.¢., those containing coal, vein material, gob, or
backfill), the elemental stresses and displacements are
related through the elemental material properties. For
simple linear elastic materials, the boundary conditions are

7, =GD,/t,
Ty=G D,/t,
2, = ED,/t, ©

where G is a shear modulus, E is Young’s modulus, and
tis the element thickness, It is straightforward to extend
these boundary conditions to more complex nonlinear
material properties.

Source element (k, !)

LL[//L///

Field element (i,j) |

Figure 10.—Sourc
seam plane.

ment (k, 1) and field element (J, j) within