Thanh-Tung Tung Nguyen

Julien Réthoré
email: julien.rethore@ec-nantes.fr

Julien Yvonnet

Marie-Christine Baietto

M-C Baietto Lamcos

Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials

Keywords: Crack propagation, Phase field, Polycrystals, Complex crack networks, Cohesive zone model

come

Introduction

Prediction of the strength and failure of materials by simulating damage and microcracking at the microscopic level is a major task to study more accurately the durability of structures. Many engineering materials such as ceramics, alloys, rocks or metals are heterogeneous, containing an assembly of many irregularly-shaped phases or crystallites arranged in a random fashion. Microstructural features in polycrystalline materials such as topology, orientation or grain boundary characteristics, strongly influence the failure mechanisms as well as the macroscopic mechanical response in these class of materials.

Intergranular and transgranular are the two main types of cracks observed in the failure of polycrystalline materials. The intergranular fracture is a fracture that follows the grains of the material, where cracks take place along the grain boundary. We refer to transgranular fracture, when the crack occurs through the grains. In the transgranular fracture the fracture often changes direction from grain to grain following the different lattice orientation, i.e the crack propagates following some preferential cleavage planes. In other words, the fracture in polycrystalline material is strongly anisotropic and dependent on the microstructure. The competition between intergranular and transgranular cracking depends on the loading conditions, bulk and grain boundary behavior,...(see: [START_REF] Akao | Mechanical properties of sintered hydroxyapatite for prosthetic applications[END_REF]). We can observe either only intergranular [START_REF] Was | The mechanism of intergranular cracking of Ni-Cr-Fe alloys in sodium tetrathionate[END_REF][START_REF] Farkas | Intergranular fracture in nanocrystalline metals[END_REF] or only transgranular crack propagation [START_REF] Munroe | Observation of < 001 > dislocations and a mechanism for transgranular fracture on {001} in FeAl[END_REF], or both of them [START_REF] Aguilar-Elguézabal | Fracture behaviour of α-Al2O3 ceramics reinforced with a mixture of single-wall and multi-wall carbon nanotubes[END_REF]. The studies in the literature have shown that the crack path strongly depends on local differences in toughness (grain interior versus grain boundaries), which significantly influences the macroscopic response of the structure. Hence, understanding microstructural effects in polycrystals is crucial to predict the damage mechanisms, but requires numerical developments to take into account the complexity of local phenomena.

Many numerical methods have been developed to investigate the fracture of polycrystalline solids. The eXtended Finite Element Method (XFEM) [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Moran | Extended finite element method for three-dimensional crack modeling[END_REF] has been used in the work of Sukumar et al [START_REF] Sukumar | Brittle fracture in polycrystalline microstructures with the extended finite element method[END_REF] to study the transition from intergranular to transgranular fracture by changing the proportion of toughness between grain boundary and grain interior. The results qualitatively agree with the observation in [START_REF] Yang | Microstructural effects in the fracture of brittle materials[END_REF]. However, simulating crack nucleation in XFEM is a difficult task. Moreover, this technique requires the construction of level-set functions to describe the crack, which can be cumbersome for multiple cracks interaction and/or 3D implementation.

Cohesive zone models (CZM) is another technique widely used in this field. The main idea of CZM is based on the concept introduced in the work of Dugdale [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and Barenblatt et al [START_REF] Barenblatt | The formation of equilibrium cracks during brittle fracture. General ideas and hypotheses. axially-symmetric cracks[END_REF] to address the stress singularity at the crack tip. In these models, all nonlinearities take place in a cohesive zone or fracture process zone of the material. Therein, the damage in this zone is described by cohesive laws (traction displacement-jump relation). The application of CZM to polycrystalline brittle materials can be found in [START_REF] Espinosa | A grain level model for the study of failure initiation and evolution in polycrystalline brittle materials. part II: numerical examples[END_REF][START_REF] Espinosa | A grain level model for the study of failure initiation and evolution in polycrystalline brittle materials. part I: Theory and numerical implementation[END_REF], where the initiation and evolution of failure under quasi-static and dynamic loading was investigated. In the study of Wei et al [START_REF] Wei | Grain-boundary sliding and separation in polycrystalline metals: application to nanocrystalline fcc metals[END_REF], CZM was used to simulate the grain boundary sliding and separation in nanocrystalline metals. The cohesive laws have also been used to model intergranular and transgranular crack propagation in [START_REF] Musienko | Simulation of inter-and transgranular crack propagation in polycrystalline aggregates due to stress corrosion cracking[END_REF] by inserting finite elements in the grain boundary areas for Zircaloy tubes; or in [START_REF] Qian | Application of multiscale cohesive zone model to simulate fracture in polycrystalline solids[END_REF] by using multiscale cohesive zone models. Other studies using CZM for polycrystals failure can be found e.g in [START_REF] Wei | Grain-boundary sliding and separation in polycrystalline metals: application to nanocrystalline fcc metals[END_REF][START_REF] Warner | Micromechanical finite element modeling of compressive fracture in confined alumina ceramic[END_REF][START_REF] Luther | Polycrystal models for the analysis of intergranular crack growth in metallic materials[END_REF][START_REF] Verhoosel | Modelling inter-and transgranular fracture in piezoelectric polycrystals[END_REF]. However, one shortcoming of CZM is that cohesive surfaces can only lie along element edges, which makes crack paths meshdependent.

There are several alternative numerical methods to study failure of polycrystalline material, such as the generalized finite element method (GFEM) [START_REF] Simone | A generalized finite element method for polycrystals with discontinuous grain boundaries[END_REF]. However, in this method the polycrystalline aggregate topology is used to define discontinuous functions. Thus, GFEM only allows the crack to stop at the element edges. Finally, the Boundary Element Method (BEM) has been used to investigate the intergranular microfracture in [START_REF] Sfantos | A boundary cohesive grain element formulation for modelling intergranular microfracture in polycrystalline brittle materials[END_REF].

More recently, the Phase Field Method (or Phase Field Model) has been proposed, based on the pioneering works of Marigo and Francfort [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]. In this method, the brittle fracture is described by using a diffuse approximation of discontinuities through an additional phase field variable [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF][START_REF] Ambrosio | Approximation of functionals depending on jumps by elliptic functionals via Γ -convergence[END_REF][START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF]. It has been shown that the solution of the associated variational problem converges to the solution of the sharp crack description implying discontinuities, in the Γ -convergence sense [START_REF] Maso | An Introduction to Γ -Convergence[END_REF][START_REF] Braides | Approximation of Free Discontinuity Problems[END_REF][START_REF] Braides | Γ -Convergence for Beginners[END_REF]. Such a method is able to simulate brittle crack initiation and propagation without any prescription of the crack geometry. Moreover, the phase field method has a low dependence to the mesh in a classical FEM framework and can handle very complex, multiple crack fronts, branching in both 2D and 3D without adhoc numerical treatment [START_REF] Miehe | A phasefield model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Bourdin | A timediscrete model for dynamic fracture based on crack regularization[END_REF][START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF]. It has been shown that the method accurately predict experimental results in quasi-britle materials, even for complex crack paths [START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF][START_REF] Nguyen | Initiation and propagation of complex 3D networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF][START_REF] Nguyen | A phase field method for modeling stress corrosion crack propagation in a nickel base alloy[END_REF]. The phase field method has been used to study failure in polycrystalline materials in [START_REF] Li | Phasefield modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy[END_REF][START_REF] Clayton | Phase field modeling of directional fracture in anisotropic polycrystals[END_REF][START_REF] Clayton | Phase field modeling and simulation of coupled fracture and twinning in single crystals and polycrystals[END_REF]. However, in the mentioned studies, the model only allows to investigate the fracture along one preferential plane and a new phase for grain boundary was defined. In [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF], we have demonstrated the potential of the method to simulate interactions between bulk brittle fracture and interfacial damage (described by cohesive laws).

This work is concerned with the development of a new phase field formulation for polycrystalline brittle materials 1 , taking into account both anisotropic brittle fracture within grains and interfacial damage at grain boundaries. First, the previous model [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF] is adapted with the new unilateral contact formulation proposed in [START_REF] Clayton | Phase field modeling of directional fracture in anisotropic polycrystals[END_REF][START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF]. This formulation is suitable for introducing elastic anisotropy. Second, anisotropic surface energy with many damage systems (see [START_REF] Nguyen | Phase field modelling of anisotropic crack propagation[END_REF]) is used. Third, CZM in the regularized frame-work is introduced to describe the decohesion of grain boundaries. In this new formulation, the phase field describes the bulk crack surface density (depending on anisotropic effect of material), as well as the interface crack density at grain boundaries. The technique allows to simulate the competition between intergranular and transgranular crack propagation without definition of a new phase for grain boundary.

The overview of the paper is as follows. In section 2 the diffuse approximation of discontinuities for anisotropic systems related to cracks and grain boundary decohesion using phase field is introduced. In section 3, a new formulation of phase field framework able to model both intergranular and transgranular failure is proposed. The finite element implementation of the model is then detailed in section 4. Finally, the potential of the proposed model is illustrated by several numerical examples in section 5.

Regularized representation of free discontinuities

Smeared approximation of cracks and grain boundary in polycrystalline material

Let Ω ⊂ R D an open domain describing a cracked solid, with D the space dimension and ∂Ω its boundary. Let Γ be 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. [START_REF] Miehe | A phasefield model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]) that such a function can be determined by solving the following boundary value problem on Ω:

   d -2 ∆d = 0 in Ω, d(x) = 1 on Γ, ∇d(x) • n = 0 on ∂Ω, (1)
where ∆(.) is the Laplacian operator, is a regularization parameter describing the actual thickness 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:

d(x) = Arg inf d∈S d Γ (d) , (2)
with S d = {d|d(x) = 1 on Γ ∀x ∈ Γ } and where

Γ (d) = Ω γ(d, ∇d) dΩ (3)
represents the total crack length. In (3), γ(d, ∇d) denotes the crack density function per unit volume, defined herein by:

γ(d, ∇d) = 1 2 d 2 + 2 ∇d • ∇d. (4)
The grain boundaries are here described in the same manner, introducing a field β(x). The field β(x) satisfies:

   β(x) -2 β (x) β(x) = 0 in Ω, β(x) = 1 on Γ B , ∇β(x) • n = 0 on ∂Ω, ( 5
)
where β is the regularization parameter describing the width of the regularized grain boundaries. Similarly, [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF] is the Euler-Lagrange equation associated with the variational problem:

β(x) = Arg inf β∈S β Γ β (β, ∇β) , (6)
where

S β = {β|β(x) = 1 ∀x ∈ Γ B } and Γ β (β) = Ω γ β (β, ∇β) dΩ, (7)
where Γ β represents the total grain boundaries length and γ β is defined by

γ β (β, ∇β) = 1 2 β β(x) 2 + β 2 ∇β(x) • ∇β(x). (8)
For β → 0 the above variational principle leads to the exact description of the sharp grain boundaries Γ B .

Smeared approximation of cracks for anisotropic material

In the anisotropic case, the formulation ( 4) is extended to anisotropic materials. An anisotropic crack surface density function is written by the following expression:

γ(d, ∇d, ω) = 1 2 d 2 + 2 ω : (∇d ⊗ ∇d), (9)
where ω is a second-order structural tensor, being invariant with respect to rotations (characterizing the material anisotropy). As considered in [START_REF] Clayton | Phase field modeling of directional fracture in anisotropic polycrystals[END_REF], to make the energy release rate orientation-dependent, the tensor ω can be defined by:

ω = 1 + α(1 -M ⊗ M), (10)
where M denotes the unit vector normal to the preferential cleavage plane (with respect to the material coordinates), and α 1 is used to prevent damage to develop on planes not normal to M. In the case of isotropic material α = 0.

In the case of cleavage planes with multiple discrete orientations, we introduce the multiple phase field d to quantify the damage accumulation on each such plane. Thus the total crack length is here rewritten as:

Γ (d, ω i ) = i Γ i (d i , ω i ) = i Ω γ i (d i , ∇d i , ω i ) dΩ (11)
The variational derivative of the crack density function for each damage system is now defined as:

δγ i (d i , ∇d i , ω i ) = d i -∇ 2 d i + α i ∇ 2 d i -M i ⊗ M i : ∇∇d i (12)
where ω i belongs to the set ω 1 , ω 2 , ...ω n is calculated from Equation [START_REF] Braides | Γ -Convergence for Beginners[END_REF] cleavage plane (depending on the unit vector normal M i and coefficient α i ). Note that the regularization parameter can be different for each plane. However, in this work we assume that is the same for all planes.

In order to estimate the anisotropy introduced in the surface energy, we compute the distribution of surface energy g c Γ (d i , ω i ) for a phase field d i (x), which is determined numerically by the following variational problem:

d i (x) = Arg inf di∈S d Ω γ(d i , ∇d i , ω i ) dΩ , (13)
A Benchmark problem as described in Fig. 2 is considered. The phase field variable d i is assigned a unit value within a disc of radius r = 0.05B at the center of a square domain of size B. Once the computation of multi-phase-field is performed, the corresponding surface energy distribution is calculated at the integration points. For the case of material containing two preferential directions M = [1 0], and M = [0 1], g c = 1 J/m 2 and α ranges between 0 and 10, polar plots of element average value of surface energy and reciprocal surface energy at a distance r = 0.05B + 2 from the domain center are depicted in Fig. 3. For α = 0, the surface energy does not depend on the orientation and isotropy is recovered. When α is increased from 0 to 10, there is a strong modification of the surface energy distribution as illustrated in Fig. 3(a). This provides the anisotropic behavior of fracture evolution. The sensitivity of the reciprocal surface energy with respect to α is analyzed in Fig. 3(b). It is clearly shown that for α ≤ 1 the reciprocal surface energy distribution is convex but it becomes non-convex for α > 1. A smooth transition from weak to strong anisotropy can thus be modeled within the proposed framework (see [START_REF] Li | Phasefield modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy[END_REF] for more details on this aspect).

Smeared displacement jump approximation (JA1)

The displacement jump [[u(x)]] created by grain boundary decohesion is approximated as a smooth transition. It can be defined by using a Taylor expansion at first order of the assumed smoothed regularized displacement field as in our previous work [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF], see Appendix A. The advantage of this approximation is that it does not require additional variables to describe the jump at the interface, and thus can drastically reduce computational times compared to what is propose in [START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF].

Another approximation, denoted here as (JA1), based on the addition of a new field to described displacement jump following the work [START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF], is also tested in this work. The additional field approach is computationally more expensive, but seems to improve the convergence when applied to polycrystalline materials with multiple junctions between grains. The details of the comparison between the two approximations is discussed in Section 5.1.2.

Let Γ B be a grain boundary. In the context of the approximation (JA1), the discrete displacement jump at Γ B is approximated by the auxiliary jump field v(x):

[[u]] Γ B = Ω v(x)γ β (β, ∇β) dΩ = v(x)Γ β . ( 14
)
It has been shown in [START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF] that in this context, the auxiliary field is required to be constant in the normal direction of the grain boundary:

∂v ∂n Γ B = 0. ( 15
)
The description of this approximation is provided in Fig. 4, wherein the displacement jump across the grain boundary is described. A one-dimensional illustration of smeared grain boundary is plotted by the dashed line. The displacement jump v(x)Γ β only exists in the smeared grain boundary where Γ β > 0. However, the auxiliary field is here not only defined on the interface but over the entire domain. In order to reduce the computation costs, we can define the distributed internal discontinuity region only in a subdomain such that:

Γ β = x ∈ Ω | β(x) > 10 -3 . ( 16
)
The normal vector n Γ B to Γ B can be also approximated as:

n Γ β (x) = n Γ B (x) = ∇β(x) ∇β(x) . (17)
The second approximation scheme (JA2) is described in the Appendix A.

3 Phase field model

Energy functional

Let us consider a polycrystalline solid as described in section 2.1. The solid contains both cracks and grain boundaries, implying strong displacement discontinuities. In a standard framework of sharp discontinuity description, the total energy is given by: 18 In [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] g c is the fracture energy, Ψ B is a grain boundary energy function depending on the displacement jump across the grain boundary Γ B and κ is a history parameter to distinguish between loading and unloading. Once, the cohesive fracture is Ψ B , the cohesive tractions are derived through its differentiation (see more details in Appendix C):

E = Ω W u (ε(u)) dΩ+ Γ g c dΓ + Γ B Ψ B ([[u]] , κ)dΓ.(
t([[u]] , κ) = ∂Ψ B ([[u]] , κ) ∂ [[u]] (19)
If a regularized description for strong discontinuities related to both cracks and grain boundaries is adopted (substituting [[u]] by v(x)), then the infinitesimal strain tensor can be decomposed into a part related to the bulk and a part induced by the smoothed jump at the grain boundaries, denoted by ε e and ε, respectively:

ε = ε e + ε. ( 20
)
This decomposition is proposed here such that ε → 0 away from the grain boundary, i.e. when β(x) → 0. The fracture energy is also regularized by the crack density function γ i (d i , ∇d i , ω i ) for each damage variable d i , the energy of grain boundary is regularized by γ β (β, ∇β) and the strain energy is replaced by energy of a damageable material W e u (ε e (u), d). The above functional ( 18) is rewritten as:

E = Ω W e u (ε e (u, β), d) dΩ +g c i Ω γ i (d i , ∇d i , ω i ) dΩ+ Ω Ψ B (v, κ)γ β (β) dΩ,( 21
)
where γ i and γ β have been defined in section (2.1). This definition implies that for β(x) → 0 (away from the interface) then γ β → 0 and ε e → ε. As a consequence of these definitions, we recover the regularized energy functional for brittle fracture without grain boundaries [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Miehe | A phasefield model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Nguyen | Phase field modelling of anisotropic crack propagation[END_REF]:

E = Ω W u (ε(u), d) dΩ+g c i Ω γ i (d i , ∇d i , ω i ) dΩ.( 22
)
In ( 21), we identify

W = W e u (ε e (u, β), d) + g c i γ i (d i , ∇d i , ω i ) + Ψ B (v, κ)γ β (β) (23)
as the total energy density.

Phase field problem

In this subsection we formulate a crack phase field evolution law that can guarantee the irreversibility of the process. Assuming isothermal process and without external mircoforces, a reduced form of the Clausius-Duhem inequality can be written as:

i - δW δd i ḋi ≥ 0, (24)
where δW δd i is the variational derivative of W with respect to the phase field d i . The total dissipation due to damage is defined by:

φ = i φ i = i - Ω δW δd i ḋi ≥ 0, (25)
where each dissipation φ i corresponds to one damage variable d i . Assuming that the damage systems are independent, then the condition (25) becomes:

       φ 1 ≥ 0, φ 2 ≥ 0, ... φ n ≥ 0. ( 26
)
The damage parameter d i is an increasing function, i.e ḋi ≥ 0. Thus for each damage variable, we have:

F i = - δW δd i = - ∂W e u ∂d i -g c δγ i (d i , ∇d i , ω i ) ≤ 0, (27)
where δγ i (d i , ∇d i , ω i ) being defined in Equation ( 12) and F i is the thermodynamic force associated with d i .

Invoking the principle of maximum dissipation, the dissipation φ i must be maximum under the constraint [START_REF] Mehl | Pressure dependence of the elastic moduli in aluminum-rich Al-Li compounds[END_REF]. It yields that for ḋi > 0, F i = 0 for each damage system, which gives (see more details in [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF], [START_REF] Nguyen | Phase field modelling of anisotropic crack propagation[END_REF]):

F i = - ∂W e u (u, d i ) ∂d i -g c δγ i (d i , ∇d i , ω i ) = 0. ( 28
)
In order to prevent the issue of cracks interpenetration in compression mode, a unilateral contact formulations is used. There are two common kinds of formulations often chosen for phase field method. (a) it can be based on the decomposition of the strain tensor into positive and negative strain components (see [START_REF] Miehe | A phasefield model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF]), (b) or based on the decomposition of strain tensor into spherical and deviatoric strain components (see [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF]). The first one shows a good agreement with experimental observation in brittle material [START_REF] Nguyen | Modeling of complex microcracking in cement based materials by combining numerical simulations based on a phase-field method and experimental 3D imaging[END_REF]. However it is difficult to apply it for anisotropic materials, because of the lack of general formulation for the elastic tensor. In this work, the second approach is adopted. The elastic strain is decomposed into deviatoric ε e dev and spherical ε e sph parts. Then, it is assumed that the damage is created by expansion only (positive spherical part) and shear:

W e u =                1 2 g(d i ) ε e : C 0 : ε e if tr ε e ≥ 0, 1 2 ε e sph : C 0 : ε e sph + 1 2 g(d i ) ε e dev : C 0 : ε e dev if tr ε e < 0, (29)
where tr ε e is the trace operator for a second-order tensor ε e and C 0 denotes the initial elastic tensor of the material, possibly anisotropic. The degradation function g(d i ) in ( 29) is assumed to have the simple form:

g(d i ) = i (1 -d i ) 2 + k. ( 30
)
The function g(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. [START_REF] Braides | Approximation of Free Discontinuity Problems[END_REF]) and g(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 chosen here (1-d i ) 2 is the simplest case to ensure the existence of a regular solution in the sense of Carfagni. Alternatively quartic function or cubic functions have been introduced in [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF], [START_REF] Borden | Isogeometric analysis of phase-field models for dynamic brittle and ductile fracture[END_REF]. The small parameter k << 1 is used to maintain the well-posedness of the problem for partially broken parts of the domain.

By introducing the bulk modulus k 0 for the undamaged material (relating the spherical part of the strain to the spherical part of the stress), the elastic tensor is now written as: [START_REF] Munroe | Observation of < 001 > dislocations and a mechanism for transgranular fracture on {001} in FeAl[END_REF] where the sign function sign -(x) = 1 if x < 0 and sign -(x) = 0 if x ≥ 0. The strain energy is now rewritten as:

C(d) = g(d)C 0 + k 0 1 ⊗ 1 1 -g(d) sign -(tr ε e ),
W u (ε e (u), d) = 1 2 ε e : C(d) : ε e . (32)
From ( 28), the thermodynamic force F i associated to ḋi satisfies:

F i = (1 -d i ) j =i (1 -d j ) 2 ε e : C h : ε e -g c δγ i (d i , ∇d i , ω i ) = 0, (33)
where

C h = C 0 (x) -k 0 (x) 1 ⊗ 1 sign -(tr ε e (x, τ )). As (1 -d i ) j =i (1 -d j )
2 ε e : C h : ε e ≥ 0, and g c > 0 then

δγ i (d i , ∇d i , ω i ) ≥ 0. ( 34
)
Expressing the variation of the crack length:

d dt Γ i = Ω δγ i (d i , ∇d i , ω i ) ḋi dΩ, (35)
we can check that due to (34)

Γ i ≥ 0, (36)
verifying the non-reversible evolution of cracks. In addition, to handle loading and unloading histories, [START_REF] Miehe | A phasefield model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] introduced the strain history functional, which is expressed in the present context by: [START_REF] Moran | Extended finite element method for three-dimensional crack modeling[END_REF] and using (1), it yields the following phase field problem to be solved to evaluate the field d i (x, t) at time t:

H i (x, t) = max τ ∈[0,t] j =i (1 -d j ) 2 ε e : C h : ε e . ( 37
)
H i (x, t) is substituted to j =i (1 -d j ) 2 ε e : C h : ε e in
   (1 -d i )H i -g c δγ i (d i , ∇d i , ω i ) = 0 in Ω, d i (x) = 1 on Γ, ∇d i (x) • n = 0 on ∂Ω. ( 38
)
The associated weak form for each damage system is obtained as:

Ω H i + g c d i δd + g c ∇d i 1 + α 1 -M i ⊗ M i (δd i ) dΩ = Ω H i δd dΩ. ( 39
)
Remark: the choice of has been discussed e.g in [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF]. In the mentioned works, a relationship has been established between and material parameters. For example, assuming uniaxial traction of a bar, the relationship:

σ + c = 3 16 3 2 Eg c , (40)
where σ + c is the critical value of stress in traction. As mentioned in [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF], the shape of the damaged zone (and hence the approximate crack path) can significantly change with the length if is not sufficiently small with respect to the size L of the domain Ω. The authors then define two cases: (a) is considered as a pure numerical parameter of the regularized model of brittle fracture; (b) is seen as a material parameter for a gradient damage model. In the first case, it is recommended to take as small as possible to better approximate brittle fracture, with regards to the size of the mesh. In the second case, should be identified by experimental data through relationships like Eq. [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]. Such analysis has been conducted with experimental validations in [START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF][START_REF] Nguyen | Initiation and propagation of complex 3D networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF].

Displacement problem

Considering the quasi-static equilibrium equations without body forces, the corresponding equations describing the problem are written as:

       ∇ • σ = 0 ∀x ∈ Ω, n Γ β • σ = t([[u]] , κ) ∀x ∈ Γ β , u = u ∀x ∈ ∂Ω u , n t • σ = t ∀x ∈ ∂Ω t , (41)
where σ is the Cauchy stress in the bulk material; t([[u]] , κ) are the tractions at grain boundaries introduced in [START_REF] Herbig | 3-d growth of a short fatigue crack within a polycrystalline microstructure studied using combined diffraction and phase-contrast x-ray tomography[END_REF], with outward normal n Γ B ; t, u are the prescribed tractions and the prescribed displacements on ∂Ω t , ∂Ω u , respectively. In the regularized framework, this boundary value problem is rewritten as:

                   ∇ • σ(u, v, d) = 0 ∀x ∈ Ω, γ β t(v, κ) -n Γ β • σ(u, v, d) = ζ ∂ 2 v ∂(x n ) 2 ∀x ∈ Γ β , ∂v ∂x n = 0 ∀x ∈ ∂Γ β , u = u ∀x ∈ ∂Ω u , n t • σ = t ∀x ∈ ∂Ω t , (42)
where

x n = (x-x Γ B )•n Γ B and x Γ B = argmin y∈Γ ( y -x ),
the positive constant ζ is introduced to ensure that the auxiliary displacement jump field is constant in the normal direction. The Cauchy stress is here defined by:

σ(u, v, d) = ∂W e u ∂ε e = C(d) : ε e (43)
where the expression of the elastic strain ε e is given by (see [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]):

ε e = ∇ s u -n Γ β ⊗ s vγ β , (44)
where

(∇ s u) ij = 1 2 (u i,j + u j,i ) and n Γ β ⊗ s v ij = 1 2 (n i v j + v i n j ).
From [START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF] we identify ε as:

ε = n Γ β ⊗ s vγ β . ( 45
)
Using the variation in the displacement for (42) 1 , the variation in the displacement jump for (42) 2 , we obtain the corresponding weak form:

                           Ω σ : ε e (δu) dΩ - ∂Ωt t • δu dS = 0 Ω γ β t(v, κ)δv -σ : ε (δv) -ζ ∂v ∂x n ∂δv ∂x n dΩ - ∂Γ β ∂v ∂x n δv dS β = 0 (46)
4 Discretization and numerical implementation

FEM discretization of the displacement problem

We introduce the following classical FEM discretization:

u = N u u e , δu = N u δu e , ∆u = N u ∆u e , v = N v v e , δv = N v δv e , ∆v = N v ∆v e , (47)
where u e , δu e and ∆u e denote nodal displacement components, nodal trial function components and nodal incremental displacement components, respectively. v e , δv e and ∆v e denote nodal displacement jump components in one element, nodal trial function jump components and nodal incremental displacement jump components, respectively. Furthermore, using usual notations, strains are obtained from:

[ε e ] (u) = B u u e , [ε e ] (δu) = B u δu e , [ε e ] (∆u) = B u ∆u e , [ε] (v) = B v v e , [ε] (δv) = B v δv e , [ε] (∆v) = B v ∆v e , (48)
with ∂v x n = G v v e . (49)
A standard Newton-Raphson method can be used to solve the non-linear problem [START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF]. Introducing the above discretization yields the following linear system is solved at each iteration:

K tan ∆ũ ∆ṽ = -R, (50)
where:

K tan = K uu K uv K vu K vv , R = F u F v , (51)
with:

K uu = Ω B T u CB u dΩ, K uv = Ω -γ β B T u CB v dΩ, K vu = Ω -γ β B T v CB u dΩ, K vv = Ω γ 2 β B T v CB v + γ β N T v ∂t(v, κ) ∂v N v +ζG T v G v dΩ,
and

F u = Ω B T u CB u u e dΩ - Ω γ β B T u CB v v e dΩ, F v = Ω -γ β B T v CB u u e + γ 2 β B T v CB v v e +γ β N T v t(v, κ)v e + ζG T v G v v e dΩ,
where C is the matrix form corresponding to the fourthorder elastic tensor C in Eq. [START_REF] Munroe | Observation of < 001 > dislocations and a mechanism for transgranular fracture on {001} in FeAl[END_REF]. In these equations, {∆ũ, ∆ṽ} T is a column vector containing the nodal values of {∆u, ∆v}, respectively.

FEM discretization of the phase field problem

In the present work, a staggered procedure is employed, i.e. we solve alternatively the phase field problem and the mechanical problem. Using the history functions, given displacements from the mechanical problem, the phase field problem is linear. The phase field and phase field gradient for each preferential system are approximated in one element by

d i (x) = N d (x)d e i and ∇d i (x) = B d (x)d e i , (52)
where N d (x) and B d (x) are vectors and matrices of shape functions and of shape functions derivatives, respectively, and d e i are the nodal values of d i . The same discretization is employed for the test function:

δd(x) = N d (x)δd e and ∇δd(x) = B d (x)δd e . ( 53
)
Introducing this FEM discretization in [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF] results into a linear system of equations:

K di di = F di , ( 54
)
where di is a column vector containing the nodal values of d i and

K di = Ω g c + H i N T d N d + g c B T d ωB d dΩ ( 55
)
and

F di = Ω N T d H i (u n ) dΩ. ( 56
)

Overall algorithm

The overall algorithm is described in the following.

-Initialization 1.1 Initialize the displacement field u 0 (x), the phase field d 0 i (x) for each preferential plane, and the strain-history functional H 0 i = 0.

END END

The present code is parallelized and has been implemented in Matlab ®. The simulations have been run using 12 cores.

Numerical examples

Model assessments

In this first example, the proposed model will be assessed. We consider a benchmark problem consisting in a square domain which contains two grains and an initial crack, as depicted in Fig. 5. The size of the square domain is

L × B = 1 × 1 mm 2 .
Ψ B [[u]] = g B c 1 -1 + [[u]] v n exp - [[u]] v n . ( 57
)
The normal traction law is obtained as:

t [[u]] = g B c [[u]] v 2 n exp - [[u]] v n , (58)
where v n = g B c /(t u exp(1)). The fracture strength is t u = 0.1 GPa, and can be dependent on the angular mismatch between grains as described in the Appendix C. However, for the sake of simplicity and due to the lack of experimental data, t u is considered here as independent for misorientation between two grains. The fracture energy within the grain boundary and in the bulk material are g B c = 1.45 × 10 -3 kN/mm, g c = 1.5 × 10 -3 kN/mm. The preferential cleavage plane is oriented at -45 • with respect to the horizontal axis within the entire domain. The computation is performed with monotonic displacement increments of ∆U = 10 -4 mm for 150 load increments. The displacements are prescribed along the y-direction for upper edge (y = B) while the displacement along x are here free. Along the lower edge (y = 0), the displacements along y are fixed to zero, while the displacements along x are free. Plane strain condition is assumed.

Influences of length scale

The proposed framework uses two internal length parameters, for smeared cracks and β for smeared grain boundaries. As mentioned in Section 3.1, these length scales are related to the material parameters and strongly influence the predicted results of the phase field solution. In this section, its effects on the fracture behavior will be investigated in more details.

Firstly, to analyze the influences of , the regularized parameter for grain boundary is fixed, β = 0.012 mm, while three values of are taken: = 0.01, 0.012 and 0.015 mm. The comparison of the crack paths is depicted in Fig. 6. Increasing of length scale results in larger cracks, but the crack path is the same. Moreover, the length scale influences the time of the crack onset, a large value of leads to a crack appearing earlier. That will change the global behavior of the structure as depicted on the stress -displacement curve in Fig. 7 Secondly, to analyze the influences of β , the crack internal length is fixed, = 0.012 mm, while three values of β are taken: β = 0.01, 0.012 and 0.015 mm. The comparison of obtained crack paths is depicted in Fig. 8. The crack path is the same for different β . However, the decreasing of regularized parameter for grain boundary increases the cohesive tractions, i.e increasing its stiffness and the material becomes more brittle. Hence, at the same prescribed displacement, reducing β will provide an increasing of the crack length. The details of this effect is depicted on the stress -displacement curve in Fig. 9

Assessment of displacement jump approximation

In this example, the regularization parameter is chosen as β = = 0.012 mm. Two approximations for the displacement jump will be compared: (JA1) is based on the additional auxiliary field, (JA2) is based on the approach described in the Appendix A.

We study the behavior of the structure at a chosen loading for the two approximations (JA1) and (JA2). The displacement field in the y-direction and the corresponding damage states are plotted at loading step U = 0.008 mm in Fig. 10. We also plot the displacement field along a line of investigation (x = 0.5L). The results from both algorithms match very well. In both cases, we observe the development of damage at the crack tip. We also capture the displacement jump across the grain boundary. This phenomenon shows a good agreement with the observation in the literature. By adding the smeared cohesive law for grain boundaries in the phase field method, we can model in this framework both intergranular and intragranular fracture, and interactions between both without defining a new phase.

The application of the second jump approximation (JA2) for polycrystalline materials induces some convergence problem when solving the displacement problem. That may be due to the effects of orientation, or oscillation of displacement jump at the vertex of grains, where the approximation do not work very well. However, we note that this method is still robust for many cases of heterogeneous materials, and that about 3 -4 times of the computational time can be saved as compared to the (JA1) method. Then, defining a method without an additional field (as in the (JA2) method) but maintaining the robustness of nonlinear computations (as in the (JA1) method) is still an open question and requires more investigations.

Orientation effects

In order to investigate the influences of preferential orientation on the fracture behavior, the preferential direction of Grain 1, denoted here as ϕ, is chosen within a range from 0 • to 90 • , while the preferential direction of Grain 2 is fixed (see Fig. 12). The material parameters are set the same as Section 5.1.2.

The obtained crack paths are depicted in Fig. 13. We observe here an interesting phenomenon: transgranular fracture is obtained for ϕ ≤ 60 while when increasing the orientation ϕ > 60, a transition to intergranular fracture is obtained. The crack propagation direction does not follow the preferential direction anymore, which may be due to the unfavorable orientation for fracture growing. This illustrates that even if the grain boundary fracture energy does not depend on the misorientation between grains, its influence on the stress distribution around the grain boundary, which lead to can promote a transition to intergranular fracture. 5.2 Crack propagation in a bi-material containing a grain boundary

Influence of grain boundaries

In this next example, we analyze the propagation of a crack when it passes across a grain boundary. Two cases are considered: (1) the grain boundary is unfavorable for crack propagation due to different angles between the preferential cleavage plane and the grain boundary; (2) the grain boundary is favorable for crack propagation. For this purpose, the structure described in the previous example (Section 5.1) is used for the first case. The details of the second structure are depicted in Fig. 14(b).

The structure is meshed with triangular elements with h max e = 0.05 mm and h min e = 0.003 mm. The material parameters are chosen the same as in the previous test. The computation is performed with monotonic displacement increments ∆U = 10 -4 mm for 150 load increments.

Figs. [START_REF] Espinosa | A grain level model for the study of failure initiation and evolution in polycrystalline brittle materials. part II: numerical examples[END_REF] show the damage evolution (phase field d(x)) for the structure ST1. When the direction of grain boundary and preferential plane are strongly different, only a small influence of grain boundary on crack trajectory is observed. After it crosses the grain boundary, then the crack keeps its preferential direction related to the grain's orientation.

On the contrary, when the grain boundary and preferential plane have similar orientations. The crack propagation has a tendency to follow the grain boundary. In other words, the crack is attracted within the grain boundary, see Figs. [START_REF] Faber | Crack deflection processes-I. theory[END_REF].

Note that the intergranular fracture can be obtained when the grain boundary fracture parameter is chosen smaller than that one of the grains. We perform another simulation on the structure ST1, with a fracture energy chosen as g B c = g c /5. The comparison of crack morphology for the three different situations are depicted in Figs. [START_REF] Farkas | Intergranular fracture in nanocrystalline metals[END_REF].

The load-displacement curve for the three configurations is plotted in Fig. 18. We observe a significant change of the overall behavior when the crack propagates across the interface (depicted by arrows in the Fig. 18). The unfavorable grain boundary orientation leads to a resistance of the structure ST2 that is higher than for the structure ST1, while the intergranular fracture exhibits the smallest resistance.

Crack deflection at grain boundary: tilt/twist

In this second example, we investigate the influence of the orientations of the preferential slip planes in neighbouring grains on the crack propagation. Difference in the orientation between two grains can result in the twist and/or tilt of the crack when crossing a grain boundary [START_REF] Faber | Crack deflection processes-I. theory[END_REF]. To analyze this phenomena, we consider a 2-grain structure in 2D and 3D, depicted in Fig. 19. Tilting is investigating with a 2D configuration whereas twisting requires a 3D computation. The length of the domain for the 2D case is L × B = 2 × 1 mm 2 . For the 3D domain, different thickness will be considered.

The two models contain an initial crack perpendicular to the y-axis of length 0.5 mm in the x direction. The details of the geometries are depicted in Fig. 19. The numerical simulations are run using the parameters as defined in Table . 1. The elastic tensor is considered of cubic symmetry. In the grain natural coordinates system ( 0 • orientation), this elastic tensor is chosen as follows:

Initial crack

C 0 g =        
280 120 120 0 0 0 120 280 120 0 0 0 120 120 280 0 0 0 0 0 0 80 0 0 0 0 0 0 80 0 0 0 0 0 0 80

        (GPa.) (59)
The bulk modulus is obtained as k 0 = (C 11 +2C 12 )/3 = 173.3 GPa. We use the same cohesive law as in the previous example. The fracture strength is t u = 0.15 GPa and is considered independent for the angular mismatch between grains. The fracture energy in the grain boundary and in the grains are g B c = 1.25 × 10 -3 kN/mm, g c = 1.5 × 10 -3 kN/mm. The ratio g B c /g c = 0.83 is chosen here to promote the transgranular crack propagation. Each grain has two preferential cleavage planes, which are described in Table 2 The boundary conditions are as follows: on the lower surface (y = 0), the y-displacements are fixed, while the x-displacements (and z-displacements in 3D case) are free. On the upper surface, the x-displacements (and z-displacements in 3D case) are free, while the y-displacements are prescribed to an increasing uniform value of ∆U during the simulation. Plane strain is assumed for the 2D case.

The microcracking distribution (damage variables d 1 (x), d 2 (x) ) for the 2D structure is depicted in Fig. 20 for several loading steps. To better visualize the final crack state, we also plot an equivalent damage field [START_REF] Nguyen | Phase field modelling of anisotropic crack propagation[END_REF]). The results show the propagation of a crack along one slip plane in each grain. The final crack path shows a good agreement with the experimental observation of tilt character of the grain boundaries in literature [START_REF] Zhai | A crystallographic mechanism for fatigue crack propagation through grain boundaries[END_REF][START_REF] King | Three-dimensional in situ observations of short fatigue crack growth in magnesium[END_REF].

d eq = d 1 + d 2 -d 1 d 2 (see
The results of crack propagation in the 3D structure with dimension L × B × H = 2 × 0.3 × 1 mm 3 are depicted in Fig. 21. The crack plane twisting is clearly captured, reproducing very well the phenomena of failure in polycrystalline materials [START_REF] Faber | Crack deflection processes-I. theory[END_REF][START_REF] Zhai | The grain boundary geometry for optimum resistance to growth of short fatigue cracks in high strength Al-alloys[END_REF]. By changing the dimension of the 3D structure to L×B ×H = 2×1.5×1 mm 3 , in Fig. 22, the crack front finally forms a M (or W) as often observed in experiments. The comparison of the results for these two different thicknesses clearly shows that there is a length scale associated with the formation of the twisting patterns. For further comparisons with experiments it would also be interesting to analyze in details the transient regime between the initial flat crack and the final rough crack surface.

Structure Grain 1 Grain 2 2D structure n 1 = [1 ; 0] n 1 = [1/ √ 2 ; 1/ √ 2] n 2 = [0 ; 1] n 2 = [1/ √ 2 ; -1/ √ 2] 3D structure n 1 = [0 ; 1 ; 0] n 1 = [0 ; 1/ √ 2 ; 1/ √ 2] n 2 = [0 ; -1 ; 0] n 2 = [0 ; 1/ √ 2 ; -1/ √ 2]
Table 2: Crack deflection at grain boundary test: description of preferential cleavage planes The numerical results of this section show that this new model can predict very well the crack deflection (including crack plan twist and tilt) and path across a grain boundary. The fracture energy in the bulk, g c , is set to 0.1 N/m. The regularization parameters are chosen as β = 0.03 µm for smeared grain boundaries and = 0.025 µm for the phase field. Plane strain condition is assumed.

A bi-linear cohesive law (see Fig. 24) has been used in this example:

t = k 1 [[u]] if [[u]] < δ 1 , t = k 1 δ 1 + k 2 ([[u]] -δ 1 ) if [[u]] ≥ δ 1 , with k 2 = - k 1 δ 1 δ 1 -δ 2 .
(61)

In that case, the cohesive fracture energy along a grain boundary is defined by:

g B c = k 1 δ 1 δ 2 2 (62)
and the fracture strength is assumed independent of the angular mismatch between grains, reads : The influence of the fracture energy of grain boundary g B c and the anisotropic factor α is analyzed by using three different sets of parameters: (C 1 ) t u = 1 MPa, δ 1 = 0.01 µm, δ 2 = 0.1 µm, or equivalently g B c = 0.05 N/m, α = 50; (C 2 ) t u = 1 MPa, δ 1 = 0.01 µm, δ 2 = 0.1 µm, or equivalently g B c = 0.05 N/m, α = 5; (C 3 ) t u = 1.5 MPa, δ 1 = 0.01 µm, δ 2 = 0.12 µm, or equivalently g B c = 0.09 N/m, α = 50. In the two first cases (C 1 ) and (C 2 ), the fracture energy in the grain boundary is chosen smaller than the one of the bulk to promote intergranular fracture or at least to mix inter and transgranular fracture. The anisotropic factor in case (C 2 ) is taken α = 5 to allow smoother crack bifurcation in cases of strong misorientation. In the last case (C 3 ), we have chosen g B c /g c = 0.9 to promote transgranular fracture.

t u = k 1 δ 1 . ( 63
) [[u]]
In the crystallographic orientation, the cleavage planes are [1 0] and [0 1]. The damage variable corresponding to these two preferential directions are denoted d 1 (x) and d 2 (x), respectively. Monotonic tensile displacement increments of ∆U = 4.10 -3 µm have been prescribed for d i < 0.8 in all elements and ∆U = 6.10 -4 µm for d i ≥ 0.8 at one integration point.

The computational time for each loading step is about 90 seconds. Each simulation took about 7 hours. The results of crack propagation for the case (C 1 ) are depicted in Fig. 25, where the different phase field variables are plotted. Both inter and transgranular fracture occur in this situation (g B c /g c = 0.5). In addition to the first crack derived from the initial crack, we also observe the initiation of other crack networks at grain boundary (see Fig. 25(b)(c)). Another interesting point in this simulation is that due to the high value of α only one cleavage plane is activated in each grain. Then depending on the misorientation of the neighbooring grains and the orientation of the grain boundary, the crack choses to keep the former cleavage active, to switch to the other cleavage plane or to run within the grain boundary (see Fig. 26).

A comparison of the final crack paths for the three different cases is proposed in Fig. 27. In the second case (C 2 ), when the anisotropy coefficient is chosen smaller α = 5, transgranular failure is less constrained to develop along one single cleavage plane. Then, in some grains the two systems are activated, resulting in a crack orientation perpendicular to the applied stress and in between the orientations of the cleavage planes. In this simulation, after running through about half of the analyzed volume, the crack finds a more favorable intergranular path.

In the last case (C 3 ), the fracture energy in grain boundary is set to an higher value g B c /g c = 0.9. Then, the anisotropic factor being set to its value for the case (C 1 ) with α = 50, transgranular fracture is promoted. These results demonstrate the ability of the proposed model to capture all three different scenarios of fracture in polycrystalline materials, namely intergranular, transgranular and a combined inter transgranular fracture.

The load-displacement curves for the three cases are depicted in Fig. 28. The elastic behavior for two first case are the same, while the resistance of material in intergranular fracture case (C 1 ) is quite smaller than case (C 2 ). The case (C 3 ) exhibits a different behavior, because the cohesive energy is chosen bigger in this case

(k C3 1 = 150 > k C1,C2 1
= 100), i.e that also changes the behavior of material even in elastic state. The elastic tensor is of cubic symmetry. In the crystalographic orientation, it is given as:

C 0 g =         280
1/ √ 3 1/ √ 3 1/ √ 3], [-1/ √ 3 -1/ √ 3 1/ √ 3], [-1/ √ 3 1/ √ 3 1/ √ 3] and [1/ √ 3 -1/ √ 3 1/ √ 3
] in crystallographic orientations are chosen, corresponding to four damage variables d 1 (x), d 2 (x), d 3 (x) and d 4 (x). The equivalent damage variable d eq (x) is here defined by:

d eq = i d i - i =j d i d j + i =j =k d i d j d k - i d i , (65)
where i, j, k = 1, 2, 3, 4. The boundary conditions are as follows: on the lower surface (z = 0), the z-displacements The computational time for each loading step is about 3700 seconds, the entire simulation took about one week. Fig. 31 shows the evolution of the crack for different loading steps, where the equivalent damage variable d eq (x) (see Eq. ( 65))is depicted. We observe here the variation of the crack surface orientation according to the cleavage planes. To see better these effects, we also plot the the crack morphology through several planes of investigations in x-directions and directions. The results are presented in Figs. [START_REF] Musienko | Simulation of inter-and transgranular crack propagation in polycrystalline aggregates due to stress corrosion cracking[END_REF].

The 3D simulation demonstrates the performance of this frame work to simulate a very large problem in 3D in a robust manner, to study failure phenomena of polycrystalline materials.

Conclusions

A new phase field framework has been introduced in this work to simulate failure in polycrystalline materials. Both elastic anisotropy and anisotropic surface energy are taken into account in this model. Moreover, the failure in grain boundary is described by a cohesive law incorporated in the regularized variational framework. Only phase field variables are used to describe both the crack density within the grains and along the grain boundaries; it can be used as an internal variable to model irreversible damage of cohesive laws type in this region. The new model provides an efficient tool to simulate crack propagation in polycrystals, including both intergranular and transgranular mechanisms with many preferential cleavage planes. Therefore, the proposed method allows for the analysis of the interactions and the competition between intergranular and transgranular fracture.

The analysis of crack growth behavior when the crack passes through the grain boundary shows a very good agreement with observations in literature. The proposed model can reproduce both crack twisting and tilting. The analysis of polycrystalline aggregates in 2D and 3D show the ability of the model to handle complex crack shapes and the combination of intergranular and transgranular fracture modes.

The model is also able to deal with crack initiation whether within the grains or at the grain boundaries. Crack merging, coalescence and bifurcation is also possible within the proposed approach. The model invokes only a few parameters, mainly the elastic tensor and the fracture energy of the grains, and of the grain boundaries. The identification of these parameters and the validation of the model might thus be possible using the rich data provided by modern 3D imaging techniques like X-ray tomography combined with displacement measurements by digital volume correlation [35, [START_REF] Herbig | 3-d growth of a short fatigue crack within a polycrystalline microstructure studied using combined diffraction and phase-contrast x-ray tomography[END_REF]. This will be the purpose of future works.

48] and Diffraction Computed Tomography

A Smeared displacement jump approximation (JA2)

In [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF] we proposed a method to approximate the displacement jump [[u(x)]] created by interface decohesion as a smooth transition without additional variable, defined as follows:

[[u(x)]] w(x) = h∇u(x)n Γ B , (66)
where w(x) denotes the smoothed displacement jump approximation, h is a small scalar parameter and n Γ B the normal vector to Γ B at the point x.

The corresponding BVP can be rewritten as:

       ∇ • σ(u, d) = 0 ∀x ∈ Ω, γ β t(w, κ) -n Γ β • σ(u, d) = 0 ∀x ∈ Γ β , u = u ∀x ∈ ∂Ω u , n t • σ = t ∀x ∈ ∂Ω t . (67)
The Cauchy stress is here defined as follows: where the expression of elastic strain ε e is given [START_REF] Sukumar | Brittle fracture in polycrystalline microstructures with the extended finite element method[END_REF]. Using the variation in the displacement for (67) 1 , the variation in the displacement jump for (67) 2 , we obtain the corresponding weak form: (69)

            
The details of linearization and numerical implementation can be found in the work [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF].

B Anisotropic elasticity

In a polycrystalline material, the elastic moduli mentioned in [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] depend on the relative orientation of different grains. Hence the rotation of the grains must be taken into account to determine the elastic stiffness tensor for each grain. This can be achieved by using the reference transformation tensor and grains coordinate system. Suppose that the components of the elastic stiffness tensor matrix C 0 g (in Voigt notation) are given in grains coordinates, then the position-dependent elastic stiffness tensor with respect to the reference coordinate system is given by:

C 0 = P T C 0 g P, ( 70
)
where P is the transformation tensor in Voigt's notation. In the general 3D case, the transformation matrix P can be defined by the expression: where L ij are the components of the rotation matrix.

P =        L 2
For most materials treated in this work we considered an anisotropic behavior with cubic symmetry, where the elastic stiffness tensor contains only three independent parameters, the corresponding matrix in 2D reads:

C 0 g =   C 11 C 12 0 C 12 C 11 0 0 0 C 44   . ( 72
)
The corresponding bulk modulus is given as k 0 = (C 11 + 2C 12 )/3 (see [START_REF] Jamal | Elastic constants of cubic crystals[END_REF][START_REF] Mehl | Pressure dependence of the elastic moduli in aluminum-rich Al-Li compounds[END_REF]). The rotation of material orientation can be simply determined following [START_REF] Li | Phasefield modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy[END_REF]. The transformation matrix is here written as:

P =   c 2 s 2 2cs s 2 c 2 -2cs -cs cs c 2 -s 2   ( 73
)
where c = cosϕ 1 , s = sinϕ 1 .

C Cohesive zone model for grain boundary

The Cohesive Zone Model (CZM) is used to model the grain boundary decohesion. An example of traction separation law is described in Fig. 33 Fracture strength Fig. 33: Separation law used to simulate grain boundary decohesion, where t u is fracture strength, t([[u]] , κ), and Ψ B are defined in [START_REF] Herbig | 3-d growth of a short fatigue crack within a polycrystalline microstructure studied using combined diffraction and phase-contrast x-ray tomography[END_REF] In the case of polycrystalline materials, the fracture strength t u may depend on the misorientation between neighboring single crystals, defined by:

∆θ = θ i 1 -θ i 2 , (74)
where θ i 1 , θ i 2 , are the orientations of two crystals along the considered grain boundary on plane i.

Fig. 1 :

 1 Fig. 1: Regularized representation of a crack: (a) sharp crack model of a polycrystalline body containing a crack; (b) smeared representation of the crack and joint grains.

Fig. 2 :

 2 Fig.2: Benchmark problem for the analysis of surface energy distribution. The phase field variable d i is prescribed to a unit value within the disc of radius 0.05B at the center of the domain.

Fig. 3 :

 3 Fig. 3: Polar plots of the surface energy and reciprocal surface energy for the case of a material containing two preferential directions M = [1 0], M = [0 1] and α varying from 0 to 10.

Fig. 4 :

 4 Fig. 4: Approximation of the displacement jump across the grain boundary.

Fig. 5 :

 5 Fig. 5: Benchmark problem for testing the approximation of displacement jump: geometry and boundary conditions

Fig. 6 :Fig. 7 :

 67 Fig. 6: Influences of internal length , crack path comparison corresponding to prescribed displacement U = 0.115 mm.

Fig. 8 :Fig. 9 :

 89 Fig. 8: Influences of internal length β , crack path comparison corresponding to prescribed displacement U = 0.105 mm.

Fig. 10 :Fig. 11 :Fig. 12 :Fig. 13 :Fig. 14 :Fig. 15 :Fig. 16 :

 10111213141516 Fig. 10: Displacement field in the y-direction u y (x) and damage d(x) corresponding to U = 0.008 mm: (a) first displacement jump approximation (JA1), (b) second displacement jump approximation (JA2)

Fig. 17 :Fig. 18 :

 1718 Fig. 17: Comparison of crack propagation when the crack passes across the grain boundary: (a) structure ST1; (b) structure ST2 with g B c = 0.95g c ; (c) structure ST2 with g B c = 0.2g c .

Fig. 19 :

 19 Fig. 19: Configurations for the twist and tilt tests, geometry and boundary conditions: (a) 2D problem; (b) 3D problem

Fig. 20 :

 20 Fig. 20: Crack propagation in the 2D structure with two grains, phase field corresponding to: (a) U = 0.01 mm; (b) U = 0.0125 mm and (c) U = 0.013 mm.

Fig. 21 :

 21 Fig. 21: Crack propagation in the 3D structure with two grains, L × B × H = 2 × 0.3 × 1 mm 3 : the equivalent phase field corresponding to: (a) U = 0.011 mm; (b) U = 0.0125 mm; (c) U = 0.0135 mm and (d) U = 0.014.

Fig. 22 :Fig. 23 :

 2223 Fig. 22: Crack propagation in the 3D structure with two grains, L × B × H = 2 × 1.5 × 1 mm 3 : the equivalent phase field d eq (x) is plotted for: (a) U = 0.0012 mm; (b) U = 0.0135 mm and (c) U = 0.014 mm (d) U = 0.0145 mm.

Fig. 24 :

 24 Fig. 24: Definition of the linear cohesive law

Fig. 25 :Fig. 26 :

 2526 Fig. 25: Crack propagation in a 2D polycrystal for the case (C 1 ), the phase field is plotted for: (a) U = 0.52 µm; (b) U = 0.6 µm and (c) U = 0.645 µm.

  (a) Crack path for (C 1 ) g B c = 0.05 N/m, α = 50. (b) Crack path for (C 2 ) g B c = 0.05 N/m, α = 10. (c) Crack path for (C 3 ) g B c = 0.09 N/m, α = 50.

Fig. 27 : 50 Fig. 28 :

 275028 Fig. 27: Crack growth behavior in a 2D polycrystal: comparison of crack paths for the three cases

  elastic direction of elastic anisotropy and preferential damage directions are chosen randomly in each grain. The bulk modulus is k 0 = (C 11 +2C 12 )/3 = 173.3 MPa. The fracture energy of grain boundary and of bulk are set to g B c = 80 N/µm, g c = 100 N/µm, respectively. The regularization parameters are chosen as β = 0.18 µm for smeared grain boundaries and = 0.12 µm for the phase field. Four cleavage planes [

I n i t i a l c r aFig. 29 :

 29 Fig. 29: Crack propagation in 3D polycrystal containing 150 grains: geometry and boundary conditions.

Fig. 30 :

 30 Fig. 30: Detailed view of the polycrystal: overall view, surface mesh and zoom on several local regions

Fig. 31 :

 31 Fig. 31: Crack propagation in the 3D structure with 150 grains, equivalent phase field d eq (x) corresponding to: (a) U = 0.125 µm; (b) U = 0.145 µm; (c) U = 0.16 µm; (d) U = 0.175 µm.

  (a) Crack morphology at plane x = 3 µm. (b) Crack morphology at plane x = 4 µm.

  (c) Crack morphology at plane y = 4 µm.

Fig. 32 :

 32 Fig. 32: Crack morphology in the 3D structure with 150 grains, equivalent phase field d eq (x) corresponding to U = 0.175 µm.

  σ(u, d) = ∂W e u ∂ε e = g(d)C 0 + k 0 1 ⊗ 1 1 -g(d) sign -(tr ε e ) : ε e (68)

Ωσ

  : ε e (δu)dΩ -∂Ω t t • δudS = 0, Ω γ β t(w, κ)δw -σ : ε (δw) dΩ = 0.

  1.3 Compute the smeared grain boundary β(x) by solving[START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF].-FOR all loading increments (pseudo time t n+1 ):

	Given u n , d n i and H n i (x):
	2.1 Compute the strain history functional for each
	damage system H n+1 i 2.2 Compute d n+1 i (x) by solving the linear phase (x) by (37).
	field problem (54).
	2.3 Compute u n+1 (x) and v n+1 (x):
	Initialize u 2.3.1 Compute	∆ũ k+1 ∆ṽ k+1	from (50).
	2.3.2 Update	u k+1 v k+1	=	u k v k	+	∆ũ k+1 ∆ṽ k+1	.
	2.3.3 (.)

k = u n , v k = v n (displacement and displacement jump of time t n ) WHILE [∆u k+1 ∆v k+1 ]

T > , << 1: k ← (.) k+1 and go to 2.3.1.

Table 1 :

 1 Crack deflection at grain boundary tests: numerical parameters.

	Structure	h min e	(mm) h max e	(mm)	beta (mm)	(mm)
	2D structure	0.008	0.02	0.016	0.012
	3D structure	0.01	0.03	0.016	0.012

  2L 23 L 13 L 12 L 11 L 21 L 22 L 31 L 32 L 11 L 22 + L 21 L 12 L 11 L 13 L 21 L 23 L 31 L 33 L 13 L 21 + L 23 L 11 L 13 L 12 L 23 L 22 L 33 L 32 L 13 L 22 + L 23 L 12 2L 31 L 11 2L 21 L 31 2L 12 L 32 2L 32 L 22 2L 13 L 33 2L 32 L 33 L 31 L 12 + L 11 L 32 L 31 L 22 + L 21 L 32 L 11 L 33 + L 31 L 13 L 33 L 21 + L 23 L 31 L 12 L 33 + L 13 L 32 L 33 L 22 + L 23 L 32

				
				    	(71)
				
	11 L 2 12 L 2 13	L 2 21 L 2 22 L 2 23	L 2 31 L 2 32 L 2 33	2L 21 L 11 2L 22 L 12

Acknowledgements

The authors would like to acknowledge the National French Research Agency (ANR) for its financial support under contract MATETPRO ANR-12-RMNP-0020 (ECCOFIC project). The authors would like also to thank their partners MATEIS, Areva, Andra, Institut de la Corrosion, ACXCOR and MISTRAS by their participation into the fruitful discussion during this work. Finally, the fnancial support of Institut Universitaire de France (IUF) is acknowledged for J. Yvonnet.

Following the work [START_REF] Luther | Polycrystal models for the analysis of intergranular crack growth in metallic materials[END_REF], the fracture strength is then expressed by the following:

where t avg u , ∆t u are the average and the maximal fracture strength deviation, respectively.

Fig. 34: 2D description of crystal orientation between two grains

An illustration for the case of 2D plane stress is depicted in Fig. 34. In that case, the grain boundary fracture strength is written by: t u (∆θ) = t avg u + ∆t u cos(4∆θ).

(76)