A stress and displacement discontinuity element method for elastic transversely anisotropic rock

The paper presents closed‐form solutions for stress and displacement influence functions for stress discontinuity (SD) and displacement discontinuity (DD) elements, for a two‐dimensional plane‐strain elastic, transversely anisotropic medium. The solutions for SD elements are based on Kelvin's problem and for DD elements on the concept of dipoles. Stress and displacement influence functions are derived for the following elements: constant SD, linear SD, constant DD, linear DD, square root DD, parabolic DD, constant DD surface, and linear DD surface elements. The formulations are incorporated into FROCK, a hybridized boundary element method code, and are validated by providing comparisons between the results from FROCK and the finite element code ABAQUS. A limited parametric analysis shows the effects of slight anisotropy on the stress field around the tip of a crack and of the orientation of the crack with respect to the axes of elastic symmetry. Copyright © 2014 John Wiley & Sons, Ltd.     

Analyses of Inclined Cracks Neighboring Two Iso-Path Cracks in Rock-Like Specimens Under Compression

Crack propagation process in pre-cracked rock like specimens has been studied experimentally and numerically considering three cracks in the middle part of each specimen. The rock-like specimens are specially prepared from Portland pozzolana cement, fine sands and water. These pre-cracked cylindrical specimens (each containing a single inclined crack in the neighborhood of two iso-path cracks) are experimentally tested under compressive loading. The same problems are numerically simulated by a modified displacement discontinuity method using higher order displacement discontinuity elements and higher order special crack tip elements for crack tip treatment to increase the accuracy of the Mode I and Mode II stress intensity factors obtained based on linear elastic fracture mechanics theory. The crack propagation and coalescence paths of the inclined crack are estimated by implementing a suitable iteration algorithm of incremental crack length extension in a direction predicted by using the maximum tangential stress criterion. The numerical and analytical crack extension analyses are compared which are in good agreement and show the validity, applicability and accuracy of the present work.     

裂纹扩展过程的动态仿真技术

岩石断裂力学在岩土工程领域应用越来越广泛,其数值模拟技术是研究的主要方向之一,也是制约断裂力学发展的瓶颈。基于扩展有限元思想,依托ABAQUS用户子单元二次开发技术,引入Heaviside函数模拟不连续位移场,采用裂尖 逼近函数模拟裂尖位移场,通过J积分求解裂尖的应力强度因子（SIF）,最后采用最大周向应力准则,预测裂纹扩展方向,实现了裂纹动态演化过程的扩展有限元数值模拟。通过3个经典断裂力学模型,对比了计算值与理论值,验证了所开发程序的适用性和有效性     

Uniform asymptotic Green's functions for efficient modeling of cracks in elastic layers with relative shear deformation controlled by linear springs

We present a uniform asymptotic solution (UAS) for a displacement discontinuity (DD) that lies within the middle layer of a three‐layer elastic medium in which relative shear deformation between parallel interfaces is controlled by linear springs. The DD is assumed to be normal to the two interfaces between the elastic media. Using the Fourier transform method we construct a leading term in the asymptotic expansion for the spectral coefficient functions for a DD in a three‐layer‐spring medium. Although a closed‐form solution will require a solution in terms of an infinite series, we demonstrate how this UAS can be used to construct highly efficient and accurate solutions even in the case in which the DD actually touches the interface. We compare the results using the Green's function UAS solution for a crack crossing a soft interface with results obtained using a multi‐layer boundary element method. We also present results from an implementation of the UAS Green's function approach in a pseudo‐3D hydraulic fracturing simulator to analyze the effect of interface shear deformation on the fracture propagation process. These results are compared with field measurements. Copyright © 2008 John Wiley & Sons, Ltd.     

Computational model coupling mode II discrete fracture propagation with continuum damage zone evolution

We propose a numerical method that couples a cohesive zone model (CZM) and a finite element‐based continuum damage mechanics (CDM) model. The CZM represents a mode II macro‐fracture, and CDM finite elements (FE) represent the damage zone of the CZM. The coupled CZM/CDM model can capture the flow of energy that takes place between the bulk material that forms the matrix and the macroscopic fracture surfaces. The CDM model, which does not account for micro‐crack interaction, is calibrated against triaxial compression tests performed on Bakken shale, so as to reproduce the stress/strain curve before the failure peak. Based on a comparison with Kachanov's micro‐mechanical model, we confirm that the critical micro‐crack density value equal to 0.3 reflects the point at which crack interaction cannot be neglected. The CZM is assigned a pure mode II cohesive law that accounts for the dependence of the shear strength and energy release rate on confining pressure. The cohesive shear strength of the CZM is calibrated by calculating the shear stress necessary to reach a CDM damage of 0.3 during a direct shear test. We find that the shear cohesive strength of the CZM depends linearly on the confining pressure. Triaxial compression tests are simulated, in which the shale sample is modeled as an FE CDM continuum that contains a predefined thin cohesive zone representing the idealized shear fracture plane. The shear energy release rate of the CZM is fitted in order to match to the post‐peak stress/strain curves obtained during experimental tests performed on Bakken shale. We find that the energy release rate depends linearly on the shear cohesive strength. We then use the calibrated shale rheology to simulate the propagation of a meter‐scale mode II fracture. Under low confining pressure, the macroscopic crack (CZM) and its damaged zone (CDM) propagate simultaneously (i.e., during the same loading increments). Under high confining pressure, the fracture propagates in slip‐friction, that is, the debonding of the cohesive zone alternates with the propagation of continuum damage. The computational method is applicable to a range of geological injection problems including hydraulic fracturing and fluid storage and should be further enhanced by the addition of mode I and mixed mode (I+II+III) propagation. Copyright © 2016 John Wiley & Sons, Ltd.     

扩展有限元法在岩石断裂力学中的应用

岩石裂纹的扩展是一个经典的不连续问题,常规有限元方法难以实现裂纹扩展过程的仿真模拟。扩展有限元法(XFEM)实现了计算网格与不连续面相互独立,因此模拟移动的不连续面时无需对网格进行重新剖分。本文介绍了XFEM基本原理和岩石断裂力学常用判据,尝试对岩石类材料单缝Ⅰ型三点弯曲、单缝剪切和双缝平板实验进行模拟。分析结果表明:扩展有限元模拟岩石类材料断裂问题不受网格划分限制,裂纹以实际应力场分布随机扩展;直观地给出岩样的微裂纹产生、演化,直至完全破坏的全过程,并与实验结果吻合。该方法能够应用到岩石断裂力学方面的研究,模拟岩石类材料的宏细观破坏过程,为解决复杂问题提供了方便的途径。     

Propagation of a hydraulic fracture parallel to a free surface

This paper analyses the plane strain problem of a fracture, driven by injection of an incompressible viscous Newtonian fluid, which propagates parallel to the free surface of an elastic half‐plane. The problem is governed by a hyper‐singular integral equation, which relates crack opening to net pressure according to elasticity, and by the lubrication equations which describe the laminar fluid flow inside the fracture. The challenge in solving this problem results from the changing nature of the elasticity operator with growth of the fracture, and from the existence of a lag zone of a priori unknown length between the crack tip and the fluid front. Scaling of the governing equations indicates that the evolution problem depends in general on two numbers, one which can be interpreted as a dimensionless toughness and the other as a dimensionless confining stress. The numerical method adopted to solve this non‐linear evolution problem combines the displacement discontinuity method and a finite difference scheme on a fixed grid, together with a technique to track both crack and fluid fronts. It is shown that the solution evolves in time between two asymptotic similarity solutions. The small time asymptotic solution corresponding to the solution of a hydraulic fracture in an infinite medium under zero confining stress, and the large time to a solution where the aperture of the fracture is similar to the transverse deflection of a beam clamped at both ends and subjected to a uniformly distributed load. It is shown that the size of the lag decreases (to eventually vanish) with increasing toughness and compressive confining stress. Copyright © 2005 John Wiley & Sons, Ltd.     

巷道裂隙围岩稳定性影响因素的数值分析

对具有单一裂隙的巷道围岩问题,用边界元分区算法沿裂隙面将围岩分区,对每个子域分别建立边界积分方程。在裂端配置1/4面力奇异单元,对张开或闭合的裂隙,既模拟了裂端位移 的变化规律,又模拟了裂端应力1/ 的变化规律。以应力强度因子为指标,依据岩石断裂理论,系统地分析了裂隙几何位置、双向加载侧压比、裂隙面摩擦系数、裂隙长度等因素对围岩稳定性的影响,得出一些有价值的结论。该数值方法适用于分层介质中及复杂几何形状的巷道裂隙围岩稳定性分析。     

Effects of non‐uniform traction and specimen height in the direct shear test on stress and deformation in a rock fracture

In the direct shear test (DST), an internal moment is distributed within the rock specimen by non‐coaxial shear loads applied to the specimen, which cause non‐uniform distributions of both the traction on the loading planes and the stress and deformation in the specimen. To examine the validity of the DST for a rock fracture and to clarify the effect of specimen height, both the stress and deformation in a fracture in the DST were analyzed for specimens with three different heights using a three‐dimensional finite element method with quadratic joint elements for a fracture model. The constitutive law of the fracture considers the dependence of the non‐linear behavior of closure on shear displacement and that of shear stiffness on normal stress and was implemented in simulation code to give a conceptional fracture with uniform mechanical properties to extract only the effect of non‐uniform traction on the stress and deformation in the fracture. The results showed that both normal and shear stresses are concentrated near the end edges of the fracture, and these stress concentrations decrease with a decrease in the specimen height according to the magnitude of the moment produced by the non‐coaxial shear loads. Furthermore, although closure is greater near the end edges of the fracture, where normal stress is concentrated, this concentration of closure is not so significant within the range of this study because of the non‐linear behavior of closure, that is, closure does not significantly increase with an increase in normal stress at large normal stresses. Copyright © 2012 John Wiley & Sons, Ltd.     

Fracture mechanics analysis of multiple cracks in anisotropic media

This paper presents a single‐domain boundary element method (BEM) for linear elastic fracture mechanics analysis in the two‐dimensional anisotropic material. In this formulation, the displacement integral equation is collocated on the un‐cracked boundary only, and the traction integral equation is collocated on one side of the crack surface only. A special crack‐tip element was introduced to capture exactly the crack‐tip behavior. A computer program with the FORTRAN language has been developed to effectively calculate the stress intensity factors of an anisotropic material. This BEM program has been verified having a good accuracy with the previous researches. Furthermore, by analyzing the different anisotropic degree cracks in a finite plate, we found that the stress intensity factors of crack tips had apparent influence by the geometry forms of cracks and media with different anisotropic degrees. Copyright © 2010 John Wiley & Sons, Ltd.     

Evaluating the stress intensity factors of anisotropic bimaterials using boundary element method

This paper presents a boundary element method (BEM) procedure for a linear elastic fracture mechanics analysis in two‐dimensional anisotropic bimaterials. In this formulation, a displacement integral equation is only collocated on the uncracked boundary, and a traction integral equation is only collocated on one side of the crack surface. A fundamental solution (Green's function) for anisotropic bimaterials is also derived and implemented into the boundary integral formulation so that except for the interfacial crack part, the discretization along the interface can be avoided. A special crack‐tip element is introduced to capture the exact crack‐tip behavior. A computer program using FORTRAN has been developed to effectively calculate the stress intensity factors of an anisotropic bimaterial. This BEM program has been verified to have a good accuracy with previous studies. In addition, a central cracked bimaterial Brazilian specimen constituting cement and gypsum is prepared to conduct the Brazilian test under diametral loading. The result shows that the numerical analysis can predict relatively well the direction of crack initiation and the path of crack propagation. Copyright © 2007 John Wiley & Sons, Ltd.     

An algorithm for the simulation of curvilinear plane-strain and axisymmetric hydraulic fractures with lag using the universal meshes

We present an algorithm to simulate curvilinear hydraulic fractures in plane strain and axisymmetry. We restrict our attention to sharp fractures propagating in an isotropic, linear elastic medium and driven by the injection of a laminar, Newtonian fluid governed by lubrication theory, and we require the existence of a finite lag region between the fluid front and the crack tip. The key novelty of our approach is in how we discretize the evolving crack and fluid domains: we utilize universal meshes (UMs), a technique to create conforming triangulations of a problem domain by only perturbing nodes of a universal background mesh in the vicinity of the boundary. In this way, we construct meshes, which conform to the crack and to the fluid front. This allows us to build standard piecewise linear finite element spaces and to monolithically solve the quasistatic hydraulic fracture problem for the displacement field in the rock and the pressure in the fluid. We demonstrate the performance of our algorithms through three examples: a convergence study in plane strain, a comparison with experiments in axisymmetry, and a novel case of a fracture in a narrow pay zone.     

Fuzzy parameters analysis of time‐dependent fracture of concrete dam models

In order to apply the mechanical properties (measured on material specimens or laboratory‐sized models) to large structures (such as concrete dams), a non‐linear theory able to predict the size‐scale effect has to be used. One of these theories was first proposed by Hillerborg and co‐workers (fictitious crack model) and is based on the earlier works by Barenblatt and Dugdale for metals (cohesive crack model). It is based on the existence of a fracture process zone (FPZ), where the material undergoes strain softening. The behaviour of the material outside the FPZ is linear elastic. A large number of short‐time laboratory tests were executed, by varying the load, under crack mouth opening displacement control. Since concrete exhibits a time‐dependent behaviour, an interaction between creep and micro‐crack growth occurs in the FPZ. Therefore, different testing conditions can be applied: rupture can be achieved by keeping the load constant before peak value (pre‐peak tests), or after peak value and after an unloading and reloading procedure (post‐peak tests). The crack propagation rate is shown to be small enough to neglect inertial forces and large enough to keep the time‐dependent behaviour of the process zone as dominant compared to the behaviour of the undamaged and viscoelastic zone. Due to the variability in material microstructure from one specimen to another, experimental data show large ranges of scatter. Well established methods in probability theory require sufficient experimental data in order to assume a probability density distribution. The objective of this study is to investigate the ranges of variation of the time response under constant load in simple structural elements associated with pre‐selected variation (fuzziness) in the main material parameters. For situations where the values of the material parameters are of a non‐stochastic nature, the fuzzy set approach to modelling variability has been proposed as a better and more natural approach. Copyright © 2002 John Wiley & Sons, Ltd.     

Numerical simulation of crack propagation in layered formations

In the present study, a boundary element method based on the higher order displacement discontinuity formulation is presented to solve the general problem of hydraulic fracture propagation in layered formations. Displacement collocation technique is employed to model the higher order displacement variation along the crack and the special crack tip element near its ends. The hydraulic fracture propagation and its interaction with the layer interface in non-homogenous rock materials are studied by the proposed semi-analytical (hybridized boundary element-boundary collocation) method. The maximum tangential stress criterion (or σ-criterion) of fracture mechanics considering different elastic constants (Young modulus and Poisson’s ratio) is used to obtain the fracture path. The fracture propagation from stiff to soft and soft to stiff media for cracks having different inclination angles is modeled, and the effects of elastic constants on the hydraulic fracture propagation is studied. The results show that if the hydraulic fracture originates in the stiffer layer, its capability to cross the layer increases and is vice versa for the softer material. The comparison of the results gained from the numerical method with those in the literature show a good performance of the method in the case of propagation of hydraulic fracture in layered formations.     

水力裂缝穿越非连续面扩展时的断裂过程研究

非连续面发育是非常规油气储层的显著地质特征之一,水力裂缝能否穿越非连续面扩展会关系到压裂的改造效果。为研究水力裂缝穿越非连续面扩展时断裂过程区（fracture process zone,简称FPZ）发育特征,采用自主设计的可视化压裂试验装置对含预制摩擦界面的砂岩平板试件开展水力压裂试验。基于数字图像相关法实时监测了水力裂缝正交穿越界面扩展过程中的位移及应变场特征。试验结果表明,水力裂缝穿越界面扩展之前,断裂过程区已经开始跨越界面发育;裂缝能否穿越界面扩展在FPZ的初始发育阶段已经注定,不受FPZ内应力软化过程影响。基于Renshaw-Pollard准则建立了考虑FPZ边界范围的裂缝穿越非连续面扩展准则,并通过前人及文中试验数据进行了可靠性验证。相比而言,改进准则更准确地考虑了裂缝前端线弹性断裂力学的适用范围。研究发现FPZ长宽比对裂缝穿越界面扩展准则有显著影响,相同条件下,FPZ长宽比越大,裂缝正交穿越界面扩展所需要的摩擦系数下限值越小。     

Modelling crack propagation in concrete structures with a two scale approach

A simplified computational technique based on a refined global–local method is applied to the failure analysis of concrete structures. The technique distinguishes the scale of the structure, modelled with large size finite elements, from the scale at which material non‐linearity occurs due to progressive cracking and macro‐crack propagation. The finite element solution is split into two parts: a linear elastic analysis on a coarse mesh over the entire structure and a non‐linear analysis over a small part of the structure where a dense finite element grid is employed. In the non‐linear calculation, a non‐local damage model is implemented. These two computations are coupled with the help of an iterative scheme. The size and location of the region where a non‐linear analysis is performed, is adapted to follow the development of the damage zone. Numerical examples of mode I fracture of concrete specimens with straight and curved cracks are presented. Copyright © 2003 John Wiley & Sons, Ltd.     

Residual strength of slip zone soils

Slip zones of ancient landslides are commonly composed of fine-grained soils with amount of coarse-grained particle. Residual strength of slip zone soil is an important parameter for evaluating reactivation potential and understanding progressive failure mechanism. In this study, the residual strength is examined by in situ direct shear tests, improved laboratory reversal shear box test, precut specimen triaxial shear test and ring shear test. Some residual shear behaviors are recognized. Field residual strength is the average operational resistance along the sliding surface not an ideal drained strength, which is less than peak and greater than residual strength measured in laboratory. Stress–displacement curves obtained from in situ shear and laboratory reversal direct shear demonstrate strain-hardening which have no significant peak, but the shear stress is decreased gradually with increasing displacement. Residual friction coefficient depends on the normal stress, and this dependence is relevant to the interaction of rolling and sliding of particles. Residual friction angle is closely related to coarse fraction and dry density, appearing a linear increase with increasing coarse fraction and a form of polynomial function with increasing dry density. The influence of shearing rate on residual strength can be negligible.     

Hydro‐mechanical modeling of cohesive crack propagation in multiphase porous media using the extended finite element method

In this paper, a numerical model is developed for the fully coupled hydro‐mechanical analysis of deformable, progressively fracturing porous media interacting with the flow of two immiscible, compressible wetting and non‐wetting pore fluids, in which the coupling between various processes is taken into account. The governing equations involving the coupled solid skeleton deformation and two‐phase fluid flow in partially saturated porous media including cohesive cracks are derived within the framework of the generalized Biot theory. The fluid flow within the crack is simulated using the Darcy law in which the permeability variation with porosity because of the cracking of the solid skeleton is accounted. The cohesive crack model is integrated into the numerical modeling by means of which the nonlinear fracture processes occurring along the fracture process zone are simulated. The solid phase displacement, the wetting phase pressure and the capillary pressure are taken as the primary variables of the three‐phase formulation. The other variables are incorporated into the model via the experimentally determined functions, which specify the relationship between the hydraulic properties of the fracturing porous medium, that is saturation, permeability and capillary pressure. The spatial discretization is implemented by employing the extended finite element method, and the time domain discretization is performed using the generalized Newmark scheme to derive the final system of fully coupled nonlinear equations of the hydro‐mechanical problem. It is illustrated that by allowing for the interaction between various processes, that is the solid skeleton deformation, the wetting and the non‐wetting pore fluid flow and the cohesive crack propagation, the effect of the presence of the geomechanical discontinuity can be completely captured. Copyright © 2012 John Wiley & Sons, Ltd.     

Numerical modeling of dynamic rock fracture with a combined 3D continuum viscodamage‐embedded discontinuity model

This paper deals with numerical modeling of dynamic failure phenomena in rate‐sensitive quasi‐brittle materials, such as rocks, with initial microcrack populations. To this end, a continuum viscodamage‐embedded discontinuity model is developed and tested in full 3D setting. The model describes the pre‐peak nonlinear and rate‐sensitive hardening response of the material behavior, representing the fracture‐process zone creation, by a rate‐dependent continuum damage model. The post‐peak response, involving the macrocrack creation accompanied by exponential softening, is formulated by using an embedded displacement discontinuity model. The finite element implementation of this model relies upon the linear tetrahedral element, which seems appropriate for explicit dynamic analyses involving stress wave propagation. The problems of crack locking and spreading typical of embedded discontinuity models are addressed in this paper. A combination of two remedies, the inclusion of viscosity in the spirit of Wang's viscoplastic consistency approach and introduction of isotropic damaging into the embedded discontinuity model, is shown to be effective in the present explicit dynamics setting. The model performance is illustrated by several numerical simulations. In particular, the dynamic Brazilian disc test and the Kalthoff–Winkler experiment show that the present model provides realistic predictions with the correct failure modes and rate‐dependent tensile strengths of rock at different loading rates. The ability of initial embedded discontinuity populations to model the initial microcrack populations in rocks is also successfully tested. Copyright © 2016 John Wiley & Sons, Ltd.     

A rate‐dependent cohesive crack model based on anisotropic damage coupled to plasticity

In quasi‐brittle material the complex process of decohesion between particles in microcracks and localization of the displacement field into macrocracks is limited to a narrow fracture zone, and it is often modelled with cohesive crack models. Since the anisotropic nature of the decohesion process in separation and sliding is essential, it is particularly focused in this paper. Moreover, for cyclic and dynamic loading the unloading, load reversal (including crack closure) and rate dependency are essential features that are included in a new model. The modelling of degradation is based on a ‘localized’ version of anisotropic continuum damage coupled to inelasticity. The concept of strain energy equivalence between the states in the effective and nominal settings is adopted in order to define the free energy of the interface. The proposed fracture criterion is of the Mohr type, with a smooth transition of the failure and kinematics (slip and dilatation) characteristics between tension and shear. The chosen potential, of the Lemaitre‐type, for evolution of the dissipative processes is additively decomposed into plastic and damage parts, and non‐associative constitutive equations are obtained. The constitutive equations are integrated by applying the backward Euler rule and by using Newton iteration. The proposed model is assessed analytically and numerically and a typical calibration procedure for concrete is proposed. Copyright © 2006 John Wiley & Sons, Ltd.