A Study on Phase-Field Models for Brittle Fracture
AA Study on Phase-Field Models for Brittle Fracture
Fei Zhang ∗ , Weizhang Huang † , Xianping Li ‡ , Shicheng Zhang § Abstract
In the phase-field modeling of brittle fracture, anisotropic constitutive assumptions for thedegradation of stored elastic energy due to fracture are crucial to preventing cracking in com-pression and obtaining physically sound numerical solutions. Three energy decomposition mod-els, the spectral decomposition, the volumetric-deviatoric split, and an improved volumetric-deviatoric split, and their effects on the performance of the phase-field modeling are studied.Meanwhile, anisotropic degradation of stiffness may lead to a small energy remaining on cracksurfaces, which violates crack boundary conditions and can cause unphysical crack openings andpropagation. A simple yet effective treatment for this is proposed: define a critically damagedzone with a threshold parameter and then degrade both the active and passive energies in thezone. A dynamic mesh adaptation finite element method is employed for the numerical solutionof the corresponding elasticity system. Four examples, including two benchmark ones, one withcomplex crack systems, and one based on an experimental setting, are considered. Numericalresults show that the spectral decomposition and improved volumetric-deviatoric split models,together with the improvement treatment of crack boundary conditions, can lead to crack prop-agation results that are comparable with the existing computational and experimental results.It is also shown that the numerical results are not very sensitive to the parameter defining thecritically damaged zone.
AMS 2010 Mathematics Subject Classification.
Key Words. brittle fracture, phase-field modeling, constitutive assumption, critically damagedzone, moving mesh, finite element method
In recent years, the phase-field model for brittle fracture based on the variational approach ofFrancfort and Marigo [7] has become a commonly used numerical simulation technique for engi-neering designs because it can handle complex cracks and crack initiation and propagation moreeasily than other methods. The basic idea of the phase-field modeling is to describe cracks by acontinuous scalar field variable d , which is used to indicate whether the material is damaged or not.This variable d depends on a parameter l describing the actual width of the smeared cracks and ∗ College of Petroleum Engineering, China University of Petroleum – Beijing, 18 Fuxue Road, Changping, Beijing102249, China ( fzhang [email protected] ) † Department of Mathematics, the University of Kansas, Lawrence, Kansas 66045, U.S.A. ( [email protected] ) ‡ Department of Mathematics and Statistics, University of Missouri – Kansas City, 5120 Rockhill Road, KansasCity, Missouri 64110, U.S.A. ( [email protected] ) § College of Petroleum Engineering, China University of Petroleum – Beijing, 18 Fuxue Road, Changping, Beijing102249, China ( [email protected] ) a r X i v : . [ m a t h . NA ] M a y as a value of 0 or close to 0 near the cracks and 1 away from the cracks. There are three majoradvantages of the phase-field modeling for brittle fracture over other methods. Firstly, the behaviorof the crack is completely determined by a coupled system of partial differential equations (PDEs)based on the energy functional. Therefore, additional calculations such as stress-intensity factorsare not required to determine the crack initiation and propagation. Secondly, complex fracturenetworks can be easily handled since crack merging and branching do not require explicitly keepingtrack of fracture interfaces. Thirdly, sharp but smooth interfaces (instead of discontinuities) canbe introduced into the displacement field.Since it was first proposed by Bourdin et al. [5, 7], the phase-field modeling for brittle fracturehas attracted considerable attention and significant progress has been made; e.g., see [1, 3, 4, 15,17, 19, 24]. However challenges still exist. In the phase-field modeling, constitutive assumptionsfor the degradation of energy due to fracture can be categorized into two groups, isotropic modelsand anisotropic models. In the former group, the degradation function acts on the whole storedbulk energy, which means that energy is released due to fracture in both tension and compression.Thus, crack propagation may also arise under compressive load state, which is physically unrealistic.On the other hand, in order to overcome this unphysical feature, the elastically stored energy isdecomposed into active and passive parts where only the former is degraded by the phase field.Two commonly used energy decomposition models have been proposed in the past. Miehe etal. [19] introduced a fully anisotropic constitutive model for the degradation of energy based onthe spectral decomposition of strains with the assumption that crack evolution is induced by thepositive principal strains. The other model is the unilateral contact model proposed by Amor etal. [1] that splits the strain into volumetric and deviatoric parts, with the expansive volumetricpart and the total deviatoric part being degraded. Since the choice of the energy splitting controlsthe energy contribution in the damage evolution, different splitting models can significantly affectnumerical approximations in the phase-field modeling of cracks. In this work we shall study thesemodels plus an improved version of the unilateral contact model.In the phase-field modeling, a pre-existing crack is used to initialize crack which can be modeledas a discrete crack in the geometry or an induced crack in the phase-field. For the former, it hasbeen successfully applied in phase-field models for single initial crack. However, it is difficult tohand complex crack boundary conditions since the location of the initial crack is mesh-dependent.For the latter, an initial strain-history field is introduced to define the location of the induced crack.One of the major advantages of this treatment is that initial cracks can be placed anywhere in thedomain without referring to the mesh, which makes it possible to deal with complex initial cracks.The induced crack model was first proposed by Borden et al. [3] and significant improvements havebeen made [3, 17, 21]. However, there is a small energy remaining in the totally damaged zonedue to the anisotropic degradation of stiffness. For the induced crack model, stress remains in theinterior of the initial crack and increases with external loads before the crack begins to propagate.This violates the vanishing stress condition on the crack surface and often results in unphysical crackpropagation. May et al. [17] have observed that with the induced crack setting for a single notchedshear test, numerical results are very different from those with discrete crack boundary conditions.To overcome this problem, Strol and Seelig [22, 23] proposed a novel treatment of crack boundaryconditions in which crack orientation is taken into account so that both the positive normal stresson the crack surfaces and the shear stress along the frictionless crack surface vanish. However,2he establishment of this constitutive assumption is not based on the variational approach, whichmakes the phase-field model more complicated to implement.Another important issue is that the phase-field modeling approximates the original discreteproblem as l → h (cid:28) l or at least h < l , where h denotes the meshspacing. A very fine mesh is needed to fulfill the condition if only a uniform mesh is used in thecomputation, which increases the computational cost significantly. Moreover, cracks can propagateunder continuous load. Thus, it is natural to use a dynamic mesh adaptation strategy to improvethe efficiency of the simulation. In this work we use the moving mesh PDE (MMPDE) method[6, 12, 13, 14] to concentrate mesh elements around evolving cracks. The reader is referred to Zhanget al. [25] for the study of the application of the MMPDE method to the phase-field modeling ofbrittle fracture.The objective of this paper is twofold. The first is to investigate two commonly used energydecomposition models, the spectral decomposition model [18] and the volumetric-deviatoric splitmodel [1], and their impacts on the performance of the phase-field modeling for brittle fracture.The former degrades the energy related to the tensile strain component while the latter degradesboth the expansive volumetric strain energy and the total deviatoric strain energy. In the lattercase, the compressive deviatoric strain energy is also accounted for contributing to the damageprocess, which may cause unphysical crack propagation. An improved volumetric-deviatoric modelis proposed and studied to avoid this difficulty.The second goal is to study the treatment of crack boundary conditions for the induced crackmodel. As mention before, the induced crack model leads to a small remaining stress in the damagedarea, which can cause unrealistic crack boundary conditions such as normal stress remaining onthe crack surface. To overcome this difficulty, we propose a simple yet effective treatment of crackboundary conditions: define a critically damaged zone with 0 < d < d cr , where d cr is a positiveparameter, and then degrade both the active and passive components of the energy in this zone.We shall consider four two-dimensional examples to verify the treatment. The first two are classicalbenchmark problems, single edge notched tension and shear tests. The third example is designed todemonstrate the ability of the models to handle complex cracks. The last example is also a singleedge notched shear test but with the physical parameters and domain geometry chosen basedon an experiment setting. Numerical results show that the spectral decomposition and improvedvolumetric-deviatoric models together with the improved treatment of crack boundary conditionslead to correct crack propagation for all of the examples.The present work is organized as follows. Section 2 is devoted to the description of the phase-fieldmodeling with three energy decomposition models and the improved treatment of crack boundaryconditions. A moving mesh finite element method based on the MMPDE moving mesh methodis described in Section 3 for the elasticity system. Numerical results for four two-dimensionalexamples are presented in Section 4. Finally, Section 5 contains conclusions. We consider small strain isotropic elasticity models with cracks. Let Ω be a two-dimensionalbounded domain filled with elastic material and having the boundary ∂ Ω = ∂ Ω t ∪ ∂ Ω u , where3he surface traction t is specified on ∂ Ω t and the displacement u is given on ∂ Ω u . Denote thedisplacement by u . Then the strain tensor is given by (cid:15) = 12 (cid:0) ∇ u + ( ∇ u ) T (cid:1) , where ∇ u is the displacement gradient tensor. For isotropic material, the strain energy density isgiven by Hooke’s law as Ψ e ( (cid:15) ) = λ (cid:15) )) + µ tr( (cid:15) ) , (1)where λ and µ are the Lam´e constants and tr( · ) denotes the trace of a tensor.For an elastic body with a given crack Γ, we take the approach of Francfort and Marigo [7] todefine the total energy as W ( (cid:15), Γ) = W e ( (cid:15), Γ) + W c (Γ) ≡ (cid:90) Ω \ Γ Ψ e ( (cid:15) ) d Ω + (cid:90) Γ g c dS, (2)where W e ( (cid:15), Γ) represents the energy stored in the bulk of the elastic body, W c (Γ) is the energyinput required to create the crack according to the Griffith criterion, and g c is the fracture energydensity that is the amount of energy needed to create a unit area of fracture surface.The model (2) treats the crack Γ as a discontinuous object, which has been a challenge in nu-merical simulation. Here we use the phase-filed approach with which Γ is smeared into a continuousobject. Specifically, a phase-field variable d ( x, t ) is used to represent Γ along with the parameter l that describes the width of the smooth approximation of the crack. The value of d ( x, t ) is zeroor close to zero near the crack and one away from the crack (see Fig. 1(a)). The fracture energy W c (Γ) is approximated by the smeared total fracture energy [5] as W lc = (cid:90) Ω g c l (cid:0) ( d − + 4 l |∇ d | (cid:1) d Ω . (3)The elastic energy needs to be modified to reflect the loss of material stiffness in the damagedzone. A degradation function g ( d ) is commonly used for this purpose, which is required to satisfythe following property g (0) = 0 : Damage occurred for d = 0 and this part should vanish; g (1) = 1 : No damage occurs for d = 1; g (cid:48) (0) = 0 : No more changes after the fully broken state; g (cid:48) (1) (cid:54) = 0 : The damage has to be initiated at the onset.In our computation, we use a common choice g ( d ) = d . A straightforward attempt for themodification of the elastic energy is to degrade the total stored energy in the damaged region, i.e.,Ψ le ( (cid:15), d ) = g ( d )Ψ e ( (cid:15) ) . (4)Unfortunately, this can result in crack opening in compressed regions, which is physically unrealistic.To avoid this, it is common to decompose the total stored energy into two components asΨ e ( (cid:15) ) = Ψ e,act ( (cid:15) ) + Ψ e,pas ( (cid:15) ) , (5)4 < 1( ψ + ( ε ) degraded , ψ - ( ε ) kept no change )d = 1 l (a) Regularized crack Γ l ( d ) by the phase-field approx-imation Critical damaged region : d < d cr ( ψ + ( ε ) and ψ - ( ε ) are both degraded ) d = 1Damaged region : d cr < d < 1( ψ + ( ε ) degraded , ψ - ( ε ) kept no change ) (b) The critically damaged region indicated by thecritical value d cr Figure 1: A sketch of the phase-field modeling of brittle crack.where Ψ e,act ( (cid:15) ) and Ψ e,pas ( (cid:15) ) represent the active and passive components, respectively, with onlythe former contributing to the damage process that results in fracture. (Two commonly useddecomposition models will be discussed in the next subsection.) Then, the damage model reads asΨ le ( (cid:15), d ) = g ( d )Ψ e,act ( (cid:15) ) + Ψ e,pas ( (cid:15) ) , (6)and the corresponding total energy is W l = W le + W lc = (cid:90) Ω (cid:16) ( d + k l )Ψ e,act ( (cid:15) ) + Ψ e,pas ( (cid:15) ) + g c l (cid:0) ( d − + 4 l |∇ d | (cid:1)(cid:17) d Ω , (7)where k l (cid:28) l is a regularization parameter added to avoid degeneracy. It is emphasized that withthis setting, only the active component of the elastic energy is degraded in damaged regions. Letting W = ( d + k l )Ψ e,act ( (cid:15) ) + Ψ e,pas ( (cid:15) ) + g c l (cid:0) ( d − + 4 l |∇ d | (cid:1) , we obtain the variation of the energy as δ W l = (cid:90) Ω ∂W∂d δd d Ω + (cid:90) Ω ∂W∂ ∇ d · ∇ δd d Ω + (cid:90) Ω ∂W∂(cid:15) : (cid:15) ( δu ) d Ω , where A : B is the inner product of tensors A and B , i.e., A : B = (cid:80) i,j A i,j B i,j . Define the functionspaces V u = (cid:8) ϕ | ϕ ∈ H (Ω) , ϕ = u on ∂ Ω u (cid:9) ,V u = (cid:8) ϕ | ϕ ∈ H (Ω) , ϕ = 0 on ∂ Ω u (cid:9) , where H (Ω) is a Sobolev space defined as H (Ω) = (cid:26) ϕ | (cid:90) Ω ϕ d Ω < + ∞ , (cid:90) Ω |∇ ϕ | d Ω < + ∞ (cid:27) . V d = H (Ω). The weak formulation for the phase-field model is to find d ∈ V d and u ∈ V u such that (cid:90) Ω (cid:18)(cid:18) d H + g c ( d − l (cid:19) δd + 2 g c l ∇ d · ∇ δd (cid:19) d Ω = 0 , ∀ δd ∈ V d (8) (cid:90) Ω σ : (cid:15) ( δu ) d Ω = (cid:90) Ω t t · δu dS + (cid:90) Ω f · δu d Ω , ∀ δu ∈ V u (9)where H = Ψ e,act ( (cid:15) ) and σ is the Cauchy stress defined as σ ≡ ∂W∂(cid:15) = ( d + k l ) ∂ Ψ e,act ∂(cid:15) + ∂ Ψ e,pas ∂(cid:15) . (10)To ensure that cracks can only grow (i.e., crack irreversibility), we replace H = Ψ e,act ( (cid:15) ) in (8) by H = max s ≤ t Ψ e,act ( (cid:15) )( s ) , (11)where t is the quasi-time corresponding to the load increments. In this subsection we discuss two commonly used models and an improved one for energy decom-position.
We first consider a commonly-used model proposed by Miehe et al. [18] based on the spectraldecomposition of the strain tensor. To this end, we define the positive-negative decomposition of ascalar function f as f = f + + f − , f + = f + | f | , f − = f − | f | . For the strain tensor, we have (cid:15) = (cid:15) + + (cid:15) − , (cid:15) + = Q diag( λ +1 , ..., λ + n ) Q T , (cid:15) − = Q diag( λ − , ..., λ − n ) Q T , where Q diag( λ , ..., λ n ) Q T is the eigen-decomposition (or the spectral decomposition) of (cid:15) . Noticethat (cid:15) + represents the tensile strain component that contributes to the damage process resulting incrack initiation and propagation while (cid:15) − represents the compression strain component that doesnot contribute to the damage process. Based on this, the active and passive components of theelastic energy are given byΨ e,act ( (cid:15) ) = λ (cid:0) (tr( (cid:15) )) + (cid:1) + µ tr(( (cid:15) + ) ) , Ψ e,pas ( (cid:15) ) = λ (cid:0) (tr( (cid:15) )) − (cid:1) + µ tr(( (cid:15) − ) ) . The elastic energy density in (6) and the Cauchy stress in (10) can be written asΨ le ( (cid:15), d ) = ( d + k l ) (cid:18) λ ((tr( (cid:15) )) + ) µ tr (cid:0) ( (cid:15) + ) (cid:1)(cid:19) + (cid:18) λ ((tr( (cid:15) )) − ) µ tr (cid:0) ( (cid:15) − ) (cid:1)(cid:19) ,σ = ( d + k l ) (cid:0) λ (tr( (cid:15) )) + I + 2 µ(cid:15) + (cid:1) + (cid:0) λ (tr( (cid:15) )) − I + 2 µ(cid:15) − (cid:1) . In this model, only the energy component related to the tensile strain component is degraded indamaged regions . 6 .2.2 Volumetric-deviatoric split
Another commonly-used model is proposed by Amor et al. [1] based on the volumetric-deviatoricsplit (v-d split), (cid:15) = (cid:15) S + (cid:15) D , (cid:15) S = 1 m tr( (cid:15) ) I, (cid:15) D = (cid:15) − m tr( (cid:15) ) I, (12)where m is the spatial dimension, (cid:15) S is the spherical component (called the volumetric strain tensor)related to the volume change, and (cid:15) D is the deviatoric component (called the strain deviator tensor)related to distortion. In this model, the expansive volumetric strain energy and the total deviatoricstrain energy are released by the creation of new cracks whereas the compressive volumetric strainenergy is not . The active and passive components of the elastic energy are given byΨ e,act ( (cid:15) ) = κ ((tr( (cid:15) )) + ) µ tr (cid:0) ( (cid:15) D ) (cid:1) , Ψ e,pas ( (cid:15) ) = κ ((tr( (cid:15) )) − ) , where κ = λ + 2 µ/m is the bulk modulus of the material. The energy density and Cauchy stressassociated with this model areΨ le ( (cid:15), d ) = ( d + k l ) (cid:18) κ ((tr( (cid:15) )) + ) µ tr (cid:0) ( (cid:15) D ) (cid:1)(cid:19) + κ ((tr( (cid:15) )) − ) ,σ = ( d + k l ) (cid:0) κ (tr( (cid:15) )) + I + 2 µ(cid:15) D (cid:1) + κ (tr( (cid:15) )) − I. According to unilateral contact conditions, crack can only open in the regions where the materialtends to expand. For the volumetric-deviatoric split model described in the previous subsection,both the positive part of the volumetric component and the total deviatoric component of the elasticenergy are degraded in damaged regions. However, the three principal strains of the deviatoriccomponent of the strain tensor can be negative, which indicates that the compressive deviatoricstrain energy is also accounted for contributing to the damage process in the v-d split model.To avoid this unphysical feature, we propose to apply the spectral decomposition to the straindeviatoric tensor and call the resulting model as the improved v-d split model. Specifically, thespectral decomposition of (cid:15) D is (cid:15) D = (cid:15) + D + (cid:15) − D , (cid:15) + D = Q diag( λ +1 , ..., λ + m ) Q T , (cid:15) − D = Q diag( λ − , ..., λ − m ) Q T , provided that Q diag( λ , ..., λ m ) Q T is the eigen-decomposition of (cid:15) D . The degradation applies onlyto the expansive volumetric and deviatoric strain energies.
Then, the active and passive part of theelastic energy are given byΨ e,act ( (cid:15) ) = κ ((tr( (cid:15) )) + ) µ tr (cid:0) ( (cid:15) + D ) (cid:1) , Ψ e,pas ( (cid:15) ) = κ ((tr( (cid:15) )) − ) µ tr (cid:0) ( (cid:15) − D ) (cid:1) , and the corresponding energy density and Cauchy stress areΨ le ( (cid:15), d ) = ( d + k l ) (cid:18) κ ((tr( (cid:15) )) + ) µ tr (cid:0) ( (cid:15) + D ) (cid:1)(cid:19) + (cid:18) κ ((tr( (cid:15) )) − ) µ tr (cid:0) ( (cid:15) − D ) (cid:1)(cid:19) ,σ = ( d + k l ) (cid:0) κ (tr( (cid:15) )) + I + 2 µ(cid:15) + D (cid:1) + (cid:0) κ (tr( (cid:15) )) − I + 2 µ(cid:15) − D (cid:1) . .3 Improved treatment of crack boundary conditions We recall that on any crack Γ, the material is totally damaged, the phase-field variable is zero (i.e., d ( x, t ) = 0), and the stress vanishes, viz., σ · n | Γ = 0 , (13)where n is the unit outward normal to Γ. However, there is no guarantee that this is satisfied inthe above described three models of energy decomposition. Indeed, we recall that Ψ e,pas ( (cid:15) ) is notdegraded. By close examination, we can see that it does not vanish on Γ in general for all of themodels described above; see Fig. 2(a) for a sketch for this remaining energy in the totally damagedzone. This can also be explained from (10). On Γ, we have d = 0 and σ = k l ∂ Ψ e,act ∂(cid:15) + ∂ Ψ e,pas ∂(cid:15) ≈ ∂ Ψ e,pas ∂(cid:15) , where we have used k l (cid:28) l (cid:28)
1. The fact that ∂ Ψ e,pas /∂(cid:15) does not vanish on Γ in general impliesthat it is not guaranteed that (13) be satisfied. The violation of the boundary condition can leadto unphysical crack propagation. To see the impacts, we consider a single edge notched shear test,with the domain and boundary conditions shown in Fig. 3(b). In this case, a relative displacementis expected and no stress should remain on the crack surface. However, the results obtained withthe spectral decomposition model show a stiff response (Fig. 4(b)) and unrealistic displacementunder loading (Fig. 4(a)).To avoid this problem, we propose here a modification which introduces a critically damagedzone with d < d cr and degrades the total stored energy in the zone; see Fig. 1(b). Here, d cr ∈ (0 , d cr = 0, the modified model reduces back to the original model. On theother hand, when d cr = 1, the total energy is degraded for the entire damaged region. For any d cr ∈ (0 , the total energy is degraded in the critically damaged zone with ≤ d < d cr and only theactive component of the energy is degraded in the damaged zone with d cr < d ≤ . More discussionon the choice of d cr and its effects will be given in the numerical result section §
4. Mathematically,the modified damage model reads asΨ le,m ( (cid:15), d ) = (cid:0) d + k l (cid:1) Ψ le,act,m ( (cid:15) ) + Ψ le,pas,m ( (cid:15) ) , (14)where Ψ le,act,m ( (cid:15) ) = (cid:40) Ψ le,act ( (cid:15) ) + Ψ le,pas ( (cid:15) ) = Ψ le , for 0 ≤ d ≤ d cr Ψ le,act ( (cid:15) ) , for d cr < d ≤ le,pas,m ( (cid:15) ) = (cid:40) , for 0 ≤ d ≤ d cr Ψ le,pas ( (cid:15) ) , for d cr < d ≤ . By construction, Ψ le,m ( (cid:15), d ) = 0 and σ = 0 on Γ (assuming that k l = 0) and thus the crack bound-ary condition (13) is satisfied. For simplicity, we will call this modification as ItCBC (Improvedtreatment of Crack Boundary Conditions).For the single edge notched shear test considered earlier, the results with ItCBC with d cr = 0 . hase-field variable d E ne r g y den s i t y (a) Original model: Ψ le Phase-field variable d E ne r g y den s i t y (b) Modified model: Ψ le,m Phase-field variable d E ne r g y den s i t y (c) Comparison of the Ψ le and Ψ le,m Figure 2: Energy stored in the bulk of the elastic body. XY . mm . mm (a) Tension test XY . mm . mm U (b) Shear test Figure 3: Domain and boundary conditions for single edge notched tests, (a) tension test forExample 1 and (b) shear test for Example 2. 9 (a) Mesh on deformed domain (b) Von Mises stress distribution on thedeformed domain
Figure 4: The single edge notched shear test for the energy spectral decomposition model. Thevon Mises stress is defined as ( σ x + σ y − σ x σ y + 3 τ xy ) . (a) Mesh on deformed domain (b) Von Mises stress distribution on thedeformed domain Figure 5: The single edge notched shear test with ItCBC with d cr = 0 . A moving mesh finite element method
In this section we briefly describe a moving mesh finite element method for solving the phase-fieldproblem (8) and (9). It was first considered for the phase-field modeling of brittle fracture by Zhanget al. in [25]. The reader is referred to the reference for more detailed discussion of the method.
Let T h be a simplicial mesh for the domain Ω and denote N and N v as the number of its elementsand vertices, respectively. The function spaces V d , V u and V u are approximated by V hd = (cid:8) ϕ h | ϕ h ∈ C (Ω); ϕ h | K ∈ P ( K ) , ∀ K ∈ T h (cid:9) ⊂ V d ,V hu = (cid:8) ϕ h | ϕ h ∈ C (Ω); ϕ h | ∂ Ω u = u ; ϕ h | K ∈ P ( K ) , ∀ K ∈ T h (cid:9) ⊂ V u ,V ,hu = (cid:8) ϕ h | ϕ h ∈ C (Ω); ϕ h | ∂ Ω u = 0; ϕ h | K ∈ P ( K ) , ∀ K ∈ T h (cid:9) ⊂ V u , where P ( K ) is the set of polynomials of degree less than or equal to 1 defined on K . For thephase-field problem (8) and (9), the linear finite element approximation is to find d h ∈ V hd and u h ∈ V hu such that (cid:90) Ω (cid:18)(cid:18) d h H + g c ( d h − l (cid:19) ϕ h + 2 g c l ∇ d h · ∇ ϕ h (cid:19) d Ω = 0 , ∀ ϕ h ∈ V hd (15) (cid:90) Ω σ ( u h ) : (cid:15) ( ϕ h ) d Ω = (cid:90) Ω t t · ϕ h dS + (cid:90) Ω f · ϕ h d Ω , ∀ ϕ h ∈ V ,hu . (16)In our computation, we solve (15) for d h and (16) for u h alternately. The procedure leads tosmaller and easier systems to solve since d h and u h are decoupled and (15) is linear about d h . Werecall that the energy density is decomposed into active and passive components with only theformer contributing to the evolution of cracks. This decomposition results in a non-smooth elasticenergy, increases nonlinearity of the displacement system, and makes Newton’s iteration often diffi-cult to converge. Three regularization methods have been proposed by Zhang et al. [25] to smoothpositive and negative eigenvalue functions via a switching technique (sonic-point regularization) orconvolution with a smoothed delta function. In this work, we use the sonic-point regularizationmethod with which the positive and negative eigenvalue functions are replaced by λ + α = λ + √ λ + α , λ − α = λ − √ λ + α , where α > t being introduced torepresent the load increments. The solution procedure from t n to t n +1 is described as follows.(i) Suppose the mesh T nh at time t n and the history field H nh in (11) (defined on T nh ) are given.(ii) Compute the phase-field variable d n +1 h and new mesh T n +1 h as follows. • Let T n +1 , h = T nh ; • For k = 1 : kk
11 Compute H on T n +1 ,kh using linear interpolation of H nh ;- Compute d n +1 ,kh using (15) and H on T n +1 ,kh ;- If k < kk , compute the new mesh T n +1 ,k +1 h by the MMPDE moving mesh methodbased on T n +1 ,kh and d n +1 ,kh ; see Section 3.2. • Let T n +1 h = T n +1 ,kkh and d n +1 h = d n +1 ,kkh .(iii) Compute the displacement field u n +1 h by solving the nonlinear system (16) based on d n +1 h and T n +1 h . Newton’s iteration is used.(iv) Compute Ψ l,n +1 e,act,h ( (cid:15) ( u n +1 h )) and set H n +1 h = max { Ψ l,n +1 e,act,h ( (cid:15) ( u n +1 h )) , ˜ H nh } , where ˜ H nh is the linearinterpolation of H nh from the old mesh T nh to the new mesh T n +1 h .The parameter kk in the above procedure determines the adaptivity of the mesh T n +1 h to thephase-field variable d n +1 h . Our experience shows that kk = 5 is sufficient for the mesh to be welladaptive to d n +1 h . We use the MMPDE moving mesh method [12, 14] for generating the new mesh T n +1 h . The methodis based on the M -uniform mesh approach with which a nonuniform adaptive mesh is viewed as auniform one in the metric specified by a tensor M . The metric tensor M is assumed to be symmetricand uniformly positive definite on Ω and determines the shape, size, and orientation of the meshelements through the so-called equidistribution and alignment conditions (see (18) and (19) below).Denote H ( d h ) as a recovered Hessian of d h and assume its eigen-decomposition is given by H ( d h ) = Q diag( λ , λ ) Q T . Then we choose the metric tensor in our computation as M = det( I + | H ( d h ) | ) − ( I + | H ( d h ) | ) , (17)where | H ( d h ) | = Q diag( | λ | , | λ | ) Q T . Since (17) is based on the Hessian of the phase-field variable d , the mesh elements are concentrated around the crack where the curvature of d is large. Theform (17) is known to be optimal for the L norm of linear interpolation error (e.g., see [14]).Denote x , ...., x N v as the coordinates of the vertices of T h . Let ˆ T c,h = { ˆ ξ , ..., ˆ ξ N v } be thereference computational mesh which is taken as the initial physical mesh. We also denote T c,h = { ξ , ..., ξ N v } as an intermediate computational mesh. We assume that T h , ˆ T c,h , and T c,h have thesame number of elements and vertices and the same connectivity. Then for any element K ∈ T h ,there exists a corresponding element K c ∈ T c,h . Let F K be the affine mapping from K c to K and F (cid:48) K be its Jacobian matrix. Denote the vertices of K and K c as x K , x K , x K and ξ K , ξ K , ξ K ,respectively. We define the edge matrices of K and K c as E K = [ x K − x K , x K − x K ] , ˆ E K = [ ξ K − ξ K , ξ K − ξ K ] . It is easy to show that F (cid:48) K = E K ˆ E K − , ( F (cid:48) K ) − = ˆ E K E − K . It is known (e.g., see [9, 14]) that an M -uniform mesh T h approximately satisfies the equidistri-bution and alignment conditions | K | (cid:112) det( M K ) = | Ω h | | K c || Ω c | , ∀ K ∈ T h (18)122 tr (cid:0) ( F (cid:48) K ) T M K F (cid:48) K (cid:1) = det (cid:0) ( F (cid:48) K ) T M K F (cid:48) K (cid:1) , K ∈ T h (19)where | K | and | K c | denote the area of K and K c , respectively, M K is the average of M over K ,det( · ) denotes the determinant of a matrix, and | Ω h | = (cid:88) K ∈T h | K | (cid:112) det( M K ) , | Ω c | = (cid:88) K c ∈T c,h | K c | . An energy function that combines the above two conditions has been proposed by Huang [8] as I h ( T h ; T c,h ) = (cid:88) K ∈T h | K | G ( J K , det( J K ) , M K ) , (20)where J K = ( F (cid:48) K ) − and G ( J K , det( J K ) , M K ) = 13 (cid:112) det( M K ) (cid:0) tr( J K M K J TK ) (cid:1) + 2 (cid:112) det( M K ) (cid:32) det( J K ) (cid:112) det( M K ) (cid:33) . In principle, a new physical mesh T n +1 h can be found by minimizing I h with respect to T h for agiven T c,h (for example, taken as ˆ T c,h ). However, this minimization can be difficult and costly since I h is generally not convex. Here, we use the ξ -formulation of the MMPDE moving mesh method.That is, we first take T h = T nh and define the moving mesh equation as a gradient system of I h with the coordinates of the computational vertices, i.e., dξ j dt = − P j τ (cid:18) ∂I h ∂ξ j (cid:19) T , j = 1 , ..., N v (21)where τ > P j = det( M ( x j )) ischosen to make (21) invariant under the scaling transformation of M . Using the analytical formulasfor the derivatives [10], we can rewrite (21) into dξ j dt = P j τ (cid:88) K ∈ ω j | K | v Kj K , (22)where ω j is the element patch associated with the j -th vertex, j K is the corresponding local indexof the vertex on K , and v Kj K is the local velocity of the vertex that is given by (cid:34) ( v K ) T ( v K ) T (cid:35) = − E − K ∂G∂ J − ∂G∂ det( J ) det( ˆ E K )det( E K ) ˆ E − K , v K = − (cid:88) i =1 v Kj . The derivatives of the function G in the above equation are given by ∂G∂ J = (cid:112) det( M ) (cid:0) tr( JM − J T ) (cid:1) M − J T ,∂G∂ det( J ) = √ M ) − det( J ) . Note that the velocities for the boundary vertices need to be modified appropriately so thatthey stay on the boundary. After that, (22) is integrated using the Matlab R (cid:13) ODE solver ode15s t n to t n +1 with ˆ T c,h as the initial mesh. The new computational mesh is denoted as T n +1 c,h .This mesh and the physical mesh T nh form a piecewise linear correspondence, i.e., T nh = Φ h ( T n +1 c,h ).The new physical mesh is then defined as T n +1 h = Φ h ( ˆ T c,h ), which can be computed using linearinterpolation.It is worth pointing out that an x -formulation of the MMPDE moving mesh method can also beobtained by taking T c,h = ˆ T c,h and using the gradient system of I h with respect to the coordinatesof the physical vertices. Although more complicated to implement (e.g., see [10]) than the ξ -formulation, the x -formulation has the advantage that its generated mesh is theoretically andnumerically guaranteed to be nonsingular if it is so initially (cf. [11]). On the other hand, a similartheoretical result has not yet been proven for the ξ -formulation although numerical experiment hasshown that it also produces no mesh tangling or crossing. In this section we present numerical results for four examples. The first two are classical benchmarkproblems, a single edge notched tension test and a shear test. The third example is designed todemonstrate the ability of our method to handle complex cracks. The last example is also a singleedge notched shear test but with the physical parameters and domain geometry chosen based on anexperiment setting. The numerical results obtained with the three energy decomposition modelsand their modified ones are presented and compared. The effects of the critical damage threshold d cr and the improved treatment of crack boundary conditions on the numerical solution are discussed.Unless stated otherwise, the following choices of the parameters are used: α = 10 − in the sonic-point regularization method, N = 6 ,
400 for the size of an adaptive mesh, and l = 0 . We first consider a single edge notched tension test from Miehe et al. [18]. The geometry andboundary conditions are show in Fig. 3(a). The bottom edge of the domain is fixed. The top edgeis fixed along the x -direction while along the y -direction a uniform displacement U is increased withtime to drive the crack propagation. The following material properties are used in our computation: λ = 121 .
15 kN/mm , µ = 80 .
77 kN/mm , and g c = 2 . × − kN/mm. Due to the brutalcharacter of the crack propagation, we choose two displacement increments for the computation,∆ U = 10 − mm for the first 500 time steps and ∆ U = 10 − mm afterwards. For comparisonpurpose, we compute the surface load vector on the top edge as F = ( F x , F y ) ≡ (cid:90) top edge σ ( (cid:15) ) · n dl, where n is the unit outward normal to the top edge. We are interested in F y for the tension testand F x for the shear test.We now investigate the effects of different energy decomposition models and the crack boundaryconditions on the numerical results. An initial triangular mesh is constructed from a rectangularmesh by subdividing each rectangle into four triangles along diagonal directions. Typical adaptivemeshes and contours of the phase field and von Mises stress distribution using three decomposition14odels are shown in Fig. 6. The load-deflection curves are shown in Fig. 7. For the spectraldecomposition and v-d split method, the resulting crack path and load-deflection curves agree wellwith the results obtained by Miehe at al. [19] using pre-refined mesh around the regions of thecrack and its expected propagation path. Moreover, the stress is always concentrated in the regionaround the crack tip during crack evolution as shown in Figs. 6(g) and 6(h). Inconsistent resultsare obtained with for the improved v-d split model. The stress growth and concentration occur inthe totally damaged area as shown in Fig. 6(i). The reason of this unphysical phenomenon is thatfor the tension test, this split model leads to a small remaining energy in the damaged zone.Next, we examine the effects of ItCBC (cf. § d cr . UsingItCBC with d cr = 0 .
4, typical adaptive meshes and contours of the phase-field variable and vonMises stress at U = 5 . × − mm are shown in Fig. 8. One can see that the effects of ItCBC forboth the spectral decomposition and v-d split models are small. On the other hand, the effects ofItCBC on the improved v-d split model are significant. Indeed, the von Mises stress distributionfor the latter is now comparable to those obtained with the spectral decomposition and v-d splitmodels; see Fig. 8(i).The load-deflection curves for three decomposition models with various values of d cr (0.2, 0.4,0.6, 0.8) are shown in Figs. 9 and 10. As we can see in Figs. 9(a) and 9(b), there is no significantdifference between the original model and the model with ItCBC with various values of d cr for thespectral decomposition and v-d split. However, for the improved v-d split, smaller d cr values leadto underestimates of the load after crack starts propagating, see 9(c). The work done by boundarytractions and body forces (external work) on an elastic solid are stored inside the body in theform of strain energy. According to the Griffith theory, the energy required to create new cracksurfaces is transformed from the stored elastic strain energy. When d cr is large, more strain energyis degraded and stored energy becomes less. In this case, more external work is needed to generatenew crack, which means that the peak of reaction force increases with d cr (see Fig. 9).Figs. 10(b), 10(c), and 10(d) show that for ItCBC with d cr ∈ [0 . , . d cr values ( d cr = 0 . We now consider a single edge notched shear test. The domain and boundary conditions are shownin Fig. 3(b). The bottom edge of the domain is fixed and the top edge is fixed along the y -directionwhile a uniform x -displacement U is increased with time to drive the crack propagation. Thematerial properties are the same as the tension test in Example 1, that is, λ = 121 .
15 kN/mm , µ = 80 .
77 kN/mm , and g c = 2 . × − kN/mm. The displacement increment is chosen as∆ U = 10 − mm for the computation.We first investigate the effects of the three decomposition models and ItCBC on the crack prop-agation and the distribution of the stress. The spectral decomposition model without ItCBC leadsto the development of a secondary crack along the left edge of the domain and the concentrationof the stress in the damaged region; see Fig. 11. This (unphysical) phenomenon has also beenobserved by May et al. [17] where the initial notch is modeled with d = 0. The results with ItCBCmodification ( d cr = 0 .
4) are shown in Fig. 12. It can be seen that a secondary crack does not occurand the stress concentrates at the turning point and crack tip only. The crack path agrees well15 a) spectral decomposition (b) v-d split (c) improved v-d split (d) spectral decomposition (e) v-d split (f) improved v-d split(g) spectral decomposition (h) v-d split (i) improved v-d split
Figure 6: Example 1. Meshes and contours of the phase-field and von Mises stress distribution areplotted at U = 5 . × − mm. Three energy decomposition models are used. displacement U [mm] -3 l oad F [ K N ] spectral decompositionv-d split Figure 7: Example 1. The load-deflection curves are obtained for spectral and v-d split decompo-sition models. 16 a) spectral decomposition (b) v-d split (c) improved v-d split (d) spectral decomposition (e) v-d split (f) improved v-d split(g) spectral decomposition (h) v-d split (i) improved v-d split
Figure 8: Example 1. Meshes and contours of the phase-field and von Mises stress distribution areplotted at U = 5 . × − mm. Three decomposition models with ItCBC ( d cr = 0 .
4) are used.17 displacement U [mm] -3 l oad F [ K N ] d cr = 0.2d cr = 0.4d cr = 0.6d cr = 0.8 (a) spectral decomposition displacement U [mm] -3 l oad F [ K N ] d cr = 0.2d cr = 0.4d cr = 0.6d cr = 0.8 (b) v-d split displacement U [mm] -3 l oad F [ K N ] d cr = 0.2d cr = 0.4d cr = 0.6d cr = 0.8 (c) improved v-d split Figure 9: Example 1. The load-deflection curves are obtained using three decomposition modelswith or without ItCBC. (a) spectral decomposition with various d cr ; (b) v-d split with various d cr ;(c) improved v-d split with various d cr . 18 displacement U [mm] -3 l oad F [ K N ] spectral decomposition d cr = 0.2v-d split d cr = 0.2improved v-d split d cr = 0.2 (a) d cr = 0 . displacement U [mm] -3 l oad F [ K N ] spectral decomposition d cr = 0.4v-d split d cr = 0.4improved v-d split d cr = 0.4 (b) d cr = 0 . displacement U [mm] -3 l oad F [ K N ] spectral decomposition d cr = 0.6v-d split d cr = 0.6improved v-d split d cr = 0.6 (c) d cr = 0 . displacement U [mm] -3 l oad F [ K N ] spectral decomposition d cr = 0.8v-d split d cr = 0.8improved v-d split d cr = 0.8 (d) d cr = 0 . Figure 10: Example 1. The load-deflection curves are compared for three energy decompositionmethods with ItCBC. (a) d cr = 0 .
2; (b) d cr = 0 .
4; (c) d cr = 0 .
6; (d) d cr = 0 . d cr = 0 . d cr = 0 . d cr . The load-deflection curves obtained usingthe spectral decomposition and improved v-d split model with various values of d cr are plotted inFigs. 17 and 18.As can be seen in Fig. 17(a), the effects of d cr on the load-deflection are small for the spectraldecomposition model when d cr is small ( d cr ≤ . d cr (e.g., d cr = 0 . d cr ; see Fig. 17(b). It is also interesting to observe that thecurves for the two decomposition methods with d cr = 0 . In this example, we test the modeling of complex crack systems. We consider a square plate ofwidth 2 mm with two or five cracks. The domain and boundary conditions are shown in Figs. 19(a)and 19(b), respectively. For both problems, the bottom edge of the domain is fixed. The top edgeis fixed along y -direction while a uniform x -displacement U is increased with time. The materialparameter are the same as in Example 1 except g c = 2 . × − kN/mm for the five-crack problem.For the two-crack problem, the lengths of Crack 1 and Crack 2 are 0 . . ◦ and 64 . ◦ , respectively.Fig. 20 shows the adaptive mesh and phase-field and von Mises stress distribution for spectraldecomposition with original crack boundary conditions. The results for all three decompositionmethods with ItCBC ( d cr = 0 .
4) are shown in Figs. 21, 22, and 23, respectively. The results usingimproved v-d split in Fig. 23 agree well with those using spectral decomposition as shown in Fig.21 whereas the v-d split (even with ItCBC) gives unphysical crack propagation. This finding isconsistent with the observations made from the previous examples.For the five-crack problem, the lengths of Crack 1 to Crack 5 are 0 . .
35 mm, 0 .
35 mm,0 . . ◦ , 45 ◦ , 17 . ◦ , 28 . ◦ and 9 ◦ , respectively. The phase-field andstress distribution for the spectral decomposition method without and with ItCBC ( d cr = 0 .
4) areshown in Fig. 24 and 25. Again, unphysical stress concentration in the damaged zone is observedfor the original treatment of crack boundary conditions but vanishes with ItCBC.20 a) U = 2 . × − mm (b) U = 2 . × − mm (c) U = 2 . × − mm (d) U = 2 . × − mm (e) U = 2 . × − mm (f) U = 2 . × − mm(g) U = 2 . × − mm (h) U = 2 . × − mm (i) U = 2 . × − mm Figure 11: Example 2. Meshes and contours of the phase-field and von Mises stress distributionare plotted at U = 2 . × − , 2 . × − , and 2 . × − mm. The spectral decompositionmodel is used. 21 a) U = 1 . × − mm (b) U = 1 . × − mm (c) U = 1 . × − mm (d) U = 1 . × − mm (e) U = 1 . × − mm (f) U = 1 . × − mm(g) U = 1 . × − mm (h) U = 1 . × − mm (i) U = 1 . × − mm Figure 12: Example 2. Meshes and contours of the phase-field and von Mises stress distributionare plotted at U = 1 . × − , 1 . × − , and 1 . × − mm. The spectral decomposition modelwith ItCBC ( d cr = 0 .
4) is used. 22 a) U = 1 . × − mm (b) U = 1 . × − mm (c) U = 1 . × − mm (d) U = 1 . × − mm (e) U = 1 . × − mm (f) U = 1 . × − mm(g) U = 1 . × − mm (h) U = 1 . × − mm (i) U = 1 . × − mm Figure 13: Example 2. Meshes and contours of the phase-field and von Mises stress distributionare plotted at U = 1 . × − , 1 . × − , and 1 . × − mm. The v-d split model is used.23 a) U = 1 . × − mm (b) U = 1 . × − mm (c) U = 1 . × − mm (d) U = 1 . × − mm (e) U = 1 . × − mm (f) U = 1 . × − mm(g) U = 1 . × − mm (h) U = 1 . × − mm (i) U = 1 . × − mm Figure 14: Example 2. Meshes and contours of the phase-field and von Mises stress distributionare plotted at U = 1 . × − , 1 . × − , and 1 . × − mm. The v-d split model with ItCBC( d cr = 0 .
4) is used. 24 a) U = 3 . × − mm (b) U = 3 . × − mm (c) U = 3 . × − mm (d) U = 3 . × − mm (e) U = 3 . × − mm (f) U = 3 . × − mm(g) U = 3 . × − mm (h) U = 3 . × − mm (i) U = 3 . × − mm Figure 15: Example 2. Meshes and contours of the phase-field and von Mises stress distributionare plotted at U = 3 . × − , 3 . × − , and 3 . × − mm. The improved v-d split model isused. 25 a) U = 1 . × − mm (b) U = 1 . × − mm (c) U = 1 . × − mm (d) U = 1 . × − mm (e) U = 1 . × − mm (f) U = 1 . × − mm(g) U = 1 . × − mm (h) U = 1 . × − mm (i) U = 1 . × − mm Figure 16: Example 2. Meshes and contours of the phase-field and von Mises stress distributionare plotted at U = 1 . × − , 1 . × − , and 1 . × − mm. The improved v-d split modelwith ItCBC ( d cr = 0 .
4) is used. 26 displacement U [mm] l oad F [ K N ] d cr = 0.2d cr = 0.4d cr = 0.6d cr = 0.8 (a) spectral decomposition displacement U [mm] l oad F [ K N ] d cr = 0.2d cr = 0.4d cr = 0.6d cr = 0.8 (b) improved v-d split Figure 17: Example 2. The load-deflection curves are plotted for two decomposition models. (a)The spectral decomposition model with ItCBC (with various d cr ); (b) The improved v-d split modelwith ItCBC (with various d cr ). displacement U [mm] l oad F [ K N ] spectral decomposition d cr = 0.2improved v-d split d cr = 0.2 (a) d cr = 0 . displacement U [mm] l oad F [ K N ] spectral decomposition d cr = 0.4improved v-d split d cr = 0.4 (b) d cr = 0 . displacement U [mm] l oad F [ K N ] spectral decomposition d cr = 0.6improved v-d split d cr = 0.6 (c) d cr = 0 . displacement U [mm] l oad F [ K N ] spectral decomposition d cr = 0.8improved v-d split d cr = 0.8 (d) d cr = 0 . Figure 18: Example 2. The load-deflection curves are compared for two decomposition models withItCBC. (a) d cr = 0 .
2; (b) d cr = 0 .
4; (c) d cr = 0 .
6; (d) d cr = 0 . Y mm . m m . mm (a) junction between two cracks . m m . m m . m m . m m . mm XY mm (b) random distribution of cracks Figure 19: Example 3. Domain and boundary conditions for the shear test with multiple cracks,(a) two cracks, (b) five cracks.
To further verify the decomposition models, we compare the simulation results with the experi-mental results obtained by Lee et al. [16]. The physical experiments were performed on a modifiedversion of the shear apparatus used in Reber et al. [20]. Gelatin was used as the material in theexperiments. The sketch of the experimental setup can be seen in Fig. 26(a). The spring blackarrow represents the direction of the shear force. The bottom side of the table moves under theconstant velocity while the fixed side is stationary. For the numerical computation, a rectangularplate with a length of 120 mm and a width of 70 mm is considered. The initial fracture is locatedat the middle of the plate with the length of 30 mm. The top edge of the domain is fixed andthe bottom edge is fixed along y -direction while a uniform x -displacement U is increased with time(∆ U = 5 × − mm). We consider a mesh size as 41 ×
41 ( N = 6 , l = 1 . . × Pa, the Poisson ratio is 0 .
45 and the fracture toughness is 1 .
96 Pa · m. (The Poisson ratiofor gelatin is 0.499. We choose the Poisson ratio to be 0.45 to avoid the locking effects in the finiteelement approximation of elasticity problems where the performance of certain commonly usedfinite elements deteriorates when the Poisson ratio is close to 0.5; e.g. see Babuˇ s ka and Suri [2].)These material properties correspond to λ = 4 . × − kN/mm , µ = 4 . × − kN/mm ,and g c = 1 . × − kN/mm. The results from the computation with three decomposition modelsand the experiment are comparable qualitatively in crack propagation. As can be seen in Fig. 27,with ItCBC ( d cr = 0.4), the spectral decomposition and the improved volumetric-deviatoric splitlead to results comparable with the physical experiment. The displacement load is almost the same( U = 3 . U = 3 . a) U = 6 . × − mm (b) U = 6 . × − mm (c) U = 6 . × − mm (d) U = 6 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (e) U = 6 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (f) U = 6 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (g) U = 6 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (h) U = 6 . × − mm(i) U = 6 . × − mm (j) U = 6 . × − mm (k) U = 6 . × − mm (l) U = 6 . × − mm Figure 20: Example 3. The mesh and contours of the phase-field and von Mises stress distributionduring crack evolution for the two-crack shear test with l = 0 . N = 10 ,
000 (51 × a) U = 2 . × − mm (b) U = 3 . × − mm (c) U = 3 . × − mm (d) U = 4 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (e) U = 2 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (f) U = 3 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (g) U = 3 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (h) U = 4 . × − mm(i) U = 2 . × − mm (j) U = 3 . × − mm (k) U = 3 . × − mm (l) U = 4 . × − mm Figure 21: Example 3. The mesh and contours of the phase-field and von Mises stress distributionduring crack evolution for the two-crack shear test with l = 0 . N = 10 ,
000 (51 × a) U = 2 . × − mm (b) U = 2 . × − mm (c) U = 2 . × − mm (d) U = 3 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (e) U = 2 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (f) U = 3 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (g) U = 2 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (h) U = 3 . × − mm(i) U = 2 . × − mm (j) U = 2 . × − mm (k) U = 2 . × − mm (l) U = 3 . × − mm Figure 22: Example 3. The mesh and contours of the phase-field and von Mises stress distributionduring crack evolution for the two-crack shear test with l = 0 . N = 10 ,
000 (51 × a) U = 2 . × − mm (b) U = 3 . × − mm (c) U = 3 . × − mm (d) U = 4 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (e) U = 2 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (f) U = 3 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (g) U = 3 . × − mm -1 -0.5 0 0.5 1-1-0.500.51 (h) U = 4 . × − mm(i) U = 2 . × − mm (j) U = 3 . × − mm (k) U = 3 . × − mm (l) U = 4 . × − mm Figure 23: Example 3. The mesh and contours of the phase-field and von Mises stress distributionduring crack evolution for the two-crack shear test with l = 0 . N = 10 ,
000 (51 × a) U = 0 .
179 mm (b) U = 0 .
182 mm (c) U = 0 .
184 mm (d) U = 0 .
186 mm (e) U = 0 . ./ (f) U = 0 .
179 mm (g) U = 0 .
182 mm (h) U = 0 .
184 mm (i) U = 0 .
186 mm (j) U = 0 . Figure 24: Example 3. The contour of the phase-field during crack evolution for the five-crackshear test with l = 0 . N = 25 ,
600 (81 × (a) U = 0 .
092 mm (b) U = 0 .
097 mm (c) U = 0 . U = 0 .
127 mm (e) U = 0 . U = 0 .
092 mm (g) U = 0 .
097 mm (h) U = 0 . U = 0 .
127 mm (j) U = 0 . Figure 25: Example 3. The contour of the phase-field during crack evolution for the five-crackshear test with l = 0 . N = 25 ,
600 (81 × i x ed s i de P u ll ed s i de (a) experiment setup XY mm
120 mm U30 mm (b) computation setup
Figure 26: Example 4. Sketch of the experiment and computation setup.model material, friction less deformation (no friction exist along the crack surface) and perfectboundary conditions may not accurately describe the experimental setting. Nevertheless, this willbe an interesting topic for future investigations. Moreover, as can be seen in Fig. 27(g), 27(h) and27(i), the volumetric-deviatoric split model leads to incorrect crack propagation../
A crucial component in the phase-field modeling of brittle fracture is the decomposition of theenergy into the active and passive components, with only the former contributing to the damageprocess that results in fracture and thus being degraded. An improper decomposition can lead tounphysical crack openings and propagation. In the previous sections, we have studied two commonlyused decomposition models, the spectral decomposition model [18] and the volumetric-deviatoric(v-d) split model [1]. The active energy consists of those related to the tensile strain component inthe former model and the expansive volumetric strain energy plus the total deviatoric strain energyin the latter model. Four examples have been used for those studies. The first two are classicalbenchmark problems, single edge notched tension and shear tests. The third example is designed todemonstrate the ability of the models to handle complex cracks. The last example is also a singleedge notched shear test but with the physical parameters and domain geometry chosen based onan experiment setting.Numerical results have shown that the v-d split model can lead to unphysical crack propagationfor some physical settings for which the spectral decomposition model gives correct crack propaga-tion; e.g., see Example 2 (a single edge notched shear test) and particularly Fig. 13. An explanationfor this is that the v-d split model includes the total deviatoric strain energy in the active energyand thus accounts both the compressive and expansive deviatoric strain energies for contributingto the damage process. To avoid this, we have proposed an improved v-d model (cf. Section 2.2.3)which degrades only the the expansive volumetric and deviatoric strain energies in the damagezone. This improved v-d model improves the v-d model and leads to results comparable with thoseobtained with the spectral decomposition model.34hysical experiment: (a) U = 3 . U = 5 . U = 9 mm Spectral decomposition: (d) U = 3 . U = 4 mm (f) U = 4 . v-d split: (g) U = 3 . U = 3 . U = 3 . Improved v-d split: (j) U = 4 . U = 4 . U = 4 . Figure 27: Example 4. Comparison of geometry of the fracture obtained from the numericalcomputation with different decomposition models (with ItCBC d cr = 0 .
4) and the experiment [16].35umerical results have also shown that all of the models, including the spectral decompositionand improved v-d models, can still lead to unphysical crack openings and propagation for variousphysical settings. A close examination on this has revealed that they do not generally satisfy thevanishing stress condition (13) and there is stress remaining on the crack surface. We have proposeda simple yet effective remedy (ItCBC or the Improved treatment of Crack Boundary Conditions, cf.Section 2.3): define a critically damaged zone with 0 < d < d cr , where d cr is a positive parameter,and then degrade both the active and passive components of the energy in this zone. It has beenshown that the spectral decomposition and improved v-d models with ItCBC lead to correct crackpropagation for all of the examples we have tested. It should be emphasized that they includeExample 3 with multiple cracks and Example 4 which is based on an experimental setting. In thelatter case, both the spectral decomposition and improved v-d models with ItCBC yield comparablecrack propagation results that also agree well qualitatively with the experiment. Moreover, it hasbeen shown that the numerical results are not very sensitive to the choice of d cr although d cr inthe range of [0.4, 0.6] seems to work best.Finally, the numerical examples have demonstrated that the MMPDE moving mesh method isable to dynamically concentrate the mesh elements around propagating cracks even for complexcrack systems. Acknowledgments.
F.Z. was supported by China Scholarship Council (CSC) and China University of Petroleum- Beijing (CUPB) for his research visit to the University of Kansas from September of 2015 toSeptember of 2017. F.Z. is thankful to Department of Mathematics of the University of Kansas forthe hospitality during his visit.
References [1] H. Amor, J.-J. Marigo, and C. Maurini. Regularized formulation of the variational brittlefracture with unilateral contact: Numerical experiments.
J. Mech. Phys. Solids. , 57:1209–1229, 2009.[2] I. Babuˇ s ka and M. Suri. Locking effects in the finite element approximation of elasticityproblems. Numer. Math. , 62:439-463, 1992.[3] M. J. Borden.
Isogeometric Analysis of Phase-field Models for Dynamic Brittle and DuctileFracture . The University of Texas at Austin, 2012. PhD thesis.[4] M. J. Borden, T. J. R. Hughes, C. M. Landis, and C. V. Verhoosel. A higher-order phase-fieldmodel for brittle fracture: formulation and analysis within the isogeometric analysis framework.
Comput. Methods Appl. Mech. Eng. , 273:100–118, 2014.[5] B. Bourdin, G. A. Francfort, and J. J. Marigo. Numerical experiments in revisited brittlefracture.
J. Mech. Phys. Solids , 48:797–826, 2000.[6] C. Budd, W. Huang, and R. Russell. Adaptivity with moving grids.
Acta Numerica , 18:111–241, 2009. 367] G. A. Francfort and J. J. Marigo. Revisiting brittle fracture as an energy minimization problem.
J. Mech. Phys. Solids , 46:1319–1342, 1998.[8] W. Huang. Variational mesh adaptation: isotropy and equidistribution.
J. Comput. Phys. ,174:903–924, 2001.[9] W. Huang. Mathematical principles of anisotropic mesh adaptation.
Comm. Comput. Phys. ,1:276–310, 2006.[10] W. Huang and L. Kamenski. A geometric discretization and a simple implementation forvariational mesh generation and adaptation.
J. Comput. Phys. , 301:322–337, 2015.[11] W. Huang and L. Kamenski. On the mesh nonsingularity of the moving mesh pde method.
Math. Comp. , 87:1887–1911, 2018.[12] W. Huang, Y. Ren, and R. D. Russell. Moving mesh methods based on moving mesh partialdifferential equations.
J. Comput. Phys. , 113:279–290, 1994.[13] W. Huang, Y. Ren, and R. D. Russell. Moving mesh partial differential equations (MMPDEs)based upon the equidistribution principle.
SIAM J. Numer. Anal. , 31:709–730, 1994.[14] W. Huang and R. D. Russell.
Adaptive Moving Mesh Methods . Springer, New York, 2011.Applied Mathematical Sciences Series, Vol. 174.[15] C. Kuhn and R. Muller. A continuum phase field model for fracture.
Eng. Frac. Mech. ,77:3625–3634, 2010.[16] S. Lee, J. E. Reber, N. W. Hayman, and M. F. Wheeler. Investigation of wing crack forma-tion with a combined phase-field and experimental approach.
Geophysical Research Letters. ,43:7946–7952, 2016.[17] S. May, J. Vignollet, and R. de Borst. A numerical assessment of phase-field models for brittleand cohesive fracture: Γ-convergence and stress oscillations.
European J. Mech. A/Solids. ,52:72–84, 2015.[18] C. Miehe, M. Hofacker, and F. Welschinger. A phase field model for rate-independent crackpropagation: Robust algorithmic implementation based on operator splits.
Comput. MethodsAppl. Mech. Eng. , 199:2765–2778, 2010.[19] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field modelsof fracture: Variational principles and multi-field FE implementations.
Int. J. Numer. Meth.Eng. , 83:1273–1311, 2010.[20] J.E. Reber, L.L. Lavier, and N.W. Hayman. Experimental demonstration of a semi-brittleorigin for crustal strain transients.
Nat. Geosci. , 8:712–715, 2015.[21] D. Schillinger, M. J. Bordenr, and H. K. Stolarski. Isogeometric collocation for phase-fieldfracture models.
Comput. Methods Appl. Mech. Eng. , 284:583–610, 2015.3722] M. Stronl and T. Seelig. A novel treatment of crack boundary conditions in phase field modelsof fracture.
Appl. Math. Mech. , 15:155–156, 2015.[23] M. Stronl and T. Seelig. On constitutive assumptions in phase field approaches to brittlefracture. , pages 3705–3712, 2016. Catania, Italy.[24] J. Vignollet, S. May, R. de Borst, and C. V. Verhoosel. Phase-field models for brittle andcohesive fracture.
Meccanica , 49:2587–2601, 2014.[25] F. Zhang, W. Huang, X. Li, and S. Zhang. Moving mesh finite element simulation for phase-field modeling of brittle fracture and convergence of Newton’s iteration.