A micromechanics-informed phase field model for brittle fracture accounting for the unilateral constraint
Yangyuanchen Liu, Cheng Cheng, Vahid Ziaei-Rad, Yongxing Shen
AA MICROMECHANICS-INFORMED PHASE FIELDMODEL FOR BRITTLE FRACTURE ACCOUNTINGFOR THE UNILATERAL CONSTRAINT
Yangyuanchen Liu a,b , Cheng Cheng a,b , Vahid Ziaei-Rad a , Yongxing Shen a, ∗ a University of Michigan – Shanghai Jiao Tong University Joint Institute, Shanghai JiaoTong University, 800 Dongchuan Road, Shanghai, 200240, China b These authors contributed equally to this work.
Abstract
We propose a new direction-dependent model for the unilateral constraint in-volved in the phase field approach to fracture and also in the continuous damagemechanics models. The construction of this phase field model is informed bymicromechanical modeling through the homogenization theory, where the repre-sentative volume element (RVE) has a planar crack in the center. The proposedmodel is made closely match the response of the RVE, including the friction-less self-contact condition. This homogenization approach allows to identify adirection-dependent phase field model with the tension-compression split ob-tained from cracked microstructures. One important feature of the proposedmodel is that unlike most other models, the material degradation is consistentlydetermined without artificial assumptions or ad hoc parameters with no physicalinterpretation, thus, a more realistic modeling is resulted. With standard testssuch as uniaxial loadings, three-point bending, simple shear, and through-cracktests, the proposed model predicts reasonable crack paths. Moreover, with theRVE response as a benchmark, the proposed model gives rise to an accuratestress-strain curve under shear loads, more accurate than most existing models.
Keywords:
Phase field approach to fracture, Variational theory of fracture,Micromechanics, Homogenization theory, Unilateral constraint
1. Introduction
The prediction of failure mechanisms due to crack initiation and propagationin solids is of great significance for engineering applications. Griffith’s theory[1] provides a criterion for crack propagation but it is unable to predict crackinitiation, merging, and branching. ∗ Corresponding author
Email address: [email protected] (Yongxing Shen)
Preprint submitted to Engineering Fracture Mechanics July 9, 2020 a r X i v : . [ c s . C E ] J u l ased on the Griffith’s theory [1], Francfort and Marigo [2] proposed a vari-ational theory of fracture and Bourdin et al. [3] regularized this theory fornumerical computation, the outcome of which is also named the phase field ap-proach to fracture. The phase field modeling of brittle fracture has shown itsadvantages on simulating complex fracture processes, see [3, 4, 5]. One of the main challenges in phase field modeling is how to account forthe unilateral constraint in material degradation. The original model pro-posed in [3] adopts an isotropic response of the cracked solid, i.e., it assumesthat both tension and compression loads contribute to cracking. This leadsto unphysical crack propagation under a compressive load. With the need of direction-dependent material degradation, a number of models have been de-veloped. Early models that do distinguish tension vs. compression were thoseproposed in [4], [5], and [6]. Amor et al. [4] built their model based on thevolumetric-deviatoric (V-D) split of the strain tensor, while Miehe et al. [5]based their model on the spectral decomposition of the strain tensor. Around roughly the same time, Freddi and Royer-Carfagni [6] proposed a V-D split formasonry-like materials that accounts for the Poisson effect. The V-D split [4]and the spectral decomposition [5] are both widely used. However, it is wellknown that in certain tests, unphysical predictions are resulted, e.g. bendingtest with the V-D split and through-crack shear test with the spectral decom- position. It is also noted that a nonvariational formulation of the V-D split isproposed in [7]. This model preserves the linearity of the symmetric mechanicalbehavior, and provides more realistic crack propagation at compressive loads,though a modified evolution equation is used for the phase field variable.Recently, more sophisticated models based on the local crack orientation were proposed. For example, in the two models by Strobl and Seelig [8, 9], thecrack orientation is assumed to coincide with the local gradient direction of thephase field. Note that both models are variationally inconsistent, thus, theyneed to introduce ad hoc criteria for the evolution of phase field, as opposedto variationally consistent models following energy minimization. In the model proposed by Steinke and Kaliske [10], a crack coordinate system is introduced torepresent the local crack orientation, and the split is based on the decompositionof the stress tensor with respect to the said crack orientation. Wu et al. [11]utilized a positive/ negative projection based on the spectral decomposition tomodel the unilateral behavior, where the crack orientation is assumed to be the principal direction of the effective stress tensor. This model is expected to sharethe same disadvantage as those employing the spectral decomposition.In summary, in these models to-date, the crack state is approximately iden-tified by the definition of an artificial condition, e.g., sign of volumetric strain,principal strain/ stress, or crack normal stress. In addition, the decomposition of strains and stresses into crack-driving and persistent portions are somehowartificially introduced into these models, which can lead to predicting unrealisticcrack propagation and mechanical behavior.In 2018, Cheng [12] and Cheng and Shen [13] proposed an innovative modelfor the unilateral constraint involved in the phase field approach for brittle fracture for a two-dimensional case. Through homogenization techniques, a2ink is established between the cracked microstructures and the correspondingmacroscopic phase field model. For this purpose the authors built the modelwith inspiration by Dascalu et al. [14]. The main feature of this approachis that as the constitutive model is based on the RVE response, no artificial condition or material parameter is introduced into the model, unlike for theaforementioned split models. In this way, the constitutive model is obtainedeither with analytical expressions or by numerical fitting.We note that independently, Storm et al. [15] constructed a homogenizationframework for their phase field model in a somewhat similar way. This paper is a continuation of [13]. We present a three-dimensional phasefield formulation, which allows deriving a variationally consistent direction-dependent material degradation. More precisely, we model the macroscopiccrack as a collection of microscopic cracks, each located in the center of a verysmall volume element, which is called a representative volume element (RVE), as in the homogenization theory. Then the behavior of each RVE is obtainedthrough a detailed analysis, in particular with the frictionless self-contact condi-tion exactly enforced. The overall behavior of the RVE is then used to constructthe phase field model.The proposed model is subjected to a number of simple tests and compared with a few existing models. In the uniaxial compression test, three-point bend-ing test, and simple shear test, the proposed model behaves similarly as thespectral decomposition model [5]. In the through-crack shear tests, it is wellknown that the spectral model responds in an unreasonable way, i.e., the ma-terial does not break completely across the crack. The proposed model avoids such phenomenon, as does the V-D split model [4].Another feature is the response of the proposed model in the case of a non-through crack subjected to shear tests. In such tests, the proposed model pre-dicts a response very close to the RVE analysis, and it is the only one thatprovides such uniformly accurate result. One reason that enables such accurate prediction is that the model has not only a tensile term and a compressive term,but also a shear term, of which both the tensile and shear terms will be degradedby the phase field, but to different degrees. While the degradation function fortensile tests is pre-selected, that for shear is obtained by fitting. Moreover, thefitted expression has an explicit dependence on the Poisson ratio. One restriction of the model is that the formulation is obtained for isotropiclinear elasticity, though generalization to anisotropy or finite elasticity is possi-ble.The content of the paper structured as follows. In Section 2, a number ofexisting models are recapitulated. In Section 3, the construction of the proposed phase field model is introduced in detail. In Section 4, several numerical exam-ples are tested and results from existing models are compared with the proposedmodel. It is found that only the proposed model gives an accurate response inall such simple tests, with a benchmark as a more detailed and costly three-dimensional RVE. Finally, a summary of the features of the proposed model is presented in Section 5. 3 . Existing phase field models In this section, essential ingredients of a few phase field models are recapit-ulated and compared, with a focus on the tension-compression decomposition.Most phase field approaches to brittle fracture are based on the variationaltheory of fracture proposed in [2], which in turn is built on the theory of brittlefracture by [1]. In the framework of [2], of central importance is the followingenergy functional: Π ( u , Γ) = (cid:90) Ω Ψ ( ε ( u )) dΩ + (cid:90) Γ g c d s, where Ω ⊂ R n , n = 2 ,
3, is the domain occupied by the material at the initialconfiguration, Γ is the crack path, g c is the energy release rate of crack propaga-tion, and Ψ denotes the strain energy density of the (pristine) material, whichis a function of the strain tensor ε ( u ) = [ ∇ u +( ∇ u ) T ] /
2. For an isotropic linearelastic material, the expression of Ψ reads Ψ ( ε ) = λ ε ) + µ ε : ε , where λ and µ are Lam´e constants such that µ > λ + 2 µ > In order to construct a numerical simulation method for the variational the-ory of fracture, [3] put forward a regularized formulation, in which the functionalΠ is regularized as Π l :Π l [ u , d ] := (cid:90) Ω Ψ [ ε ( u ) , d ] dΩ + g c (cid:90) Ω (cid:18) d l + l |∇ d | (cid:19) dΩ , which is a functional of the displacement field u and the phase field d , see Figure1. Here l is a regularization length scale parameter. Note that for d we adopt aconvention opposite to that of [3], i.e., d = 1 in this manuscript represents thecrack and d = 0 pristine material.For generality, we let t N denote the traction field applied on Γ N ⊂ ∂ Ω and b the body force field, then the expression of Π l is modified to beΠ l [ u , d ] = (cid:90) Ω Ψ [ ε ( u ) , d ] dΩ − (cid:90) Γ N t N · u d Γ − (cid:90) Ω b · u dΩ+ g c (cid:90) Ω (cid:18) d l + l |∇ d | (cid:19) dΩ . (1)We particularly focus on the strain energy density term Ψ ( ε , d ). Once Ψ isdecided, the Cauchy stress tensor is given by σ ( ε , d ) = ∂Ψ∂ ε , and the driving force for the evolution of d is given by − ∂Ψ/∂d . The requirements on Ψ are: 4 igure 1: Comparison of a sharp crack and its phase field representation. Theleft figure illustrates the sharp crack model; the right figure shows the phase fieldrepresentation of the crack through a smooth transition area, in which the phasefield variable is between 0 and 1. The phase field is one of the arguments in theregularized variational theory of fracture, also named the phase field approachto fracture.1. Ψ ( ε , d = 0) = Ψ ( ε ).2. Ψ ( ε , d ) ≥ ε and d .3. ∂Ψ ( ε , d ) /∂d ≤
0, for all ε and d . General expression of Ψ ( ε , d ) . In most models to follow, Ψ ( ε , d ) takes the fol-lowing form: Ψ ( ε , d ) = g ( d ) Ψ + ( ε ) + Ψ − ( ε ) , where Ψ + and Ψ − are the crack-driving and persistent parts of the strain energydensity, respectively, and g ( d ) is the degradation function, often taken as g ( d ) =(1 − d ) + k , with k a small positive number merely to prevent the tangent stiffnessmatrix of the system from being singular.The Cauchy stress is thus given by σ ( ε , d ) = g ( d ) σ + + σ − , σ ± := ∂Ψ ± ∂ ε . The evolution equation for the phase field variable follows from the mini-mization principle of Π l , i.e., ∂Ψ∂d + g c (cid:18) d(cid:96) − (cid:96) ∇ d (cid:19) = 0 . (2)The split models constructed on the basis of this evolution equation are varia- tionally consistent.We next proceed to examine three models.5 he isotropic model proposed by [3]. In the formulation of [3], the degradedstrain energy density Ψ is expressed as the entire strain energy density of theunbroken material Ψ multiplied by the degradation function g ( d ). More pre-cisely, Ψ ( ε , d ) = g ( d ) Ψ ( ε ) , σ ( ε , d ) = g ( d )[ λ (tr ε ) + 2 µ ε ] ,∂Ψ∂d = g (cid:48) ( d ) Ψ ( ε ) . As well known, since the degradation function is applied to the entire Ψ , thismodel does not distinguish the fracture behavior between tension and compres-sion, resulting in unphysical results when the crack is subjected to compressive or shear loads. The V-D (volumetric-deviatoric) decomposition model proposed by [4].
Thismodel splits the strain energy density into volumetric and deviatoric parts: Ψ ( ε ) = Ψ + ( ε ) + Ψ − ( ε ), where Ψ + ( ε ) = K (cid:104) tr ε (cid:105) + µ (cid:107) dev ε (cid:107) ,Ψ − ( ε ) = K (cid:104) tr ε (cid:105) − , σ ( ε , d ) = g ( d ) ( K (cid:104) tr ε (cid:105) + + 2 µ dev ε ) + K (cid:104) tr ε (cid:105) − ,∂Ψ∂d = g (cid:48) ( d ) (cid:18) K (cid:104) tr ε (cid:105) + µ (cid:107) dev ε (cid:107) (cid:19) , where K = λ + 2 µ/ ε := ε − (1 / ε ) , and (cid:104) a (cid:105) ± :=( a ±| a | ) /
2. This split is based on the assumption that the strain energy resultingin the decrease of the local volume will not be degraded.
The spectral decomposition model proposed by [16].
In this model, the split isbased on the spectral decomposition of the strain tensor. Let { ε a } a =1 be theprincipal strains and { E a } a =1 be the corresponding orthonormal principal di-rections. The expressions are as follows: Ψ ± ( ε ) = λ (cid:104) tr ε (cid:105) ± + µ (cid:88) a =1 (cid:104) ε a (cid:105) ± , σ ± ( ε ) = λ (cid:104) tr ε (cid:105) ± + 2 µ (cid:88) a =1 (cid:104) ε a (cid:105) ± E a ⊗ E a ,∂Ψ∂d = g (cid:48) ( d ) (cid:32) λ (cid:104) tr ε (cid:105) + µ (cid:88) a =1 (cid:104) ε a (cid:105) (cid:33) . Related models in damage mechanics.
The last two splits can trace their roots in the damage mechanics community. For example, the split adopted by in [4]is similar to that by [17], while the decomposition proposed in [16] is similar tothe formulations of [18, 19, 20, 21]. 6 he anisotropic model proposed by [11].
In this model, Wu et al. utilized apositive/negative projection based on the spectral decomposition. The spectraldecomposition of the effective stress is σ = λ (tr ε ) + 2 µ ε , σ = (cid:88) i =1 σ i p i ⊗ p i = σ + + σ − , σ ± = (cid:88) i =1 σ ± i p i ⊗ p i , where σ i and p i represent the i th, i = 1 , ,
3, eigenvalue and eigenvector of theeffective stress σ , respectively. The principal values σ + i are determined as σ +1 = (cid:104) σ (cid:105) ,σ +2 = (cid:104) max( σ , ˜ νσ ) (cid:105) ,σ +3 = (cid:104) max[max( σ , ν ( σ + σ )) , ˜ νσ ] (cid:105) , where ˜ ν := ν/ (1 − ν ). The SS models proposed by [8, 9]. [8, 9] proposed two models, which we willrefer to as SS1 and SS2 models. In these models, the expression of σ in termsof ε depends on the local crack orientation, which could be either the “active”case or the “passive” case. The local crack orientation is given by n = ∇ d/ |∇ d | .Also the authors define N = n ⊗ n . Both models take the following form: σ = (cid:40) σ act , if ε : N > , σ pas , if ε : N < . The SS1 model reduces the stiffness parallel to the crack in both active andpassive cases, which reads σ act = g ( d ) ( λ tr( ε ) + 2 µ ε ) , σ pas = g ( d ) λ tr( ε ) + 2 g ( d ) µ ε + (1 − g ( d ))( λ + 2 µ ) ( ε : N ) N . The SS2 model only degrades the stiffness normal to the crack. Its expressionis σ act = (cid:18) λ + ( g ( d ) − λ λ + 2 µ (cid:19) (tr ε ) + 2 µ ε + ( g ( d ) − (cid:18) λ + λ λ + 2 µ (cid:19) ((tr ε ) N + ( ε : N ) )+ 4(1 − g ( d )) (cid:18) λ + 2 µ − λ λ + 2 µ (cid:19) ( ε : N ) N + µ ( g ( d ) −
1) ( N · ε + ε · N ) , σ pas = λ tr( ε ) + 2 µ ε + 4 µ (1 − g ( d ))( ε : N ) N + µ ( g ( d ) −
1) ( N · ε + ε · N ) . ε > ε = 0, and n = e , where e i is the i th Cartesian coordinate axis,when d = 1,( σ act ) = (cid:26) λ − λ λ + 2 µ − (cid:18) λ + λ λ + 2 µ (cid:19) + 4 (cid:18) λ + 2 µ − λ λ + 2 µ (cid:19)(cid:27) ε , which is not positive for a certain range of λ/µ , giving rise to an unstablematerial. The SK model proposed by [10].
In this model, the split is based on the de-composition of the stress tensor with respect to the crack orientation. For eachpoint, a crack coordinate system is defined with mutually orthonormal vectors n , s , and t . The crack orientation n is obtained from the maximum principalstress direction in this model. The expressions of the directional split are asfollows: σ = g ( d ) σ + + σ − , σ + = H ( σ nn ) (cid:20) n ⊗ n + λλ + 2 µ ( s ⊗ s + t ⊗ t ) (cid:21) ( λ + 2 µ n ⊗ n ) : ε + µ [( n ⊗ s + s ⊗ n ) ⊗ ( n ⊗ s + s ⊗ n ) + ( t ⊗ n + n ⊗ t ) ⊗ ( t ⊗ n + n ⊗ t )] : ε , σ − = (cid:20) H ( − σ nn ) n ⊗ n − H ( σ nn ) λλ + 2 µ ( s ⊗ s + t ⊗ t ) (cid:21) ⊗ ( λ + 2 µ n ⊗ n ) : ε + ( s ⊗ s ) ⊗ ( λ + 2 µ s ⊗ s ) : ε + ( t ⊗ t ) ⊗ ( λ + 2 µ t ⊗ t ) : ε + µ [( s ⊗ t + t ⊗ s ) ⊗ ( t ⊗ s + s ⊗ t )] : ε , where the degradation function g ( d ; b ) is given by g ( d ; b ) = exp( bd ) − ( b ( d −
1) + 1) exp( b )( b −
1) exp( b ) + 1 , where b is an extra parameter to be specified. Remarks.
All the aforementioned models except the SS models are variation-ally consistent, while the SS models are variationally inconsistent. With the relaxation of variational consistency, the efficiency of the numerical calculationmay be improved. Nevertheless, non-variational models give rise to additionalcomplications, in that they need to introduce extra damage or failure criteriafor the evolution of d . In contrast, variationally consistent models are closelyrelated to Griffith’s theory and would better fit in a thermodynamic framework. Hence, the proposed model will also be made variationally consistent.
3. The proposed micromechanics-informed phase field model
In this section, the steps for constructing the proposed phase field model aredetailed. 8n order to construct our model, we borrow the basic concepts from Dascalu et al. [14], in which the authors applied asymptotic homogenization techniquesto model the overall behavior of a damaged elastic body. The main ingredient isa micromechanical energy analysis, performed on a finite-sized cell, which leads,through homogenization, to a macroscopic evolution equation for the damagevariable. Here, we follow the same idea to obtain a three-dimensional phase field model that incorporates direction-dependent material degradation.Essentially we model a possibly curved crack as a collection of fictitious smallflat cracks with a particular spatial distribution. More precisely, we assume thata portion of the solid can be divided into many cube-shaped subdomains (orsquare-shaped in 2D), at the center of each of which there exists a said small planar crack (or straight crack in 2D), similar to the configuration proposed in[14]. The behavior of the macroscopic crack is then described by the collectionof such small cracks, see Figure 2. Here we borrow a terminology of the homog-enization theory and call each cube-shaped subdomain a representative volumeelement (RVE), and will apply standard homogenization techniques to extract its response subjected to mechanical loading.
Figure 2:
Modeling a macroscopic crack as a collection of fictitious smallstraight cracks (planar crack in 3D), each in one RVE.In a nutshell, we will establish a one-to-one relation between the crack di-mension in the RVE and the local phase field d . Then for any macroscopic inputstrain tensor ε ∈ R n × n ( n = 2 ,
3) with standard homogenization technique, wecan calculate an average stress σ ∈ R n × n , and then obtain C ∈ R n × n × n × n , theeffective secant modulus, such that σ ( ε , d ) = C ( ε , d ) : ε . (3)The expression of C ( ε , d ) can be obtained by comparison with analyticalexpressions or by numerical fitting. This way Ψ ( ε , d ) can also be obtained. Consider the RVE shown in Figure 3. Let the material be isotropic elastic with Lam´e parameters λ and µ such that µ > λ + 2 µ >
0. We denote9 igure 3:
The boundary value problem of RVEby B the domain occupied by the RVE, a bounded open subdomain of R n witha smooth boundary. Without loss of generality, let Γ ⊂ B be a micro-cracknormal to the x -axis. In the sequel we will test three types of micro-cracks:(i) a straight crack at the center of a 2D RVE; (ii) a penny-shaped crack at the center of a 3D RVE; (iii) a square-shaped crack with rounded fillets ( r = 0 . L )at the center of the 3D RVE, as shown in Figure 3. Details of the numericalexperiments on the RVE are provided in Appendix A.The computational domain is then B s = B \ Γ . In B s , the equilibrium10quation and boundary conditions are div σ = , in B s , (4a) σ = C : ∇ u , in B s , (4b) u = ε · x , on ∂ B , (4c)either σ · n = , (cid:74) u · n (cid:75) > , or (cid:74) σ · n (cid:75) = , n · σ · n < , t · σ · n = , (cid:74) u · n (cid:75) = 0 , (cid:41) on Γ, (4d)where (4d) is the frictionless contact condition for Γ , n and t are the unit normaland tangent vectors of the crack Γ , respectively. Let + and − denote the twosides of Γ , and for any vector or tensor a , we let a ± denote its value evaluatedon the ± side of Γ . Also we denote n ± the unit outward normal to the ± sideof Γ . The jump operator (cid:74) · (cid:75) is defined such that (cid:74) a · n (cid:75) = a + · n + + a − · n − . Once the solution of problem (4) is obtained, the total strain energy and themacroscopic stress tensor is given by σ = 1meas( B ) (cid:90) B s σ d B , (5)which will be matched with (5) to obtain the expression of σ ( ε , d ) and Ψ ( ε , d ),if fitting is needed. Note that meas( B ) denotes the measure of B . By construction, this RVE behaves differently in tension and in compression,a desired behavior. Under the above assumptions, the effective elastic responsewill be orthotropic in both cases. In this work, the numerical solution to theproblem (4) is obtained with the commercial software Abaqus. More details aregiven in Appendix A. For convenience, we will drop the bar of σ and ε from now on within Section 3. From the setup of the RVE, it can be concluded that for a given d , if the signsof the components of ε (i.e., the loading mode being tensile or compressive) arefixed, then the secant modulus C is independent of ε . In other words, C is alsothe tangent modulus . We hereby propose an ansatz for the macroscopic strainenergy density (3): Ψ ( ε , d ) = 12 (cid:110) H ( β ( ε )) ε : C t : ε (cid:124) (cid:123)(cid:122) (cid:125) tension + [1 − H ( β ( ε ))] ε : C c : ε (cid:124) (cid:123)(cid:122) (cid:125) compression (cid:111) = 12 (cid:110) H ( β ( ε )) (cid:0) C t ε + C t ε + C t ε + 2 C t ε ε +2 C t ε ε + 2 C t ε ε + 4 C t ε + 4 C t ε + 4 C t ε (cid:1) + [1 − H ( β ( ε ))] (cid:0) C c ε + C c ε + C c ε + 2 C c ε ε +2 C c ε ε + 2 C c ε ε + 4 C c ε + 4 C c ε + 4 C c ε (cid:1)(cid:111) , (6)11here H is the Heaviside function, and β ( ε ), whose expression will be deter-mined later, switches between the cases in which the RVE is in tension or com-pression, under the given strain ε . Tensors C t = C t ( d ) and C c = C c ( d ) denote the effective tensile and compres-sive moduli, respectively, and the sense (tensile vs. compressive) is determinedby β ( ε ). Here the phase field d is closely related to the crack length ratio r a := 2 a/L , where 2 a is the crack length, whose definition depends on the shapeof Γ , and L is the side length of the RVE. For any d (or r a ) and ε , C tijkl and C cijkl can be obtained either by inspection or by numerically analyzing the RVE.Hence we construct C tijkl and C cijkl using the following ansatz C mijkl ( d ) = g ( d ) P mijkl + Q mijkl , m = t, c. Then Ψ ( ε , d ) = 12 (cid:110) H ( β ( ε )) (cid:8) g ( d ) ε : P t : ε + ε : Q t : ε (cid:9) + [1 − H ( β ( ε ))] { g ( d ) ε : P c : ε + ε : Q c : ε } (cid:111) . (7)Equations (6) and (7) then give rise to the following constitutive relations,when the mode m = t or c is known: σ = [ g ( d ) P m + Q m ] ε + [ g ( d ) P m + Q m ] ε + [ g ( d ) P m + Q m ] ε , (8a) σ = [ g ( d ) P m + Q m ] ε + [ g ( d ) P m + Q m ] ε + [ g ( d ) P m + Q m ] ε , (8b) σ = [ g ( d ) P m + Q m ] ε + [ g ( d ) P m + Q m ] ε + [ g ( d ) P m + Q m ] ε , (8c) σ = 2 [ g ( d ) P m + Q m ] ε , (8d) σ = 2 [ g ( d ) P m + Q m ] ε , (8e) σ = 2 [ g ( d ) P m + Q m ] ε . (8f) Remark 3.1.
Here we do not associate tension (or compression) with the degra-dation function g ( d ) , to make the model more general at the beginning. Remark 3.2.
At this point, the remaining task is to determine { P mijkl } and { Q mijkl } , which totals to 36 unknowns. This can be achieved in two means: (a) with prior knowledge, i.e., by making (7) reduce to usual expressions in thespecial cases of uniaxial tension, etc., and (b) by fitting, i.e., by matching (8) and (5) . Our strategy is to first deduce as many coefficients as possible withprior knowledge, and then resort to fitting. Before ending this section we list one of the constraints for { P mijkl } and { Q mijkl } , which is that Ψ ( ε , d = 0) = Ψ ( ε ) for either the tensile case H ( β ( ε )) = 112r the compressive case H ( β ( ε )) = 0. That is to say, Ψ ( ε ) = 12 (cid:110) [ P m + Q m ] ε + [ P m + Q m ] ε + [ P m + Q m ] ε + 2 [ P m + Q m ] ε ε + 2 [ P m + Q m ] ε ε + 2 [ P m + Q m ] ε ε + 4 [ P m + Q m ] ε + 4 [ P m + Q m ] ε + 4 [ P m + Q m ] ε (cid:111) ,m = t, c. (9)Here for a linear isotropic material, Ψ ( ε ) takes the following form: Ψ ( ε ) = λ ε + ε + ε ) + µ (cid:0) ε + ε + ε + 2 ε + 2 ε + 2 ε (cid:1) . (10)Equating (9) and (10) we get the following equations, for m = t, c , P m + Q m = P m + Q m = P m + Q m = λ + 2 µ, (11a) P m + Q m = P m + Q m = P m + Q m = λ, (11b) P m + Q m = P m + Q m = P m + Q m = µ. (11c) { P mijkl } and { Q mijkl } with prior knowlegde In this section we will determine most of the coefficients { P mijkl } and { Q mijkl } by specializing the model to simple tests.First, it is noted that the symmetry in the crack plane imply: P m = P m , Q m = Q m , (12a) P m = P m , Q m = Q m , (12b) P m = P m , Q m = Q m . (12c)i. Tensile tests. We first consider the uniaxial tension in the x -direction, i.e.,with macroscopic strain ε > m = t . Then(8b) yields σ = (cid:2) g ( d ) P t + Q t (cid:3) ε . We further assume, as consistent with existing models, that the degradationtakes full effect in this case, i.e., Q t = 0, which yields P t = λ + 2 µ per(11a). Next we consider the uniaxial tension loading in the x -direction, i.e., ε > σ = (cid:2) g ( d ) P t + Q t (cid:3) ε . In accordance with numerical results (see Figure A.28) we found that Q t =0 and hence (11b) can be simplified to P t = λ . Also, referring to (12b),we obtain Q t = 0, and P t = λ .13i. Compressive tests. We now examine the case of the uniaxial compressivetest in the x -direction, i.e., the only non-zero component of strain tensoris ε <
0. Then σ = [ g ( d ) P c + Q c ] ε . In this case, there is no degradation, i.e., P c = 0. Equation (11a) thengives Q c = λ + 2 µ . We then move on to the case of uniaxial compression in the x -direction,i.e., ε < σ = [ g ( d ) P c + Q c ] ε ,σ = [ g ( d ) P c + Q c ] ε . Numerical results (see Figures A.26 and A.28) confirm that there is nodegradation in this case, which yields P c = P c = 0. Then, (11a) and(11b) lead to Q c = λ + 2 µ and Q c = λ , respectively. Also, due to(12a) and (12b) we obtain P c = P c = 0, thus Q c = λ + 2 µ and Q c = λ . iii. Shear test.Consider a shear test, with ε (cid:54) = 0, while other components are zero. Bysymmetry, we can deduce P t = P c , Q t = Q c , P t = P c , Q t = Q c , P t = P c , Q t = Q c . (13)Per (13), henceforth we will simply use P and Q to denote P m and Q m , respectively. However, in this case, it is not easy to determinea simple analytic expression to separate P m from Q m , also P m from Q m . In the sequel we will resort to fitting to determine these coefficients. Summary of results so far.
Up to now we have obtained the following resultsof the coefficients:(i) Q t = Q t = Q t = P c = P c = P c = P c = P c = 0,(ii) P t = Q c = Q c = Q c = λ + 2 µ ,(iii) P t = P t = Q c = Q c = λ , (iv) P + Q = P + Q = P + Q = µ .Then the unknown coefficients at this point are: P t , Q t (cid:0) = λ + 2 µ − P t (cid:1) , P t , Q t (cid:0) = λ − P t (cid:1) , P c , Q c (= λ − P c ) , P , Q (= µ − P ) , P , Q (= µ − P ) . Ψ ( ε , d ) = 12 { H ( β ( ε )) (cid:110) g ( d ) (cid:2) P t ε + ( λ + 2 µ ) ε + P t ε + 2 λε ε + 2 P t ε ε + 2 λε ε (cid:3) + ( λ + 2 µ − P t ) ε + ( λ + 2 µ − P t ) ε + 2( λ − P t ) ε ε (cid:111) + [1 − H ( β ( ε ))] (cid:110) g ( d ) P c ε ε + ( λ + 2 µ ) (cid:0) ε + ε + ε (cid:1) + 2 λε ε + 2( λ − P c ) ε ε + 2 λε ε (cid:111) + 2 g ( d ) (cid:0) P ε + P ε + P ε (cid:1) + 2( µ − P ) (cid:0) ε + ε (cid:1) + 2( µ − P ) ε . P t , P t , and P c According to the numerical experiments (see Figures A.26 and A.29), thecoefficient P t is determined by making P t ε + ( λ + 2 µ ) ε + P t ε + 2 λε ε + 2 P t ε ε + 2 λε ε a perfect square. This allows us to determine that P t = P t = λ λ + 2 µ , and that β ( ε ) should read β ( ε ) = λε + ( λ + 2 µ ) ε + λε , which is just σ for the undamaged material. An alternative expression for β ( ε ) is given in terms of the Poisson ratio, νε + (1 − ν ) ε + νε .It is also worth mentioning that according to numerical results, we obtain P c = 0, and hence Q c = λ . P and P by fitting It turns out that P does not seem to depend on the material parameterswith a closed-form analytic expression, see Figure A.28. Hence we fit the non-dimensional expression P /µ as a function of the phase field d , for various Poisson ratios. More precisely, we fit P /µ as a quadratic function of d for d ∈ [0 , ν .Here one constraint follows. When d = 1, ε (cid:54) = 0, ε = ε = ε = 0,we need to enforce σ = 0, otherwise the through-crack shear test would notgo through, see Section 4.5. Thus we need P = µ when d = 1. With thisconstraint, we write P = µ [ a ( ν )(1 − d ) + b ( ν )(1 − d ) + 1] . The coefficients a ( ν ) and b ( ν ) obtained from least square fits are given inTable 1. Note that for d close to 1, an analysis based on [22] may be able toyield more accurate expression. Also, according to the numerical results, we obtain: P = 0 , Q = µ. able 1: Fitted coefficients for P = µ [ a ( ν )(1 − d ) + b ( ν )(1 − d ) + 1]. Eachrow is fitted by least squares with four data points. ν a ( ν ) b ( ν ) Correlation coefficients0 . − . − . . − . − . . . − . .
35 0 . − . . . − . .
45 1 . − . In summary, the strain energy density Ψ is given as follows: Ψ ( ε , d ) = H ( β ( ε )) Ψ t + (cid:0) − H ( β ( ε )) (cid:1) Ψ c + Ψ s . (14)We then provide the expressions of relevant quantities in different variables. The model expressed with the x -axis normal to the crack Γ . In this case, β ( ε ) = λε + ( λ + 2 µ ) ε + λε ,Ψ t := 12( λ + 2 µ ) (cid:104) g ( d ) ( λε + ( λ + 2 µ ) ε + λε ) + 4 µ ( λ + µ ) (cid:0) ε + ε (cid:1) + 4 λµε ε (cid:105) ,Ψ c := 12 (cid:104) λ ( ε + ε + ε ) + 2 µ (cid:0) ε + ε + ε (cid:1)(cid:105) ,Ψ s := g s ( d ; ν )2 µ (cid:0) ε + ε (cid:1) + 2 µε , (15)where g s ( d ; ν ) = 1 + [ g ( d ) − (cid:0) a ( ν )(1 − d ) + b ( ν )(1 − d ) + 1 (cid:1) is named the shear degradation function . The model expressed in terms of invariants and pseudo-invariants of the straintensor.
To obtain an expression suitable for any crack orientation, we definetwo invariants I = tr ε = ε + ε + ε ,I = 12 (cid:104) (tr ε ) − tr (cid:0) ε (cid:1)(cid:105) = ε ε + ε ε + ε ε − ε − ε − ε , and two pseudo-invariants I = n · εn , I = n · ε n , n := ∇ d (cid:107)∇ d (cid:107) . Now we can express the terms of (14) in terms of these invariants and pseudo-invariants as: β ( ε ) := λI + 2 µI ,Ψ t := 12( λ + 2 µ ) (cid:104) g ( d ) ( λI + 2 µI ) + 4 µ ( λ + µ ) (cid:0) I + I − I − I (cid:1) + 4 λµ ( I − I I + I ) (cid:105) ,Ψ c := 12 (cid:2) λI + 2 µ (cid:0) I + 2 I − I − I (cid:1)(cid:3) ,Ψ s := g s ( d ; ν )2 µ (cid:0) I − I (cid:1) . (16)The derivation is given in Appendix B.With the above development, in any coordinate system, i.e., not necessarilyrelated to the crack orientation, the constitutive relation can be obtained as, inVoigt’s notation, σ ( ε , d ) = ∂Ψ∂ ε = ∂Ψ∂I + ∂Ψ∂I ε + ε ε + ε ε + ε − ε − ε − ε + ∂Ψ∂I n n n n n n n n n + ∂Ψ∂I ε n + 2 ε n n + 2 ε n n ε n n + 2 ε n + 2 ε n n ε n n + 2 ε n n + 2 ε n ε n n + ε n + ε n n + ε n n + ε n n + ε n ε n + ε n n + ε n n + ε n n + ε n n + ε n ε n + ε n n + ε n n + ε n n + ε n + ε n n . The corresponding moduli can be written as C ( ε , d ) = ∂ Ψ∂I T + ∂Ψ∂I − −
00 0 0 0 0 − + ∂ Ψ∂I n n n n n n n n n n n n n n n n n n T + ∂Ψ∂I n n n n n n n n n n n n n n n n n n n n + n n n n n n n n n n n n + n n n n n n n n n n n n + n .
17s now Ψ depends on ∇ d , the strong form for d is modified to be ∂Ψ∂d − ∇ · ∂Ψ∂ ( ∇ d ) + g c (cid:18) d(cid:96) − (cid:96) ∆ d (cid:19) = 0 in Ω , where ∇ · ∂Ψ∂ ( ∇ d ) = H ( β ( ε )) (cid:18) ∂Ψ t ∂I ∇ · ∂I ∂ ( ∇ d ) + ∂Ψ t ∂I ∇ · ∂I ∂ ( ∇ d ) (cid:19) + ∂Ψ s ∂I ∇ · ∂I ∂ ( ∇ d ) + ∂Ψ s ∂I ∇ · ∂I ∂ ( ∇ d ) , where the partial derivatives can be obtained via ∂I ∂ n = 2 εn , ∂I ∂ n = 2 ε n , and ∂ n ∂ ( ∇ d ) = 1 (cid:107)∇ d (cid:107) − (cid:107)∇ d (cid:107) ∇ d ⊗ ∇ d. Finally, in order to regularize the abrupt change of the crack normal where ∇ d is near zero, we introduce a replacement of I and I in the constitutiveexpressions in the form of ˆ I , = tanh (cid:0) αl |∇ d | (cid:1) I , , where α is a small positive number. Here we take α = 5 . × − .
4. Numerical examples
In this section, the proposed micromechanics-informed phase field model iscompared with a few existing models using several numerical examples. Theseexamples are categorized into two types. Sections 4.1 through 4.6 solve prob-lems in the macroscopic scale, including a uniaxial tension test, a compression test, a shear test, a symmetric bending test, and two through-crack tests. Inthese examples the proposed model is compared against earlier models, i.e.,the isotropic, V-D, and spectral models. These results show that the proposedmodel gives reasonable crack paths in all cases. Sections 4.7 and 4.8 comparethe constitutive response for a fixed phase field value among a number of mod- els. In these examples, the proposed model is compared with the SK model andthe SS models, as well as the earlier models. We will see that only the proposedmodel show a very close response to the benchmark, by construction.In the sequel, the marker
Isotropic refers to the original model proposedby [3],
V-D the volumetric-deviatoric decomposition model proposed by [4],
Spectral the spectral decomposition model proposed by [5], SK the modelproposed by [10], and SS1 and
SS2 the two models proposed by [8]. Finally,
Proposed refers to the proposed micromechanics-informed phase field model.Here we use an in-house code for the macroscopic simulations, while forgetting the benchmark solution with the RVE analysis in Sections 4.7 and 4.8, the software Abaqus is adopted. 18 .1. Uniaxial tension test
The first example is a benchmark simulation of a uniaxial tension test. Con-sider a square plate with a horizontal initial crack located at the middle heightfrom the left side to the center of the specimen.
The material constants are chosen as in Table 2.
Table 2:
Material parameters used in the tension and compression simulations λ (GPa) µ (GPa) g c (mJ / mm ) l (mm)121 .
15 80 .
77 2.7 40The boundary value problem and the mesh are illustrated in Figure 4. Thedomain is Ω = [0 , and the pre-existing crack is Γ = (0 , × { } . Weset d = 1 at Γ , leave the left and right edges traction free, and minimize (1) with u D ≡ u D e on the top edge [0 , × { } and u D = − u D e on the bottom edge [0 , × { } . For each load step we update the boundary condition as u D ← u D + ∆ u , where the displacement increment ∆ u is set to 0 . u D = 0 . u = 0 . Figure 4:
Schematic and mesh of the uniaxial tension test [unit: mm].Figure 5 compares the phase field profiles among different models at a certain displacement load. This load is prior to the one corresponding to the rapidpropagation of the crack. From Figure 5 it can be concluded that the crackstarts to propagate almost at the same displacement load. In addition, theresults of the isotropic, V-D, and spectral models are almost the same, while19he crack in the proposed model propagates a little bit faster. Next, in Figure
6, the results of the V-D model and the proposed model are compared.
Figure 5:
Phase field profiles obtained from the uniaxial tension test at u D =0 . Figure 6:
Phase field profiles of the uniaxial tension test from the V-D [(a)-(d)]and proposed [(e)-(h)] models at various load levels.After increasing the displacement continuously, it can be found that thecrack in the V-D model propagate faster, while the crack in the proposed modelpropagates slower, and finally both of them will be totally broken when thedisplacement reaches u D = 0 . The force-displacement curves of the three models are shown in the Figure7. Overall speaking, the compared models behavior very similarly for a uniaxialtension load.
The second example is a simulation to show the behaviors of the models under uniaxial compression. The material constants are given in Table 2.20 igure 7:
Force-displacement curves of the uniaxial tension test resulting fromthe V-D, spectral and proposed phase field models.21he mesh given in Figure 4(b) is used. The only difference is that the dis-placement condition is opposite. As mentioned previously, the isotropic modelis not expected to distinguish between tension and compression while the othersare.
The the phase field profiles resulting from different models at u D = 0 . because in a uniaxial loading, the deviatoric part of the strain is not zero. Incontrast, in the spectral and proposed models, the phase field value does notincrease all over the domain. Figure 8:
Phase field profiles obtained from the uniaxial compression test at u D = 0 . force decreases to 0 can be clearly observed. The spectral and proposed modelsdisplay almost the same stiffness, while the V-D model is less stiff, indicat-ing that some part of the solid is degraded, in accordance with the observablephase field profile in Figure 8(b), which shows a slight damage. Again, this isattributed to the non-zero deviatoric strain under a uniaxial compression load. The deformed configurations at u D = 0 . with respect to the line x = 500, as if the crack does not exist. Based on theabove results, we can conclude that within the load range, only the spectral andproposed models perform reasonably in the uniaxial compression test. This test is a classical benchmark problem which has been analyzed in, for example, [5]. The material constants of this test are shown in Table 3, whichare the same as those used in [5].The geometry, loads, and the mesh are illustrated in Figure 11. The domain Ω is the shaded part, which is contained in the box (0 , × (0 , igure 9: Force-displacement curves of the uniaxial compression test resultingfrom different models.
Figure 10:
Deformed mesh (with the displacement field 1000 times magnified)of the uniaxial compression test at u D = 0 . Table 3:
Material parameters used in the three-point bending simulations λ (GPa) µ (GPa) g c (mJ / mm ) l (mm)8 12 0.5 0.06231) with u D = − u D e on the central of top edge [3 . , . × { } , where u D ← u D +∆ u . The displacement increment ∆ u is set to 0 . u D ∈ [0 , . ∪ [0 . , . u = 0 . u D ∈ [0 . , . Figure 11:
Configuration and mesh of the three-point bending [unit: mm]We show the crack evolution for the V-D, spectral and proposed models in
Figures 12 and 13. Figure 12(a)(e)(i) show the phase field profiles from differentmodels at the beginning of crack propagation and we observe some initial phasefield at equilibrium around the notch. Figure 12[(a)-(d)] and Figure 13[(a)-(d)]illustrate the phase field profiles of the bending tests resulting from the V-Dmodel. It can be observed that the phase field of the element around the point where the displacement is applied (the compression area) will increase. This isdue to the non-zero deviatoric strain in a uniaxial compression state, and thiswas also reported in [23]. On the other hand, Figures 12 and 13 show the phasefield results from the spectral model and the proposed model, respectively, bothof which are reasonable. And from Figure 13, we can see that when u D = 0 . mm, the plate is almost fully broken for both cases.Figure 14 shows the force-displacement curves resulting from the simulationsfor the V-D, spectral, and proposed models. The results conform to what weobserved in the phase field profiles.It can be concluded that only the spectral and proposed models give accept- able results for the three-point bending test. Now we investigate a square plate of a different size for a shear test. Thematerial constants are chosen as shown in Table 4.
Table 4:
Material parameters used in the shear simulations λ (GPa) µ (GPa) g c (mJ / mm ) l (mm)121 .
15 80 .
77 2.7 3.12524 igure 12:
Phase field profiles of the three-point bending test from the V-D[(a)-(d)], spectral [(e)-(h)], and proposed models[(i)-(l)] for u D in the range of[0 . , . igure 13: Phase field profiles of the three-point bending test from the V-D[(a)-(d)], spectral [(e)-(h)] and proposed models [(i)-(l)] for u D in the range of[0 . , . igure 14: Force-displacement curves of the three-point bending test for threemodels. 27he boundary value problem and the mesh are illustrated in Figure 15. The domain has a pre-existing crack Γ = (0 , × { } , and Ω = [0 , \ Γ , asshown in Figure 15(a). We minimize (1) with u D ≡ u D e on the top edge[0 , × { } , u D ≡ − u D e on the bottom edge [0 , × { } , u D ← u D + ∆ u and the left and right edges and Γ traction free. The displacement increment∆ u is set to 0 . u D ∈ [0 , . ∪ [0 . , . follow the crack propagation, the increment is reduced to ∆ u = 0 . u D ∈ [0 . , . Figure 15:
Configuration and mesh of the shear test [unit: mm].The crack propagation paths predicted by different models are quite differentin this test. In this case, we will only compare the results between V-D, thespectral and the proposed models since the isotropic model does not distinguish between tension and compression, resulting in unphysical branching crack paths,which has been reported by Ziaei-Rad and Shen [23].Figure 16 shows the phase field profiles of the shear test from V-D, thespectral and proposed model for u D ∈ [0 . , . the crack patterns in the spectral and the proposed model are more similar andhave a much larger angle. Also, the crack in V-D model run through the materialat u D = 0 . material response by the spectral model. [24] also reported this phenomenon.Figure 17 shows the subsequent crack propagation of the spectral and proposedmodels.The red lines in Figure 17(d)(h) show the initial crack paths given by theclassical mixed-mode loading result from the maximum circumferential stress28 igure 16: Phase field profiles of the shear test from the V-D [(a)-(d)], spectral[(e)-(h)] and proposed [(i)-(l)] models till u D = 0 . igure 17: Phase field profiles of the shear test from the spectral [(a)-(d)] andproposed [(e)-(h)] models till u D = 0 . Figure 18:
Force-displacement curves of the shear test resulting from differentmodels. 30riterion [25], which gives the crack extension direction according to θ c = 2 arctan K I K II ± (cid:115)(cid:18) K I K II (cid:19) + 8 , where θ c is the deviation angle from the extension line of the current crack. Forthis case the stress intensity factors K I = 0, K II >
0, hence θ c = 2 arctan( √ / ≈ ◦ .Figure 18 shows the force-displacement curves for the shear test. Again, inthis test, the V-D model differs from the spectral and proposed models. One of the criticisms [see [9]] that the spectral model [5] receives is its inabil- ity to output a reasonable result in the through-crack shear test. Hence here wewill subject the proposed model and other models to this test. Consider a squareplate with a horizontal through crack which is placed at the middle height ofthe specimen, then a shear load is applied at the top and bottom of the plate.As the crack is through, the domain is actually made of two independent parts with a smooth interface. Hence the correct response is that the two parts sliderelative to each other. Note that our model is designed to perform reasonablyby construction, see Section 3.5.The material constants for the simulations are provided in Table 5. Theconfiguration and mesh are shown in Figure 19. The domain is Ω = [0 , , and the pre-existing crack Γ = [0 , × { } is described by setting the phasefield d = 1 on Γ . We minimize (1) with u D = u D e on the top edge [0 , ×{ } , u D = − u D e on the bottom edge [0 , × { } where u D = 0 . ∂Ω traction free. Table 5:
Material parameters used in the through-crack shear simulations. λ (GPa) µ (GPa) g c (mJ / mm ) l (mm)121 .
15 80 .
77 2.7 3.125Figure 20 shows the deformed configuration of the V-D, spectral, and pro- posed models, with the displacement field magnified. The results from the V-Dand spectral models are consistent to those in [9]. In contrast, the spectralmodel shows an erroneous response, which can be seen by setting ε = ε = 0and d = 1, which still yields a non-zero σ , resulting in some remnant stiffness.It can be concluded that only the V-D and the proposed models give reasonable result for the through-crack shear test among the three compared. In order to make sure the result of the through-crack test is valid for differentload combinations, we will subject the proposed model and other models to the31 igure 19:
Configuration and mesh for the two through-crack tests [unit: mm].
Figure 20:
Deformed mesh (with the displacement field 1000 times magnified)of the through-crack shear test for different models.32ollowing test. Using the configuration and mesh shown in Figure 19, we fix thebottom of the sample, leave the left and right sides traction free, and imposethe following Dirichlet boundary condition on the top edge: u D = ∆ u sin θ e + ∆ u (1 + cos θ ) e , θ ∈ { , π/ , π/ , π/ , π } , (17)where ∆ u = 0 . θ ’s. The results of the V-D model are very similar, and are omitted for the sake of brevity. Figure 21:
Deformed configurations for the proposed model (with the displace-ment field 1000 times magnified) of the through-crack test with a circular loadpath (17).
In order to compare the proposed model with models involving the localcrack orientation, i.e., the SS1, SS2, and SK models, we perform a shear testfor a single RVE with the shear load. The configuration and mesh are shown in
Figure 22. The material parameters are shown in Table 6.The RVE domain is B = [0 , . Besides the response of the comparedmodels, we also provide the RVE analysis result from Abaqus for reference. Notethat the Abaqus result is based on a calculation with a much finer resolutionand which exactly takes into account the self-contact condition due to the crack, igure 22: Configuration and mesh of the RVE used for the 2D shear test.
Table 6:
Material parameters used in the 2D and 3D constitutive response ofshear tests λ (GPa) µ (GPa) r a . . ε = ε ( e ⊗ e + e ⊗ e ) where ε ∈ [0 , . × − ].To interpret the results to come, it is helpful to see the simplified expressionsof the models under this loading. The isotropic, V-D, SS1, and SS2 models aresimplified to σ = g ( d ) 2 µε . The spectral model is simplified to σ = [ g ( d ) + 1] µε . The SK model is simplified to σ = g ( d ; b ) 2 µε . The proposed model is simplified to σ = g s ( d ; ν ) 2 µε . In Figure 23 we show the curves of σ / (2 µε ) vs. d , illustrating the effective shear degradation. Two observations can be made. First, the proposed model isthe closest to the benchmark, which is no surprise, as the coefficients in g s ( d ; ν )34re obtained by fitting the benchmark solution. Second, it is clear that σ ofthe spectral model does not degrade to 0 when d = 1, which explains why thespectral model cannot pass the through crack shear test, see Figure 20(b). Figure 23:
Effective degradation curves σ / (2 µε ) vs. d for the 2D shear test. Abaqus indicates the response of the more costly benchmark solution obtainedwith Abaqus.
In the last group of tests, a 3D combined shear and compression load isapplied. The configuration and mesh are shown in Figure 24. The macroscopicstrain is ε = ε ( e ⊗ e + e ⊗ e ) + ε e ⊗ e where ε ∈ [0 , . × − ], ε = − × − . The RVE domain is B = [0 , , and the material parameters are givenin Table 6. The same as Section 4.7, the Abaqus result takes into account theself-contact condition due to the crack, while the response from the phase fieldmodels are obtained just by a simple function evaluation.The effective degradation curves σ / (2 µε ) vs. d are shown in Figure 25. The proposed model is still the closest to the benchmark Abaqus in this casewith mixed three dimensional loads, again, by construction.35 igure 24:
Configuration and mesh of the RVE used for the 3D combinedshear and compression test.
Figure 25:
Effective degradation curves σ / (2 µε ) vs. d for the 3D combinedshear and compression test. Abaqus indicates the response of the more costlybenchmark solution obtained with Abaqus.36 . Conclusions
We have proposed a micromechanics-informed phase field model whose macro-scopic constitutive relationship is fully determined by the fictitious microstruc- ture. This model possesses the following features:1. The proposed model is based on three shapes of micro-cracks in the RVE.Yet the results from these shapes coincide, once the phase fields are prop-erly calibrated.2. The proposed model is the only one that outputs an accurate stiffness response and crack path among the models compared.3. The model can be expressed in terms of invariants and pseudo-invariantsof the strain tensor, and hence easy to be applied to an arbitrary crackorientation.4. The methodology involved in developing the present model can be general- ized to more complicated material constitutive relations, such as anisotropyor inelasticty, to be combined with the phase field description of fracture.
Appendix A. Details of the numerical experiments on the RVE
In this appendix, details of the numerical solution of (4) are provided. Thisis based on Abaqus/CAE 6.14-2.
Configuration.
We test three types of RVEs: (i) 2D plane strain RVE witha straight crack at the center; (ii) 3D RVE with a penny-shaped crack at thecenter; (iii) 3D RVE with a square-shaped crack with rounded fillets ( r = 0 . L )at the center. The crack is always perpendicular to the x -axis. Mesh.
We use a uniform triangular mesh for the plane strain case for which the element type is
CPS3 (continuum, plane strain, 3-nodes). And we use auniform tetrahedral mesh for the 3D cases for which the element type is
C3D4 (continuum, 3 dimensional, 4-nodes).
Contact condition.
For the frictionless contact conditions, within the contactoption
Tangential Behavior , we set the
Friction Coeff to be 0. This inter- action condition is applied to the two sides of the crack Γ in RVE. Macroscopic stress.
The macroscopic stress on the RVE is taken as the volu-metric average: σ = 1 |B| (cid:90) B s σ d V. ne-to-one relation of the phase field d and the crack length ratio r a . We corre-late the phase field d with the crack length ratio r a by calibrating the responseof the RVE with uniaxial tensile loads. More precisely, we apply a series ofuniaxial tensile loading along the x -direction, in which case it is assumed thatthe degradation function g ( d ) = (1 − d ) is applied to the entire strain energydensity of the unbroken material, i.e., Ψ ( ε , d ) = (1 − d ) Ψ ( ε ) . Substituting ε = ε = 0 into (8b) yields12 g ( d ) ( λ + 2 µ ) ε = 12 C t ε , ∀ ε y > . (A.1)With (A.1), the phase field d can then be calculated as d = 1 − (cid:115) C t λ + 2 µ , where C t can be looked up from Figure A.27.The phase field d and crack length ratio r a at λ = 1 . µ = 0 . Table A.7:
Phase field d at different crack length ratios r a . r a d (plane strain) d (penny-shaped) d (square-shaped)0 0 0 00 . . . . . . . . Effective secant moduli.
With both macroscopic strain and stress known, theeffective secant moduli can be calculated using the following relationship: σ ij = C ijkl ε kl . The numerical results of the effective secant moduli are verified to be inde- pendent of ε for a given d and given signs of β ( ε ). The values are shown inFigures A.26, A.27, A.28, A.29, and A.30. Appendix B. Deriving the model in terms of invariants
The key step in the derivation is to express the strain components in termsof the (pseudo-)invariants. We first choose a very special coordinate system38 igure A.26:
Effective modulus C t vs. d .39 igure A.27: Effective modulus C t vs. d .40 igure A.28: Effective modulus C t vs. d .41 igure A.29: Effective modulus C t vs. d .such that the x -axis is normal to the crack Γ , and that ε = 0. In this specialcase, it is easy to deduce that ε = I , ε + ε = I − I , ε + ε = I − I , ε ε = I + I − I I . Applying these expressions in (15) yields (16).As a sanity check, we now still take the x -axis normal to Γ , and ε to have all zero components except ε = ε . Then I = I = I = 0, and I = − ε .It can be shown that Ψ t = Ψ c = 2 µε while Ψ s = 0, giving rise to the expectedexpression, Ψ = 2 µε . Acknowledgments
This work is supported by the Shanghai Municipal Science and Technology
Commission, Grant No. 19ZR1424200, and by the National Natural ScienceFoundation of China, Grant No. 11972227.
References [1] A. A. Griffith, The phenomena of rupture and flow in solids, PhilosophicalTransactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 221 (582–593) (1921) 163–198.42 igure A.30:
Effective modulus C t vs. d .[2] G. A. Francfort, J.-J. Marigo, Revisiting brittle fracture as an energy min-imization problem, Journal of the Mechanics and Physics of Solids 46 (8)(1998) 1319–1342.[3] B. Bourdin, G. A. Francfort, J.-J. Marigo, Numerical experiments in revis- ited brittle fracture, Journal of the Mechanics and Physics of Solids 48 (4)(2000) 797–826.[4] H. Amor, J.-J. Marigo, C. Maurini, Regularized formulation of the varia-tional brittle fracture with unilateral contact: numerical experiments, Jour-nal of the Mechanics and Physics of Solids 57 (8) (2009) 1209–1229. [5] C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: Robust algorithmic implementation basedon operator splits, Computer Methods in Applied Mechanics and Engineer-ing 199 (45) (2010) 2765–2778.[6] F. Freddi, G. Royer-Carfagni, Regularized variational theories of fracture: A unified approach, Journal of the Mechanics and Physics of Solids 58 (8)(2010) 1154–1174.[7] M. Ambati, T. Gerasimov, L. De Lorenzis, A review on phase-field mod-els of brittle fracture and a new fast hybrid formulation, ComputationalMechanics 55 (2) (2015) 383–405. [10] C. Steinke, M. Kaliske, A phase-field crack model based on directionalstress decomposition, Computational Mechanics 63 (5) (2019) 1019–1046.[11] J.-Y. Wu, V. P. Nguyen, H. Zhou, Y. Huang, A variationally consistentphase-field anisotropic damage model for fracture, Computer Methods inApplied Mechanics and Engineering 358 (2020) 112629. [12] C. Cheng, A micromechanics-based phase field model for fracture, Master’sthesis, Shanghai Jiao Tong University (2018).[13] C. Cheng, Y. Shen, A micromechanics-based phase field approach to frac-ture, in: The Third International Conference on Damage Mechanics, 2018.[14] C. Dascalu, G. Bilbie, E. Agiasofitou, Damage and size effects in elastic solids: a homogenization approach, International Journal of Solids andStructures 45 (2) (2008) 409–430.[15] J. Storm, D. Supriatna, M. Kaliske, The concept of representative crack ele-ments for phase-field fracture: Anisotropic elasticity and thermo-elasticity,International Journal for Numerical Methods in Engineering 121 (5) (2020) [17] J. Lemaitre, A Courese on Damage Mechanics, Springer-Verl, Berlin, 1992.[18] M. Ortiz, E. Popov, Accuracy and stability of integration algorithms forelastoplastic constitutive relations, International Journal for NumericalMethods in Engineering 21 (9) (1985) 1561–1576.[19] J. Simo, J. Ju, Strain- and stress-based continuum damage modelsi. formu- lation, International Journal of Solids and Structures 23 (7) (1987) 821–840.[20] S. Yazdani, H. Schreyer, An anisotropic damage model with dilatation forconcrete, Mechanics of Materials 7 (3) (1988) 231–244.[21] V. A. Lubarda, D. Krajcinovic, S. Mastilovic, Damage model for brittleelastic solids with unequal tensile and compressive strengths, Engineering
Fracture Mechanics 49 (5) (1994) 681–697.4422] X. Markenscoff, C. Dascalu, Asymptotic homogenization analysis for dam-age amplification due to singular interaction of micro-cracks, Journal of theMechanics and Physics of Solids 60 (8) (2012) 1478–1485.[23] V. Ziaei-Rad, Y. Shen, Massive parallelization of the phase field formulation for crack propagation with time adaptivity, Computer Methods in AppliedMechanics and Engineering 312 (2016) 224–253.[24] V.-D. Nguyen, E. B´echet, C. Geuzaine, L. Noels, Imposing periodic bound-ary condition on arbitrary meshes by polynomial interpolation, Computa-tional Materials Science 55 (2012) 390–406. [25] F. Erdogan, G. C. Sih, On the crack extension in plates under plane loadingand transverse shear, Journal of Basic Engineering 85 (4) (1963) 519–525.[26] C. Dascalu, B. Fran¸cois, O. Keita, A two-scale model for subcritical damagepropagation, International Journal of Solids and Structures 47 (3) (2010)493–502.560