FRIC2D is a computer program used for analyzing rock fracture propagation in two dimensions, based on the boundary element method for solving the elastic quasi-static boundary value problems. It is important for the user to understand the basic principles of the boundary element method.

The calculation of stresses and displacements in FRIC2D is based on the two-dimensional analytical solution for the problem of a constant discontinuity in displacement over a finite line segment in the x, y plane of an infinite elastic body (Crouch & Starfield 1983). This solution is, in turn, based on the solution for a point source of displacement discontinuity (Crouch 1979). One can divide the external boundary of the body as well as any internal cavities and fractures into different series of elements (boundary elements) and take the displacement discontinuity to be a different constant over each element. Then one can find a numerical solution to the problem by determining the displacement discontinuities for all elements which match the specified traction or displacement boundary conditions.

The user documentation for FRIC2D contains some good example applications of this family of Boundary Element Codes to the fields of faulting, hydraulic fracture and underground storage. An introduction to the mechanics of fracture propagation is also presented in the FRIC2D documentation. We recommend that new users read that document to learn more about possible applications of this code and to review the basics of fracture propagation theory.

Chapter Contents:

- 2.1 Constant Displacement Discontinuity in an Infinite Solid
- 2.2 Discretization of Boundary
- 2.3 Boundary Values
- 2.4 Gravitational Loading
- 2.5 Calculation of Stresses and Displacements
- 2.6 Fault Elements
- 2.7 Joint Initiation

A crack has two surfaces or boundaries, one effectively coinciding with
the other. In crossing from one side of the crack to the other, the
displacements undergo a specific change in value referred to as the
displacement discontinuity, `D _{i}`.
The components of displacement
discontinuity are defined as the difference in displacement between the
two sides of the crack (see Figure 2.8):

or

where `u _{x}` and

The two-dimensional solution for a line of constant displacement
discontinuity in an infinite solid was given by Crouch (1976). The
displacements, `u _{i}`, and stresses,

and

where `nu` is Poisson's ratio and `G` is the
elastic shear modulus. The function `f(x,y)` which solves this
problem for the line of displacement discontinuity of length
2`a` is:

The quantities in which commas appear as subscripts in equations
(2.3) and (2.4) are derivatives of the function `f`(x,y).
For example, f,_{xy} stands for
.

The analytical solution given in section 2.1 forms the basis of a
boundary element method for solving boundary value problems. The
numerical approximation to the solution of these problems can be
represented by `N` discrete straight line elements
(see Figure 2.10).

For the `j`^{th} element, a local
`x`bar, `y`bar
coordinate system is defined. The local
`x`bar, `y`bar
coordinate system is obtained by a transformation of the global
`x, y` system as follows (see Figure 2.11):

For the midpoint (`x _{i}`,

The displacements and stresses due to displacement discontinuity
components `Dx`bar,
`Dy`bar
over the `j`^{th} element can be written by inspection of
equations (2.3) and (2.4) after making appropriate coordinate
transformation.

and

The displacement and stress components relative to the local
`x`_{bar}',
`y`_{bar}'
coordinate system at the point `i` can be obtained by appropriate
coordinate transformation (Figure 2.12):

where `gamma`=
`beta _{i}` -

and

Now describe the (`s,n`) coordinate system. After setting
`Dsj` =
`Dx`bar`j` and
`Dnj` =
`Dy`bar`j`,
and noting that
`ux`bar`i` and
`uy`bar`i`
are equivalent to
`usi` and
`uni`,
we can write from equation (2.29)

where `Bssij`,
`Bsnij`,
`Bnsij`,
`Bnnij`
are the boundary influence coefficients for the displacements.

The coefficient, `Bnsij`,
for example, gives the normal (`n`) displacement at the midpoint of the
`i`^{th} element due to a constant unit shear (`s`)
displacement discontinuity over the `j`^{th} element.

Similarly, by substituting traction components
`tsi` and
`tni` for the stress
components
`sigmax`bar`y`bar`i`
and
`sigmay`bar`y`bar`i`
acting on the element in equations (2.12), we can obtain

where
`Assij`,
`Asnij`,
`Ansij`,
`Annij`
are the boundary influence coefficients for the tractions. The coefficient,
`Ansij`,
for example, gives the normal (`n`) traction
at the midpoint of the `i`^{th} element due to a
constant unit shear (`s`) displacement discontinuity over the
`j`^{th} element.

In general, the displacements and tractions at the ith element are
functions of the displacement discontinuity components at all `N`
elements, i.e. `j`=1 to `N`. We can write from equations
(2.13) and (2.15)

and

For the mechanical analysis of deformation in an elastic solid, the
boundary values can be divided into two types: traction boundary
conditions and displacement boundary conditions. If tractions
`t _{s}^{i}` =

Similarly, if displacements
`u _{s}^{i}` =

Mixed formulations in which
`t _{s}^{i}` =

The quantities
`b _{s}^{i}` and

The equations (2.20) can be solved by using either Gaussian
elimination or Gauss-Siedel methods. Detailed descriptions of these
methods for solving linear algebraic equations are available in many
textbooks (Fox 1965, Desai 1979).
FRIC2D adopts the Gaussian elimination method to solve the
2`N` linear equations.

When considering frictional slip along surfaces within the earth it is important to evaluate the tractions on the surface due to pressure of overlying rock. The overburden loading is modelled by superposing gravitational compression which increases linearly with depth. The total stress state at any point is the stress computed from the applied boudnary conditions and the gravitational stresses.

The vertical gravity stress is deterimined from the weight of the overburden (2.21) while the horizontal stress is a portion of the vertical stress prescribed by

(2.21)

where r is the homogeneous rock density, `g`
is acceleration of gravity, and `d` is the depth, the subscripts `x`
and `y` here refer to the global coordinate system.

For lithostatic loading* kratio* = 1

For bilateral constraint *kratio* = n/(1-n)
where v is Poisson's ratio.

After obtaining all displacement discontinuity components,
`D _{s}^{j}` and

and

The displacements and stresses at the designated points are functions of all elements' effects on these points, that is

and

Special boundary elements, termed 'fault' elements, that can accommodate
inelastic slip have been used to model frictional interfaces (Crouch 1979,
Crouch & Starfield 1990) including geologic structures such as faults.
Unlike standard boundary elements, which require prescribed loading conditions,
these frictional elements require prescription of constitutive behavior.
The elemental deform elastically and may also experience friction slip.
The elemental mechanical behavior is described by 4 parameters: shear
stiffness, `Ks`, normal stiffness, `Kn`, cohesion,
`c`, and coefficient of friction, µ .

That, is the shear and normal displacement discontinuites are proportional to the shear and compression acting on the elements.

If tensile stresses develop across the element, there is no longer contact across the the surfaces of the element; shear and normal tractions become zero.

The normal displacement discontinuties,
`Dni`,
always follow the relationships described above. However, the shear
displacement discontinuties, `Dsi`,
may be a combination of elastic shear and inelastic frictional slip.

We model each interface as a set of frictional elements and use the Coulomb
criterion to determine if interface elements slip. This criterion relates the
shear stress required for slip to the frictional resistance, `R`
(Jaeger & Cook 1979). Slip occurs under compressive (negative) normal
stress, `sigman`, if the magnitude of shear stress,
`sigmas`, exceeds the sum of the
cohesion, `c`, and the coefficient of friction, µ ,
times the normal stress (2.28).

Thus we restrict the shear stress as follows:

and

If the normal stresses are tensile, there is no longer contact across the interface; shear stresses are zero.

For compressive normal stresses, if |`sigmas`| >
`R`, the element will slip until the shear stress drops below the
frictional resistance of the element. At this stage, the stress state
for the entire system must be recalculated because slip on one element
may change the normal and shear stresses on neighboring parts of the
interface. The shear stress on elements adjacent to the slipped element
may then exceed the frictional resistance and slip as well. For example,
Figure 2.6 schematically shows how shear stresses on a bedding interface may
respond to slip along a particular element `n` through successive iterations.
After several iterations, shear stresses on the fault elements approach the
frictional resistance, and displacement discontinuities begin to converge
to the solution for that loading step.

Figure 2.6

The solution converges when the relative difference between shear stresses,
`sigmas`, of successive iterations, `k`-1
and `k` is less than the prescribed tolerance.

Equation 2.31 is a version of the relative difference formulation
(Acton 1990, p. 49) which is capable of quantifying the shear
stress differences for values near zero by adding a factor of 1 to each
term. For large shear stresses (`sigmas` >>1)
the factor of one in the denominator of equation 2 is negligible and the formula
calculates the percentage difference of shear stresses between iterations.
For small shear stresses (`sigmas`1<<1)
equation 2.31 calculates the absolute difference (Acton 1990, p. 49).
If the user prescribed a tolerance of 0.1 then the solution is considered
to be converged when the relative difference between
successive iterations (2.31) is less than 0.1. Lower tolerances require
more iterations to reach convergence. Tolerances from 0.01 to 0.001
are reccommended.

The initiation and propagation of a fracture is determined by two parameters: the stress intensity factor (a function of geometry determined analytically or numerically by solving the appropriate boundary value problems), and the fracture toughness (a material property determined experimentally). The relationship between stress intensity factor and fracture toughness is analogous to the relationship between stress and strength.

For mode I loading, the fracture will propagate along a straight path
when the stress intensity factor, `K`_{I},
exceeds the fracture toughness,
`K`_{I}^{c}. i.e.

A criterion for fracture initiation under mixed mode I/II loading was
proposed by Erdogan and Sih (1963). It is known as the maximum
circumferential stress criterion. It states that fracture initiation takes
place in a direction, `theta`0,
in which the circumferential stress,
`sigma _{theta-theta}`,
is a maximum.

The maximum circumferential stress
`sigma _{theta-theta}`
(

A new fracture will develop when the following is met:

For a uniformly loaded fracture, Olson (1991) compared the normal
and shear components of the displacement discontinuity,
`D _{n}` and

where `E` is Young's modulus of rock and `l` is
the length of the fracture tip element.

Equations (2.34) were also tested against the solution technique of Pollard & Holzhausen (1979) for the case of closely spaced, straight, echelon fractures. The margin of error was found to be less than 5% even when each fracture was divided into only two elements (Olson 1991).

If the stress intensity factors `K`_{I}
and `K`_{II} at each tip satisfy the
criterion given in equation (2.33), then a new fracture increment is
added to the end of the fracture.
The fracture will grow in the direction of the maximum
circumferential stress (2.35)

Along a frictional interface the normal and shear stresses and displacements
are calculated for the BEM solution. The normal stress tangential to the
interface is not calculated but can be found from the normal stress and
tangential displacement as described in Crouch and Starfield
(1990, pp. 100-109).
Hooke's Law expresses the normal strain, `epsilonxx`,
in terms of normal stresses `sigmaxx`,
`sigmayy` and `sigmazz`
(Timoshenko & Goodier 1934). (For the rest ofthis section `x`
and `y` refer to the coordinates local to the fault element.)

We assume that plane strain conditions dominate at depth within the earth so
that we only need to consider the stress state in the two-dimensional
`x-y` plane. The normal stress, `sigmazz`, is
related to `sigmaxx` and
`sigmayy` according to plane strain conditions
(Timoshenko & Goodier 1934).

Substituting equation 2.37 in 2.36 and rearranging 2.36 allow us to express
the normal stress, `sigmaxx`,
in terms of normal stress, `sigmayy`, and normal strain,
`epsilonxx`, as follows:

In the case of a slipping interface, the normal stresses acting tangential to the interface are discontinuous across the interface and may have different values above and below the interface. We use a superscript to distinguish the upper (positive) and lower (negative) sides of the frictional interface. The normal tangential stresses are:

The normal strain, `epsilonxx`, equals the
change in slip on the interface, (Timoshenko & Goodier 1934) .
Where the tangential displacement is constant along the
interface, the `x`-parallel strain is zero. Where the shear displacement,
`ux`, changes rapidly the `x`-parallel
strain, `epsilonxx`, will be concentrated.

Using the backward difference method (Crouch & Starfield 1990) , the displacement
gradient is determined at points between successive elements along the frictional
interface. The interface normal, `sigmayy`, and shear,
`sigmayx`, stresses are determined
between neighboring elements by taking the average of the two elements. From the
interface normal stress and the displacement gradients the tangential stresses
above and below the interface can be determined between successive elements along
the frictional interface.

The maximum principal stress is calculated along the frictional interface at points between successive elements from the normal, shear and tangential stresses. The maximum principal stress is determined above and below the interfaces from the respective normal and shear stresses (Timoshenko & Goodier 1934) .

If the maximum tensile stress exceeds the tensile strength of the rock a new
fracture will grow perpendicular to the direction of maximum tension.
The new splay crack is oriented at an angle a from the frictional interface that
is `pi`/2 from the orientation of the maximum tensile stress
(Timoshenko & Goodier 1934).

If the maximum tensile stress exceeds the tensile strength of the rock at any point along the frictional interface, a fracture is added.

The initial angle at which the new fracture grows is determined from the orientation of the maximum tensile stress at the initiation point according to the relationship: opening-mode fractures grow perpendicular to the direction of maximum tension (Jaeger & Cook 1979, Lawn 1993). The new fracture orientation is calculated and a new fracture element added to the problem. The new fracture length equals that of the frictional element from which it initiated and the fracture is considered open (no shear or normal stresses on its surfaces). The boundary value problem must be recalculated with the new fracture in place before incrementing the monotonic boundary conditions (tectonic strains).