Integrating complementarity into the 2D displacement discontinuity boundary element method to model faults and fractures with frictional contact properties

doi:10.1016/j.cageo.2011.11.017

Computers & Geosciences

Volume 45, August 2012, Pages 304-312

https://doi.org/10.1016/j.cageo.2011.11.017 Get rights and content

Abstract

We present a two-dimensional displacement discontinuity method (DDM) in combination with a complementarity solver to simulate quasi-static slip on cracks as models for faults and fractures in an otherwise homogeneous, isotropic, linear elastic material. A complementarity algorithm enforces appropriate contact boundary conditions along the cracks so that variable friction and frictional strength can be included. This method accurately computes slip and opening distributions along the cracks, displacement and stress fields within the surrounding material, and stress intensity factors at the crack tips. The DDM with complementarity is a simple yet powerful tool to investigate many aspects of the mechanical behavior of faults and fractures in Earth's brittle crust. Implementation in Excel and Matlab enables easy saving, organization, and sharing.

Introduction

Faults and fractures are commonly modeled as cracks in an otherwise homogeneous and isotropic linear elastic material. This approximation allows us to idealize the problem of faulting and fracturing in Earth's brittle crust. Although many analytical solutions for crack problems exist, such solutions are not viable for problems where the cracks have irregular geometry or the boundary conditions cannot be prescribed. Numerical tools allow us to solve new and more realistic crack problems for geoscience applications. We introduce a numerical method that merges two existing computational tools: the displacement discontinuity boundary element method (DDM) and a complementarity algorithm.

The elementary implementation of the two-dimensional DDM (Crouch and Starfield, 1983, Chapter 5) is only accurate for problems for which the contact conditions are known a priori. Incorporating a complementarity algorithm into the DDM is an effective way to numerically model quasi-static slip on cracks and the resulting deformation of the nearby material (De Bremaecker et al., 2004, Griffith et al., 2010, Mutlu and Pollard, 2008, Ritz and Pollard, 2011). Complementarity is a popular optimization method used in engineering, economics, and various scientific disciplines (e.g., Cottle et al., 1992, Harker and Pang, 1990, Murty, 1988, Pang, 1994, Pang and Trinkle, 1996). The aim of this article is to communicate this simple yet powerful tool to the geoscience community, as we and several other researchers have recently used this approach to investigate the mechanical behavior of faults and fractures in Earth's brittle crust (De Bremaecker et al., 2004, Griffith et al., 2010, Mutlu and Pollard, 2008, Ritz and Pollard, 2011). Here, the DDM with complementarity formulation is described, and examples are provided to verify the validity of this method as well as demonstrate its utility.

Section snippets

DDM formulation

DDM is a special type of boundary element method (BEM) (Crouch, 1976, Crouch and Starfield, 1983). The DDM can treat cracks embedded in an infinite homogeneous, isotropic, elastic material as curvilinear boundaries discretized into linear elements of constant displacement discontinuity (Fig. 1). A major advantage of using a BEM is that only discretization of the crack is required, whereas other popular methods, such as the finite element method, require discretization of the surrounding

Frictional contact properties

The shear and normal stress components at the center of a boundary element are related to the specified traction boundary conditions using Cauchy's Formula (Pollard and Fletcher, 2007, p. 213), and also are related to each other by the Coulomb criterion: $| σ_{n s} | \leq - μ σ_{n n} + S_{f}, σ_{n n} \leq 0 MPa .$ here μ is the coefficient of friction. Laboratory values for the coefficient of friction of rocks are generally reported to be 0≤μ≤1 (Byerlee, 1967, Hoskins et al., 1968, Jaeger, 1959), and commonly are found to be

Complementarity formulation

The DDM discussed in Section 2 does not include frictional contact properties because interpenetration is allowed, preventing contact between the two crack surfaces. By introducing a complementarity algorithm, we can correct the initial DDM solution and include the frictional contact properties discussed in Section 3.

Since the behavior of each boundary element is not known a priori, the complementarity solver utilizes the basic relationships between sticking, slipping, and opening to correct

Implementation in Matlab

We implement the DDM with complementarity in Matlab, which imports an Excel spreadsheet. All input parameters are managed in a spreadsheet so that problems can be more easily organized, saved, and shared. Alternatively, the input parameters can be prepared using other text formats, readable by Matlab. The user-defined parameters are the elastic moduli, remote Cartesian stress components, the number of fractures, the coordinates for each endpoint of every boundary element, the coefficient of

Validation

In linear elastic fracture mechanics (LEFM) theory, if the gradient of displacement at the model crack tip is infinite, the stresses are singular there. However, this singularity is physically unrealistic because natural materials cannot support infinite stress. One explanation is that real materials limit the stress concentrations at crack tips by locally deforming inelastically (Dugdale, 1960). If the size of the zone of inelastic deformation at a crack tip is negligible in comparison with

Stick or slip conditions

To illustrate the utility of the DDM with complementarity, we explore the deformation of an isolated, planar crack with various frictional contact properties. Throughout this discussion, the elastic moduli are E=50,000 MPa, and ν=0.25. The crack is divided into 500 boundary elements of equal length, lies along the x-axis, and is centered on the y-axis (Fig. 1). A remote principal stress ratio of $σ_{3}^{\infty} / σ_{1}^{\infty} = 3.5$ drives slip on the crack in all cases, with the maximum principal stress $σ_{1}^{\infty} = - 50 MPa$ and

Discussion

Faults are rarely isolated, continuous, or planar. While CEZs may warp the planar geometry of a fracture (Martel, 1997) the DDM with complementarity also can model segmented, branching, and other types of non-planar fractures. Fig. 8 schematically illustrates various non-planar fractures that have been observed in nature; segmented and stepping, ramped, intersecting and branching, splayed, and curved.

Many numerical studies addressing fracture shapes shown in Fig. 8 have developed from an

Conclusion

The DDM with complementarity accurately computes slip, opening, and stress distributions along cracks, and displacement and stress within the surrounding material, while incorporating the contact properties μ and S_f. This technique can accommodate any number of cracks, each with an arbitrary distribution of friction and frictional strength along its length. Additionally, this numerical tool can model natural fracture and fault geometries if a sufficient number of boundary elements are used to

Acknowledgments

We thank Michael Ferris for permission to use PATH as the exemplar complementarity code. We also thank two anonymous reviewers for comments that greatly improved the manuscript. Funding was provided by DOE Basic Energy Sciences, grant DE-FG02-04ER15588-5.

References (59)

G.I. Barenblatt
The mathematical theory of equilibrium cracks in brittle fracture
Advances in Applied Mechanics
(1962)
R. Bürgmann et al.
Slip distributions on faults: effects of stress gradients, inelastic deformation, heterogeneous host-rock stiffness, and fault interaction
Journal of Structural Geology
(1994)
P.A. Cowie et al.
Physical explanation for the displacement-length relationship of faults using a post-yield fracture mechanics model
Journal of Structural Geology
(1992)
J. DeBremaecker et al.
Including gravity in the displacement discontinuity method when accounting for contact interaction
International Journal of Rock Mechanics and Mining Sciences
(2004)
J.C. DeBremaecker et al.
A comparison of two algorithms for solving closed crack problems
Engineering Fracture Mechanics
(2000)
D.S. Dugdale
Yielding of steel sheets containing slits
Journal of the Mechanics and Physics of Solids
(1960)
N. Elvin et al.
A fast iterative boundary element method for solving closed crack problems
Engineering Fracture Mechanics
(1999)
E.R. Hoskins et al.
A medium-scale direct friction experiment
International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts
(1968)
L. Maerten et al.
Slip distributions on intersecting normal faults
Journal of Structural Geology
(1999)
S.J. Martel
Effects of cohesive zones on small faults and implications for secondary fracturing and fault trace geometry
Journal of Structural Geology
(1997)

J.E. Olson et al.

The initiation and growth of en-echelon veins

Journal of Structural Geology

(1991)

D.D. Pollard et al.

Theoretical displacements and stresses near fractures in rock

P.G. Resor et al.

Slip heterogeneity on a corrugated fault

Earth and Planetary Science Letters

(2009)

A.L. Thomas et al.

The geometry of echelon fractures in rock—implications from laboratory and numerical experiments

Journal of Structural Geology

(1993)

J. Wang et al.

An iterative algorithm for modeling crack closure and sliding

Engineering Fracture Mechanics

(2008)

R. Bilham et al.

The morphology of strike-slip faults: examples from the San Andreas Fault, California

Journal of Geophysical Research

(1989)

J. Byerlee

Friction of rocks

Pure and Applied Geophysics

(1978)

J.D. Byerlee

Frictional characteristics of granite under high confining pressure

Journal of Geophysical Research

(1967)

M. Cooke

Frictional Slip and Fractures Associated with Faults and Folds. PhD dissertation

(1996)

M.L. Cooke

Fracture localization along faults with spatially varying friction

Journal of Geophysical Research

(1997)

R. Cottle et al.

The Linear Complementarity Problem

(1992)

C.A. Coulomb

Sur une application des règles de maximis & minimus a quelques problèmes de statique, relatifs a l'architecture

Académie Royale des Sciences Mémoires de Mathématique & de Physique

(1773)

S.L. Crouch

Solution of plane elasticity problems by the displacement discontinuity method. I. Infinite body solution

International Journal for Numerical Methods in Engineering

(1976)

S.L. Crouch et al.

Boundary Element Methods in Solid Mechanics

(1983)

Dirkse, S.P., Ferris, M., Munson, T.S., 〈http://pages.cs.wisc.edu/∼ferris/path/〉. The PATH Solver. University of...

S.P. Dirkse et al.

The PATH solver: a non-monotone stabilization scheme for mixed complementarity problems

Optimization Methods and Software

(1995)

Erdogan, F., 1962. On the stress distribution in plates with collinear cuts under arbitrary loads, In: Proceedings of...

M.C. Ferris et al.

Feasible descent algorithms for mixed complementarity problems

Mathematical Programming

(1999)

Ferris, M.C., Munson, T.S., 1998. Complementarity Problems in GAMS and the PATH Solver, Mathematical Programming...

Cited by (23)

Investigating stress shadowing effects and fracture propagation patterns: Implications for enhanced geothermal reservoirs
2021, International Journal of Rock Mechanics and Mining Sciences
During hydraulic stimulation treatment in an enhanced geothermal (EGS) reservoir, it has been suggested that new fractures are initiated from stimulated preexisting fractures and propagate through the reservoir. In this stimulation mechanism, a fracture propagation from a preexisting natural fracture and the interaction of newly formed fractures and preexisting natural fractures play an important role in the creation of a complex fracture network. In this study, we developed a physics-based numerical model that couples fluid flow between fracture surfaces, fracture deformation, and fracture propagation driven by fluid injection to better understand how a fracture network is created by a hydraulic stimulation treatment. The results showed that the created fracture network complexity is affected by the fracture intersection angle, stress state, and injection rates. Having flow path complexity is advantageous in EGS because it increases the heat exchange areas, the fracture surface to rock volume ratio, and the general reservoir permeability. Therefore, the implication from the results of this study to make a better EGS reservoir is to stimulate an EGS reservoir with 1) well-oriented fractures, 2) a high stress ratio, and 3) a low injection rate. These reservoir conditions and stimulation design will likely make more flow paths and create a more complex fracture network.
Contact between rough rock surfaces using a dual mortar method
2020, International Journal of Rock Mechanics and Mining Sciences
Citation Excerpt :
Simulation tools that resolve small-scale roughness, while allowing computation of contact stresses and stress variations in the rock mass, are rare. Most analytical or numerical studies at the laboratory or reservoir scale rely on simplifying assumptions and regularizations, such as averaging asperity scale processes with an empirical relationship, reduction to two dimensions, using variants of boundary element methods, not modeling explicitly the mechanical contact between two surfaces, treating the surfaces as two parallel plates, using penalty methods or augmented Lagrangian methods, or only matching triangulations at the surface.4,14–23 In this work, we present a new parallel implementation of an approach that does not require any of the previously mentioned assumptions.
The mechanical behavior of fractures in rocks has strong implications for reservoir engineering applications. Deformations, and the corresponding change in contact area and aperture field, impact rock fracture stiffness and permeability, thus altering the reservoir properties significantly. Simulating contact between fractures is numerically difficult as the non-penetration constraints lead to a nonlinear problem and the surface meshes of the solid bodies on the opposing fracture sides may be non-matching. Furthermore, due to the complex geometry, the non-penetration constraints must be updated throughout the solution procedure. Here we present a novel implementation of a dual mortar method for contact. It uses a non-smooth sequential quadratic programming method as solver, and is suitable for parallel computing. We apply it to a two body contact problem consisting of realistic rock fracture geometries from the Grimsel underground laboratory in Switzerland. The contributions of this article are: (1) a novel, parallel implementation of a dual mortar method with a non-smooth sequential quadratic programming method, (2) realistic rock geometries with rough surfaces, and (3) numerical examples, which prove that the dual mortar method is capable of replicating the nonlinear closure behavior of fractures observed in laboratory experiments.
3D scene and geological modeling using integrated multi-source spatial data: Methodology, challenges, and suggestions
2020, Tunnelling and Underground Space Technology
Citation Excerpt :
The model can be used for groundwater seepage and grouting diffusion analysis, as shown in Fig. 17(a). Compared with a traditional plan and section, the 3D geological model can quickly and clearly display the three-dimensional shape and distribution of a geological body (Tertois and Mallet, 2007; Caumon et al., 2009; Ritz et al., 2012). Users can directly see the crosscutting relationship between the fault and the stratum and the coal bed, change the observation angle at will, close the secondary horizon, and highlight the expression of the spatial relationship of the three-dimensional model through various means, such as interface coloring to reflect the elevation changes (Caumon et al., 2009).
The visualization and analysis of massive, multi-source, multi-dimensional, heterogeneous, dynamic, and complex spatial data is playing an increasingly important role in 3D scene and geological modeling. There are great differences in various research fields regarding data acquisition techniques, morphological descriptions, and application purposes of spatial objects. Therefore, based on an existing three-dimensional visualization technology and analysis program of spatial data, an integrated processing method for the effective fusion and analysis of multi-source complex spatial data is proposed. This method combines a massive point cloud model with a digital elevation model, in addition to combining the complex stratum with fault modeling technology. The structural grid and geological grid are effectively combined to realize the unified modeling and visualization of natural and artificial objects. This method was then applied to a 3D modeling scene of deeply buried mines and a geological construction model for open pit mines. The comprehensive and effective utilization of multi-source geological data is realized by using this modeling method. Finally, the challenges of spatial data visualization from the aspects of the difficulties in spatial data acquisition, the complexities of the spatial structure, the constraint relations of geological bodies, and the complexity of the 3D modeling process are discussed in this study. Three specific suggestions are proposed. This study enhances the application efficiency of multi-source spatial data modeling technologies already in practice and provides reference for future applications and developments.
Slip on wavy frictional faults: Is the 3rd dimension a sticking point?
2019, Journal of Structural Geology
Citation Excerpt :
This is deemed sufficiently accurate in capturing the boundary condition set. In these equations, h is the distance from the mid-point (geometric incenter) of the boundary element to the fracture's tip (Ritz et al., 2012), Dn is the normal displacement of this element, DII is the displacement perpendicular to the crack edge, and DIII is the displacement parallel with the edge. The correction factor c is used to correct for the errors due to the numerical approximation.
The formulation for the 3D triangular displacement discontinuity boundary element method with frictional constraints is described in detail. Its accuracy in comparison to analytical solutions is then quantified. We show how this can be used to approximate stress intensity factors at the crack tips. Using this method, we go on to quantify how slip is reduced on fault surfaces with topography, where the asperities are approximated as a sinusoidal waveform, i.e. corrugations. We use stress boundary conditions (compressive) orientated such that frictional contacts shear. We show that slip reductions relative to planar faults for 2D line and 3D penny-shaped crack models are comparable within 10% when slip is perpendicular to the corrugations. Using the 3D model, we then show how slip is reduced more when corrugation wavelengths are doubled compared to the reduction due to corrugation misalignment with the slip direction. When slip is parallel with the corrugation alignment we show that reducing the out-of-plane stress, from the normal traction acting on the fault when planar to that resolved on a perpendicular plane, has the same effect as halving the length of the corrugation waveform in terms of slip reduction for a given amplitude.
The distinct element method (DEM) and the extended finite element method (XFEM) application for analysis of interaction between hydraulic and natural fractures
2018, Journal of Petroleum Science and Engineering
Hydraulic fracturing can be considered one of the most effective stimulation methods for fractured network formation. The interaction between induced and natural fractures greatly affects the success of the treatment. In this paper, the extended finite element method (XFEM) and distinct element method (DEM) have been combined to identify the propagation of hydraulic fractures in a porous medium with naturally fractured blocks. DEM is used to investigate hydraulic fractures in a natural-fracture network and evaluate the production flow rate under different apertures and lengths of hydraulic fracture. XFEM is used to simulate the deformation of the natural-fracture interface during the approaching phase of induced fractures as a no-remeshing tool. The results of DEM reveal that length of hydraulic fracturing has a significantly greater effect on production than aperture. Results from XFEM show that tensile and shear debonding of natural fractures change as a function of angle and distance from an induced fracture.
Mechanical evolution of transpression zones affected by fault interactions: Insights from 3D elasto-plastic finite element models
2018, Journal of Structural Geology
Citation Excerpt :
In this sense, fault steps are zones of slip transfer between discontinuous sub-parallel fault segments in which segments interact through their associated stress field and any hard-linkages. Various non-planar fractures and fault geometries observed in nature and studied both theoretically and numerically (e.g., Rodgers, 1980; Segall and Pollard, 1980; Aydin and Schultz, 1990) are segmented/stepping, ramped, intersecting/branching, splayed, and curved (Ritz et al., 2012, 2015; Ritz, 2013). Fault steps localize either contraction or extension as a function of their manual geometries (right- or left-stepping) and fault kinematics (left- or right-lateral) (Christie-Blick and Biddle, 1985; Woodcock and Fischer, 1986; Ramsay and Huber, 1987; Reches, 1987; Crider, 2001; Storti et al., 2003; Cunningham and Mann, 2007; Mann, 2007).
The mechanical evolution of transpression zones affected by fault interactions is investigated by a 3D elasto-plastic mechanical model solved with the finite-element method. Ductile transpression between non-rigid walls implies an upward and lateral extrusion. The model results demonstrate that a, transpression zone evolves in a 3D strain field along non-coaxial strain paths. Distributed plastic strain, slip transfer, and maximum plastic strain occur within the transpression zone. Outside the transpression zone, fault slip is reduced because deformation is accommodated by distributed plastic shear. With progressive deformation, the σ₃ axis (the minimum compressive stress) rotates within the transpression zone to form an oblique angle to the regional transport direction (∼9°–10°). The magnitude of displacement increases faster within the transpression zone than outside it. Rotation of the displacement vectors of oblique convergence with time suggests that transpression zone evolves toward an overall non-plane strain deformation. The slip decreases along fault segments and with increasing depth. This can be attributed to the accommodation of bulk shortening over adjacent fault segments. The model result shows an almost symmetrical domal uplift due to off-fault deformation, generating a doubly plunging fold and a ‘positive flower’ structure. Outside the overlap zone, expanding asymmetric basins subside to ‘negative flower’ structures on both sides of the transpression zone and are called ‘transpressional basins’. Deflection at fault segments causes the fault dip fall to less than 90° (∼86–89°) near the surface (∼1.5 km). This results in a pure-shear-dominated, triclinic, and discontinuous heterogeneous flow of the transpression zone.

View all citing articles on Scopus

View full text

Integrating complementarity into the 2D displacement discontinuity boundary element method to model faults and fractures with frictional contact properties

Abstract

Introduction

Section snippets

DDM formulation

Frictional contact properties

Complementarity formulation

Implementation in Matlab

Validation

Stick or slip conditions

Discussion

Conclusion

Acknowledgments

Advances in Applied Mechanics

Journal of Structural Geology

Journal of Structural Geology

International Journal of Rock Mechanics and Mining Sciences

Engineering Fracture Mechanics

Journal of the Mechanics and Physics of Solids

Engineering Fracture Mechanics

International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts

Journal of Structural Geology

Journal of Structural Geology

Journal of Structural Geology

Earth and Planetary Science Letters

Journal of Structural Geology

Engineering Fracture Mechanics

The morphology of strike-slip faults: examples from the San Andreas Fault, California

Journal of Geophysical Research

Friction of rocks

Pure and Applied Geophysics

Frictional characteristics of granite under high confining pressure

Journal of Geophysical Research

Frictional Slip and Fractures Associated with Faults and Folds. PhD dissertation

Fracture localization along faults with spatially varying friction

Journal of Geophysical Research

The Linear Complementarity Problem

Sur une application des règles de maximis & minimus a quelques problèmes de statique, relatifs a l'architecture

Académie Royale des Sciences Mémoires de Mathématique & de Physique

Solution of plane elasticity problems by the displacement discontinuity method. I. Infinite body solution

International Journal for Numerical Methods in Engineering

Boundary Element Methods in Solid Mechanics

The PATH solver: a non-monotone stabilization scheme for mixed complementarity problems

Optimization Methods and Software

Feasible descent algorithms for mixed complementarity problems

Mathematical Programming