Chapter
2: BACKGROUNDChapter
2: BACKGROUNDFRIC2D 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:
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, Di. The components of displacement discontinuity are defined as the difference in displacement between the two sides of the crack (see Figure 2.8):
(2.1) or
(2.2) where ux and uy are displacement components in x, y coordinate directions, respectively. The superscripts "+" and "-" denote the positive and negative sides of the crack with respect to local coordinates.
The two-dimensional solution for a line of constant displacement discontinuity in an infinite solid was given by Crouch (1976). The displacements, ui, and stresses, sigmaij, in an elastic solid surrounding the line of displacement discontinuity can be written as (see Figure 2.9) (Crouch & Starfield 1983):
(2.3) and
(2.4) 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 2a is:
(2.5) 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 jth element, a local xbar, ybar coordinate system is defined. The local xbar, ybar coordinate system is obtained by a transformation of the global x, y system as follows (see Figure 2.11):
(2.6) For the midpoint (xi, yi) of the ith element, the local coordinate with respect to the jth element's local coordinate system are
(2.7) The displacements and stresses due to displacement discontinuity components Dxbar, Dybar over the jth element can be written by inspection of equations (2.3) and (2.4) after making appropriate coordinate transformation.
(2.8) and
(2.9) The displacement and stress components relative to the local xbar', ybar' coordinate system at the point i can be obtained by appropriate coordinate transformation (Figure 2.12):
(2.10) where gamma= betai - betaj
(2.11) and
(2.12) Now describe the (s,n) coordinate system. After setting Dsj = Dxbarj and Dnj = Dybarj, and noting that uxbari and uybari are equivalent to usi and uni, we can write from equation (2.29)
}i=1,N (2.13) 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 ith element due to a constant unit shear (s) displacement discontinuity over the jth element.
(2.14) Similarly, by substituting traction components tsi and tni for the stress components sigmaxbarybari and sigmaybarybari acting on the element in equations (2.12), we can obtain
}i=1,N (2.15) 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 ith element due to a constant unit shear (s) displacement discontinuity over the jth element.
(2.16) 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)
}i=1,N (2.16) and
}i=1,N (2.17) 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 tsi = tsi(0) and tni = tni(0) are prescribed at the midpoint of the ith boundary element, two equations for the ith element are given as follows:
}i=1,N (2.18) Similarly, if displacements usi = usi(0) and uni = uni(0) are prescribed, the equations for the ith element are given as follows:
}i=1,N(2.19) Mixed formulations in which tsi = tsi(0) and uni = uni(0) or usi = usi(0) and tni = tni(0) are prescribed and handled by selecting the appropriate ones of equations (2.37) and (2.38). Processing in this way for all of the N boundary elements, one can obtain 2N algebraic equations in 2N unknown displacement discontinuities components: Dsj and Dnj (j=1, N), as follows:
}i=1,N (2.20) The quantities bsi and bni stand for the known boundary values of stress or displacement, and Cssij, etc., are the corresponding influence coefficients from equations (2.18) and (2.19).
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 2N 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 kratio. This parameter describes the ratio of hoizontal stress to vertical stress.
(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, Dsj and Dnj (j=1, N), one can use equations (2.8) and (2.9) (where Dxbar and Dybar are denoted as Dsj and Dnj, uxbar and uybar as uxbar(j) and uybar(j), and sigmaxbarxbar, sigmaybarybar and sigmaxbarybar as sigmaxbarxbar(j), sigmaybarybar(j) and sigmaxbarybar(j) ) to calculate the displacements and stresses at designated points due to the displacement discontinuity components, Dsj and Dnj, over the jth element (j=1, N). The displacement and stress components in the x, y coordinate system can be obtained by appropriate coordinate transformation:
(2.22) and
(2.23) The displacements and stresses at the designated points are functions of all elements' effects on these points, that is
(2.24) and
(2.25) 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, µ .
(2.26) 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.
(2.27) 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).
(2.28) Thus we restrict the shear stress as follows:
(2.29) and
(2.30) 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.
(2.31) 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 (sigmas1<<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, KI, exceeds the fracture toughness, KIc. 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, theta0, in which the circumferential stress, sigmatheta-theta, is a maximum.
The maximum circumferential stress sigmatheta-theta (theta0) at a fracture tip can be expressed in terms of the modes I and II stress intensity factors:
(2.32) A new fracture will develop when the following is met:
(2.33) For a uniformly loaded fracture, Olson (1991) compared the normal and shear components of the displacement discontinuity, Dn and Ds, for the fracture tip element to the analytical solution for fracture-wall opening and shearing displacements (Pollard & Segall 1987), and derived the following expressions for KI and KII
(2.34) 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 KI and KII 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)
(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.)
(2.36) 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).
(2.37) 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:
(2.38) 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:
(2.39) 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.
(2.40) 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) .
(2.41) 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).
(2.42) 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).