# Modeling of inter- and transgranular stress corrosion crack propagation in polycrystalline material by using phase field method

Thanh-Tung Nguyen, Julien Réthoré, Marie-Christine Baietto, Jose Bolivar, Marion Fregonese, Stéphane P.A. Bordas

► To cite this version:

Thanh-Tung Nguyen, Julien Réthoré, Marie-Christine Baietto, Jose Bolivar, Marion Fregonese, et al. Modeling of inter- and transgranular stress corrosion crack propagation in polycrystalline material by using phase field method. Journal of the mechanical behavior of biomedical materials, 2018, 10.1515/jmbm-2017-0024. hal-02115952

HAL Id: hal-02115952

https://hal.science/hal-02115952

Submitted on 27 May 2020

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
Thanh-Tung Nguyen, Julien Réthoré, Marie-Christine Baietto, José Bolivar, Marion Fregonese, and Stéphane P. A. Bordas

**Modelling of inter- and transgranular stress corrosion crack propagation in polycrystalline material by using phase field method**

**Abstract:** A coupled multiphysics phase field framework is proposed to model anodic dissolution induced by stress corrosion fracture growth at microstructural level. The effects of electrochemical-mechanical processes (including crystal anisotropy) are all taken into account. This new model is based upon: (i) an anisotropic phase transformation model based on a variational formulation to describe material dissolution along preferential directions; (ii) an efficient description of grain boundaries as a smeared cohesive zone; (iii) an explicit approximation to model the different electrochemical behaviors between grain boundary and grain interior. Both intergranular and transgranular stress corrosion cracking is simulated in an efficient manner. The abilities of the proposed model are illustrated through several numerical examples involving a full prediction of complex crack network growth induced by stress corrosion cracking within 2D polycrystalline models.

**Keywords:** Stress corrosion cracking, anodic dissolution, crack propagation, polycrystals, phase field, cohesive zone model

**Communicated by:** J. Réthoré

Thanh-Tung Nguyen, Stéphane P. A. Bordas, Institute of Computational Engineering, University of Luxembourg, 6 Avenue de la Fonte, 4362 Esch-sur-Alzette, Luxembourg

Marie-Christine Baietto, Université de Lyon, CNRS, INSA-Lyon, LaMCoS UMR5259, France

José Bolivar, Marion Fregonese, Université de Lyon, CNRS, INSA-Lyon, Université Lyon 1, MatelS UMR5510, France

Julien Réthoré, Université de Nantes, Ecole Centrale de Nantes, GEM, CNRS UMR 6183 CNRS, France
# 1 Introduction

Stress corrosion cracking (SCC) is a very common cause of failures for the engineering components and structures. At the microstructural scale, cracks often nucleate from corrosion pits and, depending on the kind of material and environment, it could be inter- or transgranular (IG-TG) in nature [1]. Especially, crystals orientation, grain boundary properties and mechanical loading strongly affect the switching between the inter- and transgranular stress corrosion cracking (IGSCC-TGSCC) [2, 3]. This competition provides a significant difference in the material resistance [4, 5, 6]. Therefore, understanding these aspects plays an important role for improvements in the development and application of advanced structural materials.

Many works in the literature have been proposed to deal with the simulation of microstructural effects on the stress corrosion assisted fracture growth in polycrystals. Cohesive Zone Models (CZM), based on the original idea in the work of Dugdale [7] and Barenblatt et al [8] to avoid the singularity of stress at the crack tip, is a popular technique used in this field. One of its application can be found in the study [9], where a grain-scale finite element (FE) model based on the CZM is proposed to study IGSCC in stainless steel wire. In the work [10], a 3D grain-level cohesive finite element model is proposed to study the IGSCC characteristics of AZ31 magnesium alloy due to the hydrogen embrittlement process. Another contribution based on cohesive zone domains can be found in [11] by using multiscale model to simulate the intergranular hydrogen-assisted cracking. The investigation of the microstructural effects (especially the grain boundaries) on the behaviour of the IG stress corrosion assisted short cracks can be found in [12] for the alloy 600 by using the body force method, and [13] in the fully sensitised austenitic stainless steel by using the polycrystalline aggregate model, or in [14] by using a new model based on Markov Chain theory and Monte Carlo simulations. The polycrystalline model is also used in [15, 16] to simulate the IGSCC in stainless steel. Another technique, namely dislocation-based boundary element model is proposed in [17] to microstructurally study intergranular stress corrosion short cracks in an elastoplastic medium.

Development of the numerical models at the grain size scales has so far been mostly limited to simulate IGSCC. There is a very few studies dealing with both IGSCC and TGSCC such as the work of Musienko et al [6]. The authors proposed a coupled approach (environment and plasticity) based on the CZM to model the IG-TG fracture transition. More recently, a numerical multiphysics peridynamic framework [18] based on the adsorption induced decohesion mechanism is proposed to this field. However, a detailed description of the electrochemical and
mechanical effects on stress corrosion assisted damage is still missing, especially for the common situation of anodic dissolution induced SCC.

In this work, the phase field method augmented by a smeared description of interfaces is used to develop a fully coupled multiphysics model of stress corrosion assisted damage growth [19, 20, 21]. Starting from the ideas in [22, 23], we extend them to the case of polycrystalline materials, where the anisotropic effects for both elastic energy and fracture surface energy are taken into account. The behavior of grain boundaries is described by the smeared cohesive model following the work [24, 25]. Moreover, to account for the different behavior of grain boundary and grain interior, we propose an interpolation formulation to present the smoothed transition of the electrochemical-mechanical properties between the two regions. It is expected that the new model can model both IGSCC and TGSCC depending on the electrochemical-mechanical loading.

The overview of the paper is as follows. In section 2 the fundamentals of the proposed model are presented. In section 3, an extented formulation of phase field method able to model the stress corrosion induced crack propagation accounted the crystal effects is proposed. The variational principle for the proposed model is then detailed in section 4. Finally, the potential of the proposed model is illustrated by several numerical examples in section 5.

## 2 Fundamentals of the proposed method

Let $\Omega \subset \mathbb{R}^D$ be an open domain describing a corroding system at micro scale, which generally contains four domains: (i) metal or solid phase described by temporally and spatially constant metal atom concentration $c_s^{\text{alloy}}$; (ii) electrolyte solution or liquid phase described by temporally and spatially constant electrolyte atom concentration $c_l^{\text{alloy}}$; (iii) an interfacial region where the corrosion occurs with the atom concentration taken within the range $(c_l^{\text{alloy}}, c_s^{\text{alloy}})$; (iv) and the grain boundaries region described by metal atom concentration $c_s^{\text{alloy}}$ as well, but with another behavior.

The normalization procedure is then adopted following the work [22]. We define $c_s$ as the normalized solid concentration and $c_l$ as the normalized electrolyte concentration, reading:

$$c_s = \frac{c_s^{\text{alloy}}}{c_s^{\text{alloy}}} \quad \text{and} \quad c_l = \frac{c_l^{\text{solution}}}{c_s^{\text{alloy}}}. \qquad (1)$$

The corrosion process is modeled by using the phase field method. A scalar parameter $\phi(\mathbf{x})$ is considered to describe the corroding system. $\phi(\mathbf{x}) = 1$ is
solid phase, $\phi(\mathbf{x}) = 0$ is liquid phase and the interfacial region is within the range $0 < \phi(\mathbf{x}) < 1$. Moreover, the grain boundaries are also approximated by a smeared representation defined by a temporally constant $\beta(\mathbf{x})$, taking a unit value on grain boundaries and vanishing away from it<fnref n="1" />. An illustration of a two-dimensional corroding system is depicted in Fig. 1

<figure><img src="image_2.png" /><figcaption>**Fig. 1:** 2D representation of a typical corroding system due to the anodic dissolution of a polycrystalline body</figcaption></figure>

The dissolution of material (phase transformation: solid $\phi = 1$ to liquid $\phi = 0$) is provided by the diffusional transport of metal ions from solid domain to liquid domain at interfacial domain. Note that, the interfacial region (transition domain) is represented here as a mixture of both phases (solid - liquid) with different chemical compositions by using the assumption of (Kim Kim Suzuki) KKS model [26]. Then, by introducing the interpolation function, satisfying $h(\phi = 0) = 0$,

<footnote n="1">1 The smeared grain boundaries are here described by introducing a grain boundaries density function $\gamma_{\beta}(\beta, \nabla\beta) = \frac{1}{2\ell_{\beta}}\beta(\mathbf{x})^2 + \frac{\ell_{\beta}}{2}\nabla\beta(\mathbf{x}) \cdot \nabla\beta(\mathbf{x})$, where $\ell_{\beta}$ is the regularization parameter. Then the total grain boundaries length represents by $\Gamma_{\beta}(\beta) = \int_{\Omega} \gamma_{\beta}(\beta, \nabla\beta) \, d\Omega$. Note that for $\ell_{\beta} \to 0$, this description leads to the exact description of the sharp cohesive interface, i.e $\Gamma_{\beta} \to \Gamma_B$ (see e.g. [25, 24] for more details)</footnote>
$h(\phi = 1) = 1$ and $\frac{\partial h(\phi)}{\partial \phi}|_{\phi=0}, \phi=1 = 0$; one choice: $h(\phi) = -2\phi^3 + 3\phi^2$, the ions concentration at any material point $c$ can be described by the following:

$$c = h(\phi)c_s + [1 - h(\phi)]c_l, \qquad (2)$$

As mentioned in the work [23], the phase transformation occurs when the ions concentration achieves the saturation value, denoted $c_{se}$ and $c_{le}$ for solid metal and electrolyte respectively, in the normalized framework.

## 3 Phase field model

This work focuses on the study of polycrystalline effects to the anodic dissolution induced stress corrosion crack growth. To this aim, we extend the formulation proposed in [23] to the case of polycrystalline materials. The new energy formulation is rewritten as follows:

$$E(\mathbf{u}, \llbracket \mathbf{u} \rrbracket, \phi, c) = \int_{\Omega \setminus \Gamma_B, \Gamma_\phi} \psi_u(\boldsymbol{\varepsilon}(\mathbf{u})) \, d\Omega + \int_{\Gamma_B} \psi_v(\llbracket \mathbf{u} \rrbracket, \kappa) \, d\Gamma \\ + \int_{\Gamma_\phi} \psi_\phi \, d\Gamma + \int_{\Omega} \psi_c(c, \phi) \, d\Omega, \qquad (3)$$

where $\psi_u(\boldsymbol{\varepsilon}(\mathbf{u}))$ is the strain energy density, $\psi_c(c, \phi)$ is the electrochemical energy density and $\psi_\phi$ describes the fracture energy. The term $\psi_v(\llbracket \mathbf{u} \rrbracket, \kappa)$ is added to account for the energy of grain boundaries, which is a function depending on the displacement jump across grain boundaries $\Gamma_B$ and an history parameter $\kappa$.

In the phase field framework, the regularized description for strong discontinuities related to both fracture and grain boundaries is adopted (substituting $\llbracket \mathbf{u} \rrbracket$ by $\mathbf{v}(\mathbf{x})$). A sharp grain boundary is regularized by a smeared grain boundary (substituting $\Gamma_B$ by $\Gamma_\beta = \int_\Omega \gamma_\beta(\beta, \nabla\beta) \, d\Omega$). Then its decohesion will be described by the smeared cohesive law proposed in [25]. This approach is based on a smeared displacement jump approximation (see [27, 24]), in which the grain boundary decohesion is approximated as a smooth transition $\mathbf{v}(\mathbf{x})$. That provides a decomposition of infinitesimal strain tensor into a part related to the bulk and a part induced by smoothed jumps, denoted by $\boldsymbol{\varepsilon}^e$ and $\tilde{\boldsymbol{\varepsilon}}$, respectively.

$$\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}^e + \tilde{\boldsymbol{\varepsilon}}. \qquad (4)$$
with

$$ \tilde{\varepsilon} = \mathbf{n}^{\Gamma\beta} \otimes^s \mathbf{v}^{\gamma\beta}, \qquad (5) $$

where $\mathbf{n}^{\Gamma\beta}$ is the normal vector to $\Gamma\beta$.

The strain energy density of damageable material $\psi_u(\varepsilon^e(\mathbf{u}), \phi)$ is then constructed from elastic part $\varepsilon^e$ by using unilateral contact assumptions introduced in [24]:

$$ \psi_u(\varepsilon^e(\mathbf{u}), \phi) = \frac{1}{2} [\varepsilon^e : \mathcal{C}(\phi) : \varepsilon^e], \qquad (6) $$

with the elastic tensor $\mathcal{C}(\phi)$ is defined by the following:

$$ \mathcal{C}(\phi) = h(\phi)\mathcal{C}^0 + k_0 \mathbf{1} \otimes \mathbf{1} [1 - h(\phi)] \mathrm{sign}^{-}(\mathrm{tr} \, \varepsilon^e), \qquad (7) $$

where $\mathcal{C}^0$ denotes the initial elastic tensor of the material, possibly anisotropic; $k_0$ is the bulk modulus; the sign function $\mathrm{sign}^{-}(x) = 1$ if $x < 0$ and $\mathrm{sign}^{-}(x) = 0$ if $x \ge 0$.

The energy formulation in the phase field framework is now defined by:

$$ \begin{aligned} E(\mathbf{u}, \mathbf{v}, \phi, c) = & \int_{\Omega} \psi_u(\varepsilon^e(\mathbf{u}), \phi) \, d\Omega + \int_{\Omega} \psi_v(\mathbf{v}, k) \gamma_{\beta} \, d\Omega \\ & + \int_{\Omega} \psi_{\phi}(\phi, \nabla\phi) \, d\Omega + \int_{\Omega} \psi_c(c, \phi) \, d\Omega. \end{aligned} \qquad (8) $$

The third term in Eq. (8), $\psi_{\phi}(\phi, \nabla\phi)$ describes the surface energy or the amount of energy released upon the creation of new fracture surfaces. It is chosen following the work [23]:

$$ \psi_{\phi}(\phi, \nabla\phi) = \omega_{\phi}\phi^2(1 - \phi^2) + \alpha_{\phi} : (\nabla\phi \otimes \nabla\phi), \qquad (9) $$

with $\omega_{\phi}$ is the height of the imposed double-well energy barrier, $\alpha_{\phi}$ is the second-order tensor describing gradient energy (included the anisotropic effects), which is an invariant with respect to rotations:

$$ \alpha_{\phi} = \alpha_0 [\mathbf{1} + \xi(\mathbf{1} - \mathbf{N} \otimes \mathbf{N})], \qquad (10) $$

where $\alpha_0$ is the classical gradient energy coefficient for the case of isotropic material; $\mathbf{N}$ denotes the unit vector normal to the preferential cleavage plane (with respect to the material coordinates), and $\xi \gg 0$ is used to penalize the
damage on planes not normal to N. This formulation can be extended for the case of many preferential cleavage planes. The detailed description can be found in the work [28].

Note that $\omega_{\phi}$ and $\alpha_0$ can be determined from the relation with the transition length $\ell$ and the interface energy $\sigma_{\phi}$ or the critical fracture energy $g_c$:

$$ \sigma_{\phi} = \frac{\sqrt{\alpha_0 \omega_{\phi}}}{3\sqrt{2}} \quad \text{and} \quad \ell = \alpha^* \sqrt{\frac{2\alpha_0}{\omega_{\phi}}}, \quad \text{or} \quad \omega_{\phi} = \frac{g_c}{2\ell}, \quad \alpha_0 = \frac{g_c \ell}{2}, \qquad (11) $$

where $\alpha^* \approx 2.94$ corresponds to the transition zone taken within the range $0.05 < \phi < 0.95$.

The electrochemical energy density $\psi_c(c_s, c_l, \phi)$ in Eq. 8 is chosen based on the KKS model [26], wherein the transition region is postulated as a mixture of the two phases with different chemical compositions, but equal chemical potentials:

$$ \psi_c(c_s, c_l, \phi) = h(\phi) [A(c_s - c_{se})^2] + [1 - h(\phi)] [A(c_l - c_{le})^2], \qquad (12) $$

where A is the free energy density curvature which is often identified from thermodynamical databases [29, 30]; $c_{se}$ and $c_{le}$ are the equilibrium concentrations of the coexisting phases; the concentration of solid $c_s$ and liquid $c_l$ phases can be expressed from the equilibrium condition at the transition zone as follows:

$$ \left\{ \begin{array}{l} c_s = c + [h(\phi) - 1](c_{le} - c_{se}) \\ c_l = c + h(\phi)(c_{le} - c_{se}). \end{array} \right. \qquad (13) $$

The total energy is then rewritten as $E(\mathbf{u}, \mathbf{v}, \phi, c) = \int_{\Omega} \psi \, d\Omega$ in which

$$ \psi = \psi_u + \psi_v + \psi_{\phi} + \psi_c \qquad (14) $$

is identified total energy density.

# 4 Variational principle

## 4.1 Variational principle for mechanical problem

The weak form of mechanical problem is obtained by taking an arbitrary variation in the displacement field $\delta u$ and in the auxiliary displacement jump field $\delta v$,
while the phase field is kept constant. In the absence of body forces, it is defined as:

$$ \int_{\Omega} \left( \boldsymbol{\sigma} : \delta\boldsymbol{\varepsilon}^e + \gamma_{\beta} \mathbf{t}(\mathbf{v}, \kappa) \delta v + \zeta \frac{\partial \mathbf{v}}{\partial x_n} \frac{\partial \delta v}{\partial x_n} \right) d\Omega - \int_{\partial\Omega_t} \bar{\mathbf{t}} \cdot \delta u \, dS = 0 \quad (15) $$

where $\mathbf{t}(\mathbf{v}, \kappa)$ are the tractions at grain boundaries with outward normal $\mathbf{n}^{\Gamma_{\beta}}$; $\bar{\mathbf{t}}$, $\bar{\mathbf{u}}$ are the prescribed tractions and the prescribed displacements on $\partial\Omega_t$, $\partial\Omega_u$, respectively; $x_n = (\mathbf{x} - \mathbf{x}_{\Gamma_{\beta}}) \cdot \mathbf{n}^{\Gamma_{\beta}}$ and $\mathbf{x}_{\Gamma_{\beta}} = \underset{\mathbf{y} \in \Gamma_{\beta}}{\operatorname{argmin}} (\|\mathbf{y} - \mathbf{x}\|)$, the positive constant $\zeta$ is introduced to ensure that the auxiliary displacement jump field is constant in the normal direction [27]. $\boldsymbol{\sigma}$ is the Cauchy stress in the bulk material, defined by:

$$ \sigma(\mathbf{u}, \mathbf{v}, \phi) = \frac{\partial \psi_u^e}{\partial \varepsilon^e} = \mathcal{C}(\phi) : \varepsilon^e \quad (16) $$

where the expression of the elastic strain $\varepsilon^e$ is given by (see [25]):

$$ \boldsymbol{\varepsilon}^e = \nabla^s \mathbf{u} - \mathbf{n}^{\Gamma_{\beta}} \otimes^s \mathbf{v} \gamma_{\beta}, \quad (17) $$

with $(\nabla^s \mathbf{u})_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i})$ and $(\mathbf{n}^{\Gamma_{\beta}} \otimes^s \mathbf{v})_{ij} = \frac{1}{2}(n_i v_j + v_i n_j)$. From Eqs. 15, 17, the weak form can be rewritten as:

$$ \int_{\Omega} \left( (-\nabla \cdot \boldsymbol{\sigma}) \delta u - \left[ \gamma_{\beta} (\boldsymbol{\sigma} \mathbf{n}^{\Gamma_{\beta}} - \mathbf{t}(\mathbf{v}, \kappa)) + \zeta \frac{\partial^2 \mathbf{v}}{\partial x_n^2} \right] \delta v \right) d\Omega \\ - \int_{\partial\Omega_t} \left( \bar{\mathbf{t}} - \boldsymbol{\sigma} \mathbf{n}_t \right) \delta u \, dS - \int_{\partial\Gamma_{\beta}} \frac{\partial \mathbf{v}}{\partial x_n} \delta v \, dS_{\beta} = 0 \quad (18) $$

In the case the crack curvature is assumed to be small, by considering a variation in the displacement and displacement jump, yields the BVP:

$$ \left\{ \begin{array}{l} \nabla \cdot \boldsymbol{\sigma}(\mathbf{u}, \mathbf{v}, \phi) = 0 \quad \forall \mathbf{x} \in \Omega, \\ \gamma_{\beta} [\mathbf{t}(\mathbf{v}, \kappa) - \mathbf{n}^{\Gamma_{\beta}} \cdot \boldsymbol{\sigma}(\mathbf{u}, \mathbf{v}, \phi)] = \zeta \frac{\partial^2 \mathbf{v}}{\partial(x_n)^2} \quad \forall \mathbf{x} \in \Gamma_{\beta}, \\ \frac{\partial \mathbf{v}}{\partial x_n} = 0 \quad \forall \mathbf{x} \in \partial\Gamma_{\beta}, \\ \mathbf{u} = \bar{\mathbf{u}} \quad \forall \mathbf{x} \in \partial\Omega_u, \\ \mathbf{n}_t \cdot \boldsymbol{\sigma} = \bar{\mathbf{t}} \quad \forall \mathbf{x} \in \partial\Omega_t, \end{array} \right. \quad (19) $$
The corresponding mechanical problem will be solved:

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \boldsymbol{\sigma} : \boldsymbol{\varepsilon}^e (\delta u) d\Omega - \int_{\partial\Omega_t} \bar{\boldsymbol{t}} \cdot \delta u dS = 0 \\[1.5em] \displaystyle \int_{\Omega} \left\{ \gamma_{\beta} \left[ \boldsymbol{t}(\boldsymbol{v}, \kappa) \delta v - \boldsymbol{\sigma} : \tilde{\boldsymbol{\varepsilon}} (\delta v) \right] - \zeta \frac{\partial \boldsymbol{v}}{\partial x_n} \frac{\partial \delta v}{\partial x_n} \right\} d\Omega - \int_{\partial\Gamma_{\beta}} \frac{\partial \boldsymbol{v}}{\partial x_n} \delta v dS_{\beta} = 0 \end{array} \right. \quad (20) $$

## 4.2 Variational principle for electrochemical problem

The Cahn-Hilliard and Allen-Cahn equation is used to describe the electrochemical processes of phase transformation. The diffusion of ions from solid volume to electrolyte, or the spatial evolution of conserved field $c$ is governed by a conservation law (Cahn-Hilliard). The material dissolution or the evolution of non conserved field $\phi$ is governed by Allen-Cahn equation:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial c}{\partial t}(\mathbf{x}, t) = \nabla \cdot M \left( \nabla \frac{\partial \psi}{\partial c} \right), \\[1.5em] \displaystyle \frac{\partial \phi}{\partial t}(\mathbf{x}, t) = -L_{\phi} \left( \frac{\partial \psi}{\partial \phi} - \nabla \cdot \frac{\partial \psi}{\partial \nabla \phi} \right), \end{array} \right. \quad (21) $$

where $L_{\phi}$ is the interface kinetics parameter, that describes the transformation rate; $M$ is the diffusion mobility, which can be evaluated based on an analogy between the Cahn-Hilliard equation and the Fick's second law (see [23, 22], for more details), being written as:

$$ M = \frac{D}{2A}, \quad (22) $$

with $D$ is the diffusion coefficient. An important property of polycrystalline materials needs to be considered is the different SCC phenomena between grain boundary and grain interior. The grain boundary atoms are more easily and rapidly dissolved, or corroded, than the atoms within the grains. So that, the grain boundaries oxidize or corrode more rapidly. To model this phenomenon, the diffusion coefficient $D$ is here defined as a function of the smeared description of grain boundary $\beta$ by the following:

$$ D = \frac{D_G}{(1 - \beta)^n (1 - k_D) + k_D} \quad (23) $$
where $D_G$ is the diffusion coefficient for grain interior; $n$ is an exponent; $k_D$ is the diffusion mismatch ratio between grain interior and grain boundary. This definition ensures $D = D_G$ in the grain interior, $D = D_B = D_G/k_D$ in the center of grain boundary and $D_G < D < D_B$ in the smeared region. An illustration of the variation of the diffusion coefficient crossing the grain boundary for the case of $n = 2$ and $k_D = 0.1$ is depicted in Fig. 2. This approximation can be understood as an interpolation procedure for two different regions. By using the same principle, the fracture energy for whole domain $g_c$ can be written:

$$g_c = [(1 - \beta)^n (1 - k_g) + k_g] g_c^G \quad (24)$$

where $g_c^G$ is the fracture energy in the grain interior, $k_g$ is fracture energy mismatch between the grain boundary $g_c^B$ and the grain interior $g_c^G$, defined by $k_g = g_c^B/g_c^G$.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Explicit description of the different electrochemical behavior between grain boundary and grain interior</figcaption></figure>
Using the variation in the concentration for Eq. (21)$_1$, and the variation in the phase field for Eq. (21)$_2$, the corresponding weak form was obtained:

$$ \left\{ \begin{aligned} & \int_{\Omega} \frac{\partial c}{\partial t} \delta c \, d\Omega + \int_{\Omega} M \nabla \frac{\partial \psi_c}{\partial c} \cdot \nabla \delta c \, d\Omega = 0 \\ & \int_{\Omega} \frac{\partial \phi}{\partial t} \delta \phi \, d\Omega + \int_{\Omega} L_{\phi} \delta \phi \psi \, d\Omega = 0 \end{aligned} \right. \qquad (25) $$

where $\delta_{\phi}\psi$ reads

$$ \begin{aligned} \delta_{\phi}\psi &= \left( \frac{\partial\psi}{\partial\phi} - \nabla \cdot \frac{\partial\psi}{\partial\nabla\phi} \right) \delta\phi \\ &= \left[ \frac{\partial\psi_c}{\partial\phi} + [\boldsymbol{\epsilon}^e : \frac{\partial\mathcal{C}(\phi)}{\partial\phi} : \boldsymbol{\epsilon}^e] + \omega_{\phi}g'(\phi) \right] \delta\phi + \nabla\phi \left[ \alpha_0 (\mathbf{1} + \xi(\mathbf{1} - \mathbf{N} \otimes \mathbf{N})) \right] \nabla\delta\phi. \end{aligned} \qquad (26) $$

In order to handle loading and unloading histories, we introduce the history function following the work of [31, 32]

$$ \mathcal{H}(\mathbf{x}, t) = \max_{\tau \in [0, t]} \left\{ \boldsymbol{\epsilon}^e(\mathbf{x}, \tau) : [\mathbb{C}^0(\mathbf{x}) - k_0(\mathbf{x}) \mathbf{1} \otimes \mathbf{1} \operatorname{sign}^{-}(\operatorname{tr} \boldsymbol{\epsilon}^e(\mathbf{x}, \tau))] : \boldsymbol{\epsilon}^e(\mathbf{x}, \tau) \right\} \quad (27) $$

From (26), (27), the final weak form for electrochemical processes is obtained by the following equation:

$$ \left\{ \begin{aligned} & \int_{\Omega} \frac{\partial c}{\partial t} \delta c \, d\Omega &&+ \int_{\Omega} M \nabla \frac{\partial \psi_c}{\partial c} \cdot \nabla \delta c \, d\Omega = 0 \\ & \int_{\Omega} \frac{\partial \phi}{\partial t} \delta \phi \, d\Omega &&+ \int_{\Omega} L_{\phi} \left[ \frac{\partial \psi_c}{\partial \phi} + h'(\phi) \mathcal{H}(\mathbf{x}, t) + \omega_{\phi} g'(\phi) \right] \delta \phi \, d\Omega \quad (28) \\ & + \int_{\Omega} L_{\phi} \nabla \phi \left[ \alpha_0 (\mathbf{1} + \xi (\mathbf{1} - \mathbf{N} \otimes \mathbf{N})) \right] \nabla \delta \phi \, d\Omega = 0 \end{aligned} \right. $$
## 5 Representative numerical examples

In this section, we present several numerical examples using the proposed method to study SCC of a zirconium alloy. Thanks to a tightly adherent and protective oxide film, zirconium is well known for its excellent corrosion resistance to a lot of acids (hydrochloric acid, hot organic acids, sulfuric acid,...) even when exceeding the normal boiling point temperature. Moreover, this material exhibits very good mechanical properties. However, zirconium can be corroded in the situation of hydrochloric acid containing oxidizing species such as cupric chloride, ferric chloride, or wet chlorine. Both Intergranular and Transgranular Stress Corrosion Cracking (IGSCC - TGSCC) were experimentally observed in the literature [33], [6].

In this study, anodic dissolution induced stress corrosion crack propagation of this material is investigated. The considered zirconium alloy is composed mostly of zirconium, tin, iron and chromium with a typical mass density of $\rho_{\text{alloy}} = 6.56 \text{ g/cm}^3$. The complete chemical composition is given in Table. 1. The metal concentration is finally obtained $c_s^{\text{alloy}} = 72.23 \text{ mol/l}$ (see the work [23] for more details).

<table><caption>**Table 1:** Chemical composition of the considered zirconium alloy (Zircaloy-4)</caption><thead><tr><th>Alloy</th><th>Element</th><th>Mass fraction, f<sub>m,i</sub>[%]</th><th>Molar mass, M[g/mol]</th><th>Mole fraction f<sub>M,i</sub>[%]</th></tr></thead><tbody><tr><td rowspan="6">zirconium alloy</td><td>Cr</td><td>0.10</td><td>52.00</td><td>0.175</td></tr><tr><td>Fe</td><td>0.20</td><td>55.85</td><td>0.325</td></tr><tr><td>O</td><td>0.12</td><td>16.00</td><td>0.681</td></tr><tr><td>Sn</td><td>1.40</td><td>118.71</td><td>1.071</td></tr><tr><td>Zr</td><td>98.50</td><td>91.22</td><td>97.748</td></tr></tbody></table>

The diffusion coefficient is assumed to be homogeneous for all chemical compositions, and taken as the value for pure Zr, with $D_G = 1.3 \times 10^{-8} \text{m}^2/\text{s}$ [34, 35]. The saturated concentration is taken as in the study [36]: $c_{\text{sat}} = 3 \text{ mol/L}$. Hence, the saturated concentration of solid and liquid phase are $c_{\text{se}} = 1$, $c_{\text{te}} = 3/72.23$. The free energy density curvature $A$ is identified by fitting the chemical free energies obtained from thermodynamic databases. That gives $A = 8.02 \times 10^8 \text{J/mol}$. The size of transition zone and the interface energy are
chosen $\ell = 0.5$ µm and $\sigma_{\phi} = 0.5$ kN/m (giving the equivalent fracture energy $g_c^G = 10^{-3}$ kN/mm). The interface kinetics energy is $L_{\phi} = 0.15 l/(J.s)$.

The grain boundary cohesive energy is taken following the work of Xu et al [37]:

$$ \Psi^B(\llbracket u \rrbracket) = g_c^B \left[ 1 - \left( 1 + \frac{\llbracket u \rrbracket}{v_n} \right) \exp\left(-\frac{\llbracket u \rrbracket}{v_n}\right) \right]. \quad (29) $$

The normal traction separation law is obtained as:

$$ t(\llbracket u \rrbracket) = g_c^B \frac{\llbracket u \rrbracket}{v_n^2} \exp\left(-\frac{\llbracket u \rrbracket}{v_n}\right), \quad (30) $$

with $v_n = g_c^B/(t_u \exp(1))$, and $t_u$ the fracture strength, which can be dependent on the angular mismatch between grains (see eg. [38]) as follows:

$$ t_u(\Delta\theta^i) = t_u^{\text{avg}} + \frac{1}{3}\Delta t_u \sum_i \cos(4\Delta\theta^i), \quad (31) $$

where $t_u^{\text{avg}}$, $\Delta t_u$ are the average and the maximal fracture strength deviation, respectively; the angular mismatch $\Delta\theta$ is defined as the difference between the orientations of two crystals $\theta_1^i, \theta_2^i$, along the considered grain boundary on plane i:

$$ \Delta\theta = \theta_1^i - \theta_2^i. \quad (32) $$

However, for the sake of simplicity and due to the lack of experimental data, the grain boundary fracture strength is considered here as independent for misorientation between two grains and taken as $t_u = 0.1$ GPa. The fracture energy within the grain boundary is $g_c^B = 2.85 \times 10^{-4}$ kN/mm. The elastic tensor of hcp symmetry is used. For 0° orientation, it is written as follows [39]:

$$ \mathbf{C}^0 = \begin{bmatrix} 142 & 73.4 & 0 \\ 73.4 & 142 & 0 \\ 0 & 0 & 34.3 \end{bmatrix} \text{(GPa)}. \quad (33) $$

The bulk modulus is defined as $k_0 = (C_{11}^0 + 2C_{12}^0)/3 = 96.26$ GPa. For easier numerical implementation, all parameters are normalized following the work by [29, 30].

We first consider a rectangular domain with dimension $40 \times 30$ µm² containing an initial crack and 10 grains. The geometry and mesh of this microstructure is
generated using the open source software Neper [40] (developed by Romain Quey at CNRS and MINES Saint-Etienne, France). A refined mesh was constructed using triangular elements, with element size ~ 0.15 µm, to have 6-8 elements in the transition zone. The average grain size is about $\sqrt{\frac{40 \times 30}{10}} \approx 11$ µm. The details of the geometry and of boundary conditions are described in Fig. 3. In each grain, both direction of anisotropy and preferential damage directions are generated randomly.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Geometry and boundary conditions for the polycrystalline structure containing 10 grains: (a) grain boundary and loading description; (b) description of crystals with angle of the direction of anisotropy.</figcaption></figure>

In the crystallographic orientation, the cleavage plane $N = [1 \ 0]$ is chosen. The structure will be numerically assigned to both electrochemical degradation and mechanical loading. The displacements are prescribed along the y-direction for upper edge ($y = 30$ µm) and lower edge ($y = 0$) while the displacement along x is here free. The rupture of passive film allowing electrochemical processes is manually activated from the beginning at the initial crack location. Plane strain condition is assumed.

The competition between IGSCC and TGSCC is analyzed by considering two situations:

(C1) SCC with low mechanical loading.

(C2) SCC with high mechanical loading.
The first case (C1) is set up to expect the IGSCC, because at low mechanical loading, SCC process, taking place preferentially at grain boundaries, prevails on mechanical process. Therefore IGSCC will be dominant in this situation. In the second case (C2), when the structure is subjected to higher mechanical loading, the mechanical process becomes dominant. So that, the cracks will propagate not only on grain boundaries but also within the grains depending on the direction which is suitable for fracture growth. The details of boundary conditions are given in table 2.

<table><caption>**Table 2:** Loading details of two situations: (C1) and (C2)</caption><thead><tr><th rowspan="2">Case</th><th colspan="3">V (μm/min)</th></tr><tr><th>Period [0 : 3 min]</th><th>Period [3 : 120 min]</th><th>Period [120 : 300 min]</th></tr></thead><tbody><tr><td>(C1)</td><td>0.03</td><td>1.35 × 10<sup>-3</sup></td><td>6 × 10<sup>-4</sup></td></tr><tr><td>(C2)</td><td>0.06</td><td>1.80 × 10<sup>-3</sup></td><td>6 × 10<sup>-4</sup></td></tr></tbody></table>

The resulting crack propagation for the first case (C1) is depicted in Fig. 4. The initial crack propagates towards the grain boundaries, according to the preferential direction (Fig. 4 with $U = 0.21$ μm). The strong influence of SCC is captured in this loading condition. The fracture width near the edge $x = 0$ is much larger than the region nearly crack tip, because of material dissolution induced by SCC had more time to develop in this zone. After reaching the grain boundary, the crack has a tendency to propagate along the grain boundaries (Fig. 4 with $U = 0.27, 0.30$ and $0.32$ μm). IGSCC is clearly observed.

A different fracture scenario is obtained in the second case (C2) (see Fig. 5). With the higher loading rate, the mechanical process becomes predominant. The fracture morphology is shaper, but the effect of SCC is still captured at the edge $x = 0$. The fracture growth in the grain boundary is always observed (see Fig. 5 with $U = 0.28$ μm), but then the crack in this situation mostly propagates within the grains (see Fig. 5 with $U = 0.35$ and $0.39$ μm). Under the higher mechanical loading, fracture growth is essentially due to mechanical process, that supports for both IG and TG cracks (not like only IG fracture in loading condition C1). Thus cracks choose either IG or TG, depending on which one is more favorable energetically. This phenomenon has been proved in Fig. 5 for $U = 0.41$ μm, where the orientation of the last grain is not suitable for crack propagation, the fracture has switched from TG to IG. The obtained numerical simulation results
<figure><img src="image_5.png" /><figcaption>(a) $U = 0.21 \ \mu\text{m}$ ($t = 115 \ \text{min}$)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) $U = 0.27 \ \mu\text{m}$ ($t = 205 \ \text{min}$)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) $U = 0.30 \ \mu\text{m}$ ($t = 255 \ \text{min}$)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) $U = 0.32 \ \mu\text{m}$ ($t = 285 \ \text{min}$)</figcaption></figure>

<caption>Fig. 4: Crack propagation for the first situation (C1) corresponding to the low mechanical loading. The IGSCC is clearly captured</caption>

show a very good qualitative agreement with experimental observations in the literature [41, 42, 6].

Note that the transition between IG and TG fracture is strongly affected by the interfacial properties at grain boundaries as well. By using the same mechanical loading condition as the case (C2) but taking a new grain boundary energy $g_c^B = 8.3 \times 10^{-4} \ \text{kN/mm}$ (denoted here the case (C3)), TG crack propagation is obtained as in Fig. 6(c). The crack trajectory of this last case (C3) is quite similar to one obtained from the loading case (C2), but with lower grain boundary effects. Especially at the last loading step, the fracture switches from TG to IG for the case (C2), while it is still TG for the case (C3). This result confirms the influence of the fracture energy mismatch between grain boundaries and grain interior on the transition between IGSCC and TGSCC.

The stress-displacement curves for all situations are plotted in Fig. 7. The elastic period for the three cases exhibits the same behavior. Then, the IGSCC in the case (C1) shows the weakest resistance, while the fully TGSCC in the case (C3) presents the strongest resistance. These observations demonstrate the
<figure><img src="image_9.png" /></figure>

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

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

<figure><img src="image_12.png" /><figcaption>Fig. 5: Crack propagation for the second situation (C2) corresponding to the high mechanical loading. Both IGSCC and TGSCC are observed</figcaption></figure>

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

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

<figure><img src="image_15.png" /><figcaption>Fig. 6: Comparison of crack path for different cases: (a) case study (C1) with the low mechanical loading and $g_c/g_c^B = 3.5$; (b) case study (C2) with the high mechanical loading and $g_c/g_c^B = 3.5$; (c) case study (C3) with the high mechanical loading and $g_c/g_c^B = 1.2$</figcaption></figure>

strong influence of the fracture features (IG or TG) on the macroscopic material properties, which plays a major role in the material design.

In the next example, we perform the numerical simulation on a polycrystalline structure containing 50 grains (dimension: 160 × 120 µm²) is performed in order to demonstrate the efficiency of the proposed framework. The details of geometry and boundary conditions are depicted in Fig. 8. The same values of the
<figure><img src="image_16.png" /><figcaption>Fig. 7: Comparison of the stress-displacement curve for both different mechanical loading conditions</figcaption></figure>

material parameters are considered as for the previous example. The boundary conditions are chosen here to capture both IGSCC and TGSCC as follows: the prescribed displacement rate $V = 0.1 \ \mu\text{m/min}$ for the period $t = [0 : 5 \ \text{min}]$ and $V = 0.002 \ \mu\text{m/min}$ for $t = (5 : 500 \ \text{min}]$. Only the cleavage plane $N = [1 \ 0]$ is considered.

<figure><img src="image_17.png" /><figcaption>Fig. 8: Geometry and boundary conditions for the polycrystalline structure containing 50 grains: (a) grain boundary and loading description; (b) angle of the direction of anisotropy</figcaption></figure>

The results of fracture evolution corresponding to different loading times are plotted in Fig. 9. Both crack in the grain interior and along the grain boundaries are observed. Note that the advantage of applying the phase field model is demonstrated, the model is able to capture the nucleation of several cracks when
the structure is subjected to predominant mechanical loading (see Fig. 9 with U = 1.40 µm). These cracks often nucleate at the grain boundaries (the weakest regions) and then either propagate along grain boundaries or onset and penetrate the grain (see Fig. 9 with U = 1.42 µm), so that the process requires minimal energy. Then these new cracks may merge with the initial one.

A similar behavior has been observed experimentally in the SCC crack growth in stainless steels exposed to high temperature water [43]. The process is commonly cited as *Initiation Dominant Growth* [44] and has been related to the nucleation of adjacent cracks in the vicinity of the primary crack tip (stress concentration zone), so the main defect becomes larger by coalescence with the recently nucleated cracks rather than by the direct growth of the primary crack.

<figure><img src="image_18.png" /><figcaption>Fig. 9: Crack propagation for polycrystalline structure containing 50 grains. The phase field is plotted corresponding to several loading steps</figcaption></figure>
# 6 Conclusion

In this contribution, a new multiphysics model is proposed to describe the transition between IGSCC and TGSCC at the microstructural level. The anodic dissolution induced SCC is modeled by using the anisotropic phase field method based on KKS model [26] for the electrochemical process and the idea in the work [23] for coupling with mechanical processes. The anisotropic effects of oriented crystals are taken into account for both elasticity and fracture energies. Moreover, the behavior of grain boundaries are described by the smeared cohesive model in the phase field framework to avoid defining of new material phases for this region. More specifically, we propose an approximation to account the different SCC phenomena between grain boundary and grain interior by using the diffusion coefficient and fracture energy of phase concentration depending on a smeared grain boundary field. The new method is able to simulate the stress corrosion crack propagation accounting for polycrystaline effects in a very efficient manner.

The applications to Zirconium alloy demonstrate a very good qualitative agreement between numerical predictions and experimental observations for both situation: IGSCC and TGSCC when SCC system is subjected to different loading conditions. The transition from IG to TG fracture due to the distribution of crystallographic and fracture energy mismatch between grain boundaries and grain interior is also captured. The proposed method opens a new efficient way to model the interaction and competition between IGSCC and TGSCC.

**Acknowledgment:** 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 Institut de la Corrosion, AREVA NP, Andra, MISTRAS Group and Pierre Combrade for their participation in the fruitful discussions during this work.

# References

[1] R.C. Newman. 2001 WR Whitney Award Lecture: Understanding the corrosion of stainless steel. *Corrosion*, 57(12):1030-1041, 2001.

[2] J.S. Chen, M. Salmeron, and T.M. Devine. Intergranular vs transgranular stress corrosion cracking of Cu-30Au. *Scripta metallurgica et materialia*, 26(5):739-742, 1992.

[3] J.S. Chen, M. Salmeron, and T.M. Devine. Intergranular and transgranular
stress corrosion cracking of Cu-30Au. *Corrosion science*, 34(12):2071-2097, 1993.

[4] V.Y. Gertsman and S.M. Bruemmer. Study of grain boundary character along intergranular stress corrosion crack paths in austenitic alloys. *Acta Materialia*, 49(9):1589-1598, 2001.

[5] M.A. Arafin and J.A. Szpunar. A new understanding of intergranular stress corrosion cracking resistance of pipeline steel through grain boundary character and crystallographic texture studies. *Corrosion Science*, 51(1):119-128, 2009.

[6] A. Musienko and G. Cailletaud. Simulation of inter-and transgranular crack propagation in polycrystalline aggregates due to stress corrosion cracking. *Acta materialia*, 57(13):3840-3855, 2009.

[7] D.S. Dugdale. Yielding of steel sheets containing slits. *Journal of the Mechanics and Physics of Solids*, 8:100-104, 1960.

[8] G.I. Barenblatt. The formation of equilibrium cracks during brittle fracture. General ideas and hypotheses. axially-symmetric cracks. *Journal of Applied Mathematics and Mechanics*, 23(3):622-636, 1959.

[9] I. Simonovski and L. Cizelj. Towards modeling intergranular stress corrosion cracks on grain size scales. *Nuclear Engineering and Design*, 246:107-114, 2012.

[10] J. Zhang, Z.H. Chen, and C.F. Dong. Simulating intergranular stress corrosion cracking in AZ31 using three-dimensional cohesive elements for grain structure. *Journal of Materials Engineering and Performance*, 24(12):4908-4918, 2015.

[11] J.J. Rimoli and M. Ortiz. A three-dimensional multiscale model of intergranular hydrogen-assisted cracking. *Philosophical Magazine*, 90(21):2939-2963, 2010.

[12] M. Kamaya. Influence of grain boundaries on short crack growth behaviour of igscc. *Fatigue & Fracture of Engineering Materials & Structures*, 27(6):513-521, 2004.

[13] A.P. Jivkov and T.J. Marrow. Rates of intergranular environment assisted cracking in three-dimensional model microstructures. *Theoretical and Applied Fracture Mechanics*, 48(3):187-202, 2007.

[14] M.A. Arafin and J.A. Szpunar. A novel microstructure-grain boundary character based integrated modeling approach of intergranular stress corrosion crack propagation in polycrystalline materials. *Computational Materials Science*, 47(4):890-900, 2010.

[15] M. Kamaya and M. Itakura. Simulation for intergranular stress corrosion cracking based on a three-dimensional polycrystalline model. *Engineering Fracture Mechanics*, 76(3):386-401, 2009.
[16] M. Kamaya, K. Fukuya, and T. Kitamura. Angular distribution of slip steps by three-dimensional polycrystalline model for stainless steel. *Journal of nuclear science and technology*, 46(9):901-906, 2009.

[17] A. Stoll and A.J. Wilkinson. Use of a dislocation-based boundary element model to extract crack growth rates from depth distributions of intergranular stress corrosion cracks. *Acta Materialia*, 60(13):5101-5108, 2012.

[18] D. De Meo, C. Diyaroglu, N. Zhu, E. Oterkus, and M.A. Siddiq. Modelling of stress-corrosion cracking by using peridynamics. *International Journal of Hydrogen Energy*, 41(15):6593-6609, 2016.

[19] Q. Chen, N. Ma, K. Wu, and Y. Wang. Quantitative phase field modeling of diffusion-controlled precipitate growth and dissolution in Ti-Al-V. *Scripta Materialia*, 50(4):471-476, 2004.

[20] I. Kovačević and B. Šarler. Solution of a phase-field model for dissolution of primary particles in binary aluminum alloys by an r-adaptive mesh-free method. *Materials Science and Engineering: A*, 413:423-428, 2005.

[21] Z. Xu and P. Meakin. Phase-field modeling of solute precipitation and dissolution. *The Journal of chemical physics*, 129(1):014705, 2008.

[22] W. Mai, S. Soghrati, and R.G. Buchheit. A phase field model for simulating the pitting corrosion. *Corrosion Science*, 110:157-166, 2016.

[23] T-T. Nguyen, J. Bolivar, Y. Shi, J. Réthoré, A. King, M. Fregonese, J. Adrien, J-Y. Buffiere, and M-C. Baietto. A phase field method for modeling anodic dissolution induced stress corrosion crack propagation. *Submitted*, 2017.

[24] T-T. Nguyen, J. Réthoré, J. Yvonnet, and M-C. Baietto. Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials. *Computational Mechanics*, 60:289-314, 2017.

[25] T-T. Nguyen, J. Yvonnet, Q.-Z. Zhu, M. Bornert, and C. Chateau. 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*, 312:567-595, 2016.

[26] S.G. Kim, W.T. Kim, and T. Suzuki. Phase-field model for binary alloys. *Physical review e*, 60(6):7186, 1999.

[27] C.V. Verhoosel and R. De Borst. A phase-field model for cohesive fracture. *International Journal for numerical methods in Engineering*, 96(1):43-62, 2013.

[28] T-T. Nguyen, J. Réthoré, and M-C. Baietto. Phase field modelling of anisotropic crack propagation. *European Journal of Mechanics-A/Solids*, 65:279-288, 2017.

[29] S.Y. Hu, J. Murray, H. Weiland, Z.K. Liu, and L.Q. Chen. Thermodynamic
description and growth kinetics of stoichiometric precipitates in the phase-field approach. *Calphad*, 31(2):303-312, 2007.

[30] A.A. Abubakar, S.S. Akhtar, and A.F.M. Arif. Phase field modeling of V2O5 hot corrosion kinetics in thermal barrier coatings. *Computational Materials Science*, 99:105-116, 2015.

[31] C. Miehe, M. Hofacker, and F. Welschinger. 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, 2010.

[32] T-T. Nguyen, J. Yvonnet, Q-Z. Zhu, M. Bornert, and C. Chateau. 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, 2015.

[33] B. Cox. Environmentally-induced cracking of zirconium alloys—a review. *Journal of Nuclear Materials*, 170(1):1-23, 1990.

[34] J. Prong, T. Jaszay, A. Caprani, and J.P. Frayret. Anodic behaviour of zirconium and its alloys in fluorinated nitric media dissolution-passivation model. *Journal of applied electrochemistry*, 25(11):1031-1037, 1995.

[35] D.R. Baker, A. Conte, C. Freda, and L. Ottolini. The effect of halogens on zr diffusion and zircon dissolution in hydrous metaluminous granitic melts. *Contributions to Mineralogy and Petrology*, 142(6):666-678, 2002.

[36] F.J. Mompean, J. Perrone, M. Illemassène, et al. *Chemical thermodynamics of zirconium*, volume 8. Gulf Professional Publishing, 2005.

[37] X-P Xu and A Needleman. Void nucleation by inclusion debonding in a crystal matrix. *Modelling and Simulation in Materials Science and Engineering*, 1(2):111, 1993.

[38] T. Luther and C. Könke. Polycrystal models for the analysis of intergranular crack growth in metallic materials. *Engineering Fracture Mechanics*, 76(15):2332-2343, 2009.

[39] E.S. Fisher and C.J. Renken. Single-crystal elastic moduli and the hcp? bcc transformation in ti, zr, and hf. *Physical Review*, 135(2A):A482, 1964.

[40] R. Quey. Neper: software package for polycrystal generation and meshing, 2016.

[41] M. Fregonese. Mécanismes de corrosion sous contrainte par l'iode dans le Zirconium et le Zircaloy-4 transposition aux conditions d'interaction pastillergaine dans les réacteursa l'eau pressurisée. *PhD thesisEcole des Mines de Saint-Etienne*, 1997.

[42] S.B. Farina, G.S. Duffo, and J.R. Galvele. Stress corrosion cracking of zirconium and zircaloy-4 in halide aqueous solutions. *Corrosion Science*, 45(11):2497-2512, 2003.
[43] M. Kamaya and T. Haruna. Crack initiation model for sensitized 304 stainless steel in high temperature water. *Corrosion science*, 48(9):2442-2456, 2006.

[44] M. Kamaya and T. Haruna. Influence of local stress on initiation behavior of stress corrosion cracking for sensitized 304 stainless steel. *Corrosion Science*, 49(8):3303-3324, 2007.