Phase field modelling of anisotropic crack propagation

Anisotropy is inherent to crystalline materials (among others) due to the symmetry of the atomic lattice. However, failure anisotropy is questioning the foundations of brittle failure as the equivalence between the principle of local symmetry and the maximum energy release rate criterion is no longer valid. Many experimental observations have been reported in the literature but anisotropic failure is thus still an open path for fun-damental research. The aim of the paper is to propose a phase ﬁeld model that could reproduce (energetically) non-free anisotropic crack bifurcation within a framework allowing for robust and fast numerical simulations. After the model and its ﬁnite element implementation have been detailed, its ability to capture the thought phenomenon is illustrated through several examples.


Introduction
Anisotropy is inherent to crystalline materials (among others) due to the symmetry of the atomic lattice. Depending on the level of symmetry, one can derive the symmetry class for elasticity or other physical properties. Concerning failure, anisotropy is often higher than for elasticity for a given symmetry class. However, failure anisotropy is questioning the foundations of brittle failure as the equivalence between the principle of local symmetry and the maximum energy release rate criterion is no longer valid.
Email address: julien.rethore@ec-nantes.fr (Marie-Christine Baietto) 1 Now at GEM, CNRS UMR 6183 CNRS / Ecole Centrale de Nantes / Université de Nantes crack paths which are in qualitative agreement with experimental observations. In their approach that uses a second order formulation, crack curvature is penalized. However, there is no direct link between the model parameters and the energy cost for crack bending. This is an important point as experimental evidences were reported in Takei et al. (2013).
The aim of the paper is to propose a phase field model that could reproduce (energetically) non-free anisotropic crack bifurcation. For this purpose, we choose to extend the formalism proposed by Clayton and Knap (2015) within a framework allowing for robust and fast numerical implementation. The outline of the paper is as follows: first we recall in Section 2 the useful regularized representation of cracks that is then used in the phase field model we develop in Section 3. The finite element implementation of the model is then detailed in Section 4 before its ability for modelling anisotropic failure are illustrated in Section 5.

General statements
Let Ω ⊂ R D an open domain describing a cracked solid, with D the space dimension and ∂Ω its boundary. Let Γ a curve of dimension D − 1 within Ω (see Fig. 1). In a regularized framework, the crack geometry is approximated by a smeared representation defined by a scalar parameter d(x), x ∈ Ω, taking a unit value on Γ and vanishing away from it. It can be shown (see e.g. Miehe et al. (2010)) that such a function can 2 be determined by solving the following boundary value problem on Ω: where ∆(.) is the Laplacian, is a regularization parameter describing the actual width of the smeared crack, and n the outward normal to ∂Ω. A two-dimensional illustration of this concept is depicted in Fig. 1(b). It can be shown that (1) is the Euler-Lagrange equation associated with the variational problem: represents the total crack length. In (3), γ(d, ∇d) denotes the crack density function per unit volume, defined by:

Anisotropy
In the anisotropic case, we extend this formulation to a class of anisotropic materials. An anisotropic crack surface density function is written by the following expression: Where ω ω ω is a second order structural tensor, being invariant with respect to rotations (characterizing the kind of the material's anisotropy). The terminology anisotropy is used herein to be consistent with the literature. However, this terminology may be 3 confusing as one has to distinguish between elastic anisotropy and failure anisotropy.
In the latter case, the term directionality would be more appropriate. One usually distinguishes between weak and strong anisotropy to failure. This transition is obtained when the angular distribution of the reciprocal surface energy (1/γ) becomes non-convex. This non-convexity induces forbidden direction for the crack to propagate whereas convex reciprocal surface energy does not.
As considered in Clayton and Knap (2015), to make the energy release rate orientation dependent, the tensor ω ω ω can be defined by: where M denotes the unit vector normal to the preferential cleavage plane (with respect to the material coordinates), β 0 is used to penalize the damage on planes not normal to M. Hence, in the case of isotropic material we set β = 0.
In order to estimate the anisotropy introduce in the surface energy by this formulation, we compute the distribution of surface energy g c γ (d, ω ω ω) for a phase field d(x) is determined numerically by the following variational problem: A Benchmark problem described in Fig. 2 is considered. The phase field variable d is assigned a unit value within a disc of radius r = 0.025B at the center of a square domain of size B. Once the computation of d is performed, the corresponding surface energy distribution is calculated at the integration points. For M = [1 0], g c = 1 J/m 2 and β between 0 and 20, polar plots of element average value of γ and 1/γ at a distance r = 0.025B + 2 from the domain center are depicted in Fig. 3. The fluctuation of the surface energy on these plots is due to the position of the integration element center which is not exactly located at a distance r from the center. For β = 0, the surface 4 energy does not depend on the orientation and isotropy is recovered. When β is increased from 0 to 20, there is a strong modification of the surface energy distribution as illustrated in Fig. 3(a). Namely, when beta increases the surface energy in the [1 0] orientation decreases whereas it remains constant in the [1 0] orientation. In Fig. 3(b), the sensitivity of the reciprocal surface energy to β is analyzed. It is clearly shown that for β ≤ 1 the reciprocal surface energy distribution is convex but it becomes nonconvex for β > 1. A smooth transition from weak to strong anisotropy can thus be modeled within the proposed framework.
In the case of the cleavage planes of multiple discrete orientations, we introduce the multiple phase field d d d i to quantify the damage accumulation on each such plane d i .
Thus the total crack length is here rewritten by the following: with γ i (d i , ∇d i , ω ω ω i ) is the crack density function of damage system d i corresponding to orientation dependent ω ω ω i .
The variational derivative of the crack density function for each damage system is now defined as: where ω ω ω i = [ω ω ω 1 , ω ω ω 2 , ...ω ω ω n ] is calculated from Equation 6 corresponding to each preferential cleavage plane (depending on the unit vector normal M i and coefficient β i ). Note that, the regularization parameter can be taken different for each plane.
However, in this work we assume is homogeneous.

Regularized variational framework
The variational approach to fracture mechanics provided by Francfort and Marigo (1998) introduces the following energy functional for the cracked body: where W u is the strain energy density function, ε ε ε = 1 2 ∇u + ∇u T , u is the displacement field, g c is the fracture toughness, and H D−1 is the Hausdorff surface measure giving the crack length (D = 2) or surface (D = 3). The term E u (u, Γ) represents the elastic energy stored in the cracked body, and E s (Γ) is the energy required to create the crack according to the Griffith criterion. Then, the state variables are the displacement field u and the geometry of the crack Γ. In a regularized framework (phase field method), fracture energy is regularized by the crack density function γ i (d i , ∇d i , ω ω ω i ) for each damage variable d i and the strain energy is replaced by energy of damageable material W u (ε ε ε(u), d d d i ), the above functional is substituted by the functional: The total energy is then rewritten as E = Ω W dΩ in which

Unilateral contact formulations and strain criterion with threshold for anisotropic elastic energy
In order to prevent the issue of cracks interpenetration in compression mode, many unilateral contact formulations have been proposed in the literature. For more details and practical implementation aspects, the interested reader may refer to Freddi and Royer-Carfagni (2010).
In what follow, we will present some points of view about the unilateral contact formulations for anisotropic elastic energy proposed in the work by Clayton and Knap (2015). This model has been proposed for the first time in the work of Ramtani et al. (1992), and applied to predict the behavior of rocks up to failure in the work of Comi (2001). This idea has been used in the regularized framework for brittle fracture of Amor et al. (2009) for isotropic case. The strain is decomposed into deviatoric ε ε ε dev and spheric ε ε ε sph parts. Then, it is assumed that damage is created by expansion (positive part of trace of the strain) and shear.
Where C 0 denote the initial elastic tensor of the material. The degradation function (13) is assumed to have the simple form: The function g(d d d i ) has been chosen such that g (d i = 1) = 0 to guarantee that the strain energy density function takes a finite value as the domain is locally cracked (see e.g. Braides (1998)) and g(d d d i = 0) = 1 to guarantee that the material is initially undamaged, g(d i = 1) = 0 is the limit for a fully damaged material. The quadratic function is chosen here (1 − d i ) 2 , that is the simplest case to ensure the existence of a solution regular in the sense of Carfagni. Alternatively the quartic function, and the cubic function are introduced in the work of Karma et al. (2001), Borden (2012). The small parameter k << 1 is introduced to maintain the well-posedness of the system for partially broken parts of the domain.
By introducing bulk modulus k 0 for the undamaged material (relating the spherical 7 of the strain to the spherical part of the stress), the elastic tensor now written as: Where the sign function sign − (x) = 1 if x < 0 and sign − (x) = 0 if x ≥ 0. The strain energy is now rewritten by: In order to prevent the problem of a damage-type occurring at low stress levels, a thresholded energy function is used. We seek to rewrite (12), so that the damage will only appear when the strain energy rather than an initial damage threshold ψ c , but the elastic behavior of material is still valid in classical way. Without loss of generality, we chosen ψ c = g c 2 , the total density energy defined in (12) can be re-written as: This formulation verifies that for d d d i → 0 (elastic state) which includes g(d d d i ) → 1,

Basics of thermodynamics and evolution of phase field
In this subsection we will formulate a crack phase field evolution law, that can guarantee the irreversibility of the process. Assuming isothermal process and without the external mircoforces, a reduced form of the Clausius-Duhem inequality can be written as: where A i is the variational derivative of W with respect to the phase field d i , reads: 8 Assuming that, the damage systems are independent of each other, then the condition (18) becomes: At this stage, a threshold F(A i ) such that no damage occurs satisfied the following condition: Assuming the principle of maximum dissipation then requires the dissipation Aḋ to be maximum under the constraint (21). Using the method of Lagrange multipliers and the following Lagrangian: yields the Kuhn-Tucker equations: The first equality in (23) gives:ḋ Without loss of generality, the threshold function F i (A i ) is assumed in the form From (24) and using the second inequality in (23), we obtain: 9 For each damage system, whenḋ i > 0, and from (20), (25) and the third equality in (23) , which give F i = 0 implying: where (9).
From (15), (16) and (26), we have with Expressing the variation of crack length: we can check that due to (29) a non-reversible evolution of cracks is satisfied. In addition, to ensure the positive condition of (29) and to handle loading and unloading histories, Miehe et al. (2010) introduced the strain history functional: Note that, this functional is initiated H i (x,t = 0) = 0. Therefore, following (32) the history functional will be equal zero when the strain energy bellow threshold ψ c , that prevent damage occurring at low stress levels. By substituting H i (x,t) to first term in (27), and using (1), (9) It yields the following phase field problem to be solved to evaluate the field d i (x,t) at time t:

Weak form of phase field problem
Starting from (33) 1 , multiplying by a test function δd i and integrating over Ω, we obtain: Using the property of the product rule used to find the derivatives of products: and the divergence theorem, we finally obtain the formulation of phase field problem for each preferential cleavage plane:

Weak form of displacement problem
The weak form associated with the displacement problem is found by solving the variational problem: where S u = u|u(x) =ū on ∂Ω u , u ∈ H 1 (Ω) and W ext = Ω f · udΩ + ∂Ω F F · udΓ with f and F body forces and prescribed traction over the boundary ∂Ω F . We obtain the classical weak form for u(x) ∈ S u : where the Cauchy stress σ σ σ = ∂W u ∂ε ε ε is given using (16) and (15), by:

Time-Discrete field variables in an incremental setting
In the present work, the computations are performed in quasi-static conditions. Then, the time steps introduced in the following actually refer to load increments. Introducing a time stepping, the phase field problem for preferential plane i to be solved at time t n+1 is expressed by seeking d i (x) ∈ S d i , such that: The Cauchy stress at time t n+1 is now defined by: The sign function sign − (tr ε ε ε(u n+1 )) in Equation 41 introduces non-linearity when solving the variational problem related to displacement field. In order to recast the latter as a linear problem, we use the shifted algorithm proposed in Nguyen et al.
(2015b) by introducing the prediction R − , which is determined from results of previous time step t n :

Overall algorithm
The overall algorithm is described as follows: Initialization. Initialize the displacement field u 0 (x), all phase field systems 3. Compute displacement field u n+1 (x): • Determine the prediction R − by Equation (42) • Solving the linear displacement problem (38).

Numerical model
In the simulations, the effective continuum material is considered. The effective elastic tensor is obtained by using a periodic homogenization scheme. For the 0 o case 14 the elastic tensor is: to induce a mode I opening solicitation at the crack tip.

Results
Two orientations of the material are tested as in the experiments (see Figure 6. The results are illustrated in Figure 8. In this Figure, the phase field variables are plotted at the end of the simulation. For the two orientations, during the nucleation stage both cleavage planes absorb part of the elastic energy available in the system. Then, crack initiation occurs in one of the two planes. As soon as a crack initiates, the dissipation is concentrated in the activated cleavage plane only. The phase field variable for the unactivated plane remains at the same level during the rest of the failure process. As expected the failure patterns follow one of the orientations allowed in the model (see Figure 5). For the two orientations, the selected plane in the simulation gives the same crack orientation as in the experiments.

Zigzag crack
In this example, a benchmark from the literature is reproduced. The idea of this example is to guide the crack in a direction that is forbidden for the material as observed experimentally by Wu et al. (1995) and Takei et al. (2013). In Figure 9, the evolution of the equivalent phase field variable d eq is plotted for β = 20, d eq being defined as . In this Figure, it is observed that the crack initiates in one of the direction allowed by the model and then bifurcates to the other direction before its reaches the fully constrained zone. Propagation occurs in this direction until the crack almost reached the other constrained region and then it bifurcates again. This leads to a zigzag crack pattern as depicted in Figure 9. In the formulation proposed herein, the crack results from the superposition of branches activating one of the two allowed cleavage planes , others orientations being forbidden due to strong anisotropy. This is illustrated in Figure 10, where the phase field variable of the two systems are plotted at the end of the simulation. One clearly observes how each of the cleavage plane is successively activated and unactivated to finally obtain the zigzag pattern. Subsequently, changing orientation means nucleating and initiating a new crack in an other plane what implies higher energy dissipation than propagating along the same plane. This is further illustrated by the global load v.s. displacement curve plotted in Figure 11. The red dots on this curve correspond to the snapshots in Figure 9 that are taken at the times when the crack starts a zigzag. It is clearly observed that after each of this dots a phase of lower softening (higher dissipation) starts before the propagation occurs along the new orientation.
The simulation is now carried out for a value of β = 1 inducing weak anisotropy.
In this case there is no forbidden direction for the crack to propagate. As illustrated in Figure 12, the crack runs along the direction expected from the vertical loading.
However, the crack line is slightly curved. Its initial orientation results from the com-petition between the external load and the weak anisotropy that induces a lower cost for the crack to grow in other orientation than those penalized by the model. After going along this initial orientation, the crack bends and tends to grow along the opposite orientation. It is observed in Figure 12 that this curvature is related to the progressive activation of the second cleavage plane.
From these results, it can be concluded that the model allows for capturing the effect of weak and strong anisotropy. Further, it seems that energy cost for the crack to bend is related to the initiation of damage along a different cleavage plane. In the case of weak anisotropy, this activation is progressive thus leading to low curvature whereas in the case of strong anisotropy a kink is obtained.

Conclusions
To handle failure anisotropy, the proposed model is based on the consideration of several cleavage planes, each of these planes having its own damage variable. A directionality effect is obtained by adding a penalty term in the phase field equations that prevent damage to develop along the normal to the considered cleavage planes. It is shown that the penalty coefficient allows to control the level of anisotropy and makes the model switch from weak anisotropy to strong anisotropy. An efficient numerical implementation is derived for this formulation by using history variables associated to each of considered cleavage planes. Note that compared to the formulation proposed by Li et al. (2015), a standard finite element implementation is used herein.