An efficient phase-field model of shear fractures using deviatoric stress split

We propose a phase-field model of shear fractures using the deviatoric stress decomposition. This choice allows us to use general three-dimensional Mohr–Coulomb’s failure function for formulating the relations and evaluating peak and residual stresses. We apply the model to a few benchmark problems of shear fracture and strain localization and report remarkable performance. Our model is able to capture conjugate failure modes under biaxial compression test and for the slope stability problem, a challenging task for most models of geomechanics.


Introduction
The shear failure of brittle materials in compression, also known as shear bands or localized strains, are one of the dominant modes of failure in geo-structures.It has recently emerged as an active research topic due to its interest in structural geology and engineering.The growing interest stems from its engineering applications in subsurface energy technologies, including enhanced geothermal energy systems where the hydro-shearing technique is aimed to reactivate and slide the preexisting fracture network to increase the rock mass permeability [1][2][3][4], large-scale CO 2 sequestration in deep saline aquifers [5][6][7], impoundment and level changes of artificial water reservoirs of hydropower plants [8][9][10][11] and underground natural gas storage facilities [13], where their mechanics are crucial to understanding the stability of faults and hence earthquake mechanisms [12,[14][15][16].Other engineering applications include fault and slope stability assessment [17,18], or the stability of faults during the groundwater injection and production operations [19,20].
The simulation of shear fracturing processes is a challenging task.The finite element method (FEM) has been the dominant numerical method for modeling solids and continua.Classically, two fundamentally different perspectives are proposed to study compressive fractures using FEM: -Discrete fracture models (DFM) that are based on the classical theory of Linear Elastic Fracture Mechanics (LEFM) founded by Griffith [21], Irwin [22].-Smeared fracture models (SFM) that are based on the classical theory of Continuum Damage Mechanics (CDM) proposed initially by Kachanov [23].Each class includes extensive literature dating back to the 1960s that is out of the scope of this text to cover comprehensively.Therefore we only point the interested reader to a few primary studies of each class.
Phase-field models have been extensively used for modeling brittle, cohesive, and ductile Mode-I fracture patterns, in elastic or poroelastic materials and homogeneous or heterogeneous domains [63][64][65][66][67][68][69][70][71][72][73][74][75][76][77][78] [see 79, for a detailed review].Although Lancioni and Royer-Carfagni proposed a simple extension for shear fractures, the applicability of phase-field for modeling shear failure remained virtually untouched until very recently [81][82][83].In a detailed study, Fei and Choo presented a phase-field formulation of frictional fracture based on Palmer and Rice theory and using a similar stress decomposition approach to the one proposed by Hu et al. for tensile cracks.The authors validated their model on a set of classical problems as well as various experimental setups [86].
In the present study, we propose a phase field model of shear failure that adapts the cohesive model of shear fractures proposed by Fei and Choo for deviatoric stress decomposition (DSD) instead of the proposed contact stress decomposition (CSD).Hence, we arrive at an alternative descriptor for the shear fracture orientation (i.e., the α tensor) which is solely based on the deviatoric strain.We adapt the crack driving force to be consistent with the DSD decomposition.The resulting formulation simplifies the damage criterion since it results in damaging the shear modulus.Lastly, the proposed model allows us to use the general forms of the failure functions from the classical plasticity theory and therefore is not limited to Mohr-Coulomb failure model.
In what follows, we first briefly describe the original framework based on CSD.We then discuss our generalization proposal.Lastly, we use both frameworks to model a set of benchmark problems.

Phase-field method
In this section, we first describe the general phase-field framework for modeling crack propagation in materials.We then summarize the most recent CSD shear model [83].Finally, we discuss our proposed adjustment for better stability.

Phase-field governing equations
Consider the continua Ω ∈ R D in D-dimensional space, depicted in Figure 1, with its boundary represented as Γ.The boundary Γ is subjected to Neumann boundary conditions on Γ t and Dirichlet boundary conditions on Γ u , where Γ u ∪ Γ t = Γ and Γ u ∩ Γ t = ∅.The set of discontinuities in the domain is represented by a discrete surface Γ d .
According to the phase-field formulation, the fracture's discrete surface Γ d is approximated implicitly as a continuous function over a width l using the Allen-Cahn fracture surface density function γ(d) as where d is the phase-field variable, with d = 0 presenting the intact part of the domain while d = 1 expressing a point on Γ d .w(d) is the transition function, also known as the dissipation function, defined for cohesive cracks as w(d) = d [87,88], hence c 0 = 8 3 .Accordingly, a surface integral ds is approximated using a volume integral as ds ≈ γ(d) dv.
Given the displacement field u, the small-deformation strain measure ε = (∇u+∇u T )/2, and the crack surface density function γ(d), the total energy of a fractured continua, occupying the domain Ω and bounded by the boundary Γ, shown in Figure 1, is expressed as where, Ψ external is the work done by the external traction stress τ and body force b, and expressed as The fracture energy, i.e., Ψ f racture , is the energy dissipated from the system to create a fracture surface Γ d .Given the energy release rate G c (per unit fracture length), Ψ f racture is expressed as The stored internal energy of the system Ψ internal consists of the elastic stored energy in the intact part of the domain and stored energy in the damaged part of the domain, expressed as The internal energy density function ψ(ε, d) is defined as ψ(ε, d) = 1 2 σ : ε, which consists of both inactive and damaged counterparts.For the intact part of continuum, i.e., where d = 0, the Cauchy stress tensor σ(ε, d = 0) is expressed using Hook's law as where, κ and µ are bulk and shear moduli of the intact material, respectively, and ε v is the volumetric strain, expressed as ε v = tr(ε).For the parts of the domain where d > 0, the Cauchy stress tensor is decomposed into inactive part σ I and active part σ A as The active part of the stress tensor undergoes the damage process, and g(d) is a degradation function that expresses the stress transition from bulk ( σ) to fracture ( σ).We will discuss these in more details in the next sections.Therefore, there are two solution variables associated with the phase-field formulation, the standard displacement field u and the additional phase-field variable d.Taking the variation of Ψ with respect to u and d, and following the standard weak to strong form steps of the FEM [89,90] and phase-field [73,83,85], we can arrive at the following governing relations: The irreversibility of the fracture process is guaranteed with the local history field of maximum stored shear energy H + (ε) that allows us to solve the constrained minimization of eq. ( 9) in a straightforward way [67] and avoids unphysical self-healing.H + (ε) is defined as follows: where t is time.Equation ( 9) is then rewritten as follows: Since Ḣ ≥ 0, non-negative ḋ is guaranteed and, consequently, the irreversibility of the fracture growth.We define H(ε) after describing the stress decomposition approach.
In this work, we use the Lorenz degradation function g(d) defined as [91,92]: where, M = G c /(c 0 l) and ψ c is the critical crack driving force at the material's peak strength, evaluated as ψ c = −M w (0)/g (0).The damage begins to accumulate as soon as elastic stored energy exceeds this critical threshold.

Stress decomposition. Introduction
The split of the strain energy density into crack driving and intact components defines the damage mode and fracture pattern.Up to date, two fundamental approaches are available.The approaches of the first class do not take into account the local fracture orientation, whereas the second approaches take into consideration the local crack orientation.The first group of models includes the isotropic model, the volumetric and deviatoric decomposition model, the spectral decomposition model, or the anisotropic models.The isotropic model proposed by Bourdin et al.where the entire strain energy density is degraded.The volumetric and deviatoric decomposition model proposed by Amor et al. splits the strain tensor into its volumetric and deviatoric components.This approach avoids crack inter-penetration in composites and masonry structures.The fracture is then assumed to be driven by volumetric expansion and deviatoric strains.The spectral decomposition model proposed by Miehe et al. splits the strain tensor into its principal components and only tensile components drive the fracture propagation.The anisotropic models are based on the spectral decomposition of the strain tensor using other projections, such as the eigenvalue and eigenvector of the effective stress tensor [94].
The second group of approaches take into consideration the local crack orientation.The directional model proposed by Steinke and Kaliske splits the stress tensor into the crack driving and persistent components using the fracture orientation.For each point, a fracture coordinate system is defined and the fracture orientation is obtained from the maximum principal stress direction.Strobl and Seelig and Strobl and Seelig computed the fracture orientation from the phase-field gradients.Following this way to compute the fracture direction, Liu et al. developed a phase field model based on micromechanical modeling, i.e., the macroscopic fracture is modeled as a collection of microscopic fractures.
In the following subsections, we describe the contact stress decomposition (CSD), used satisfactorily to simulate shear fractures under confining pressures, and lastly we present our proposal based on the deviatoric stress decomposition (DSD).Both models do not take into account the local fracture orientation.

Contact stress decomposition (CSD)
Since a compressive fracture behaves like a contact problem, Fei and Choo proposed a stress decomposition approach that is closely related to the contact formulation, which we refer here as CSD.It starts by considering a corotational coordinate system on the fracture surface with m and n as tangential and normal vectors to the crack surface, and m along the direction of sliding.Additionally, let us define α = (mn + nm)/2.
According to this approach and under the assumption that the fracture remains closed, i.e., no tensile fracture, the only stress component that should undergo damage is the shear stress, and other stress components remain inactive.The bulk shear stress can be expressed as where, Consider the contact shear stress as τ .Then, the inactive stress tensor is expressed as and the active stress tensor as Here, τ is the residual contact stress while the fracture is fully developed, i.e., d=1.
Remark 1 Given the Mohr-Coulomb's failure function as, with σn = n•σ •n as to normal stress on the fracture surface, and c and φ as cohesion and friction angle of the intact material, the peak and residual shear stresses are expressed as τp = c + σn tan(φ), τr = cr + σn tan(φr), where cr and φr are residual friction and cohesion at the fully developed failure state.
Remark 2 Based on the Mohr-Coulomb failure criterion, the critical plane for the failure is evaluated at two conjugate angles θ = ±(45 • − φr/2) [27] with respect to the direction of the maximum principal stress.However, the authors only consider θ = +(45 • − φr/2) [see 83, eq.56].This restriction is required otherwise m, n is not uniquely defined.

Our proposal: deviatoric stress decomposition (DSD)
The total strain tensor can be decomposed into volumetric and deviatoric parts, as ε = ε v 1 + e.We can also express the Cauchy tensor in terms of the mean confining stress p and the deviatoric stress tensor s as σ = −p1 + s.Therefore, we can re-write Hook's law for the intact part as [80,99] Given the equivalent deviatoric (Mises) stress q = ( 3 2 s : s) 1/2 and the equivalent deviatoric strain ε q = ( 2 3 e : e) 1/2 and with some algebra, we can write that Let us now define the Unit Deviator Tensor α q as α q = 2 3 e ε q , where α q = √ α q : α q = 1, ( Hook's law can therefore be expressed as qα q , where p = κε v , q = 3µε q .
Equivalent to the CSD, we can describe the compressive failure in a material as damage in the deviatoric stress component.Therefore, the compressive pressure becomes the inactive part of the stress tensor, i.e., and active stress is described as where the bulk deviatoric stress is q = 3µε q .
Remark 3 This deviatoric stress decomposition allows us to leverage the general form of virtually any failure surface that are described in the classical plasticity theory, including the Mohr-Coulomb failure function.Given the friction angle φ and cohesion coefficient c, the general form of the Mohr-Coulomb's failure criterion is expressed as Here, R M C defines the shape of the Mohr-Coulomb's failure surface and is expressed as where Θ is the Lodè angle, evaluated as cos(3Θ) = (r/q) 3 .The invariant r is the third invariant of the deviatoric stress tensor, and is defined as r = ( 9 2 tr(s 3 )) 1/3 .Based on this criterion, we can find the peak and residual Mises stresses as with φr and cr as the residual values for friction angle and cohesion at the fully damaged state.

Remark 4
We can easily replace the non-smooth Mohr-Coulomb surface R M C with some alternatives [100,101].In fact, we can potentially pick any alternative failure function available for different materials.

Crack driving force
Given τ = µε γ and τ p = p tan φ + c = µε p γ , the crack driving force relations for CSD is derived as [83] where and they showed that this model is consistent with Palmer and Rice model.Now, for the deviatoric stress decomposition discussed above, we can revise the crack driving force, given q = 3µε q and qp = (p tan φ + c)/R M C = 3µε p q , as More details on the derivation of H t and H slip for CSD approach are provided in Appendix A.

Boundary conditions
To have a complete mathematical description of the problem, we lastly need to describe the boundary conditions.Considering Figure 1, the boundary conditions are described as where ū and τ are prescribed displacement and traction forces, respectively.The steps used to solve the problem are detailed in Algorithm 1.

Applications to compressive strain localization
Here, we consider three reference problems of shear fractures, including direct shear test, biaxial compression test, and slope failure analysis.We show that our model can effectively capture multiple modes of failure concurrently.

Direct shear test
Our first example is the direct shear test.We simulate the propagation of a fracture in a long shear apparatus and we compare our results with analytical solutions and Fei and Choo's numerical simulations [83].The setup of the experiment is plotted in Figure 2. The domain is 500 mm long, 100 mm tall, and an initial 10-mm horizontal fracture is carved in the middle of the left boundary.The boundary conditions are: the bottom boundary is fixed, the top boundary is displaced horizontally, and the two lateral boundaries are fixed vertically.We neglect gravity.
The material properties are: shear modulus G = 10 MPa, Poisson's ratio ν = 0.3, cohesion strength c = 40 kPa, peak and residual friction angle φ = φ r = 15 • , shear fracture energy G c = 30 J/m 2 , and fracture's length-scale l = 2 mm.As in the previous works of Palmer and Rice and Fei and Choo, we impose the fracture propagation to be horizontal.Following Fei and Choo's simulations [83], we initialize vertical compressive normal stress to 149 kPa, which results in τ p = 80 kPa and τ r = 40 kPa.We mesh the domain near the fracture path with a mapped squared mesh of size l/4 = 0.5 mm and the remaining domain with a 1-mm free triangular mesh.
The horizontal force-displacement curve is shown in Figure 3.The agreement of the peak and residual forces provided by our numerical simulation is very satisfactory.Theoretically, the peak load, i.e., the peak shear stress times the width of the specimen, is 40 kN, and the output of our simulation is 40.387kN.In the same way, the theoretical residual load is 20 kN and the output of our simulation is 19.978 kN.We estimate the fracture energy from the forcedisplacement curve, the shaded area in Figure 3.The output of our model provides a fracture energy equal to 14.6914 J, while the theoretical value is 15 J. Therefore, we report a remarkable agreement between our simulations and expected theoretical values.We analyze the sensitivity of our model to the phase-field length parameter, l.We run several simulations of the direct shear test problem for several values of l, ranging from 1 mm to 10 mm.Results are depicted in Figure 4(a).The force-displacement curves for the four values of l confirm that the model is virtually insensitive to the phase-field length parameter.We check the mesh dependency of our model by running three problems of the long-shear apparatus problem.We fix the ratio length scale parameter to mesh size, l/h, to 20 and we run three simulations for three l-and h-values.Results are plot in Figure 4(b).The curves confirm that the model is insensible to the mesh size.We plot the phase-field distribution at three time steps in Figure 5.The peak load is given for U x = 0.8083 mm, after this value is reached the phasefield has already emerged and propagate along the whole fracture, Figure 5 (a).Afterward, the phase-field value intensifies during the softening stage, Figure 5 (b), up to the time the fracture is completely developed, Figure 5 (c).At this time, the domain is split into two parts.The upper part slips over the bottom one, and the shear stress between both parts is constant and equal to the residual shear stress, τ r = 40 kPa, resulting in a theoretical horizontal force of 20 kN.

Biaxial compression test
Our next example is a biaxial compression test.We simulate a laboratory-size specimen under plane strain, different confining pressures and with different residual friction angles.This example allows us to show the ability of the model to simulate the pressure dependence of the peak and residual strengths.We compare our numerical results with peak and residual strengths computed with a mechanical equilibrium model before and after the rupture.strength is reached.The peak strength is the same in both models since they have the same p c , c, and φ.Afterward the fracture propagates suddenly across the domain, reaching both lateral boundaries, and the vertical force suddenly sinks.Our numerical model is able to capture the fracture propagation during the transition from the peak to the residual strengths due to the adaptive time step.Moreover, the curves evidence that the phase-field model is able to simulate the residual strength, which depends on the confining pressure, the residual friction angle, and the fracture path.We run several simulations of the biaxial compression problem for several values of l, ranging from 1 mm to 10 mm.The force-displacement curves for the four values of l are included in Figure 7.As in the previous problem, the curves for the values of l confirm that the model is virtually insensitive to the phase-field length parameter.
The evolution of the phase-field variable for p c = 200 kPa, φ = 20 • , and φ r = 20 • , at three time steps is shown in Figure 8.The phase-field is almost zero when the peak strength is reached, Figure 8(a).In fact, due to the isotropic material model and homogeneous stress conditions of the biaxial test, two equally like fracture paths nucleate.This is consistent with the Mohr-Coulomb model.Nevertheless, only of the trajectories evolves and result in the final fracture pattern during the sudden decrease in the peak strength, Figure 8(b).Later, the phase-field variable increases its value along the fracture path up to the residual peak strength is reached, Figure 8(c).We simulate nine cases with several combinations of p c , φ, and φ r values.We also compute the peak and residual strengths applying mechanical equilibrium prior and after the fracture propagation.Given the fracture path, the mechanical equilibrium is illustrated in Figure 9.The total vertical force applied on the top boundary is F V , the total horizontal force on the left lateral boundary is F H , and the tangential and normal forces on the fracture path are T and N respectively.We suppose the nucleation and fracture propagation is instantaneous and the fracture path is a straight line.the angle between the fracture path and the vertical axis is θ.Then, at the onset of the fracture propagation, the tangential force on the fracture is: and once the fracture is fully developed, the tangential force on the fracture is: The mechanical equilibrium in the vertical direction is given by: and in the horizontal direction: where H is: Solving V from Eq. ( 37), substituting V in Eq. ( 38) and operating, the vertical force at the onset of the fracture propagation V p -peak strength-is: and the vertical force once the fracture is fully propagated V r -residual strength-is: We compute V p and V r for the nine simulated cases.The results are listed in Table 1.The agreement between both models is remarkable.   1 Triaxial experiment.We list the peak and residual strengths for nine cases in kN.Both strengths are computed with our phase-field model -denoted as Sim.-and using a mechanical equilibrium -denoted as M. Eq.-prior and after the fracture propagation.
< l a t e x i t s h a 1 _ b a s e 6 4 = " a F V R 5 g p V T G l q F a X 3 x w b S 7 < l a t e x i t s h a 1 _ b a s e 6 4 = " q / a y + D T v U G q U t X y p J / m O Z 1 l c S W I = " > A A A B 7 H i c b V B N S w M x E J 3 U r 1 q / q h 6 9 B I v g q e y K q M e C I B 4 r u G 2 h X U o 2 z b a h S X Z J s k J Z + h u 8 e F D E q z / I m / / G t N 2 D t j 4 Y e L w 3 w 8 g s e z P x P 6 + b 2 f g m z L l K M 8 s U X S y K M 4 F t g m e f 4 w H X j F o x c Y R Q z d 2 t m I 6 I J t S 6 f C o u B H / 5 5 V X S u q j 7 V / X L h 8 t a 4 6 6 I o w w n c A r n 4 M M 1 N O A e m h A A B Q 7 P 8 A p v S K E X 9 I 4 + F q 0 l V M w c w x + g z x 9 j b o 3 H < / l a t e x i t > 10m < l a t e x i t s h a 1 _ b a s e 6 4 = " G Q k g u q 5 I g y  < l a t e x i t s h a 1 _ b a s e 6 4 = " V j q u K M 0 N U M G F N g q 0 5 F N q j 5 B C F

Slope failure analysis
As the last example, we consider the problem of slope failure analysis reported in [33].Consider the soil slope shown in Figure 10.The domain is 20 m wide and 10 m tall, with a slope 1:1 on the left side.A 4 m wide rigid footing is placed on the crest of the slope.The slope is first subjected to a body force b = 20 kN/m 3 , and then these body-force stresses are used as the initial state for the footing loading step.Displacement at the bottom edge is fixed in both directions, while for the right edge, only horizontal displacement is fixed.As the main loading step, a displacement U y = 0.3 m is prescribed in the middle of a rigid foundation, which simulates the effect of a building imposing a stress on the slope.
The elastic parameters of the soil include E = 10 MPa and ν = 0.4.The initial friction angle and cohesion are φ = 16.7 • and c = 40 kPa, with φ r = 10 • and c r = 0 kPa as their respective residual values.The phase-field length scale parameter is set to l = 200 mm, and the domain is discretized using a free triangular mesh with mesh-size 20 mm.The resulting mesh roughly has 1M triangles and 500K vertices.The computational time takes about 12 hours in our desktop machine with i9-10900 processor with 10 cores and 20 threads.
Due to the relatively high cohesion and low friction angle, the shear-band formation for this problem is particularly interesting.If we plot the evolution of the Mohr-Coulomb's failure envelope right before the onset of fractures, as shown in Figure 11-(a), we observe that the failure should onset from both ends of the footing.This fact has also been reported by Haghighat and Pietruszczak [38], however, due to the pre-specification of only one orientation angle (θ), the crack formation from the left side was not captured by Fei and Choo [83].Therefore, to perform a comparison, we consider two cases: 8 2 u a E P w x 1 / + S 6 5 3 K / 5 + Z e 9 y r 3 R 0 n M c x S 9 b J J i k T n x y Q I 3 J G L k i V c H J P H s k z e X E e n C f n 1 X n 7 b i 0 4 + c w a + Q X n / Q v t 9 Z l R < / l a t e x i t > (a)Mohr-Coulomb Failure I. Shear band formation only from the right corner of the footing by suppressing the phase field variable to zero (d = 0) in the gray region (see Figure 10).II.Free shear-band formation, which results in two patterns from each side followed by coalesces.Additionally, we consider two critical fracture energies of G c = 10kJ/m 2 and G c = 5kJ/m 2 .The final fracture patterns of these two cases are shown in Figure 12.The evolution of phase-field variable for case I, with G c = 5kJ/m 2 , are plotted for different loading steps in Figure 13.The force-displacement response is plotted in Figure 10-(b).As the reader can find, the proposed formulation captures the peak and residual loads as well as the crack patterns accurately, and the results are consistent with those reported by Fei and Choo.The failure surface evaluated using phase-field method and the peak-load is wellaligned with potential failure surfaces and critical load F = 0.8MN/m resulting from limit-equilibrium analysis of the slope using the GeoStudio software (see Figure 11-b-d).
Lastly, we run a new set of simulations for case II.The results are plotted in Figure 14.As we find, here the model captures first a shear band formation from the left corner of the footing.This is in fact expected because of the stress-free surface of the slope creates a more critical failure condition on the left corner.The propagation of the mode, however, stops because it is directing to Mohr-Coulomb stable regions of the domain.Later, the main failure mode initiates and propagates from the right corner, and collides with the first mode somewhere underneath the footing, which is also consistent with the results of the limit state theories.A final branch is then generated and causes the ultimate failure of the slope.The pick stress, however, does not seem to be very different from those of Case I, as plotted in Figure 10-(b).
< l a t e x i t s h a 1 _ b a s e 6 4 = " 5 N 9 v t t r a q X l 6 0 8 R a H s 5 j Q A e m x p q E h k U y 3 0 8 k D I 3 x k l A 7 u R s p U C H i i / p 5 I i d R 6 K A P T K Q n 0 9 a w 3 F v / z m g l 0 L 9 o p D + M E W E i n i 7 q J w B D h c R q 4 w x W j I I a G E K q 4 u R X T P l G E g s k s b 0 J w Z 1 + e J 7 W T k n t W K t + W C 5 W r L I 4 c O k C H q I h c d I 4 q 6 A Z V k Y c o G q F n 9 I r e r C f r x X q 3 P q a t C 1 Y 2 s 4 f + w P r 8 A Y Q z l R 4 = < / l a t e x i t >

Concluding remarks
We presented a phase-field model of shear fractures using deviatoric stress decomposition (DSD).We validated the model by solving reference problems of shear fractures in geotechnical engineering.Our model has excellent performance.The main advantages of our phase-field approach are: (1) the model does not require re-meshing, (2) nucleation, propagation, and fracture path are automatically computed without the need to track fractures or pre-specify orientations, and (3) fracture joining and branching do not need additional algorithms.
For an isotropic Mohr-Coulomb material under homogeneous loading, it has been shown that there are two conjugate surfaces having the same likelihood for shear band formation.In fact, our model captures this for the biaxial compression problem without any intervention.This is the same for the slope stability problem, where our model was able to capture crack initiation from both corners of the foundation.While accurate in peak and residual force calculations, we found that the CSD model of shear fractures is more accurate in capturing such a transition.
The study was limited to modeling two-dimensional problems of compressive fracture.However, the proposed formulation is not limited to any dimensions.Therefore, we plan to explore three-dimensional models as a follow-up study.Additionally, pore-fluid consideration is critically important for modeling failure in geomaterials.This is also an area that will be considered next.Additional paths include incorporating rate-and-state friction models that are best suited for modeling geologic systems, and thermal coupling that is important for modeling geothermal systems.

A Crack driving force for deviatoric stress decomposition
Reminding that q, qp , qr denote the bulk and fractured deviatoric stresses at peak and residual stages, respectively, and q = 3µε q , with ε q and the deviatoric strain, the crack driving force during a plastic dissipation process as a result of frictional sliding can be expressed as ( Since q = g(d)q + (1 − g(d))q r , we will have, We observe that the relations are quite similar to those reported by Fei and Choo using shear stress split, except that shear stresses and strains are replaced now with deviatoric ones and therefore division by 3µ instead of µ.
Noting that total driving energy is expressed as H = H t + H slip , re-arranging eq. ( 43), we can write Now, one can substitute this relation into the phase field PDE eq. ( 9), and with 1D simplifications, integrate the phase field relation, as detailed in [83], to arrive at approximate relations for the evolution of deviatoric stress q as a function damage.Again, since the phase field PDE eq. ( 9) and driving force eq. ( 44) are very similar to those in [83], all the derivations hold identical and true for the deviatoric stress decomposition.Finally, by imposing length-scale independency to the deviatoric stress evolution, one obtains that This completes the derivation of crack driving force relations introduced in eqs.(30) and (31).

Figure 1
Figure 1 Domain Ω with boundary Γ, Dirichlet boundary Γu, and Neumann boundary Γt.The discontinuity surface is represented by Γ d with its phase-field diffused representation as γ(d).

Figure 2 Figure 3
Figure2Direct shear test setup.The domain is 500 mm long, 100 mm tall, and an initial 10-mm horizontal fracture is carved in the middle of the left boundary -red fracture-.The boundary conditions are: the bottom boundary is fixed, the top boundary is displaced horizontally, and the two lateral boundaries are fixed vertically.

lFigure 4
Figure4Force-displacement curves for the direct shear test with several phase-field length parameters, l, and mesh sizes, h.(a) Here the phase-field length parameter ranges from 1 mm to 10 mm and the mesh size is set to h = 0.2 mm.(b) Here the ratio phase-field length parameter to mesh size is set to 20.

Figure 5 Figure 6
Figure 5 The evolution of the phase-field variable at three time steps for the direct shear test.The imposed horizontal displacements, Ux, are: (a) 1 mm, (b) 2 mm, and (c) 3 mm.

Figure 7 Figure 8
Figure 7 Force-displacement curves with several phase-field length parameters, l.
w F Y w u p n 6 r U d U m s f y w Y w T 9 C M 6 k D z k j B o r t d L e E 7 k m b q 9 c c a v u D G S Z e D m p Q I 5 6 r /z V 7 c c s j V A a J q j W H c 9 N j J 9 R Z T g T O C l 1 U 4 0 J Z S M 6 w I 6 l k k a o / W x 2 7 o S c W K V P w l j Z k o b M 1 N 8 T G Y 2 0 H k e B 7 Y y o G e p F b y r + 5 3 V S E 1 7 5 G Z d J a l C y + a I w F c T E Z P o 7 6 X O F z I i x J Z Q p b m 8 l b E g V Z c Y m V L I h e I s v L 5 P m W d W 7 q J 7 f n 1 d q t 3 k c R T i C Y z g F D y 6 h B n d Q h w Y w G M E z v M K b k z g v z r v z M W8 t O P n M I f y B 8 / k D F w S O x g = = < / l a t e x i t > ux = 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " b Y 8 t / n 5 V 7 M X w p P P 1 4 8 Y Y R G 4 x e 6 M = "

r 9 8 3 < l a t e x i t s h a 1 _ b a s e 6 4 =
x p X n g 3 8 A w Z C 1 J e j 7 v c k r A U O 1 c 3 v U C 0 d F D a b 7 Y S / A Z L t o P L r B 7 U D I e X B 3 K 5 P b I u O y C P S o 8 C 5 w J yK N J V d q 5 L 7 c T 0 E g y H 6 g g W j c d O 4 R W T B R w K l i S d S P N Q k I H p M e a B v p E M t 2 K R 8 c k e N 8 w H d w N l H k + 4 B H 7 u y M m U q f 7 G q c k 0 N f T W k r + p z U j 6 J 6 2 Y u 6 H E T C f j g d 1 I 4 E h w G k y u M M V o y C G B h C q u N k V 0 z 5 R h I L J L 2 t C c K Z P n g W 1 Y s E 5 L p S uS / n y x S S O D N p F e + g A O e g E l d E l q q A q o u g R P a N X 9 G Y 9 W S / W u / U x t s 5 Z k 5 4 d 9 K e s z x 9 v 1 p s + < / l a t e x i t > b = 20 kN/m " 6 f r F

y v 8
I b G 6 A W 9 o 4 9 F a w H l M 8 f w B + j z B 0 2 m j T w = < / l a t e x i t > 4m < l a t e x i t s h a 1 _ b a s e 6 4 = " j 5 g 5N Y A I Y R h 4 J y d g t 8 1 V x 9 U h l Y Y = " > A A A B 6 n i c b V B N S 8 N A E J 2 t X 7 V + V T 1 6 W S y C p 5 J I U Y 9 F L x 4 r m r b Q h r L Z b t q l m 0 3 Y 3 Q g h 9 C d 4 8 a C I V 3 + R N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Q S K 4 N o 7 z j U p r 6 x u b W + X t y s 7 u 3 v 5 B 9 f C o r e N U U e b R W M S q G x D N B J f M M 9 w I 1 k 0 U I 1 E g W C e Y 3 M 7 8 z h N T m s f y 0 W Q J 8 y M y k j z k l B g r P X i D b F C t O X V n D r x K 3 I L U o E B r U P 3 q D 2 O a R k w a K o j W P d d J j J 8 T Z T g V b F r p p 5 o l h E 7 I i P U s l S R i 2 s / n p 0 7 x m V W G O I y V L W n w X P 0 9 k Z N I 6 y w K b G d E z F g v e z P x P 6 + X m v D a z 7 l M U s M k X S w K U 4 F N j G d / 4 y F X j B q R W U K o 4 v Z W T M d E E W p s O h U b g r v 8 8 i p p X 9 T d y 3 r j v l F r 3 h R x l O E E T u E c X L i C J t x B Cz y g M I J n e I U 3 J N A L e k c f i 9 Y S K m a O 4 Q / Q 5 w 9 G b I 3 O < / l a t e x i t >

e 3 I
X i / X / 5 L G n t l b 7 9 c O a 8 U q 0 e T O H J k k 2 y T E v H I A a m S U 1 I j d c L J P X k k z + T F e X C e n F f n b d w 6 5 U x m N s g P O O 9 f + 6 m Y x w = = < / l a t e x i t > (a)Slope Failure Model < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 4 s 6 I

Figure 10
Figure 10 Slope failure analysis.(a) The model setup for the slope failure analysis.The domain is 20 m wide and 10 m tall, with a slope 1:1.A 4 m wide rigid footing is located on the slope's crest, and is subjected to a vertical displacement Uy.The boundary conditions are: bottom edge is fixed in both directions, right edge is fixed horizontally and other faces are traction free.The grey region highlights the damage-inactive (Case I) and damage-active (case II) problems.(b) Force displacement results measured at the point of loading in the middle of the footing.The plots show absolute values.
< l a t e x i t s h a 1 _ b a s e 6 4 = " O y X E f p F o I F B D 5 p 1 h 9 L z Y O 0 5 O W r 4 = " > A A A B / X i c b V D J S g N B E O 2 J W 4 x b X G 5 e G o M Q L + O M B J e D E B T E i x L B L J A M o a f T S Z p 0 z w z d N W I c g r / i x Y M i X v 0 P b / 6 N n e W g 0 Q e r P e x 6 0 p a z K z i X 7 B + v g G F X C T u w = = < / l a t e x i t > (d)F = 0.9MN/m < l a t e x i t s h a 1 _ b a s e 6 4 = " w V 9 G D S 5 8 7 m v d o q i z m I Z Y c f c 7 z 4 s = " > A A A B / X i c b V D J S g N B E O 2 J W 4 z b u N y 8 N A Y h X s Y Z C Z q L E B T E i x L B L J A M o a f T S Z r 0 L H T X i H E I / o o X D 4 p 4 9 T + 8 + T d 2 k j l o 9 E H B 4 7 0 q q u p 5 k e A K b P v L y M z N L y w u Z Z d z K 6 t r 6 x v m 5 l Z N h b G k r E p D E c q G R x Q T P G B V 4 C B Y I 5 K M + J 5 g d W 9 w P v b r d 0 w q H g a 3 M I y Y 6 5 N e w L u c E t B S 2 9 w p 0 A N 8 c W p b p R a w e 0 i u r g / 9 U d v M 2 5 Y 9 A f 5 L n J T k U Y p K 2 / x s d U I a + y w A K o h S T c e O w E 2 I B E 4 F G + V a s W I R o Q P S Y 0 1 N A + I z 5 S a T 6 0 d 4 X y s d 3 A 2 l r g D w R P 0 5 k R B f q a H v 6 U 6 f Q F / N e m P x P 6 8 Z Q 7 f k J j y I Y m A B n S 7 q x g J D i M d R 4 A 6 X j I I Y a k K o 5 P p W T P t E E g o 6 s J w O w Z l 9 + S + p H V n O s V W 8 K e b L Z 2 k c W b S L 9 l A B O e g E l d E l q q A q o u g B P a E X 9 G o 8 G s / G m / E + b c 0 Y 6 c w 2 + g X j 4 x s S S 5 O 5 < / l a t e x i t > (c)F = 0.8MN/m < l a t e x i t s h a 1 _ b a s e 6 4 = " X 8 n 4 m x m X a w S f A R h i X J 6 U U a z j A 8 4 = " > A A A B / X i c b V D J S g N B E O 2 J W 4 z b u N y 8 N A Y h X s Y Z C c a L E B T E i x L B L J A M o a f T S Z r 0 L H T X i H E I / o o X D 4 p 4 9 T + 8 + T d 2 k j l o 9 E H B 4 7 0 q q u p 5 k e A K b P v L y M z N L y w u Z Z d z K 6 t r 6 x v m 5 l Z N h b G k r E p D E c q G R x Q T P G B V 4 C B Y I 5 K M + J 5 g d W 9 w P v b r d 0 w q H g a 3 M I y Y 6 5 N e w L u c E t B S 2 9 w p e A f 4 4 t S 2 S i 1 g 9 5 B c X R / 6 o 7 a Z t y 1 6 N V 4 N J 6 N N + N 9 2 p o x 0 p l t 9 A v G x z c P J p O 3 < / l a t e x i t > (b)F = 0.7MN/m < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 6 d H D / v l 7 G 9 P o s 3 v 0 X 2 7 6 F 5 K V k U = " > A A A C C H i c b V A 9 S w N B E N 2 L 3 / H r 1 N L C x S D E w n A n o p a i I D a C g j G B 5 A h 7 m 0 m y u H d 7 7 M 6 J 4 U h p 4 1 + x s V D E 1 p 9 g 5 7 9 x E 6 / Q 6 I O B x 3 s z z M w L E y k M e t 6 n U 5 i Y n J q e m Z 0 r z i 8 s

Figure 11 (
Figure 11 (a) Mohr-Coulomb failure function plotted right before the onset of localization.The stress around both corners of the footing appear as near-failure critical.(b-d) Mohr-Coulomb's critical surfaces using limit-equilibrium method at different load.F ≈ 0.8MN/m is found as the critical load, where factor of safety of the slope reaches 1.

2 Figure 12 Figure 13
Figure 12 Damage evolution for the case I slope stability analysis, where the left side is suppressed to have damage development.The plots show the evolution damage at Uy = 300mm for (a) Gc = 10kJ/m 2 and (b) Gc = 5kJ/m 2 .
r H f r Y 9 a a s 7 K Z f f Q H 1 u c P e s + V G A = = < / l a t e x i t > (a)U y = 100 mm < l a t e x i t s h a 1 _ b a s e 6 4 = " t D I j w H m d + M r i j u 4 a S u g a 6 a h D i 2 0

Figure 14
Figure14Damage evolution for the case II slope analysis (Gc = 2 ).In this case, the damage can initiate from either side and is primarily driven by the crack energy.Subplots a-d show the evolution of the damage parameter at different loading steps.