A finite element framework for modeling internal frictional contact in three-dimensional fractured media using unstructured tetrahedral meshes

Highlights


Introduction
Understanding the mechanical behavior of heavily cracked materials under different mechanical and thermal loads is of vital importance and great interest to a variety of fields, including material science, geothermal energy production, mining engineering, oil and gas reservoir engineering, and structural and earthquake engineering. Geological formations are examples of fractured media at large scales, where rock joints have been shown to extend to lengths ranging from hundreds to thousands of meters [1]. Pre-existing natural fractures in rock masses act as local points of mechanical weakness, and as main flow pathways, and therefore control not only the deformation and strength behavior of the rock mass, but also its flow [2][3][4] and transport [5] properties. Experimental and numerical investigations show that normal closure and shear dilation can significantly change the fracture transmissivity [6,7]. Fluid flow in fractured rock masses is therefore strongly stress-dependent, with regards to the magnitude and orientation of principal permeabilities [8,9]. Accurate prediction of fluid pressure and solid deformation in fractured rocks, therefore, requires hydro-mechanically coupled models with the ability to resolve normal and shear components of contact tractions acting on the fractures [10].
An exact geometric representation of naturally fractured media is challenging for two main reasons. The first is related to the matter of scales and fracture size distribution. Observations suggest that fracture size is governed by power-law scaling models, spanning orders of magnitudes of length scales [1]. The second issue is fracture characterization, for which non-invasive methods to map fractures in situ have yet to be developed for more accurate fracture representations. Stochastic models are therefore required to investigate deformation/flow characteristics of fractured media. Stochastic models often use idealized fracture shapes [11,12], based on a statistical description of parameters such as distributions of size and orientation [13,14]. Models that rely on an explicit representation of fractures, as opposed to a continuum formulation, have been termed Discrete Fracture Network (DFN) models [15,16]. The concept of DFN was first introduced by Long et al. [17] for homogenizing complex fracture networks, and has been extensively used for flow/transport applications [8,[18][19][20][21]. Nevertheless, despite great geometrical simplifications, this type of modeling approach is routinely applied to estimate effective values of engineering parameters relevant to fluid flow, e.g. permeability, [22] as well as mechanical deformation, e.g. Young's modulus [23][24][25][26].
The majority of numerical simulations of fracture networks have been conducted using discrete element method (DEM), whereas the use of the finite element method (FEM) has been limited to a few studies. DEM has also been very popular in simulating fracture growth and fragmentation of brittle solids, such as granular materials and rock and concrete [27][28][29][30]. DEM generally treats the fractured medium as the assemblage of separated blocks formed by connected fractures, solves the equation of motion for the blocks, and updates the contact between the block as a consequence of the motion and deformation of the blocks [31]. The distinct element method introduced by [32], with the commercial computer codes UDEC and 3DEC for 2-and 3-D problems [33,34], and the discontinuous deformation analysis (DDA) proposed by [35], have been the main approaches for analyzing the deformation and permeability of fractured rock masses [36,37,8,18]. DDA uses standard FEM meshes over blocks, and employs the penalty method for enforcing the contact constraint between blocks. A similar development to DDA is combined FEM/DEM, introduced by [38], which considers not only the block deformation but also fracturing and fragmentation of the rocks [39]. The application of DEM in modeling fracture growth and fragmentation entails the following difficulties: (1) Time-consuming and error-prone calibration of micro-to macro-properties must be performed for each material individually. Thus, elastic mechanical properties such as the Young's modulus and Poisson's ratio cannot be directly used to model elastic deformation. Moreover, the calibrated properties are scale-and mesh-size dependent [40]. (2) Fractures are not explicitly defined; in fact, they are modeled as the lack of cohesion between the particles in the material. (3) For fragmentation purposes, the materials often artificially behave as particulates or agglomerates. Therefore, 3D fragmentation simulations and qualitative pattern evaluation are scarce in relation to the maturity of DEM, possibly due to the lack of realism caused by the absence of a fracture mechanics-based crack growth models. The property calibration can be avoided by using an impulse-based method, which can be enriched with energy-conservative tracking to ensure energy conservation during contact [41,42]. Regarding the application of DEM in deformation and flow response of fractured networks, the following drawbacks are highlighted: (1) Isolated fractures are ignored when using DEM in modeling the fracture networks, and fractures are only modeled as the boundaries of isolated blocks; (2) The deformation inside blocks and the contact forces between the block are roughly approximated, as explicit methods are generally used to solve the balance equations. (3) The high stress gradients near the cracks cannot be captured accurately, and the variation of contact tractions over the contact surfaces is estimated roughly. The deformation of fracture surfaces, which controls the aperture change within fractures, is also only roughly estimated.
In contrast, the finite element method is able to capture the high gradient stress state near the crack, and can provide very accurate contact tractions based on implicit methods. Advantages of FEM for modeling fracture networks include: (1) Meso-scale definition of material properties such as elastic constants and material toughness is directly used in the method; (2) Fractures are explicitly modeled as local discontinuities in the continuum medium, where singular stress fields are modeled through adapting appropriate element size and type in the crack front region in a FE model, or appropriate enrichment functions in an XFEM formulation; (3) Fracture mechanics-based parameters such as stress intensity factor and J -integral can be used to study the onset of fracture growth and fragmentation in fractured media; (4) Crack interactions are captured accurately. The use of the finite element method to model contact between crack surfaces has been mainly limited to the XFEM formulation which allows modeling of the entire crack geometry independently of the finite element mesh [43][44][45][46][47][48][49][50][51][52][53]. Both LATIN, LArge Time INcrement, [43,44,[50][51][52] and Newton-Raphson [45,46,49] iterative strategies have been employed in dealing with the nonlinearity of the contact problem. However, Liu and Borja [47] demonstrated the superior convergence performance of Newton-Raphson method as compared to the LATIN strategy. The previous works mainly focus on XFEM modeling of two-dimensional cracks and interfaces which are initially closed, yielding a low contact precision model. Moreover, the accuracy of the contact tractions near the crack tip/front has not been investigated. This accuracy directly influences the computation of stress intensity factors when using an energy-based method such as the interaction integral. No validation of the accuracy of the stress intensity factors has been reported in previous work. The present work proposes a finite element formulation based on tetrahedral elements for enforcing contact constraint on 3D initially open cracks in high density fractured media, with special attention placed on the accurate resolution of contact tractions near the crack front, and accurate computation of stress intensity factors.
Due to the geometrical complexity of the explicit representation of fractures in a continuum medium, unstructured meshes using tetrahedral elements are preferred to model fractured media. Meshing procedures using tetrahedra are much simpler, and these elements are well suited to automatically mesh complicated geometries. Unstructured meshes have been successfully used in the context of FE computation of fracture parameters [54,55], and FE simulation of crack propagation [56][57][58][59][60][61] as well as fragmentation [62].
Frictional contact in fracture analysis has also been performed to study the role of stick and slip in the problem of the closure of a crack represented by a plane of dislocations in a bending infinite plane [63]. Methods to model frictional sliding cracks have been proposed, evaluating the role of the tip singularity on near-tip slip [64] and studying stress intensity factors of a compressive crack in equilibrium [65]. The Augmented Lagrangian (AL) method has been successful for enforcing the contact constraint accurately when computing high contact precisions, by combining the Lagrange multiplier and penalty methods to exploit the merits of both approaches [66,67]. The AL method, first applied to contact problems by Alart and Curnier [68], renders the objective function smooth across contact-no-contact conditions [69]. The advantages of AL over the penalty method include decreased ill-conditioning of the system, and essentially exact satisfaction of constraints with finite penalties. The advantage over Lagrange multipliers is that it avoids introducing unknowns to the problem by using current fixed estimates of the Lagrange multipliers. Applying AL often yields a double loop algorithm in which constant Lagrange multipliers are used with penalty terms during the inner loop to enforce the contact constraints. Then, within an outer loop, the Lagrange multipliers are updated to new values based on the computed tractions in the inner loop. This procedure increases the number of iterations, but allows enforcing contact constraints accurately by using small penalty parameters (see [70][71][72]). This type of algorithm has been referred to as nested AL [70], or Uzawa-type, algorithm [73].
The present paper presents a finite element formulation that uses unstructured quadratic tetrahedral meshes to model internal contact in fractured media. A sophisticated algorithm is developed for the treatment of frictional contact between the fracture surfaces, based on isoparametric integration-point-to-integration-point discretization of the contact contribution. The contact constraints are also enforced by using a gap-based AL. The frictional contact algorithm proposed here is novel because, firstly, a singular variation of penalty parameter is suggested near the crack front to circumvent the difficulty of zero gaps on the crack front nodes. Secondly, a gap-based AL is introduced for updating the contact tractions obtained from the penalty solution to obtain new better estimates. As opposed to the conventional traction-based AL, where Lagrange multipliers are augmented, in this proposed methodology gaps are augmented, which allows one to circumvent the difficulty of defining and augmenting Lagrange multipliers at the crack front nodes. In order to model the strain singularity along the crack front of fractures, quarter-point tetrahedral elements are used at the crack front region. Displacement correlation and domain integral methods are used to compute the point-wise stress intensity factors. The proposed FE formulation is able to compute fracture contact tractions and high stress gradients near the crack front very accurately.

Problem description
Consider a body containing randomly distributed interacting and intersecting cracks embedded in an elastic medium (see Fig. 1). Each crack is represented by two smooth surfaces attached at the tips. They can be initially closed or open with a specific aperture distribution. When open, crack surfaces are assumed traction-and cohesionfree. When in contact, it is assumed that a Coulomb frictional law governs the boundary conditions over the crack surfaces. The two surfaces of each crack intersect at a curve called the crack front or tip, at which a strain singularity occurs. The volumetric body quasi-statically compressed and friction on the crack surfaces is computed.
Consider a solid in its equilibrium configuration due to deformation u : , also deformed by an arbitrary virtual infinitesimal displacement field δu : Ω × [0, T ] → R 3 which satisfies the displacement boundary condition: δu =ū on Γ u × [0, T ]. This virtual displacement must also be admissible with regard to the displacement constraints of the contact condition. By applying the principle of virtual work to the material points of Ω , the weak form of the equilibrium equation is formed as an integral equation: subjected to the classical Coulomb friction condition and the Kuhn-Tucker inequality condition [73] for normal and tangential tractions, p and τ , on the contact master surfaces Γ m . The admissible variation of the deformation δu is also constrained by the contact condition. Here, δε = 1/2[∇δu + (∇δu) T ] is the strain tensor of the virtual displacement field δu. The virtual normal and tangential gaps are also developed based on the virtual displacement field as δg N = (δu s − δu m ) · n and δg T = (I − n ⊗ n) (δu s − δu m ). The term δΠ Ω indicates the contribution of internal stresses and body forces to the total virtual work δΠ , and the terms δΠ σ and δΠ c include the virtual work of pre-defined tractions and contact forces, respectively. In contrast to the strong form, which is a set of differential equations along with prescribed boundary conditions, the weak form is in the form of an integral equation which requires a weaker continuity on the displacement. The weak form also honors the pre-defined boundary conditions of the boundary value problem, and is best suited to be the basis of numerical approximations. More details on the problem description can be found in [74]. Remark 1. The boundary of the fractured body is divided into three sets of Dirichlet, Neumann and contact boundaries, Γ = Γ u ∪Γ σ ∪Γ c , as shown in Fig. 1. The boundaries Γ σ and Γ c can generally have overlap regions at the intersection of cracks with external boundaries. However, special care has to be taken when Γ u and Γ c overlap, as the kinematic constraints of one region can prevent enforcing the constraints by the other. In fact, the Dirichlet boundary conditions on contact surfaces can prevent p and τ from satisfying Coulomb inequality conditions. Therefore, in addition to Γ u ∩ Γ σ = ∅, the restriction Γ u ∩ Γ c = ∅ must be applied.
Remark 2. The crack front, Γ f , where master and slave surfaces meet, can only provide the boundaries of the contact region for each embedded crack, and is excluded from the contact area, Γ f ̸ ∈ Γ c . Therefore, no contact traction is defined or applied on the crack front. However, the limit of the contact traction as the crack front is approached, r → 0 where r is the normal distance from the crack front, can be non-zero: (2) Remark 3. The state of strain is singular along the crack front. This indicates that considerable stresses occur adjacent to the crack front even when very small penetration occurs between the fracture surfaces. Therefore, penetration has to be strictly penalized in a penalty treatment, so that accurate contact tractions can be obtained.

Augmented Lagrangian method with Uzawa algorithm
The numerical solutions suggest that near the crack front of an arbitrary 3D crack configuration, a plane-strain condition prevails locally (see Fig. 2), so that the three-dimensional deformation fields approach the two-dimensional plane strain fields [75,76]. According to these fields, the displacement adjacent to the crack front and over a perpendicular plane to the crack front are given as square-root functions of distance from the crack front [77,55]. Therefore, the distribution of relative displacement of the top crack surface with respect to the bottom surface also follows a square-root variation, giving the total gap function of a point located at normal distance r from the crack front as where K i , i = I, II, III is the mode i point-wise stress intensity factor, µ = E/2(1 + ν) is the shear modulus, E and ν are the Young's modulus and Poisson's ratio, the Kolosov constant κ is equal to 3−4ν for the plane strain condition, e i is the unit vector along the x i axis of the local coordinate system located on the crack front, andḡ is the gap vector at the distance r 0 from the crack front. Numerical results from quarter-point tetrahedral elements also capture this type of gap variation near the crack front [55]. This gap variation significantly influences the contact tractions obtained from a penalty formulation near the crack front. To avoid gap values influencing the contact reactions, it is suggested that an inverse square-root variation of the penalty parameter with r is used. Let ϵ 0 be a constant nominal penalty parameter at a distance r 0 from the crack front. The distribution of the penalty parameter near the crack front can then be defined as ϵ = ϵ 0 / √ r/r 0 , generating a square-root singularity of the penalty parameter near the crack front. This regularization of the penalty parameter cancels out the influence of gap variation on the contact traction variation, ensuring a finite value of stick traction very close to the crack front, as shown in Fig. 2: It can also be shown for the slip condition that finite values of the traction are obtained very close to the crack front using this treatment. Although this proposed regularization of the penalty parameter improves the contact tractions near the crack front, the difficulty of defining a Lagrange multipliers distribution over the element attached to the crack front remains. Contact tractions cannot be defined and updated over the crack front, and contact tractions over the master elements attached to the crack front cannot be directly obtained from the tractions at nodal values. This difficulty can be circumvented by augmenting gaps instead of Lagrangian multipliers in AL. In this type of treatment, the normal and tangential Lagrange multipliers are defined as where g * N and g * T are augmented normal and tangential gaps. The augmented gaps are constant during the penalty solution, inner loop, and are updated in the augmentation process, outer loop, based on the penalty solution. In this strategy, thanks to the singularity of the penalty parameter, the augmented gaps are able to determine very accurate contact tractions near the crack front. Algorithm 1 demonstrates the application of this gap-based augmentation Lagrangian treatment in the combination of return mapping strategy for updating the contact tractions in a time increment.
An active strategy is often used to update the regions in contact in every iteration in the inner loop. This strategy identifies the regions in contact by using the value of the normal traction p n+1 = ϵ(g * k N n+1 ), where p n+1 ≤ 0 denotes active contact zone. The slip or stick condition is also determined based on the value of f tr s in Algorithm 1. Algorithm 1 AL approach for the evolution of frictional contact in the time step [t n , t n+1 ] 1. Initialization: 2. Return Mapping and Solution: otherwise (slip).

Convergence check:
Consider Γ m st and Γ m sl being respectively stick and slip zones of the contact master surface. Γ m st ∪ Γ m sl therefore constitutes the active contact zone. According to the AL treatment in Algorithm 1, the normal, tangential stick and tangential slip tractions are defined by p = ϵ(g * N + g N ), τ st = ϵ(g * T + g T ), and τ sl = (µϵ|g * N + g N | + τ c )n T , respectively. Here, the superscript k and subscript n + 1 are removed for simplicity. Using these expressions, Eq. (1) is rewritten as where δΠ Ω +σ = δΠ Ω − δΠ σ and µ and τ c are the friction coefficient and cohesive stress, respectively. Eq. (6) has to be solved in an iterative manner using Newton's method, which constitutes an inner loop in Algorithm 1. Since the augmented gaps are fixed within the inner loop, the linearization of Eq. (6) depends only on the normal and tangential gaps g N and g T , and the direction of tangential traction n T . The finite element discretization, and the linearization procedure, will be described in Section 4. The augmentation process makes it possible to strictly penalize any violation of contact constraints by using a small penalty parameter. Due to the singularity, even a slight violation of contact constraints can cause considerable stresses near the cracks. Therefore, the contact constraints have to be applied as accurately as possible. This requires the augmentation to continue until the normal gap g N and the tangential stick gap g T become smaller than reasonable values of thresholds Tol N and Tol T , which are dependent on the crack size.

Finite element formulation
The entire domain of the problem is discretized with quadratic ten-noded tetrahedral elements, as shown in Fig. 3. Quarter-point tetrahedral elements are employed in the immediate neighborhood(s) of the crack front(s), while the remainder of the domain is discretized with the standard tetrahedral elements [55]. The introduction of quarter-point elements to the fractured medium is straightforward. Tetrahedral elements attached to the crack front are identified, and mid-side nodes are shifted to the quarter-point position near the crack front. The internal node numbering of the elements is then modified to become consistent with the numbering used in [55]. This allows one to use simple relations for obtaining the local coordinates of a given point inside a quarter-point element. The use of quarter-point tetrahedrals also introduces quarter-point triangular elements over the crack surfaces near the crack front, as shown in Fig. 3, whereas isoparametric quadratic triangles are used elsewhere on the fracture surfaces.

The contribution of internal/external forces
In order to describe the kinematics of tetrahedral and triangular elements, the vectors of displacements, u i , virtual displacements, δu i , and incremental displacements, ∆u i are introduced for every node i. The virtual work due to internal stresses and body forces throughout the domain is given by where N i is the shape function corresponding to node i. The summations over domains Ω and Γ σ include all tetrahedral elements, and triangular elements which are subjected to pre-defined external tractions, respectively. Ω e and Γ e denote the domain of tetrahedral and triangular elements. The sum over p includes element integration points where the bracketed quantities {} p and [] p are evaluated and multiplied by the corresponding weight w p . |J| also denotes the determinant of the coordinate Jacobian matrix of the elements. Matrix B i contains the derivative of the shape functions associated with node i in local coordinate, and D is the elasticity matrix containing the material properties [78]. The contribution is linearly dependent on the displacement, and the associated tangent matrix is the so-called stiffness matrix given by

Contact kinematics and contribution
Consider two conformed quadratic triangular elements on the crack surfaces. One is referred to as master element, m, and the other is denoted as slave element, s, as shown in Fig. 3. Considering the global coordinate system in Fig. 3, all the top and bottom crack surfaces are assigned as slave and master surfaces, respectively. Each node i on the master element is paired with the matched node i on the slave surface, constructing the nodal pair i. Assuming small deformation, the geometry of the contact surface can be represented by the reference configuration of master surface X m . The expression of the unit normal to the master surface is given by where ∂ X m /∂α =  6 i=1 X m i ∂ N i /∂α, α = ξ, η are the tangent vectors to the master surface. Depending on the internal node numbering of the master element, two normals with opposite directions are obtained from Eq. (9). A normal vector, outward to the master surface, shall be used (see Fig. 3). The two elements are separated by a small initial aperture, leading to an initial normal gap (−ĝ N ). Since the initial distance between the contact elements is small compared to the size of contact elements, it is not necessary to explicitly apply this gap between the elements in the geometrical specification. In fact, the implicit presence of the gap in the gap formulation would suffice, leaving the geometry of the domain unchanged during the entire contact analysis. Considering a certain fracture aperture distribution, the initial normal gap is distributed over the master nodes of the fracture surface. Considering the slave node i with the initial normal gapĝ i N , the initial normal gap distribution over the contact element is given bŷ On the other hand, the initial tangential gap has to be computed by taking into account the displacement field at the time that the master and slave surfaces first come into contact. Consider that the slave element penetrates the master surface during the current time increment t c = [t n , t n+1 ], when the relative displacement at the penetrating nodal pair i in the first iteration is (u s i − u m i ) c . The initial tangential gap distribution over the contact element is then approximated bŷ at the end of the first iteration and is used in the following iterations of this increment to apply the contact constraints. This initial tangential gap then remains constant, until the time at which the nodes in the nodal pair lose their contact again. The distribution of total initial gap over the master element will therefore bê In order to describe the kinematics of each nodal pair, the vectors of displacements, virtual displacements and incremental displacements are introduced as In addition, the following vector and matrix are introduced, based on the local normal at the master element: Based on these definitions and employing Eq. (13), the discretized version of gap functions in normal and tangential directions are given by The variation of the gap functions due to a virtual displacement, δg, and the variation of gap function due to incremental displacement, ∆g, in the normal and tangential directions are also given as The accurate computation of contact tractions near the crack front requires a square-root singular variation of the penalty parameter near the crack front. This type of variation can be applied simply over the entire fracture surface, in the case of well-defined crack shapes such as penny-shaped and elliptical cracks. For example, for an elliptical crack defined by the equation where a and b are the minor and major semi-axes of the ellipse in the local x ′ y ′ coordinate system, the parameter ϵ r = and corresponding values are applied on nodes over the master surface. The distribution of the penalty parameters over any master element is then obtained by where ϵ i r holds the value of ϵ r at the node i. In the cases of complex fracture geometries, the singular square-root variation of the penalty parameter can be applied only adjacent to the crack front, perhaps over the quarter-point triangular elements only. In this case, the application of a constant penalty parameter ϵ = ϵ 0 would suffice for the remainder of the fracture surface.
Augmented gaps are also obtained based on the gap discretization in Eq. (15). The (k + 1)th augmented gaps are given by where where u c i k is the displacement associated with the nodal pair i at kth iteration. n k+1 T is also the direction of tangential traction at the (k + 1)th augmentation iteration. It is obtained through the return mapping process by using the trial tangential gap of the kth augmentation iteration. Consider n k T and g k T to be respectively the direction of tangential slip traction and the tangential gap at the kth augmentation iteration. The direction of tangential traction for the next iteration is then obtained as where As gaps are evaluated at the integration points of the contact elements, the direction of the trial tangential traction is also evaluated there and updated during the augmentation process. For planar cracks, however, as the normals to all master elements are identical for each fracture, the directions of tangential traction can be stored and updated at the nodal pairs and then interpolated to the integration points as Here gaps and augmented gaps corresponding to the nodal pair i are defined as In this case, the augmented nodal gaps are stored for each nodal pair i in g * k N i and g * k T i and updated in the augmentation process.
The contribution of the normal, stick and slip contact tractions are approximated by numerically integrating over the master elements. Using Eqs. (13) and (16), the contact contributions in Eq. (6) are given by where the residual vectors are defined as Summations over area Γ m st and Γ m sl include all the elements domains Γ e in stick and slip conditions, respectively. The sum over p includes element integration points, 'gpts', of the master triangular elements, 'elems', where the bracketed quantities {} p and [] p are evaluated and multiplied by the corresponding weight w p . |J| denotes the determinant of the coordinate Jacobian matrix of the triangular elements. Eqs. (14), (15), (17) and (18) are used to evaluate the parameters in Eq. (24) at integration points. The direction of slip at integration point is also evaluated using Eq. (20). Linearization of Eq. (24) gives the tangent matrices as Remark 4. Some of the tangent matrix components in Eq. (25) are associated with the virtual or actual displacement variation of the gap at the crack front nodes, which are essentially zero. Therefore, the rows and columns corresponding to the crack front nodes must be eliminated from the tangent matrices of the contact elements attached to the crack front.

Contact algorithm and implementation
The Newton-Raphson method is often used to solve the system of nonlinear equations associated with the nonlinear characteristic of contact problems. Once the element residual vectors and tangent matrices are obtained, the residual vector and tangent matrix of the entire system of elements is developed through an assembling process. Let u = {u 1 , u 2 , . . . , u N } be the solution vector containing the displacements of all the nodes of the system, N . Equivalently, δu T = {δu 1 , δu 2 , . . . , δu N } can be defined to include the virtual displacement of the nodes in the system. By substituting Eqs. (7), (8), (24) and (25) into (6), the virtual work of the entire system is developed as Here G is the residual vector of the entire system, where δu T G = 0 indicates that a zero residual vector defines the solution vector u, constituting a set of 3 × N nonlinear equations. Here the application of Newton-Raphson method involves the convergence of a trial solution vector iteratively as where ∆u = {∆u 1 , ∆u 2 , . . . , ∆u N } T is vector of displacement corrections, and K is the tangent matrix of the entire system: K and G are updated in each iteration to include the contribution of all the regions in contact. The Newton-Raphson iteration continues until the norm of ∆u becomes less than some tolerance value. Algorithm 2 delineates the steps of the contact algorithm. An initial fracture aperture is modeled by applying an initial fluid (normal) pressure p 0 onto the fracture surfaces, and solving for the deformation. The induced initial normal gaps are then saved as nodal values over the crack surfaces. Afterwards, the fluid pressure is removed, the solid deformation due to fluid pressure is discarded, and the simulation begins by applying external compressive stress to the cube, where, as a result, the crack surfaces might go into contact. This process of introducing initial gap agrees well with the opening process of the natural fractures in geomechanical systems.
Remark 5. The elements of the tangent matrix in Eq. (25) require integration of polynomials of a maximum order of four. The Jacobian determinants of the straight-sided standard and quarter-point triangular elements are polynomials of order zero and two, respectively. Therefore, full integration of the elements of the tangent matrix requires the integration of sixth-order polynomials, which can be achieved by a seven-point integration rule [78]. Note that an integration rule with integration points on the sides cannot be employed, as the Jacobian determinant is zero along the crack front, due to the nonlinear mapping in quarter-point tetrahedrals. A four-point Gauss rule computes exactly the elements of the stiffness matrix of straight-sided standard tetrahedral elements. However, quarter-point tetrahedrals introduce a Jacobian determinant in the form of polynomials of order three. Therefore, a four-point Gauss rule provides a reduced integration scheme for these elements. It has been demonstrated that a five-point Gauss rule can integrate the elements of the stiffness matrix with a level of accuracy that is compatible with the full integration of the components of the tangent matrices of quarter-point triangles. In the present research, five-point and seven-point Gauss rules are used for tetrahedral and triangular elements, respectively.
Algorithm 2 A frictional contact algorithm for analyzing 3D fractured media using gap-based AL

Computation of fracture parameters
In the absence of body forces, the disk-shaped domain representation of point-wise J -integral and interaction integral at any point along the crack front of a planar 3D crack is given by Nejati et al. [54] where A is a disk-shaped area in the plane orthogonal to the crack front at point s, and C + and C − are the contours on the top (slave) and bottom (master) crack surfaces with the outward unit normal m = (0, −1, 0) and m = (0, 1, 0), respectively as shown in Fig. 4. σ kl , ε kl and u k are the Cartesian components of the stress tensor, strain tensor and displacement vector in the local x 1 x 2 x 3 coordinate system, respectively. σ aux kl , ε aux kl , and u aux k are the components of stress tensor, strain tensor and displacement vector due to an auxiliary state which includes the dominant terms in the linear elastic solution of a crack problem. The auxiliary fields are therefore the first terms of the Williams series expansion of crack tip fields [79]. W =  ε 0 σ kl ε kl dε and W I = 1/2(ε aux kl σ kl + ε kl σ aux kl ) are the actual and mutual strain energy densities, respectively. δ kl is the Kronecker delta, and q is a smooth scalar function in the area A, taking the value of unity on the disk's circumference, and vanishing on the crack front. After the interaction integrals corresponding to different auxiliary modes are computed, the stress intensity factors are computed from simple relations (see [54] for details).
Once the boundary value problem subjected to the contact constraints is solved, and the required augmentation steps are performed, Lagrange multipliers provide the contact tractions, and the penalty terms are essentially zero. In this case, the normal and tangential tractions are obtained by p = ϵg * N and τ = ϵg * T , respectively. The contact traction over the master (bottom) surface is therefore given by t m = p · n + τ = ϵg * , where g * = g * N n + g * T . The traction over the slave (top) surface is oriented opposite to the one applied on the master surface (t s = −t m ). The stress components σ 2l applied on both master and slave surfaces are therefore equal to σ 2l = ϵg * · e l where e l is the unit vector of local axis x l in local coordinate system as shown in Fig. 4. The domain integrals in Eqs. (30) and (31) are rewritten as The line integral involves the evaluation of singular integrands defined by the displacement gradients multiplied by the contact tractions. This implies that any small inaccuracy of contact tractions very close to the crack front could potentially significantly influence the value of line integral via the singular terms. Since the contributions of line integrals are significant to the entire J -and interaction integrals, these inaccuracies are most likely to influence the accuracy of the total value of J -and interaction integrals. Therefore, the accuracy of the results of the J -integral and stress intensity factors are heavily dependent on the accuracy of the contact traction near the fracture front. Although the choice of singular penalty parameter makes it possible to obtain accurate contact tractions very close to the crack front, local inaccuracies near the crack front may still be observed. This is because contact tractions in a penalty solution are directly affected by the local displacement inaccuracies near the crack front, which are mainly due to the low quality of randomly placed elements at a region with high stress gradients. These small inaccuracies can then greatly influence the results of J -and interaction integrals, due to the presence of singular displacement gradients. This vulnerability of the domain integrals to potential inaccuracies of the contact traction can however be circumvented by recasting the line integral. As explained earlier, the crack boundary condition dictates a square-root displacement variation of the crack surfaces near the crack front. Employment of a singular square-root penalty parameter ensures that the contact traction is held constant at the region very close to the crack front. Since the domain integrals are evaluated over small region near the crack front, a constant contact traction over C − + C + is therefore expected. Define ∆u l = u l | θ =π − u l | θ =−π and ∆u aux l = u aux l | θ =π − u aux l | θ =−π as the relative actual and auxiliary displacement of slave crack surface with respect to the master crack surface. ∆u and q vanish at the beginning and the end of C − , which helps to recast integrals in Eqs. (32) and (33) using integration by parts as Here, it is assumed that the contour C − is a straight line opposite to the x 1 direction (dC − = −dx 1 ), and the contact traction is constant, in value and direction, over the small contour C − . Recasting the contour integral in Eq. (34) is advantageous for numerical purposes, as the local contact inaccuracies are no longer able to influence the Jand interaction integrals through the singular displacement gradients. Using these new formulations of the contour integrals, the domain integrals are rewritten as The area and contour integrals in Eqs. (35) and (36) are evaluated using a set of virtual quadratic triangular and line elements [54]. These elements are referred to as virtual, since they are not used while performing the finite element solution of the boundary value problem. Consider a point s along the crack front with the local coordinate system x 1 x 2 x 3 . Due to the domain symmetry, only one-quarter of the disk of radius R d is discretized with virtual triangular elements, and the contour C − is discretized by line elements. The integration over the other three quarters is readily evaluated by the reflection of integrating points of the generated virtual elements as shown in Fig. 4. Using the virtual elements, evaluation of the domain integrals in Eqs. (35) and (36) follows the same standard Gauss-quadrature integration scheme available in any FE code:

Numerical examples
In order to demonstrate the accuracy of the proposed contact algorithm and the SIF computation procedure, the deformation behavior of the following fractured body configurations is analyzed under uniaxial uniform compression: (i) Single penny-shaped and elliptical cracks embedded in large cubes; (ii) Two interacting/intersecting penny-shaped cracks embedded in a cube; and (iii) Multiple randomly oriented, randomly placed, penny-shaped cracks in a cube. All these cracked bodies are subjected to a uniaxial compression. All procedures employed in this work are implemented into the Imperial College Geomechanics Toolkit, a geomechanics module [56,58] of the Complex System Modeling Platform (CSMP++), an object-oriented finite element based API developed for the simulation of complex geological processes [80]. The asymmetric matrix of the system of equations resulting from the finite element accumulation during contact resolution is solved using Fraunhofer SAMG Solver version 27a1 [81].
For the cases for which analytical values are available, the numerical error in the computation of the contact tractions, e c , and the SIFs, e t , are respectively evaluated by where t m A and t m N are, respectively, the analytical and numerical contact tractions on master surfaces, K A i and K N i are the pointwise analytical and numerical mode i SIFs, and integrations are performed over the master surface Γ m and the crack front L f . Single and double vertical bars indicate the absolute value of a scalar and the norm of a vector, respectively. Wherever closed form integration was not available, a trapezoidal rule has been employed to evaluate the integrals numerically. The analytical solutions for the SIFs of embedded initially-open inclined penny-shaped and elliptical cracks in infinite solids under uniaxial compression have been derived in [82,74].

Experimental setup
Consider a cube of length 2w containing single or multiple cracks as shown in Figs. 5(a) and 13(a, c). The cube is subjected to a uniform uniaxial compression in the X 2 direction over the top and bottom surfaces. The cracks lie in the plane X 2 = X 1 cot β which generates the angle of β with the direction of applied load. A horizontal single crack configuration (β = 90 • ) produces pure mode I crack deformation in the case of initially open cracks, while the inclined one (0 • < β < 90 • ) provokes a mixed-mode condition. In these configurations, a denotes the crack radius for the penny-shaped crack, and semi-major axis for the elliptical crack. The semi-minor axis b of the elliptical crack is perpendicular to the X 1 X 2 plane. Young's modulus and Poisson's ratio values of E = 10 GPa and ν = 0.3 are used in all models. The penalty parameter ϵ 0 over each fracture is determined individually as the value of Young's modulus over the average size of elements at the crack front region ϵ 0 = E/L n . The average length of the elements L n is defined as the crack front length L f over the number of crack front segments N f (L n = L f /N f ). This choice of penalty parameter generates a well-conditioned system of equations, where the values of the members corresponding to the fracture nodes in the global stiffness matrix are comparable to the value of the members corresponding to nearby nodes. In all cases, the augmentation process continues until no further improvement is seen in the contact tractions, which means the convergence of Lagrange multipliers is reached. Also at each augmentation step, the Newton-Raphson iteration continues until the norm of residual vector becomes smaller than a reasonable threshold.

Mesh
An octree-based mesh generation software was employed to generate arbitrary meshes for all geometries, using ten-noded isoparametric tetrahedral and six-noded triangular elements. This mesh generator is able to split the walls of fractures, and generate matched surface elements over the two surfaces of the cracks. For elements attached to the crack front, the nodes near the front are moved from the mid-side point to the quarter-point position to produce inverse square-root singular fields near the front. The curved crack fronts impose one curved edge for the tetrahedral elements sharing an edge with the crack front. When using quarter-point elements, the Jacobian determinant over small volumes, near the curved edges becomes negative [55]. To avoid this, the curved edges are straightened by moving the mid-side nodes to the center. The refinement of the mesh near the crack front is controlled by assigning the number of segments along the crack front. Assume that the crack front of length L f is discretized by N f segments. A parameter called the nominal length (size) of the elements in the crack front region can be defined as L n = L f /N f . The nominal element length L n represents the approximate length of the element sides near the crack front and therefore gives an approximate for the average size of the quarter-point tetrahedral elements in the crack front region. The degree of mesh refinement in the crack front region is controlled by keeping the nominal crack front element size about one twentieth of the crack length (L n ≈ a/20). Since estimations suggest that the size of the singular dominant zone depends mainly on the crack length, ranging between a/10 and a/50 [83], keeping L n ≈ a/20 ensures that the quarter-point elements at the crack front predominantly remain in the singular dominant zone, where the fields have the inverse square-root singularity. Five-, seven-, and two-point Gaussian quadrature rules are employed for the numerical integration over tetrahedral, triangular, and line elements, respectively.

Details of the SIF computation
For all crack configurations, the mesh-dependent domain radius of R d = L n has been used to generate the virtual domains and compute the fracture parameters. Domains are built at the locations of both corner and mid-side nodes of the segments along the crack front. A similar virtual mesh structure as the one proposed in [54], with four elements in the radial direction (k = 4), was used to compute the SIFs. This choice yields 112 quadratic triangular elements, containing 112 × 3 integration points, together with eight quadratic line elements, containing 8 × 2 integration points (as in [54]). In order to compute the fracture parameters, a smooth function q must be defined over the integration domain. All numerical results in this research are determined by using q = 1 − r/R d , where r = (x 2 1 + x 2 2 ) 1/2 is the distance from the disk center, and R d is the domain radius. The derivatives of this function (∂q/∂ x 1 = −x 1 /r R d and ∂q/∂ x 2 = −x 2 /r R d ) are directly evaluated at the integration points of the virtual triangular elements. For the SIF computation from the DC method, the displacements are correlated at points located at the fixed distance of r m = L n from the crack front [55].

Single penny-shaped/elliptical crack
Consider a single penny-shaped/elliptical crack in a large cube subjected to uniform compression, as shown in Fig. 5(a). The crack-length to body-width ratio of a/w = 0.1 was used in order to eliminate any influence of the cube boundaries on the crack fields. Fig. 5(b) shows the finite element mesh of the penny-shaped crack together with two close-up pictures of the mesh structure over the penny-shaped and elliptical (b/a = 0.4) cracks. The following boundary conditions are applied for this configuration: u 1 = u 2 = u 3 = 0 at the point X 1 = X 2 = X 3 = −w, u 2 = 0 over the plane X 2 = −w, and σ = 1 over the plane X 2 = w. Fig. 5(c) and (d) also show the distribution of ϵ r over the surfaces of the penny-shaped and elliptical cracks. This variation generates a square-root singular penalty variation near the crack front (ϵ = ϵ 0 /ϵ r ), which makes it possible to compute accurate contact tractions very close to the crack front. Fig. 6 shows the distribution of normalized augmented gap and contact traction over the slave surfaces of a slipping penny-shaped crack. The augmented gap maintains zero magnitude along the crack front, increasing towards the center of the crack, where it attains its maximum. The application of an inverse singular penalty parameter variation near the crack front allows computing very accurate contact tractions for the elements attached to the crack front. Application of a constant penalty parameter over the crack surface, however, would not compute accurate contact tractions near the crack front, where the contact traction tends to zero when the crack front is approached. The more accurate the displacements over the fracture surface near crack front, the more accurate will be the contact tractions. Therefore, for more accurate contact tractions, a more accurate FE displacement solution near the crack front is required. Although the accuracy of displacements increases by using quarter-point tetrahedrals at the crack front region, some inaccuracies may still remain, due to the low quality of the element resulting from the random size, shape and orientation of the elements in the crack front region [55]. A more structured mesh, in which significant variation of the size of quarterpoint tetrahedral elements along the crack front is avoided, is therefore expected to produce more accurate contact tractions. Fig. 7 also shows the variation of the analytical and numerical point-wise mixed-mode SIFs along the crack fronts of the penny-shaped crack under two different contact conditions. Analytical solutions for 3D penny-shaped and elliptical cracks embedded in infinite solids [74] are also plotted. Here, φ and ω are the polar angle of the circle, and the parametric angle of the ellipse, respectively. These results demonstrate the accuracy of the disk-shaped domain integral and displacement correlation to compute the SIFs from arbitrary meshes, even when crack surfaces are in contact. Fig. 8 shows the distribution of normalized augmented gap and contact traction over the slave surfaces of a slipping elliptical crack. The application of an inverse singular penalty parameter variation near the crack front has resulted in very accurate contact tractions. Fig. 9 shows the variation of the numerical point-wise mixed-mode SIFs along the crack fronts of the elliptical crack subjected to two different contact conditions. One important feature of Figs. 7 and 9 is the mode I crack deformation for the initially open cracks. It should be noted that open cracks are ubiquitous in the subsurface, even at great depths. The significant negative K I , which is due to the crack deformation before crack closure, can therefore considerably influence the growth behavior of cracks under compression. Fig. 10(a) shows the convergence of the contact tractions through the augmentation process. A very low penalty parameter (ϵ 0 = E/a) has been considered, where the role of the augmentation procedure in enforcing the contact constraints is significant. The errors in the contact traction and the SIFs drop considerably in the first two augmentations, indicating that very efficient enforcement of the contact constraints is achieved by only two augmentations. More than four augmentations enforce the contact constrains almost exactly, by strictly penalizing any penetration on the contact surfaces. It is evident that a larger value of the penalty parameter yields a faster convergence of the contact tractions to the exact values. In order to make the augmentation procedure more efficient, the penalty parameter has to be chosen as large as possible, without making the system of equations ill-conditioned. The contact precision for heavily fractured media can attain very high values, and therefore a lower bound of the penalty parameter, for which no ill-conditioned behavior occurs in the system, must be defined irrespective of the contact precision. A penalty parameter based on the ratio of Young's modulus to the average crack surface element size introduces penalty terms into the system that are comparable in magnitude to the values of the members in stiffness matrix corresponding to the nodes near the contact surfaces. L n implies the average size of the elements at the crack front region. A penalty parameter defined as the ratio of Young's modulus to this average size (ϵ 0 = E/L n ) can therefore be recommended as the lower bound for the value of the penalty parameter for individual cracks. For a reasonably fine mesh, suitable for the crack problems, this proposed value for the penalty parameter yields less contact traction errors and faster convergence than the one shown in Fig. 10(a). Therefore, a penalty parameter ϵ 0 = E/L n defined individually for each fracture, together with two to three augmentations, enforces the contact constraints accurately and efficiently.
Consider a ray emanating from an arbitrary point on the crack front, lying on the surface of the penny-shaped crack, and extending in a direction normal to the crack front, as shown in Fig. 5(b). Fig. 10(b) presents the variation of the magnitude of the augmented gap, the penalty parameter, and the magnitude of the contact traction, after three augmentations along this ray, plotted against the normalized distance from the crack front. A singular square-root variation of the penalty parameter forces the contact traction to have finite values very close to the crack front. The magnitude of the gap vector tends to zero when the crack front is approached (lim r →0 ∥g * ∥ = 0), and a square-root singular penalty parameter ϵ = ϵ 0 /ϵ r = ϵ 0 / √ r/a makes it possible to obtain a non-zero limit for the magnitude of the contact traction (lim r →0 ∥t m ∥ = lim r →0 ϵ∥g * ∥ ̸ = 0). Such a variation of contact traction cannot be reproduced when a constant penalty parameter is considered. In fact, the contact traction approaches zero in the case of constant penalty, as a result of the influence of zero magnitude gap along the crack front.
In order to show the significance of the singular variation of the penalty parameter near the crack front, the contact solution obtained from the application of a uniform penalty parameter is shown in Fig. 11. Fig. 11(a) shows the variation of contact tractions over the fracture surface when a uniform penalty of 2E/a is applied, and three augmentations are conducted. Fig. 11(b) presents the variation of the magnitude of the augmented gap, penalty parameter, and the magnitude of the contact traction along the ray shown in Fig. 5(b). Significant errors are introduced in the contact tractions near the crack front, due to the decay of the gaps at the crack front. With a uniform penalty parameter, a square-root decay of gaps near the crack front forces the contact tractions to follow the same variation, resulting in considerable errors in contact tractions. The significance of these errors is highlighted when one compares Figs. 10(b) and 11(b). These errors also significantly influence the accuracy of the SIFs obtained from the domain integral approach. A singular variation of the penalty parameter, however, compensates for the decay in the gaps near the crack front, and allows computation of very accurate contact tractions and stress intensity factors.
A mesh sensitivity analysis was conducted to demonstrate the accuracy of the proposed contact and stress intensity factor computation algorithms. Fig. 12(a) shows five crack surface meshes with different degrees of refinement, ranging from a coarse mesh (n = 1) to a fine one (n = 5). Fig. 12(b) and (c) present the variation of the contact traction error e c in slip and stick conditions against the average size of contact elements for different numbers of augmentations. According to these results, and those presented in Fig. 10(a), three conclusions can be reached: (i) the contact tractions converge with respect to the number of augmentations, and the rate of convergence is quadratic; (ii) for a constant number of augmentations, the contact tractions converge with respect to the mesh size L n ; (iii) the contact traction error approaches zero as the number of augmentations tends to infinity and the element size approaches zero (A → ∞, L n → 0). Fig. 12(d) and (e) also show the variation of SIF computation error e t against the normalized domain size R d /L n and normalized distance of correlation point from the crack front r m /L n , using the DI and DC methods, respectively. The SIF computation for fixed values of domain radius and correlation distance (R d = r m = L n ) are plotted against the average size of the elements (L n ), in Fig. 12(f). These results demonstrate that the SIF computation error obtained from both the DI and DC methods converges to zero with mesh refinement. The main features of these plots are as follows: (i) For coarse meshes, domain integral yields more accurate SIFs than does displacement correlation, provided that the domain integral is small enough with respect to the mesh size. (ii) For fine meshes, both methods compute very accurate SIFs, with an error in the range of e t ≈ 1%-2%.  the mesh size near the crack front is that the quarter-point element must lie entirely inside the singular-dominant zone. The size of the singular-dominant zone mainly depends on the characteristic crack length, ranging between a/10 and a/50 [83]. Therefore, a suitable mesh size for a penny-shaped crack requires L n < a/10 for which r m = R d = L n is recommended for computing accurate SIFs using the DC and DI methods. Similar values have been suggested for r m and R d in the case of cracks under tensile loadings [55,54]. It is noteworthy that for the purpose of accurate SIF computation, a local refinement at the crack front region suffices. Therefore, coarse meshes can be used for the regions far from the crack front, to avoid unnecessary computational cost (see mesh structures in Fig. 5(b)). Fig. 11. The results of the contact tractions when a uniform penalty parameter is applied over the crack surfaces (ϵ r = 0.5). (a) The distribution of normalized contact traction over the slave surface of the penny-shaped crack shown in Fig. 5(b) in slip condition. (b) The variation of the master contact tractions, t m , augmented gaps, g * and penalty parameter along a radial ray emanating from the crack front of the penny-shaped crack shown in Fig. 5(b). Geometrical details: a/w = 0.1, β = 45 • ; Mesh details: L n ≈ a/20; Contact details: ϵ 0 = E/a; Stick condition: p 0 /σ = 0, τ c /σ = 0, µ = 1.2, Slip condition: p 0 /σ = 0, τ c /σ = 0, µ = 0.2. Results in (b) are obtained after three augmentations.

Two penny-shaped cracks
Consider two cubes, one with a pair of interacting penny-shaped cracks ( Fig. 13(a)), another one with a pair of intersecting penny-shaped cracks (Fig. 13(c)). Each cube is subjected to uniform compression in the X 2 direction over the top and bottom surfaces. The geometrical details of the cracks in both configurations are given in Fig. 13. Both cracks lie in the plane X 2 = X 1 cot β which generates the angle β with the direction of applied load. Fig. 13(b) and (d) show the finite element mesh of these crack configurations. The same boundary conditions as in the single crack problem are applied. Fig. 14 shows the distribution of contact traction over the master surfaces of these penny-shaped cracks in the slip condition. The main features are as follows: (i) Fig. 14(a) shows how the singular field of one crack can influence the contact tractions enforced on the surfaces of another crack, when two close cracks are slipping next to each other. This interaction cannot be captured well, unless the singular stress state along the crack front is modeled accurately. (ii) Fig. 14(b) demonstrates the accuracy of the proposed contact algorithm in enforcing the contact constraints, even for very complex configurations involving crack intersections. The main difficulty is dealing with the corner singularity at the points where crack 2 intersects the surfaces of crack 1. (iii) Some limited inaccuracies in contact traction are visible near the crack fronts in both configurations. These inaccuracies result from the low quality of some elements, due the random placement of quarter-point tetrahedrals near the crack front [55]. Nevertheless, these slight inaccuracies are limited to a few quarter-point nodes, and the contact tractions are obtained with a high level of accuracy elsewhere. Overall, these results demonstrate the applicability of the proposed methodology to the analysis of very complex contact configurations. Fig. 15 shows the variation of the point-wise mixed-mode SIFs along the crack fronts of the penny-shaped cracks of the interacting and intersecting configurations shown in Fig. 13. The results are obtained from the displacement correlation and domain integral methods. Here, φ is the polar angle of the circle, as shown in Fig. 13. A few characteristics of the SIF variation along the crack front must be noted: (i) The results from displacement correlation and domain integral methods are in good agreement everywhere, except very close to the corner points, i.e., except near φ = −90 • and φ = 90 • in Fig. 15(d). The reason for this discrepancy is the complex stress state that exists near these points. In fact, the stress intensity factor loses its meaning at the exact corner points, since the order of singularity at this point is different from the order of the singularity for a crack; the reader is referred to [54] for more details. (ii) The interaction of the singular fields of one crack with the singular field or enforced contact conditions in the other crack, significantly influences the SIF variation along the crack front. For example, one can see these considerable interactions about φ = 0 • in Fig. 15(b), due to the interaction of the singular field of crack 1 with boundary conditions over surfaces of crack 2, and near φ = −90 • and φ = 90 • in Fig. 15  Overall, these results demonstrate the applicability of both disk-shaped domain integral, and displacement correlation methods, to accurately compute the SIFs of complex crack configurations.

Multiple planar cracks
Figs. 16 and 17 present the results of the contact simulation of two networks of fractures, where solid cubes are filled with randomly oriented, randomly placed, penny-shaped cracks. For both cases, the cube is subjected to a uniform compression in the X 2 direction over the top surface, and the same boundary conditions as in Section 5.2 are applied. This boundary condition requires that no crack intersects with the bottom face of the cube, where a Dirichlet boundary condition is applied. Fig. 16(a) shows the finite element discretization of a network of twenty four cracks with a power law size distribution, while Fig. 17(a) presents the FE mesh of a network of seventy five cracks   14. (a) The distribution of contact tractions over the master surfaces of two interacting cracks described in Fig. 13(a); Contact details: ϵ 0 = E/L n with three augmentations, crack 1: p 0 /σ = 0, τ c /σ = 0, µ = 0.2. crack 2: p 0 /σ = 0, τ c /σ = 0, µ = 0.4; (b) The distribution of contact tractions over the master surface for two intersecting fractures described in Fig. 13(c); Contact details for both fractures: ϵ 0 = E/L n with three augmentations, p 0 /σ = 0, τ c /σ = 0, µ = 0.2. of the same size. Figs. 16(b) and 17(b) show the distribution of normalized stick contact tractions over the slave surfaces, where the average contact traction error for both cases remains less than e c = 0.0003. This illustrates the ability of the proposed contact algorithm to enforce accurate contact constraints with small penalty parameters in high   contact precision problems. In fact, assigning an individual penalty parameter based on the local mesh refinement of each fracture avoids the system to become ill-conditioned, yet performing augmentation ensures high accuracy enforcement of the contact tractions. Figs. 16(c) and 17(c) also show the contact tractions over the slave surfaces when contact conditions τ c = 0.1, µ = 0.4 and τ c = 0, µ = 0.6 are applied for networks with power law and uniform size distributions, respectively. For these cases, the distributions of normalized tangential gap (slip) over the crack surfaces are presented in Figs. 16(d) and 17(d). The ratio of external load to the Young's modulus (σ/E) is very small in the context of infinitesimal strain theory. Therefore, according to the normalized slip values in Figs. 16(d) and 17(d), the values of tangential slip remain very small compared to the size of the cracks. This indicates the applicability of isoparametric integration point-to-integration point contact discretizations in geometrically linear applications, such as the linear elastic simulation of fractured media.

Conclusions
A tetrahedral-based finite element formulation is presented for the treatment of contact between fracture surfaces in high density fractured media. In this framework, the application of a singular square-root penalty parameter near the crack front ensures the accurate enforcement of contact constraint close to the crack front. The introduced gap-based AL approach also circumvents the difficulty of not being able to define contact tractions at the nodes located along the crack front. The proposed contact algorithm is able to enforce the contact constraint accurately over the crack surfaces of heavily fractured media, even with small penalty parameters. Results from numerical experiments on a cube containing single and multiple cracks indicate the accuracy and reliability of the FE framework introduced in this paper. Accurately computed SIFs, obtained using the displacement correlation and disk-shaped domain integral methods, indicate the accuracy of these methods using unstructured meshes. The presented method is compatible with coupled flow-mechanics and discrete fracture growth approaches that can be readily applied to investigate the effects of friction on fracture and fault permeability and reactivation.