Numerical modeling of wing crack propagation accounting for fracture contact mechanics

As a consequence of shearing, wing cracks can emerge from pre-existing fractures. The process involves the interaction of sliding of the existing fracture surfaces and the tensile material failure that creates wing cracks. This work devises a numerical model to investigate how wing cracks emerge, propagate and connect pre-existing fractures under shear processes. A mathematical and numerical model for wing crack propagation based on linear elastic fracture mechanics that also accounts for fracture contact mechanics is presented. Computational efficiency is ensured by an adaptive remeshing technique. The numerical model is verified and validated through a comparison of the analytical and experimental results. Additional numerical examples illustrate the performance of the method for complex test cases where wing-cracks develop for multiple pre-existing and interacting fractures.


Introduction
Wing cracks can develop from a pre-existing fracture when the fracture is subjected to shear processes. This occurs for many applications where fractured media are subjected to anisotropic stress regimes. For example, in fractured subsurface systems, fractures will slip if shear forces overcome the cohesion and frictional strength of the contact between the fracture surfaces. This can occur due to natural changes in tectonic stresses, but the process can also be induced by fluid injection, such as in situations of geothermal reservoirs. In the latter case, elevated pressures reduce the effective normal stress on the fracture, ultimately causing slip if the reduction in the normal stress is sufficient for the shear forces to overcome the cohesion and frictional resistance of the fracture. The slip of the fracture surfaces in opposite directions can cause the fracture to propagate in the form of wing cracks, possibly creating enhanced reservoir connectivity (Cheng et al., 2019;Jung, 2013;McClure and Horne, 2014;Norbeck et al., 2018) .Understanding this mechanism is, thus, crucial in the simulation of fractured subsurface formations. 2 Many experimental studies have been published that consider the formation, growth and connection of wing cracks caused by external compressive loading in specimens made of rock or rock-like materials (Haeri et al., 2014a(Haeri et al., , 2014bHorii and Nemat-Nasser, 1985;Ingraffea and Manu, 1980). In these experiments, if the pre-existing fracture is not perpendicular to the external load, wing cracks emerge at the tip and tend to align with the direction of the maximum compressive stress. The same conclusion is drawn from mathematical modeling. Based on the finite element method (FEM), Ingraffea and Heuze (Ingraffea and Heuze, 1980) predicted the propagation of wing cracks in rock structures by using all three stress, energy and strain criteria. Primary crack trajectories predicted by the stress and energy criteria are in good agreement with the observed trajectories. Based on the phase-field model (Bryant and Sun, 2018) and a modified phase-field model (Zhang et al., 2017), wing crack propagation was modeled using energy criteria that divided the active energy density into distinct parts corresponding to different crack modes (mode I and mode II). Sharafisafa and Nazem (Sharafisafa and Nazem, 2014) used the vector level set with both the discrete element method and the extended finite element method (XFEM) to model the wing crack propagation and coalescence in fractured rock masses. It was recognized that using the FEM with stress criteria (the maximum tangential stress -MTS) seems to be the best choice for wing crack modeling because of the simplification and agreement with the observed trajectories (Gonçalves da Silva and Einstein, 2013; Ingraffea and Heuze, 1980). While wing cracks develop as tensile fractures, the pre-existing fractures that wing cracks emerge from may be either open or in contact. This necessitates the inclusion of fracture contact mechanics in the wing crack models (Oden and Pires, 1983;Ueber et al., 2008). There are several ways to formulate the contact mechanics corresponding to the different types of discretization. For example, based on the semianalytical displacement discontinuity method, Kamali and Ghassemi (Kamali and Ghassemi, 2018) developed a simulation model in which the closed natural fractures were represented by so-called contact displacement discontinuity elements (Asgian, 1988), approximating the contact mechanics condition. However, the approach has limitations in dealing with the interaction between multiple fractures due to inherent limitations of the semianalytical displacement discontinuity method.
In FEM, equilibrium of linear elasticity is described as a displacement field minimizes the potential energy. Then, the contact should be considered as an inequality constraint of the optimization formulation. This means that the potential energy is minimized while satisfying a contact constraint assumed to be a nonpenetration condition between the surfaces of the 3 fracture. The inequality constraint can be solved by some methods, such as the active-set, Frank-Wolfe, penalty or barrier methods (Hüeber and Wohlmuth, 2005;Ueber et al., 2008).
An inherent problem in the simulation of fracture propagation is the disparate length scales. While the simulation domain can be quite large, the fracturing processes occur on a scale that is several orders of magnitude smaller. Moreover, most numerical methods for fracture propagation are dependent on resolving the fracture in a grid; however, the fracture path is not known a priori. A possible remedy for both of these issues is to apply adaptive remeshing (ARM) techniques to refine and adjust the mesh around an advancing fracture path.
This paper presents a mathematical model and corresponding numerical solution approach to simulate the development of wing cracks while accounting for fracture contact mechanics. First, in Section 2, the mathematical model for wing crack propagation is formulated based on the linear elasticity theory, in combination with the criteria for fracture propagation. Fracture surfaces are allowed to be in contact or fully open, modeled by contact mechanics. Section 3 presents the numerical solution approach. The governing equations are discretized using a finite element method with collapsed quarter-point elements at the fracture tips. This is combined with an adaptive remeshing technique based on error estimates and Laplacian smoothing. The contact mechanics are implemented by using an active set method.
Section 4 presents several numerical test cases. The obtained results are compared with both the analytical and experimental data to verify, validate and show the accuracy of the proposed model and procedure. Finally, more complex test cases where wing cracks develop for multiple pre-existing and interacting fractures show the capability of the proposed approach in modeling the development of wing cracks under shear processes.

Governing equations
A mathematical model for wing crack propagation based on linear elastic fracture mechanics is presented in the following section. Emphasis is placed on the conditions on the boundaries of existing and newly formed fracture paths. We also describe the criterion used to decide when, where and how far a fracture will propagate.

Elasticity
Consider a domain 2


with an outward unit normal vector n on its boundary and a pre-existing fracture with boundaries denoted by C   as shown in Fig. 1 where ( ) ( ) ( ) = fx σ x n x is the contact force, which vanishes in the case of an open fracture.
A nonpenetration condition is enforced in the normal direction of the fracture segments. This condition is governed by the Karush-Kuhn-Tucker (KKT) condition for the normal displacement jump and the contact force, which reads in which CC : +−   →  is a mapping that projects a point from side C For the tangential direction of C  , two types of conditions are considered: Either the fracture surfaces are modeled as frictionless, or the displacement jump in the tangential direction is specified. For frictionless fracture surfaces, the tangential traction is zero at C  ; i.e.
where ( ) tx is the tangential vector initiating from x at side C +  . The assumption of zero friction leads to an exaggeration of the slip but is acceptable herein, as the trajectory of the fracture is the primary quantity of interest. For the friction-free case, the deformation of the elastic medium and the fractures contained within are driven by the body force, external boundary conditions on D  and N  , or displacements on other fractures.
The second type of condition is a specified displacement jump in the tangential direction of the fracture, i.e.
where the total slip at C  , 0 u , is considered as known. This type of condition is relevant to mimic the slip along an existing fracture, which in applications, may be triggered by effects not considered in the present model.
The wing cracks emerging due to the shear force on existing fractures are tensile cracks (Bobet and Einstein, 1998;Wong and Einstein, 2009). This means that their surfaces are in not 6 contact and both normal and tangential tractions at the corresponding fracture faces are zero; i.e. 0  =  = σ n σ t (8) The wing cracks are not present in the computational domain at the start of the simulations. Indeed, the computation of the point of failure and the paths of the wing crack that develop are the main challenges that are addressed in this work.

Failure and propagation
The wing crack growth processes are governed by a fracture criterion. From the mathematical model for elastic deformation, the stress at an arbitrary point can be directly calculated for a certain problem. In this work, we chose to adapt the fracture criterion based on the maximum tangential stress (MTS) (Erdogan and Sih, 1963), which is simple and sufficiently accurate (Gonçalves da Silva and Einstein, 2013; Ingraffea and Heuze, 1980), to predict the initiation and propagation angle of wing cracks. This criterion states that a crack grows when the maximum average tangential stress in the fracture process zone ahead of the crack tip reaches a critical value. Moreover, the crack growth direction coincides with the direction of the maximum average tangential stress along a constant radius around the crack tip. In polar coordinates ( ) , r  with the origin at the crack tip, the tangential stress has the following form (Erdogan and Sih, 1963) When the wing crack emerges, i.e., the criterion shown in Eq. (10) is satisfied, the increment of each fracture needs to be determined. For a single crack propagation, the increment is defined by a fixed distance such as the crack tip rosette radius h. In the case where more than one crack grows simultaneously, the tips with the highest energy in the fracture set advance significantly further than the others (Paluszny and Matthäi, 2009). The increment for each tip is defined by the Paris-type law (Paris and Erdogan, 1963;Renshaw and Pollard, 1994) where adv i L and Gi are the propagation length and the energy release rate for the i th propagation crack, respectively, Lmax is the maximum length increase at any propagation step, and the exponent  is a numerical parameter, which is set to 0.35 in this work (Renshaw and Pollard, 1994). For a general fracture in a two-dimensional domain, the energy release around the fracture tip is given by Here, I k and II k are the local mode I and mode II stress intensity factors at the tip obtained by summing the normal and shear stresses (Anderson, 2017), respectively 0 0 0 0  I  I  II   33  11  2  3cos  cos  3sin  3sin  4 2 2 4 2 2

Discretization
This section presents the finite element discretization of the governing equations presented in Section 2, together with an adaptive remeshing technique. The propagation of wing cracks is complicated, and their trajectories are difficult to achieve by analytical or semianalytical approaches, particularly when multiple fractures interact. In this case, numerical solutions by means of the finite element method (FEM) are a common approach. The finite 8 element formulation is based on the weak formulation established from the governing equation and states: Find where T V and S V are the test space and the solution space satisfying the inhomogeneous Dirichlet boundary conditions (so-called essential boundary conditions), respectivelly. T V and S V are defined by where 1 () H  is the Sobolev space of functions that are square integrable and have a square integrable first derivative. In Eq. (17) L is the differential operator, and D is the material matrix modified from C and defined by , for plane stress (20) , for plane strain (21)

Deformation and contact mechanics
The approximate solution of Eq. (1) In this work, e  are chosen as triangular elements.
The numerical approximation employed on the grid differs between the interior elements that are close to the fracture tip and the interior elements that are not. On elements 9 that are not connected to the crack tip, the displacement field is approximated as a quadratic function, which is expressed in terms of the displaced values, h d , at the three vertices and the midpoints of the three edges such that where Ni are the shape functions of a 6-node triangular plane isoparametric element defined by Eq. A(1).
To represent the stress singularity at the fracture tip, quarter-point elements (QPE) (Barsoum, 1977) are employed. Each element around the crack tip, as shown in Fig. 2 (a), is mapped by an 8-node plane isoparametric quadrilateral element, as shown in Fig. 2 where Ni are the shape functions defined by Eq. A(2). Then, the displacement field is approximated through the displacements at 6 nodes, h d , as a quadratic function such that By using the approximation given in Eq. (25), the numerical stress is singular at the crack tip, similar to the analytical formula shown in Eq. (9). More details are shown in Appendix A. By substituting Eqs. (23) and (25) into Eq. (17), the discretized system can be written where K and Fb are the global stiffness matrix and global body load vector, respectively, and are obtained by the assembly of the stiffness matrix and body load vector of each element (Ke and Fe) that are expressed as where B is the gradient matrix defined as The contact mechanics relation defined in Eq.
(3) introduces a nonlinearity in the system in that the boundary condition on the fracture depends on whether the fracture is in contact or not.
This nonlinearity is treated at each pair of contact points by the active set strategy (Hüeber and Wohlmuth, 2005). The details of the active set algorithm are shown in Fig. 3 and explained as follows: (1) Set k = 1, initialize d as an initial solution, predict a set of possible contact points P V and assume the actual contact zone 1 where c is a positive constant depending on the material. If Eq. (31) is satisfied, either x . Therefore, the pair ( )

 
, xx should be considered as the contact points for the calculation in the next step.
(4) Check if the contact zone at step k, k , is the same as step k + 1, 1 k + . If yes then stop, else, the nonpenetration condition x at the contact points is counted for the system by using the Lagrangian multiplier, then go to step (2).

Fracture propagation
The modeling of wing crack propagation is based on two assumptions. First, the wing crack emerges from the tip of the fracture, and second, a crack stops growing whenever its tip reaches a domain boundary or another fracture. That is, we do not consider the fracture propagation that crosses other fractures.
As detailed in Section 2.2, the fracture computation is based on the stress intensity factor evaluation. In this work, we compute SIFs by using the nodal displacement correlation technique (Parks, 1974) in conjunction with QPE (Barsoum, 1976;Henshell and Shaw, 1975) that not only captures the singularity of the stresses but also considerably improves the displacement near the crack tip, resulting in a more accurate computation of the SIFs (Khoei et al., 2008). Through the displacement of the QPE around a crack tip, these SIFs can be calculated as (Chen and Kuang, 1992;Kuang and Chen, 1993) ( )  Fig. 4 (b). This process also entails that the grid geometry is updated in the vicinity of the crack, as detailed in the next subsection.

Adaptive remeshing
The accuracy of the numerical simulation depends on the quality of the mesh that is affected by the geometric discretization errors and the gradients of the solution within the individual elements. In this work, we use the adaptive mesh refinement to obtain a solution that satisfies a given mesh discretization error while minimizing the number of elements. The adaptive remeshing (ARM) process involves two techniques: first, mesh refinement based on the error estimator (Zienkiewicz and Zhu, 1987) is used to improve the accuracy of the numerical solution, and second, Laplacian smoothing (Buell and Bush, 1973;Field, 1988) is used to improve the quality of the mesh. ,

Error estimator and refinements
The error estimator is based on the comparison between the numerical stress computed directly from the computed displacement field and a recovered stress with higher regularity.
The numerical stress is directly computed by Eq. (1): The quadratic approximation for the displacement renders a numerical stress that is a piecewise linear function on the elements and discontinuous across the interelement boundaries. To recover a globally continuous stress, we first define a nodal stress i The refinement is then performed based on a calculated error estimator. The essence of this process is to balance the errors between the elements. This means that the elements in regions of high error are locally refined. This process is repeated until the desired accuracy is   problem, the element i  that needs to be refined ( Fig. 6 (a)) is divided into four subelements by the connection between the midpoints of the edges (Fig. 6 (b)). Three hanging nodes appear.
These nodes are removed by connecting it to an opposite vertex, as shown in Fig. 6 (c).
a) Element needing to be refined b) Element divided into four subelements c) Removal of hanging nodes Fig. 6. The refinement processes.

The mesh smoothing process
The mesh refinement process proposed above is local and therefore has a low implementation cost. However, the locality sometimes causes triangles with undesirable properties, such as overlapping elements. We improve the quality of the mesh by using the Laplacian smoothing process that is defined as follows: Let triangles , 1,..., A precaution is taken to guarantee that the new coordinate assigned to  x will define valid triangles. The new coordinate for  x is immediately used for all subsequent Laplacian smoothing of other coordinates.
The general algorithm for the fracture propagation simulation in conjunction with the adaptive remeshing and accounting for the fracture contact mechanics is presented in Fig. 7. The item "numerical solver" requires the solution of the contact mechanics problem, as shown in Fig. 3.

Method verification
To evaluate the new approach, we consider a benchmark problem with the propagation of an isolated crack in a medium that undergoes tensile or shear stress. For this problem, the performance of the second order finite element method and quarter point elements (FEM-QPE) with and without adaptive mesh refinement is compared to that of conventional finite elements. The performance is measured in terms of the accuracy of the strain energy and SIF computation under grid refinement.
The model problem is a thin rectangular plate (length L, width b, and thickness t) including a pre-existing edge fracture, which is subject to tensile (mode I) or shear (mode II) stress as illustrated in Fig. 8. To make a fair comparison, as shown in Fig. 9, the FEM and FEM-QPE use a unique mesh, while the variant of the latter method that includes ARM (FEM-ARM-QPE) uses a multi-size-mesh controlled by the error estimator. The strain energy is given by ( ) For the plane stress singularity problem, the rate of convergence of the numerical solution is bounded satisfying (Pin and Pian, 1973a) ( ) where u h and u exact denote the solution from the numerical method and the exact solution, and n and  are the spatial dimension of the domain and the singularity degree of solution near the point of singularity, respectively. In the current case, n = 2 and 12  = (by Eq. (9)); hence, the convergence of the strain energy is linear with h.
a) Edge-cracked plate subjected to a mode I loading condition b) Edge-cracked plate subjected to a mode II loading condition For the mode I study, the comparisons of the convergence of the strain energy between the three different methods are shown in Fig. 10 (a, b). A linear convergence rate for the strain energy can be observed, in accordance with Eq. (39) and the conclusions by previous published studies (A. Mirza and D. Olson, 1978;Pin and Pian, 1973b). However, the FEM-QPE is significantly more accurate than the FEM. The convergence rate of the FEM-ARM-QPE is better than that of the FEM, and its accuracy approaches that of the FEM-QPE if the mesh refinement is sufficiently good. The comparison with the analytical solution (Tada et al., 2000) for the stress intensity factor is shown in Fig. 10. (c). The QPEs considerably improve the solution near and ahead of the crack tip and result in a more accurate computation of the SIF.
The ARM technique reduces the computational cost while still ensures the accuracy of the computation of the strain energy and stress intensity factor. This is confirmed by Table 1, which shows the total number of degrees of freedom (DOF) and total number of elements for the three methods under grid refinement.  For the mode II study, the reference solutions are obtained by ANSYS for the strain energy and by ABAQUS (Treifi et al., 2008) for KII. As shown in Fig. 11, as in the case of tensile stress, the QPE improves the accuracy of the solutions close to the stress singularity, and the ARM technique preserves the accuracy with fewer degrees of freedom.

Model validation
To further validate the presented numerical model, FEM-ARM-QPE, the initiation and propagation of wing cracks from the ends of a pre-existing fracture under uniaxial compression loading are investigated. The test case focuses on the accuracy in the computation of the SIFs and fracture propagation paths for a case where both analytical (Atkinson et al., 1982) and experimental (Haeri et al., 2014b) data are available.
The computational domain is a disc-shaped rock specimen containing a central single pre-existing fracture, as shown in Fig. 12. Here, R and t denote the radius and thickness of the disc, and 2a is the length of the fracture. The fracture is inclined at an angle  to the vertical direction at the center of the specimen. The specimen is compressed by two line loads MPa·m 1/2 (Haeri et al., 2014b). On the existing fracture, a no-friction condition is assigned in the tangential direction. For this problem, the analytical solution for the SIFs is given by (Atkinson et al., 1982) Table 2 and Fig. 13. The computation of the SIFs by the FEM-ARM-QPE model is in good agreement with both the analytical solution (Atkinson et al., 1982) and the experimental data (Haeri et al., 2014b) for various inclination angles of the pre-existing crack. numerically. We do not consider this a concern for accuracy of the numerical method and note that the simulations consistently predict propagation towards the locations of the point loads, independent of the fracture rotation angle.    Fig. 16 and Fig. 17. The number of DOFs is approximately doubled at the end of the simulation. When the pre-existing fracture (2) experiences slip u0 = 0.012 mm in the tangential displacements, wing cracks emerge at its tips.
They form an angle of approximately 70 degrees to the main fracture. By increasing the slip until u0 = 0.147 mm, wing cracks from fracture (2) gradually turn in the direction perpendicular to the minimum principle stress and connect to fractures (1) and (3). By increasing the slip until u0 = 0.3105 mm, two wing cracks newly emerge from fracture (1) and (3) and propagate away.

Propagation of multiple fractures driven by the shearing boundary conditions
Finally, we consider a more complex case with multiple closed pre-existing fractures arbitrarily appearing in a specimen, as illustrated in Fig. 18. The proposed FEM-AMR-QPE technique is ideally suited to address the complexity of this problem in an efficient and accurate manner.  Step 1, DOF = 15268 0.0746 MPa

Conclusions
This work presented a numerical model for wing crack initiation and propagation due to shear slip. The governing mathematical model is based on linear elastic fracture mechanics and contact mechanics, along with failure and propagation criteria for multiple mixed-mode fracture propagation. The numerical solution approach is based on a combination of the finite element method combined with quarter point elements to handle the singularity at the fracture tips. The fracture contact mechanics are solved by using the active set strategy. In addition, an adaptive remeshing based on an error estimator and Laplacian smoothing for implementation is utilized for computational efficiency.  A (2) • The details of the QPE formulation.
Substituting specific coordinates of 6 nodes ( ) , ii xy as shown in Fig. 2. Definition of the elements around the crack tip (a) and an 8-node plane isoparametric element (b).