# Phase field modelling of anisotropic crack propagation

Thanh Tung Nguyen, Julien Réthoré, Marie-Christine Baietto

► To cite this version:

Thanh Tung Nguyen, Julien Réthoré, Marie-Christine Baietto. Phase field modelling of anisotropic crack propagation. European Journal of Mechanics - A/Solids, 2017, 10.1016/j.euromechsol.2017.05.002. hal-01519161

HAL Id: hal-01519161

https://hal.science/hal-01519161

Submitted on 5 May 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Phase field modelling of anisotropic crack propagation

Thanh Tung Nguyen, Julien Réthoré<fnref n="1" />, Marie-Christine Baietto

LaMCoS, Université de Lyon / INSA Lyon / CNRS UMR 5259
Bat. Jacquard, 27 Avenue Jean Capelle, F-69621 Villeurbanne, Cedex, France

## Abstract

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 fundamental research. The aim of the paper is to propose a phase field 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 finite element implementation have been detailed, its ability to capture the thought phenomenon is illustrated through several examples.

**Keywords:** Crack propagation, Phase field, Anisotropy

## 1. 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)

<footnote n="1">¹Now at GEM, CNRS UMR 6183 CNRS / Ecole Centrale de Nantes / Université de Nantes</footnote>
Many experimental observations have been reported in the literature as for example in Wu et al. (1995). More recent theoretical investigations proposed by Hakim and Karma (2009) have motivated model experiments on thin films by Takei et al. (2013). Even more recently, Li et al. (2014) developed a numerical implementation of a phase field model with strong anisotropy including the explicit dependency of fracture energy to the local crack orientation. In their paper, they have been able to reproduce zigzag 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.

## 2. Regularized representation of free discontinuities

### 2.1. General statements

Let $\Omega \subset \mathbb{R}^D$ an open domain describing a cracked solid, with $D$ the space dimension and $\partial\Omega$ its boundary. Let $\Gamma$ a curve of dimension $D-1$ within $\Omega$ (see Fig. 1). In a regularized framework, the crack geometry is approximated by a smeared representation defined by a scalar parameter $d(\mathbf{x})$, $\mathbf{x} \in \Omega$, taking a unit value on $\Gamma$ and vanishing away from it. It can be shown (see e.g. Miehe et al. (2010)) that such a function can
be determined by solving the following boundary value problem on $\Omega$:

$$ egin{cases} d - \ell^2 \Delta d = 0 & \text{in } \Omega, \\ d(\mathbf{x}) = 1 & \text{on } \Gamma, \\ \nabla d(\mathbf{x}) \cdot \mathbf{n} = 0 & \text{on } \partial\Omega, \end{cases} \qquad (1) $$

where $\Delta(.)$ is the Laplacian, $\ell$ is a regularization parameter describing the actual width of the smeared crack, and $\mathbf{n}$ the outward normal to $\partial\Omega$. 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:

$$ d(\mathbf{x}) = \operatorname{Arg} \left\{ \inf_{d \in S_d} \Gamma_{\ell}(d) \right\}, \qquad (2) $$

with $S_d = \{d | d(\mathbf{x}) = 1 \text{ on } \Gamma \ \forall \mathbf{x} \in \Gamma\}$ and where

$$ \Gamma_{\ell}(d) = \int_{\Omega} \gamma(d, \nabla d) d\Omega \qquad (3) $$

represents the total crack length. In (3), $\gamma(d, \nabla d)$ denotes the crack density function per unit volume, defined by:

$$ \gamma(d, \nabla d) = \frac{1}{2\ell} d^2 + \frac{\ell}{2} \nabla d \cdot \nabla d. \qquad (4) $$

## 2.2. 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:

$$ \gamma(d, \nabla d, \mathbf{\omega}) = \frac{1}{2\ell} d^2 + \frac{\ell}{2} \mathbf{\omega} : (\nabla d \otimes \nabla d). \qquad (5) $$

Where $\mathbf{\omega}$ 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
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 $\omega$ can be defined by:

$$ \omega = 1 + \beta(1 - M \otimes M) \quad (6) $$

where $M$ denotes the unit vector normal to the preferential cleavage plane (with respect to the material coordinates), $\beta \gg 0$ is used to penalize the damage on planes not normal to $M$. Hence, in the case of isotropic material we set $\beta = 0$.

In order to estimate the anisotropy introduce in the surface energy by this formulation, we compute the distribution of surface energy $g_c\gamma_l(d, \omega)$ for a phase field $d(x)$ is determined numerically by the following variational problem:

$$ d(x) = \operatorname{Arg} \left\{ \inf_{d \in S_d} \int_{\Omega} \gamma(d, \nabla d, \omega) d\Omega \right\}, \quad (7) $$

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 \text{ J/m}^2$ and $\beta$ between 0 and 20, polar plots of element average value of $\gamma$ and $1/\gamma$ at a distance $r = 0.025B + 2\ell$ 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 $\beta = 0$, the surface
energy does not depend on the orientation and isotropy is recovered. When $\beta$ 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 $\beta$ is analyzed. It is clearly shown that for $\beta \le 1$ the reciprocal surface energy distribution is convex but it becomes non-convex for $\beta > 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_i$ to quantify the damage accumulation on each such plane $d_i$. Thus the total crack length is here rewritten by the following:

$$ \Gamma_{\ell}(d_i, \omega_i) = \sum_i \Gamma_{\ell}^i(d_i, \omega_i) = \sum_i \int_{\Omega} \gamma_i(d_i, \nabla d_i, \omega_i) d\Omega, \quad (8) $$

with $\gamma_i(d_i, \nabla d_i, \omega_i)$ is the crack density function of damage system $d_i$ corresponding to orientation dependent $\omega_i$.

The variational derivative of the crack density function for each damage system is now defined as:

$$ \delta\gamma_i(d_i, \nabla d_i, \omega_i) = \frac{d_i}{\ell} - \ell \left[ \nabla^2 d_i + \beta_i \left( \nabla^2 d_i - \mathbf{M}_i \otimes \mathbf{M}_i : \nabla \nabla d_i \right) \right] \quad (9) $$

where $\omega_i = [\omega_1, \omega_2, \dots, \omega_n]$ is calculated from Equation 6 corresponding to each preferential cleavage plane (depending on the unit vector normal $\mathbf{M}_i$ and coefficient $\beta_i$). Note that, the regularization parameter $\ell$ can be taken different for each plane. However, in this work we assume $\ell$ is homogeneous.
### 3. Phase field method

#### 3.1. Regularized variational framework

The variational approach to fracture mechanics provided by Francfort and Marigo (1998) introduces the following energy functional for the cracked body:

$$E(\mathbf{u}, \Gamma) = E_u(\mathbf{u}, \Gamma) + E_s(\Gamma) = \int_{\Omega \setminus \Gamma} W_u(\boldsymbol{\epsilon}(\mathbf{u})) d\Omega + g_c \mathcal{H}^{D-1}(\Gamma) \quad (10)$$

where $W_u$ is the strain energy density function, $\boldsymbol{\epsilon} = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$, $\mathbf{u}$ is the displacement field, $g_c$ is the fracture toughness, and $\mathcal{H}^{D-1}$ is the Hausdorff surface measure giving the crack length ($D = 2$) or surface ($D = 3$). The term $E_u(\mathbf{u}, \Gamma)$ represents the elastic energy stored in the cracked body, and $E_s(\Gamma)$ is the energy required to create the crack according to the Griffith criterion. Then, the state variables are the displacement field $\mathbf{u}$ and the geometry of the crack $\Gamma$. In a regularized framework (phase field method), fracture energy is regularized by the crack density function $\gamma_i(d_i, \nabla d_i, \omega_i)$ for each damage variable $d_i$ and the strain energy is replaced by energy of damageable material $W_u(\boldsymbol{\epsilon}(\mathbf{u}), \mathbf{d}_i)$, the above functional is substituted by the functional:

$$E(\mathbf{u}, \mathbf{d}_i) = \int_{\Omega} W_u(\boldsymbol{\epsilon}(\mathbf{u}), \mathbf{d}_i) d\Omega + g_c \sum_i \int_{\Omega} \gamma_i(d_i, \nabla d_i, \omega_i) d\Omega, \quad (11)$$

The total energy is then rewritten as $E = \int_{\Omega} W d\Omega$ in which

$$W = W_u(\boldsymbol{\epsilon}(\mathbf{u}), \mathbf{d}_i) + g_c \sum_i \gamma_i(d_i, \nabla d_i, \omega_i) \quad (12)$$

#### 3.2. 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 $\boldsymbol{\epsilon}_{dev}$ and spheric $\boldsymbol{\epsilon}_{sph}$ parts. Then, it is assumed that damage is created by expansion (positive part of trace of the strain) and shear.

$$ W_u(\boldsymbol{\epsilon}(\mathbf{u}), d) = \begin{cases} \frac{1}{2}g(\boldsymbol{d}_i)[\boldsymbol{\epsilon} : \mathbb{C}^0 : \boldsymbol{\epsilon}] & \text{if } \operatorname{tr} \boldsymbol{\epsilon} \ge 0 \\ \frac{1}{2}[\boldsymbol{\epsilon}_{sph} : \mathbb{C}^0 : \boldsymbol{\epsilon}_{sph}] + g(\boldsymbol{d}_i)[\boldsymbol{\epsilon}_{dev} : \mathbb{C}^0 : \boldsymbol{\epsilon}_{dev}] & \text{if } \operatorname{tr} \boldsymbol{\epsilon} < 0 \end{cases} \quad (13) $$

Where $\mathbb{C}^0$ denote the initial elastic tensor of the material. The degradation function $g(\boldsymbol{d}_i)$ in (13) is assumed to have the simple form:

$$ g(\boldsymbol{d}_i) = (1-k) \prod_i (1-d_i)^2 + k. \quad (14) $$

The function $g(\boldsymbol{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(\boldsymbol{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 \ll 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
of the strain to the spherical part of the stress), the elastic tensor now written as:

$$ \mathbb{C} = g(\mathbf{d}_i)\mathbb{C}^0 + k_0\mathbf{1} \otimes \mathbf{1} [1 - g(\mathbf{d}_i)] \text{sign}^{-}(\text{tr}\boldsymbol{\varepsilon}) \quad (15) $$

Where the sign function $\text{sign}^{-}(x) = 1$ if $x < 0$ and $\text{sign}^{-}(x) = 0$ if $x \ge 0$. The strain energy is now rewritten by:

$$ W_u(\boldsymbol{\varepsilon}(\mathbf{u}), \mathbf{d}_i) = \frac{1}{2} [\boldsymbol{\varepsilon} : \mathbb{C} : \boldsymbol{\varepsilon}] \quad (16) $$

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 $\psi_c$, but the elastic behavior of material is still valid in classical way. Without loss of generality, we chosen $\psi_c = \frac{g_c}{2\ell}$, the total density energy defined in (12) can be re-written as:

$$ W = W_u(\boldsymbol{\varepsilon}(\mathbf{u}), \mathbf{d}_i) - g(\mathbf{d}_i)\psi_c + \psi_c + 2\psi_c\ell \sum_i \gamma_i(d_i, \nabla d_i, \mathbf{\omega}_i) \quad (17) $$

This formulation verifies that for $\mathbf{d}_i \to 0$ (elastic state) which includes $g(\mathbf{d}_i) \to 1$, $\gamma_i(d_i, \nabla d_i, \mathbf{\omega}_i) \to 0$, we recover the classical form of elastic behavior of material. Note that the similar idea of this method can be found in the work by Frémond and Nedjar (1996); Pham et al. (2011); Miehe et al. (2015).

### 3.3. 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 microforces, a reduced form of the Clausius-Duhem inequality can be written as:

$$ \sum_i \mathcal{A}_i \dot{d}_i \ge 0 \quad (18) $$

where $\mathcal{A}_i$ is the variational derivative of $W$ with respect to the phase field $d_i$, reads:

$$ \mathcal{A}_i = -\frac{\delta W}{\delta d_i} = -\frac{\partial W}{\partial d_i} + \nabla \cdot \left( \frac{\partial W}{\partial \nabla d_i} \right) \quad (19) $$
Assuming that, the damage systems are independent of each other, then the condition (18) becomes:

$$ \left\{ \begin{array}{l} \mathcal{A}_1 \dot{d}_1 \ge 0 \\ \mathcal{A}_2 \dot{d}_2 \ge 0 \\ \dots \\ \mathcal{A}_n \dot{d}_n \ge 0 \end{array} \right. \qquad (20) $$

At this stage, a threshold $F(\mathcal{A}_i)$ such that no damage occurs satisfied the following condition:

$$ F(\mathcal{A}_i) \le 0. \qquad (21) $$

Assuming the principle of maximum dissipation then requires the dissipation $\mathcal{A}\dot{d}$ to be maximum under the constraint (21). Using the method of Lagrange multipliers and the following Lagrangian:

$$ \mathcal{L}_i = -\mathcal{A}_i \dot{d}_i + \lambda_i F_i(\mathcal{A}_i) \qquad (22) $$

yields the Kuhn-Tucker equations:

$$ \frac{\partial \mathcal{L}_i}{\partial \mathcal{A}_i} = 0, \quad \lambda_i \ge 0, F_i \le 0, \lambda_i F_i = 0. \qquad (23) $$

The first equality in (23) gives:

$$ \dot{d}_i = \lambda_i \frac{\partial F_i(\mathcal{A}_i)}{\partial \mathcal{A}_i}. \qquad (24) $$

Without loss of generality, the threshold function $F_i(\mathcal{A}_i)$ is assumed in the form $F_i(\mathcal{A}_i) = \mathcal{A}_i$. From (24) and using the second inequality in (23), we obtain:

$$ \dot{d}_i = \lambda_i \ge 0 \qquad (25) $$
For each damage system, when $\dot{d}_i > 0$, and from (20), (25) and the third equality in (23), which give $F_i = 0$ implying:

$$F_i = -\frac{\delta W}{\delta d_i} = -\frac{\partial W_u}{\partial d_i} - \frac{\partial (g(\mathbf{d}_i) \psi_c)}{\partial d_i} - 2\psi_c \ell \delta\gamma_i(d_i, \nabla d_i, \mathbf{\omega}_i) = 0, \quad (26)$$

where $\delta\gamma_i(d_i, \nabla d_i, \mathbf{\omega}_i)$ is defined in (9).

From (15), (16) and (26), we have

$$-\frac{\delta W}{\delta d_i} = 2(1-d_i)(1-k) \prod_{j \neq i} (1-d_j)^2 \left( \varepsilon : \mathbb{C}_h : \varepsilon - \psi_c \right) - 2\psi_c \ell \delta\gamma_i(d_i, \nabla d_i, \mathbf{\omega}_i) = 0, \quad (27)$$

with

$$\mathbb{C}_h = \frac{\mathbb{C}^0 - k_0 \mathbf{1} \otimes \mathbf{1} \operatorname{sign}^{-}(\operatorname{tr} \boldsymbol{\varepsilon})}{2} \quad (28)$$

As $(1-d_i)(1-k) \prod_j (1-d_j)^2 (\boldsymbol{\varepsilon} : \mathbb{C}_h : \boldsymbol{\varepsilon} - \psi_c) \ge 0$ and $\psi_c \ell > 0$, then

$$\delta\gamma_i(d_i, \nabla d_i, \mathbf{\omega}_i) \ge 0 \quad (29)$$

Expressing the variation of crack length:

$$\frac{d}{dt} \Gamma_{\ell}^{i} = \int_{\Omega} \delta\gamma_{i}(d_{i}, \nabla d_{i}, \mathbf{\omega}_{i}) \, \dot{d}_{i} d\Omega, \quad (30)$$

we can check that due to (29)

$$\dot{\Gamma}_{\ell}^{i} \geq 0, \quad (31)$$

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:

$$\mathcal{H}_i(\mathbf{x}, t) = \max_{\tau \in [0, t]} \left\{ (1-k) \prod_j (1-d_j)^2 \left( \boldsymbol{\varepsilon}(\mathbf{x}, \tau) : \mathbb{C}_h(\mathbf{x}, \tau) : \boldsymbol{\varepsilon}(\mathbf{x}, \tau) - \psi_c \right) \right\}. \quad (32)$$
Note that, this functional is initiated $\mathcal{H}_i(\mathbf{x}, t = 0) = 0$. Therefore, following (32) the history functional will be equal zero when the strain energy below threshold $\psi_c$, that prevent damage occurring at low stress levels. By substituting $\mathcal{H}_i(\mathbf{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(\mathbf{x}, t)$ at time $t$:

$$ \left\{ \begin{array}{l} (1 - d_i) \mathcal{H}_i - \psi_c \left( d_i - \ell^2 \left[ \nabla^2 d_i + \beta_i (\nabla^2 d_i - \mathbf{M}_i \otimes \mathbf{M}_i : \nabla \nabla d_i) \right] \right) = 0 \quad \text{in } \Omega, \\ d_i(\mathbf{x}) = 1 \quad \text{on } \Gamma, \\ \nabla d_i(\mathbf{x}) \cdot \mathbf{n} = 0 \quad \text{on } \partial \Omega. \end{array} \right. \qquad (33) $$

# 4. Finite element discretization

## 4.1. Weak form of phase field problem

Starting from (33)$_1$, multiplying by a test function $\delta d_i$ and integrating over $\Omega$, we obtain:

$$ \int_{\Omega} \left\{ (1 - d_i) \mathcal{H}_i \delta d_i - \psi_c \left( d_i - \ell^2 \left[ \nabla^2 d_i + \beta_i (\nabla^2 d_i - \mathbf{M}_i \otimes \mathbf{M}_i : \nabla \nabla d_i) \right] \right) \delta d_i \right\} d\Omega = 0. \qquad (34) $$

Using the property of the product rule used to find the derivatives of products:

$$ (\Delta d_i) \delta d_i = \nabla \cdot (\nabla d_i \delta d_i) - \nabla d_i \cdot \nabla (\delta d_i), \qquad (35) $$

and the divergence theorem, we finally obtain the formulation of phase field problem for each preferential cleavage plane:

$$ \int_{\Omega} \left\{ (\mathcal{H}_i + \psi_c) d_i \delta d_i + \psi_c \ell^2 \nabla d_i \left[ \mathbf{1} + \beta_i (\mathbf{1} - \mathbf{M}_i \otimes \mathbf{M}_i) \right] \nabla (\delta d_i) \right\} d\Omega = \int_{\Omega} \mathcal{H}_i \delta d_i d\Omega. \qquad (36) $$
## 4.2. Weak form of displacement problem

The weak form associated with the displacement problem is found by solving the variational problem:

$$ \mathbf{u}(\mathbf{x}) = \operatorname{Arg} \left\{ \inf_{\mathbf{u} \in S_u} (E(\mathbf{u}, \mathbf{d}_i) - W^{\text{ext}}) \right\} \quad (37) $$

where $S_u = \{\mathbf{u}|\mathbf{u}(\mathbf{x}) = \bar{\mathbf{u}} \text{ on } \partial\Omega_u, \mathbf{u} \in H^1(\Omega)\}$ and $W^{\text{ext}} = \int_{\Omega} \mathbf{f} \cdot \mathbf{u} d\Omega + \int_{\partial\Omega_F} \bar{\mathbf{F}} \cdot \mathbf{u} d\Gamma$ with $\mathbf{f}$ and $\bar{\mathbf{F}}$ body forces and prescribed traction over the boundary $\partial\Omega_F$. We obtain the classical weak form for $\mathbf{u}(\mathbf{x}) \in S_u$:

$$ \int_{\Omega} \boldsymbol{\sigma} : \boldsymbol{\varepsilon}(\delta \mathbf{u}) d\Omega = \int_{\Omega} \mathbf{f} \cdot \delta \mathbf{u} d\Omega + \int_{\partial \Omega_F} \bar{\mathbf{F}} \cdot \delta \mathbf{u} d\Gamma \quad \forall \delta \mathbf{u} \in H_0^1(\Omega), \quad (38) $$

where the Cauchy stress $\boldsymbol{\sigma} = \frac{\partial W_u}{\partial \boldsymbol{\varepsilon}}$ is given using (16) and (15), by:

$$ \boldsymbol{\sigma} = \left[ g(\boldsymbol{d}_i)\mathbb{C}^0 + k_0 \mathbf{1} \otimes \mathbf{1} \left[ 1 - g(\boldsymbol{d}_i) \right] \operatorname{sign}^{-}(\operatorname{tr} \boldsymbol{\varepsilon}(\mathbf{u})) \right] : \boldsymbol{\varepsilon}(\mathbf{u}) \quad (39) $$

## 4.3. 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(\mathbf{x}) \in S_{d_i}$, such that:

$$ \begin{aligned} & \int_{\Omega} \left\{ (\mathcal{H}_i^n + \psi_c) d_i^{n+1} \delta d_i^{n+1} + \psi_c l^2 \nabla d_i^{n+1} \left[ \mathbf{1} + \beta_i (\mathbf{1} - \mathbf{M}_i \otimes \mathbf{M}_i) \right] \nabla (\delta d_i^{n+1}) \right\} d\Omega \\ & = \int_{\Omega} \mathcal{H}_i^n \delta d_i^{n+1} d\Omega. \end{aligned} \quad (40) $$

The Cauchy stress at time $t^{n+1}$ is now defined by:

$$ \boldsymbol{\sigma}^{n+1} = \left[ g(\boldsymbol{d}_i^{n+1})\mathbb{C}^0 + k_0 \mathbf{1} \otimes \mathbf{1} \left[ 1 - g(\boldsymbol{d}_i^{n+1}) \right] \operatorname{sign}^{-}(\operatorname{tr} \boldsymbol{\varepsilon}(\mathbf{u}^{n+1})) \right] : \boldsymbol{\varepsilon}(\mathbf{u}^{n+1}) \quad (41) $$
The sign function $\text{sign}^{-}(\text{tr}\boldsymbol{\varepsilon}(\mathbf{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 $\mathcal{R}_-$, which is determined from results of previous time step $t^n$:

$$ \mathcal{R}_- = \text{sign}^{-}(\text{tr}\boldsymbol{\varepsilon}(\mathbf{u}^n)) = \begin{cases} 1 & \text{if } \text{tr}\boldsymbol{\varepsilon}(\mathbf{u}^n) < 0 \\ 0 & \text{if } \text{tr}\boldsymbol{\varepsilon}(\mathbf{u}^n) \ge 0 \end{cases} \quad (42) $$

### 4.4. Overall algorithm

The overall algorithm is described as follows:

**Initialization.** Initialize the displacement field $\mathbf{u}^0(\mathbf{x})$, all phase field systems $\mathbf{d}_i^0(\mathbf{x})$, and the strain-history functional for each systems $\mathcal{H}_i^0 = 0$.

**WHILE** $t^{n+1} \le T$, given $\mathbf{u}^n$, $\mathbf{d}_i^n$ and $\mathcal{H}_i^n$,

1. Compute the strain history functional for each damage systems $\mathcal{H}_i^{n+1}$ by (32).

2. Compute damage fields $\mathbf{d}_i^{n+1}(\mathbf{x})$ by solving the linear phase field problem (36).

3. Compute displacement field $\mathbf{u}_{n+1}(\mathbf{x})$:

• Determine the prediction $\mathcal{R}_-$ by Equation (42)

• Solving the linear displacement problem (38).

4. $(.)_n \leftarrow (.)_{n+1}$ and go to (1).

## 5. Examples

### 5.1. Crack propagation in a model material

In this first example, a comparison with experiments on a model honeycomb material obtained by additive manufacturing is proposed.
### 5.1.1. *Experimental observation*

Failure experiments on honeycomb materials have been performed. The specimen geometry is defined in Figure 4. The left arms of the specimen are fixed in the grips of a standard tensile device and submitted to a constant vertical displacement speed of 0.1mm/min. The specimens were obtained by a 3D printer. They are made from photo-sensitive ABS-type polymer powder. The bulk material obtained from this process is isotropic. Its elastic behavior is defined by a Young's modulus of 1.4 GPa and a Poisson's ratio of 0.4. A pre-crack has been designed by simply not printing some cells along the symmetry plane of the specimen. The results presented in this paper are based on the analysis of a material generated by a periodic tilling of the $D_2$ cell presented in Figure 5. This cell is invariant by rotation of $2\pi/2$. It is also invariant by symmetry with respect to horizontal vertical planes. The edges are 140 µm thick and 1100 µm long. At the macroscopic scale, when the honeycomb is considered as a continuum, the material is orthotropic for elasticity. For this unit cell geometry, experiments until rupture have been carried out. Two orientations ($0^\circ$ and $90^\circ$) of the material structure with respect to the specimen axis are tested. The crack paths are shown in Figure 6 and they clearly show how the material orientation influences failure in such material. A closer view of the broken cells is presented in Figure 7. From these experimental observations, it appears that failure occurs at the corners of the unit cells. This observation holds for all the unit cell geometries tested experimentally. Consequently, the cleavage planes for this material are defined by the red lines in Figure 5. In the model proposed herein, two possible orientations of cleavage planes have to be considered.

### 5.1.2. *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^\circ$ case
the elastic tensor is:

$$ C^0 = \begin{bmatrix} 420 & 40 & 0 \\ 40 & 180 & 0 \\ 0 & 0 & 30 \end{bmatrix} \text{MPa} \qquad (43) $$

and the initial bulk modulus $k_0$ is 247 MPa. The other model parameters are $g_c$ which is set to 20 J/m², $\ell = 0.01$ mm and $\beta = 20$ to induce strong anisotropy. For simulating the experiments, a rectangular domain of 60 × 120 mm is considered. An initial crack of length $a = 60$ mm is defined by simply removing the elements of the rectangular finite element mesh cut by this crack. A regular mesh of triangular elements of edge 0.01 mm is used in the simulations, it consists in about 30 000 elements. The vertical displacement of the nodes along the left edge of the rectangular domain is prescribed to induce a mode I opening solicitation at the crack tip.

### 5.1.3. 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.

### 5.2. 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). A square domain of 2 × 2 mm² is analyzed. The displacement is fully constrained (zero horizontal displacement
and remote vertical displacement) in the upper and lower bands of the domain to force
the crack to propagate within the rest of the domain, a band of 0.8 mm remaining free
to deform. The elastic behavior of the material is isotropic. The young modulus is set
to 10 GPa and the Poisson's ratio to 0.2. The fracture energy $g_c$ is 250 J/m$^2$ and the
regularization parameter $\ell$ is 32 $\mu$m. The preferential cleavage planes are oriented at
$\pm 45^\circ$ with respect to the horizontal axis. The mesh consists in about 36 000 triangular
elements of 16 $\mu$m size.

In Figure 9, the evolution of the equivalent phase field variable $d_{eq}$ is plotted for
$\beta = 20$, $d_{eq}$ being defined as $d_1 + d_2 - d_1d_2$ i.e. such that $1 - d_{eq} = (1 - d_1)(1 - d_2)$. 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 it 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 simu-
lation. One clearly observes how each of the cleavage plane is successively activated
and unactivated to finally obtain the zigzag pattern. Subsequently, changing orienta-
tion 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 $\beta = 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 go-
ing 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.

**6. 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 direc-
tionality 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 pro-
posed by Li et al. (2015), a standard finite element implementation is used herein.
Further, and as illustrated by the examples, the proposed model induces additional en-
ergy dissipation when the crack switches from one cleavage plane to an other one. The
definition of the implicitly considered fracture energy distribution as a function of the
kinking angle can be easily extended to 3D as only the normal to the cleavage plane
as to be defined. Also, with the use of the history function, the initial non-linear phase
filed problem recasts in two linear problems which further increases the applicability
of the proposed model to model 3D complex microstructures.
Among the prospects of this work, one can mention: further validation by compar-
ison to model experiments as performed by Takei et al. (2013) and 3D implementation
and application to polycrystalline structures. For alloys, grain boundaries have a con-
stitutive behavior that differs from the bulk. It would thus be useful to couple the
proposed anisotropic phase field model to a cohesive model for modeling failure at the
grain boundaries. For this purpose, Nguyen et al. (2015a) recently proposed a frame-
work to integrate cohesive interfaces within a phase field model. This approach may
be considered to improve the relevance of the model to model the failure behaviour of
polycrystalline systems.

Acknowledgements

The support of the french nationnal research agency (ANR) is gratefully acknowl-
edged through the grant ANR-12-RMNP-0020-02 ECCOFIC.

References

Amor, H., Marigo, J.-J., Maurini, C., 2009. Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments. Journal of the Mechanics and Physics of Solids 57 (8), 1209–1229.

Borden, M., 2012. Isogeometric analysis of phase-field models for dynamic brittle and ductile fracture.

Braides, D., 1998. Approximation of Free Discontinuity Problems. Springer Verlag: Berlin.

Clayton, J., Knap, J., 2015. Phase field modeling of directional fracture in anisotropic polycrystals. Computational Materials Science 98, 158–169.

Comi, C., 2001. A non-local model with tension and compression damage mechanisms. European Journal of Mechanics-A/Solids 20 (1), 1–22.
Francfort, G., Marigo, J., 1998. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids 46(8), 1319–1342.

Freddi, F., Royer-Carfagni, G., 2010. Regularized variational theories of fracture: a unified approach. Journal of the Mechanics and Physics of Solids 58 (8), 1154–1174.

Frémond, M., Nedjar, B., 1996. Damage, gradient of damage and principle of virtual power. International Journal of Solids and Structures 33 (8), 1083–1103.

Hakim, V., Karma, A., 2009. Laws of crack motion and phase-field models of fracture. Journal of the Mechanics and Physics of Solids 57 (2), 342–368.

Karma, A., Kessler, D., Levine, H., 2001. Phase-field model of mode iii dynamic fracture. Physical Review Letters 87 (4), 45501.

Li, B., Peco, C., Millán, D., Arias, I., Arroyo, M., 2015. Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy. International Journal for Numerical Methods in Engineering 102 (3-4), 711–727.

Li, J., Proudhon, H., Roos, A., Chiaruttini, V., Forest, S., 2014. Crystal plasticity finite element simulation of crack growth in single crystals. Computational Materials Science 94, 191–197.

Miehe, C., Hofacker, M., Welschinger, F., 2010. A phasefield model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering 199, 2765–2778.

Miehe, C., Schänzel, L.-M., Ulmer, H., 2015. Phase field modeling of fracture in multi-physics problems. part i. balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids. Computer Methods in Applied Mechanics and Engineering 294, 449–485.
Nguyen, T., Yvonnet, J., Zhu, Q.-Z., Bornert, M., Chateau, C., 2015a. A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography. Computer Methods in Applied Mechanics and Engineering.

Nguyen, T., Yvonnet, J., Zhu, Q.-Z., Bornert, M., Chateau, C., 2015b. A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure. Engineering Fracture Mechanics 139, 18–39.

Pham, K., Amor, H., Marigo, J.-J., Maurini, C., 2011. Gradient damage models and their use to approximate brittle fracture. International Journal of Damage Mechanics 20 (4), 618–652.

Ramtani, S., Berthaud, Y., Mazars, J., 1992. Orthotropic behavior of concrete with directional aspects: modelling and experiments. Nuclear Engineering and Design 133 (1), 97–111.

Takei, A., Roman, B., Bico, J., Hamm, E., Melo, F., 2013. Forbidden directions for the fracture of thin anisotropic sheets: an analogy with the wulff plot. Physical review letters 110 (14), 144301.

Wu, C., McKinney, K., Rice, R., 1995. Zig-zag crack propagation in mgal 2 o 4 crystals. Journal of materials science letters 14 (7), 474–477.
Figures

<figure><img src="image_2.png" /><figcaption>Figure 1: Regularized representation of a crack: two-dimensional case: (a) sharp crack model; (b) regularized representation through phase field.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Benchmark problem for the analysis of surface energy distribution. The phase field variable $d$ is prescribed a unit value within the disc of radius 0.025B at the centre of the domain.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Polar plots of the surface energy $\Gamma_{\ell}$ (a) and the reciprocal surface energy $1/\Gamma_{\ell}$ (b) for for $\mathbf{M} = [1 \ 0]$ and $\beta$ varying from 0 to 20.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Specimen geometry for the failure experiments</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: *D*₂ unit cell geometry for 0° orientation and its admissible cleavage planes (red lines).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Experimental crack paths for two orientation of the material architecture: (a) 0°, (b) 90°</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Close view of the failure mechanisms in the analyzed architecture materials.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Phase field variable in the activated (a,b) and unactivated (c,d) cleavage planes for the two orientations of the model material.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Evolution of the equivalent phase field $d_{eq}$ during the propagation of the guided crack.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Distribution of the two phase field variables at the end of the propagation of the guided crack for strong anisotropy ($\beta = 20$).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Load displacement curve for the example of the guided crack. The red dots corresponds to the states depicted in Figure 9.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Distribution of the two phase field variables at the end of the propagation of the guided crack for weak anisotropy ($\beta = 1$).</figcaption></figure>