Ethel Djeumen
email: etheldjeumnen@gmail.com

Gergely Molnár
email: gergely.molnar@insa-lyon.fr

Nicolas Tardif

Michel Coret

Jean Desquines

Tatiana Taurines

Marie-Christine Baietto

Modeling diffusive phase transformation and fracture in viscoplastic materials

Keywords: multiphase-field model, multiphysical couplings, viscoplasticity, diffusion, phase transformation, fracture

One of the most feared failure modes in engineering is fracture as it can occur suddenly and potentially lead to catastrophic consequences. The combined action of mechanical stresses and phase transformations can significantly affect the response of structures. For example, in case of a loss of coolant accident in a pressurized water nuclear reactor, an interaction between oxygen diffusion, oxidation, creep and cracking phenomena is observed [START_REF] Campello | Identification of the steady-state creep behavior of Zircaloy-4 claddings under simulated Loss-Of-Coolant Accident conditions based on a coupled experimental/numerical approach[END_REF][START_REF] Chosson | Creep flow and fracture behavior of the oxygen-enriched alpha phase in zirconium alloys[END_REF][START_REF] Desquines | High temperature oxidation and room temperature axial strength of pre-oxidized zircaloy-4 cladding after a simulated loca[END_REF][START_REF] Dominguez | Effect of pre-oxide and hydrogen on creep of Zircaloy-4 at 1123K[END_REF][START_REF] Thieurmel | Contribution to the understanding of brittle fracture conditions of zirconium alloy fuel cladding tubes during LOCA transient[END_REF]. Modeling of the mechanical behavior coupled with the phase transformation of a multilayer material is therefore an important issue. This work aims at developing a numerical platform coupling nonlinear mechanical behavior to phase transformation and fracture.

The phase-field approach is used to represent fracture and phase transformation. This approach based on energy minimization lends itself well to numerical implementation and seems to be the simplest tool for modeling multi-physics problems in a single framework [START_REF] Wu | Chapter One -Phasefield modeling of fracture[END_REF]. The phase-field approach is well suited for phase transformation and easily allows the consideration of the diffusion of chemical species.

Regarding fracture, the phase-field approach is suitable for brittle fracture. However, it is well known that ductile and brittle fracture involve different mechanisms. In this work, the fracture phase-field approach is chosen to represent both brittle and ductile fracture responses. The description of the material aspects related to the ductile fracture is not treated and it is assumed that the phases are homogeneous.

The phase-field approach is employed in order to deal with singularities such as fracture and phase transformation (Ambati et al., 2015b;[START_REF] Singer-Loginova | The phase field technique for modeling multiphase materials[END_REF]. It has become increasingly popular in recent years due to its simplicity and flexibility. The main idea behind the phase-field approach is to smooth singularities in order to describe them with auxiliary scalar fields.

In the phase-field approach for phase transformation, the sharp interface is replaced by a diffuse interface of non-zero thickness (δ). With this regularization, the phases are described by a continuous order parameter (φ) governing the local interface kinetics. This phase-field approach has a decisive advantage compared to sharp interface models as explicit interface tracking is not required. Sharp interface models are not suitable when dealing with complex morphological evolutions because the interface tracking at each time step becomes computationally quite unfeasible.

The phase-field approach for phase transformation can be used to model both microstructures and mesoscale complex morphological evolutions such as grain-boundary evolution [START_REF] Bhattacharya | A phase-field investigation of the effect of grain-boundary diffusion on austenite to ferrite transformation[END_REF][START_REF] Moelans | Quantitative analysis of grain boundary properties in a generalized phase field model for grain growth in anisotropic systems[END_REF][START_REF] Suwa | Three-dimensional phase field simulation of the effect of anisotropy in grain-boundary mobility on growth kinetics and morphology of grain structure[END_REF] and solidification [START_REF] Zaeem | Modeling dendritic solidification of al-3%cu using cellular automaton and phase-field methods[END_REF][START_REF] Labergère | Thermomechanical model for solidification and cooling simulation of Ni-based superalloy components[END_REF].

It can also be coupled with other physical ingredients such as elasticity [START_REF] Dreyer | A study of the coarsening in tin/lead solders[END_REF][START_REF] Düsing | A coupled phase field/diffusional/mechanical framework for simulation of upper and lower bainitic transformation[END_REF]Wang et al., 1993b), plasticity (Ammar et al., 2009a;[START_REF] Zhao | Elasto-plastic phase field modelling of oxidation of zirconium alloys[END_REF] and viscoplasticity [START_REF] Cottura | A phase field model incorporating strain gradient viscoplasticity: Application to rafting in ni-base superalloys[END_REF][START_REF] De Rancourt | Homogenization of viscoplastic constitutive laws within a phase field approach[END_REF].

Similarly, in the fracture phase-field approach, the sharp crack is replaced by a regularized diffused damage. In order to mathematically tackle the problem, a length scale (l c ) is introduced, which controls the magnitude of the damage penetration into a solid phase. A diffuse damage variable (d) connects broken and unbroken regions. The fracture phase-field approach is based on the original description of Griffith for homogeneous brittle materials. [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF][START_REF] Griffith | The theory of rupture[END_REF] suggested that a pre-existing crack can only propagate if the released potential energy per unit surface is greater than a critical value. The critical energy release rate is postulated to be a material constant that indicates its fracture toughness. The main advantage of the frac-ture phase-field approach is its simplicity. Compared to discrete methods such as the cohesive zone model [START_REF] Falk | A critical evaluation of cohesive zone models of dynamic fractur[END_REF][START_REF] Zhou | Dynamic crack propagation with cohesive elements: a methodology to address mesh dependency[END_REF], the so-called strong discontinuity technique [START_REF] Huespe | STRONG DISCONTINUITY APPROACH IN DYNAMIC FRACTURE SIMULATIONS[END_REF] or the extended finite element method [START_REF] Elguedj | Appropriate extended functions for X-FEM simulation of plastic fracture mechanics[END_REF][START_REF] Gravouil | Stabilized global-local X-FEM for 3D non-planar frictional crack using relevant meshes[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Réthoré | An energy-conserving scheme for dynamic crack growth using the eXtended finite element method[END_REF], the diffused damage models (such as the phase-field or the thick level set approach [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF]) provide a suitable framework for solving the problem with classical numerical methods without requiring special criteria or numerical tools. Ambati et al. (2015b) gives a detailed overview of the existing brittle fracture phase-field models. The fracture phase-field provides a simulation framework to adequately describe crack initiation, propagation, curvilinear crack paths, branching or coalescence [START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF]. The effect of plasticity (Ambati et al., 2015a;[START_REF] Duda | A phase-field/gradient damage model for brittle fracture in elastic-plastic solids[END_REF]Molnár et al., 2020a;[START_REF] Zhang | Ductile damage modelling with locking-free regularised gtn model[END_REF], viscoplasticity [START_REF] Shanthraj | A phase field model for damage in elasto-viscoplastic materials[END_REF], hyperelasticity [START_REF] Mandal | A length scale insensitive anisotropic phase field fracture model for hyperelastic composites[END_REF][START_REF] Peng | Phase field simulation for fracture behavior of hyperelastic material at large deformation based on edge-based smoothed finite element method[END_REF][START_REF] Yin | Hyperelastic phase-field fracture mechanics modeling of the toughening induced by bouligand structures in natural materials[END_REF], fatigue [START_REF] Lo | A phase-field model for fatigue crack growth[END_REF][START_REF] Mesgarnejad | Phase-field models for fatigue crack growth[END_REF], and chemical diffusion (Miehe et al., 2016b;[START_REF] Nguyen | A phase field method for modeling anodic dissolution induced stress corrosion crack propagation[END_REF][START_REF] Schneider | Phase-Field Modeling of Diffusion Coupled Crack Propagation Processes[END_REF] can also be included in a thermodynamically consistent manner.

The diffusive phase transformations that heterogeneous materials undergo can significantly affect their fracture behavior. However, a general approach to study the fracture behavior of these non-static heterogeneous materials is still missing. [START_REF] Schneider | Phase-field modeling of crack propagation in multiphase systems[END_REF] proposed a multiphase-field model that combines the fracture phase-field and the phase-field for the phase transformation process, but the application of this model is limited to diffusionless transformation such as martensitic transformation as the diffusion process is not modeled.

In this work, we combine the fracture phase-field approach with the phase-field approach for diffusive phase transformation processes to study the fracture behavior of solid-state materials undergoing mesoscale diffusive phase transformations. Compared to previous works of [START_REF] Schneider | Phase-field modeling of crack propagation in multiphase systems[END_REF], the proposed model is more generalised by including in the formulation the chemical diffusion which is the driving force of the diffusive phase transformation. Four physical phenomena are therefore considered including fracture, diffusion, phase transformation and viscoplasticity.The simultaneous consideration of these four physics is important for the study of the combined action of a mechanical stress and an aggressive environment encountered in several applications in industry. For example, the phenomenon of loss of coolant accident in a nuclear reactor [START_REF] Campello | Identification of the steady-state creep behavior of Zircaloy-4 claddings under simulated Loss-Of-Coolant Accident conditions based on a coupled experimental/numerical approach[END_REF][START_REF] Chosson | Creep flow and fracture behavior of the oxygen-enriched alpha phase in zirconium alloys[END_REF][START_REF] Desquines | High temperature oxidation and room temperature axial strength of pre-oxidized zircaloy-4 cladding after a simulated loca[END_REF][START_REF] Dominguez | Effect of pre-oxide and hydrogen on creep of Zircaloy-4 at 1123K[END_REF][START_REF] Thieurmel | Contribution to the understanding of brittle fracture conditions of zirconium alloy fuel cladding tubes during LOCA transient[END_REF] or the well-known phenomenon of stress corrosion cracking [START_REF] Sieradzki | Stress-corrosion cracking[END_REF]. The formulation of the model is carried out in an appropriate framework that easily allows its extension with new physical ingredients or additional couplings.

The present paper is structured as follows. Section 2 presents the consistent framework combining the various physical phenomena being modeled (fracture, diffusion, phase transformation and viscoplasticity) into a single multiphase-field model. The theoretical background and modeled couplings are specified. In Section 3, through several benchmark examples of increasing complexity, the proposed model is used to investigate the effect of viscoplastic activity and diffuse phase transformation on the fracture behavior. The impact on fracture behavior is assessed by examining the critical loading, energy density until failure and crack topology. Finally, Section 4 concludes this work.

Methods

The phase-field approach for fracture and phase transformation are based on the diffuse representation of the discontinuity in the volume as illustrated in Fig. 1. For fracture, we define a diffuse damage variable with viscoplastic activity, the mass transport or diffusion, the fracture and the phase transformation problems.

In general, the total energy of the multi-physics problem composed of a mechanical (L u ), chemical (L ch ) and fracture (L d ) part can be written as :

(1)

L (∇u, φ, c, d) = L u (∇u, φ, c, d) + L ch (∇u, φ, c, d) + L d (∇u, φ, c, d)
where u, φ, c, d are respectively the displacement vector, the order parameter, the concentration of the diffusing chemical species and the damage variable. Eq. ( 1) is written in the general case considering all possible interaction between the modelled physics. In the following, only the couplings summarized in Fig. 2 are considered. The damage reduces the stiffness of the material, while the elastic (undamaged) strain energy contributes to the crack surface opening. The phase transformation is driven by the free energy which is composed of the stored elastic strain energy and chemical energy. In the interface where the two phases coexist, a mixing or interpolation method is used to include the energy contributions of both phases.

Several couplings are taken into account but the list is not exhaustive. However, the formulation is proposed in a framework that allows the introduction of new couplings in a relatively easy way. The effect of the diffusing species on the mechanical and fracture properties (Le [START_REF] Saux | An approach to study oxidation-induced stresses in Zr alloys oxidized at high temperature[END_REF][START_REF] Nguyen | A phase field method for modeling stress corrosion crack propagation in a nickel base alloy[END_REF] and the coupling between diffusion and viscoplasticity [START_REF] Le Claire | A mathematical analysis of diffusion in dislocations. I. Application to concentration 'tails[END_REF][START_REF] Le Claire | A mathematical analysis of diffusion in dislocations. II. Influence at low densities on measured diffusion coefficients[END_REF][START_REF] Le Claire | A mathematical analysis of diffusion in dislocations. III. Diffusion in a dislocation array with diffusion zone overlap[END_REF] are not taken into account. We also neglect the explicit coupling between phase transformation and fracture, thus their constitutive laws can be formulated separately. Griffith's theory [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF] is employed for the fracture phase-field approach, while the Allen-Cahn theory [START_REF] Allen | A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening 27[END_REF] is used for phase transformation phase-field.

In fracture phase-field, the mechanical energy is used to control crack surface opening. The consideration of the contribution of plastic or viscoplastic dissipated energy to fracture is controversial in the literature (Molnár et al., 2020a). Some authors do not include the effect of plastic (Ambati et al., 2015a) or viscoplastic [START_REF] Shanthraj | A phase field model for damage in elasto-viscoplastic materials[END_REF] energy on the crack surface opening while others do [START_REF] Badnava | A Phase Field Model for Rate-Dependent Ductile Fracture[END_REF][START_REF] Fang | Phase field fracture in elasto-plastic solids: Abaqus implementation and case studies[END_REF]. This is a question of scale. Usually, it is assumed that very small-scale plasticity can cause void formation, leading to macroscopic cracks, for example, when the characteristic size of the model approaches the grain size. In this work, only the stored elastic strain energy contributes to the crack surface opening while viscoplasticity remains an independent dissipation.

In phase transformation phase-field, the so-called free energy, which contains both stored elastic strain energy and chemical energy, is used to govern the phase transformation. The motivation for adding a mechanical effect in the phase transformation problem comes from the idea that a strain accommodation leads to a development of certain optimal microstructures [START_REF] Khachaturyan | Theory of macroscopic periodicity ofr a phase transition in the solid state 29[END_REF].

In the case of an unstable crack, the monolithic solution of the fully coupled displacement and fracture problem becomes unstable. Therefore, the displacement problem is weakly coupled to the evolution of fracture topology using the staggered scheme proposed by Miehe et al. (2010a). As shown in Fig. 3, at each increment the displacement problem is solved on the basis of the last converged fracture topology and vice versa.

Displacement problem

The displacement problem is solved under the following assumptions:

• damage only reduces the tensile stiffness, the compressive stiffness is not affected,

• the newly formed phase inherits the total viscoplastic strain of the parent phase,

• there is no elastic limit in viscoplasticity, • the problem is restricted to small deformations.

The energy of the displacement problem is formulated as follows:

(2)

L u (∇u, φ, d) = E ext -E el (∇u, φ, d) -Ω u (∇u, φ, d),
where u is the displacement field. E ext , E el and Ω u are respectively the external work, the stored elastic strain energy and the viscoplastic strain energy.

The external work done by the body (f v ) and boundary (t) forces illustrated in Fig. 3 is given by:

(3)

E ext = � V (f v • u)dV + � ∂V (t • u)dS .
The stored elastic energy is split into tensile (ψ + 0 ) and compression (ψ - 0 ) components [START_REF] Bernard | Damage growth modeling using the Thick Level Set (TLS) approach: Efficient discretization for quasi-static loadings[END_REF]:

(4) E el (∇u, φ, d) = � V ψ el (∇u, φ, d)dV (5) ψ el (∇u, φ, d) = g (d) ψ + 0 � ε el (∇u, φ) � + ψ - 0 � ε el (∇u, φ) � .
The decomposition of the elastic strain energy allows a degradation only in tension. g(d) is the degradation function:

(6) g(d) = (1 -d) 2 + p,
where p is a small number responsible for numerical stability.

The degraded Cauchy stress tensor is then given by :

(7) σ = g(d) ∂ψ + 0 ∂ε el + ∂ψ - 0 ∂ε el ,
where ε el is the elastic strain tensor. The degradation function directly acts on the tensile component of the stress tensor.

The total strain field (ε) is divided into elastic (ε el ), viscoplastic (ε vp ) and eigenstrain (ε * ) components:

(8) ε = ε el + ε vp + ε * .
In this work, ε * represents the volumetric change due to the phase transformation. Assuming small deformations, the strain tensor is defined by the symmetric part of the displacement gradient:

(9) ε = ∇ s u.
In the interface, a mixing scheme is applied to determine the effective behavior. Several mixing schemes exist in the literature (A. Reuss, 1929;Ammar et al., 2009a;[START_REF] Durga | Evaluation of interfacial excess contributions in different phasefield models for elastically inhomogeneous systems[END_REF][START_REF] Mosler | A novel homogenization method for phase field approaches based on partial rank-one relaxation[END_REF][START_REF] Voigt | Thermodynamically-consistent phase-field models for solidification[END_REF][START_REF] Wang | Shape instability during precipitate growth in coherent solids[END_REF]). Voigt's, Reuss's and Khachaturyan's schemes are the most used. For the two former ones, the equivalent elastic strain is calculated from the equivalent compliance tensor. In a damaged region (d = 1), this compliance tensor is no longer well defined and tends towards infinity which can therefore induce numerical convergence issues.

Herein, we have retained Khachaturyan's scheme, which does not need this compliance tensor to compute equivalent elastic strain. Therefore, the total strain is assumed to be identical in both coexisting phases while the mixture rules are applied to viscoplastic strain components, eigenstrains and to Lamé parameters (λ, µ):

(10) ε = ε a = ε b , (11) ε vp (∇u, φ) = h(φ)ε vp a (∇u) + h(φ)ε vp b (∇u), (12) ε * (φ) = h(φ)ε * a + h(φ)ε * b , (13.1) λ(φ) = h(φ)λ a + h(φ)λ b (13.2) µ(φ) = h(φ)µ a + h(φ)µ b .
The subscripts a and b indicate the two coexisting phases. The standard choice of the homogenization function h(φ) = 1 -h(φ) is (Wang et al., 1993a):

(14) h(φ) = φ 2 (3 -2φ).
This mixing scheme, referred to as interpolated scheme by [START_REF] Ammar | Modelling inheritance of plastic deformation during migration of phase boundaries using a phase field method[END_REF], assumes that the newly formed phase inherits the total viscoplastic strain of the parent phase in contrast to the homogenization scheme of [START_REF] Voigt | Thermodynamically-consistent phase-field models for solidification[END_REF] and A. Reuss (1929).

The decomposition of elastic energy is given by [START_REF] Bernard | Damage growth modeling using the Thick Level Set (TLS) approach: Efficient discretization for quasi-static loadings[END_REF]:

(15) ψ el � ε el , φ, d � = µ(φ) 3 � i=1 � �ε i � 2 -+ g(d)�ε i � 2 + � + λ(φ) 2 � g(d)�tr(ε)� 2 + + �tr(ε)� 2 - � .
For the practical implementation of this decomposition, readers are referred to previous work by Molnár et al. (2020a).

As for the elastic strain energy in Eq. ( 4), the viscoplastic strain energy is expressed as:

(16) Ω u (∇u, φ, d) = � V � t � ψvp (∇u, φ, d) � dtdV, (17) ψvp (∇u, φ, d) = σ(∇u, d) : εvp (∇u, φ).
The viscoplastic dissipation is calculated from the degraded stress tensor (see Eq. ( 7)). Since viscoplastic strain does not induce a change in volume, the viscoplastic energy can be simplified according to:

(18) ψvp (∇u, φ) = σ vm (d) εvp vm (∇u, φ),
where σ vm and εvp vm are respectively the von Mises stress and equivalent viscoplastic strain rate. An elastic limit is not integrated into this work. The equivalent viscoplastic strain rate in each of the phases a, b is described by Norton's law:

(19) εvp vm,i = A i σ n i vm ,
where A i and n i are Norton's material parameters. Note that in Eq. ( 19) the equivalent viscoplastic strain rate in each of the phases is calculated from the equivalent von Mises stress. The rule of mixture used in this work does not give access to the stress tensor of each phase.

By taking the variation of Eq. ( 2), the corresponding strong form of the displacement problem can be obtained:

(20.1)

δL u = 0 ∀δu → ∇ • σ -f v = 0 in V, (20.2) σ • n = t on Γ N , (20.3) u = u on Γ u .

Fracture problem

The fracture problem is formulated under the following main assumptions:

• the crack surface opening is driven by the tensile elastic (undamaged) strain energy,

• the length scale (l c see Fig. 1(b)) is constant in the whole domain,

• the critical energy release rate per unit surface (g c ) is treated as a material parameter,

• the damage is irreversible.

The energy functional to solve the crack topology is formulated as follows:

(21)

L d (∇u, φ, d) = Ω d (φ, d, ∇d) + T (φ, d) -E el (∇u, φ, d) ,
where Ω d is the fracture energy and T is the threshold function. The regularized phase-field Griffith fracture energy can be expressed as [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF]:

(22) Ω d (φ, d) = � Γ g c (φ)dΓ ≈ � V g c (φ)γ(d)dV, (23) γ(d, ∇d) = 1 2l c d 2 + l c 2 |∇d| 2 ,
where γ is the fracture energy density and g c the equivalent surface energy needed to create a unit fracture surface and l c is the length scale parameter which controls the scale of damage penetration. For details on this approximation, readers are invited to consult previous works of Miehe et al. (2010a) or [START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF].

The work of Molnár et al. (2020a) shows a unique correlation between l c and the material's strength, based on the maximum load and the crack topology, which allows us to give the length scale parameter physical meaning. However initially, l c is used as a numerical aid. Therefore, the mixture rule is applied only to g c but not to l c :

(24) g c (φ) = h(φ)g ca + h(φ)g cb
The threshold function in Eq. ( 21) is chosen as (Miehe et al., 2016b;Molnár et al., 2020a):

(25) T (d) = � V � (1 -g (d)) ψ c (φ) � dV, with ψ c (φ) = g c (φ) 2l c .
In practice, to solve the fracture topology, the elastic strain energy in Eq. ( 21) is replaced by a so-called history field H (Miehe et al., 2010a) to prevent crack healing:

H 0 = 0, (26.1)
H n+1 = max            ψ + 0 (∇u, φ) -ψ c (φ) H n (26.2)
This history field (H) helps to satisfy the Karush-Kuhn-Tucker conditions [START_REF] Singh | A fracture-controlled path-following technique for phase-field modeling of brittle fracture[END_REF]. In addition to the introduction of the history field, the damage irreversibility is explicitly enforced in the displacement problem ( ḋ ≥ 0) as there is no mathematical proof that the history field alone is sufficient to ensure damage irreversibility. Therefore, the energy functional to solve the crack topology can be rewritten as:

(27) L d (∇u, φ, d) = Ω d - � V � g (d) H (∇u, φ) -ψ c (φ) � dV.
As for the displacement problem the strong form of fracture topology can be deduced as:

(28.1)

δL d = 0 ∀δd → g c l c (d -l 2 c Δd) = 2(1 -d)H in V, (28.2) ∇d • n = 0 on Γ,

Coupled diffusion-phase transformation problem

The coupled problem of diffusion and phase transformation is formulated by assuming that the diffusion of chemical species and the effect of the elastic energy on the phase transformation are isotropic. The starting point for deriving the evolution laws of the diffusion and phase transformation problem is the total free energy functional defined by :

(29) L ch (∇u, c, φ, d) = � V ψ(∇u, c, φ, d)dV,
where ψ is the Helmholtz free energy density. ψ is the amount of energy stored in the system per unit volume. Thus, ψ is made up of the stored elastic strain energy (ψ el ), the chemical energy (ψ ch ) and interface energy (ψ inter f ace ) contributions [START_REF] De Rancourt | Homogenization of viscoplastic constitutive laws within a phase field approach[END_REF]:

(30) ψ(∇u, c, φ, d) = ψ el (∇u, φ, d) + ψ ch (c, φ) + ψ inter f ace (φ, ∇φ).
The viscoplastic energy in Eq. ( 16) as well as the fracture energy in Eq. ( 22) are dissipated energies and therefore do not form part of the free energy in Eq. ( 30). Although the effects of fracture and viscoplasticity do not appear explicitly in Eq. ( 30), they have an indirect effect through the stored elastic strain energy.

Indeed, the energy dissipated by fracture and viscoplasticity comes from the stored elastic energy. By dissipating the stored elastic energy that is part of the free energy, fracture and viscoplasticity contribute indirectly to the coupled diffusion-phase transformation problem.

The equivalent chemical free energy density (ψ ch ) is obtained by a mixing law on the homogeneous chemical free energy densities of both phases (ψ ch i ) [START_REF] Kim | Interfacial compositions of solid and liquid in a phase-field model with finite interface thickness for isothermal solidification in binary alloys[END_REF]:

(31) ψ ch (c, φ) = h(φ)ψ ch a (c) + h(φ)ψ ch b .
For the sake of simplicity, the chemical free energy density ψ ch i of each phase i = {a, b} is chosen as a square function of the concentration c [START_REF] Echebarria | Quantitative phase-field model of alloy solidification[END_REF]:

(32) ψ ch i (c) = 1 2 k i (c -c i ) 2 ,
where k i , and c i are respectively the curvature of the free energy density and the equilibrium concentration of each phase.

The interface contribution ψ inter f ace is given by [START_REF] Kim | Interfacial compositions of solid and liquid in a phase-field model with finite interface thickness for isothermal solidification in binary alloys[END_REF]:

(33) ψ inter f ace (φ, ∇φ) = Wg ch (φ) + α 2 ∇φ • ∇φ,
where W is the height of the energy barrier accounting for the free energy penalty of the interface, α is the coefficient of the gradient term of the interface energy in Eq. ( 33) and g ch (φ) is the double well potential [START_REF] Raabe | Continuum scale simulation of engineering materials:Fundamentals microstructures process applications[END_REF]: The diffusion and phase transformation problem are governed by the corresponding driving force defined as the functional derivative of the total free energy. We obtain the Cahn-Hilliard equation by deriving the total free energy in Eq. ( 29) with respect to c:

(34) g ch (φ) = φ 2 (1 -φ) 2 .
(35) ∂c(x, t) ∂t = ∇ • � -L(φ)∇ • � δF δc �� = ∇ • � -L(φ)∇ � ∂ψ ∂c �� .
The concentration (c) is a conserved field for which the spatial average is constant over time. Boundary conditions are used to take into account the exchanges on the edges. In the same way, the Allen-Cahn equation is obtained by deriving the total free energy in Eq. ( 29) with respect to φ. However, the order parameter (φ) is not a conserved field:

(36) ∂φ(x, t) ∂t = M δF δφ = M ∂ψ ∂φ -∇ • ∂ψ ∂∇φ ,
where M is the interface mobility parameter which must be as large as possible in the case of diffusioncontrolled phase transformation [START_REF] De Rancourt | Homogenization of viscoplastic constitutive laws within a phase field approach[END_REF], L(φ) is the homogenized diffusion coefficient (Onsager coefficient) given by (Ammar et al., 2009a): where D i is the diffusion coefficient of the phases. L(φ) is chosen so as to recover Fick's law in a homogeneous case. For a homogeneous material (φ = 0), Eq. ( 35) is equivalent to Fick's law:

(37) L(φ) = h(φ) D a k a + h(φ) D b k b ,
∂c ∂t = ∇ • � -L(0)∇ � ∂ψ ∂c (0) �� (38.1) = ∇ • � - D b k b ∇ (k b (c -c b )) � (38.2) = -D b Δc. (38.3)
The interface energy (ω) and the interface thickness (δ) (see Fig. 1) can be deduced from α and W by (Ammar et al., 2009a;[START_REF] Kim | Interfacial compositions of solid and liquid in a phase-field model with finite interface thickness for isothermal solidification in binary alloys[END_REF]:

(39) ω = √ αW 3 √ 2 , (40.1) δ = θ � 2α W , (40.2) θ = ln � 1 -η η � .
Eq. (40.1) is established assuming that the interface region ranges from φ = κ to φ = 1κ. In this work, κ = 0.05 was chosen.

To be consistent, the model must verify the laws of thermodynamics. A variation in g c during a fracture process can compromise the energy balance and thus violate the first principle of thermodynamics.

This problem is observed in phenomenological phase-field fatigue models [START_REF] Carrara | A framework to model the fatigue behavior of brittle materials based on a variational phase-field approach[END_REF] and ratedependent fracture [START_REF] Yin | Formulation and implementation of strain rate-dependent fracture toughness in context of the phase-field method[END_REF], where the authors redefined the dissipation rate. Thus the overall fracture surface is no more the measure of dissipated energy in fracture. The criterion ḋ > 0 explicitly imposed can also affect the energy balance if g c varies. It is also important to recall that the fracture threshold energy introduced in Eq. ( 21), which has no real physical basis, can also be considered as a violation of thermodynamic consistency. The second law of thermodynamics is discussed in detail in Appendix A.

Staggered time-integration algorithm

The proposed model has been implemented in the commercial Abaqus finite element code trough a UEL routine. The flowchart in Fig. 5 shows the basic incrementation process. At each increment, the displacement problem, diffusion problem and phase transformation problem are solved together on the basis of the last converged fracture topology and vice versa. A Newton-Raphson method can be used to solve the global problem iteratively:

(41) Cracking and viscoplasticity do not occur on the same time scale. Therefore, an adaptive time step control is required to make the simulation computationally feasible. The time step was automatically controlled by the local tensile elastic strain energy increment:

                                                       K uu K uφ 0 0 K φu K φφ K φc 0 0 K cφ K cc 0 0 0 0 K dd                                                                                                               u n+1 φ n+1 c n+1 d n+1                                                        = -                                                        r u n r φ n r c n r d n                                                        ,
(42) dH ≤ ηψ c .
This allows the time step to be increased until crack initiation is reached. In the case of unstable crack propagation, the local tensile elastic strain energy at the crack tip becomes independent of the time step which is then reduced to a minimum value of dt = 10 -9 s. This way, the unstable crack propagation is detected precisely. More details on this automatic time step control based on the energy constraint can be found in (Molnár et al., 2020b).

Model testing strategy

The implemented model is tested by examining quantities that describe the fracture response of a structure. As shown in Fig. 6, the critical load (σ c ) is defined as the maximum amount of tensile stress that the material can withstand before failure. The energy density until failure (U t ) gives an indication of the toughness of a structure. U t represents the area below the stress-strain curve:

Examples and results

This section analyzes the effect of viscoplasticity and diffusive phase transformation on the critical load, energy density until failure and crack topology using examples of increasing complexity. In Section 3.1, it is shown that the proposed model is able to simulate some well-known behaviors for viscoplastic solids, such as ductile failure at low strain rates and brittle failure at high strain rates. The role of viscoplasticity at the crack tip is also highlighted in the case of the simple extension (mode I). From Section 3.2 to Section 3.4, the model is exploited to study the fracture behavior in a heterogeneous solid. The case of a rigid layer growing on a soft matrix is discussed in Section 3.2. The effect of phase transformations and the internal stresses they generate are illustrated. Section 3.3 deals with the fracture behaviour of a heterogeneous solid consisting of an inclusion evolving in a matrix. The effect of the stiffness ratio between the inclusion and the matrix on the crack topology is presented. The case of several randomly distributed inclusions in the matrix is treated in Section 3.4.

Parameters used for the simulations are dimensionless so that the results can be extrapolated to other sets of parameters. The length, time and energy density scales are defined respectively by the size of the system L, the characteristic diffusion time τ = L 2 /max(D i ) and the chemical free energy density curvature

E d = β min(k i ).
Here, β is an arbitrary positive constant. The concentrations are also scaled by the difference between the equilibrium concentrations (c a -c b ) of the two phases [START_REF] De Rancourt | Homogenization of viscoplastic constitutive laws within a phase field approach[END_REF]. Parameters to be scaled, denoted with tildes, are summarized in Table 1.

Table 1

Definitions of the dimensionless parameters.

Physical

Dimensionless

Length: l[m] l = l L Times: t[s] t = t τ
Mechanical properties (n is the Norton's exponent)

Young's modulus:

E[Pa] Ẽ = E E d Norton's constant: A[J -n • m 3n • s -1 ] Ã = A • E n d • τ Fracture surface energy: g c [J • m -2 ] gc = g c E d • L Fracture length scale: l c [m] lc = l c L Chemical proporties concentration: c c = c c a -c b Interface energy: ω[J • m -2 ] ω = ω E d • L interface thickness: δ[m] δ = δ L Interface mobility: M[J -1 • m 3 • s -1 ] M = M • E d • τ Diffusivity: D[m 2 • s -1 ] D = D • τ L 2 Chemical free energy curvature: k[J • m -3 ] k = k E d
All presented results are dimensionless but for the sake of clarity the tilde symbol has been omitted. The simulations are carried out under plane strain conditions.

Competition between fracture and viscoplasticity

The proposed phase-field model does not explicitely predict a coupling between fracture and viscoplasticity. However, these two phenomena are indirectly linked through the elastic strain energy in Eq. ( 4) that can be dissipated either by cracking or by viscoplasticity.

Homogeneous solution

A single 2D homogeneous plane strain element is the simplest case, where the interaction between fracture and viscoplasticity can be understood. A square plate with dimensions of 1 × 1 is subjected to a simple shear strain rate by moving its top side in the x direction while its bottom side is fixed in both x and y directions. Young's modulus was set to E = 210 with ν = 0.3, g c = 2 × 10 -3 and l c = 2 × 10 -4 . The energy increment constraint in Eq. ( 42) was set to η = 1%. Only cracking and viscoplasticity are considered in this example.

Regardless of cracking, the viscoplastic problem can be solved analytically. By combining Eq. ( 7) and Eq. ( 15) which allow to calculate the stress with Norton's law in Eq. ( 18) in the case of a pure shear in a homogeneous solid, we obtain the following differential scalar equation:

(44) τxy = µ � γxy -A( √ 3) n+1 τ n xy � , where µ = E 2(1 + ν)
is the shear modulus. The steady shear stress (τ ∞ xy ) described by τxy = 0 can be deduced from Eq. ( 44) by:

(45) τ ∞ xy = 3 -n + 1 2n • � γxy A � 1 n ,
the corresponding steady stored tensile elastic strain energy (ψ ∞ ) is given by:

(46) ψ ∞ = � τ ∞ xy � 2 4µ
.

It becomes apparent from Eq. ( 45) that the viscoplastic behavior is only dependent on the ratio [γ xy /A] , and not on absolute values. This remains the case when fracture is taken into account, regardless of the loading conditions.

To reach the failure of the element, the steady tensile elastic strain energy (ψ ∞ ) must be greater than the fracture threshold (ψ c ). Therefore, the threshold ratio ([γ xy /A] c ) required to reach this fracture threshold is given by: (47) For n = 1 , the global solution of Eq. ( 44) is trivial and the tensile (ψ el 0,+ ) component of elastic energy can be calculated as:

ψ ∞ ≥ ψ c ⇒ � γxy /A � c = √ 3 � 6µg c l c � n 2 .
(48)

ψ el 0,+ = ψ ∞ � 1 -exp � - µγ xy τ ∞ xy �� 2 .
In Fig. 7, it can be seen that for [γ xy /A] ≤ [γ xy /A] c the purely viscoplastic analytical solution (red continuous line) given by Eq. ( 48) is consistent with the numerical phase-field solution. The steady tensile elastic energy does not exceed the fracture threshold (ψ c ). The viscoplastic dissipation rate is high enough to keep the tensile elastic energy below the fracture threshold. For [γ xy /A] = 1.3[γ xy /A] c , the viscoplastic dissipation rate is no longer high enough. The material cracks as soon as the fracture threshold is reached (ψ el + /ψ c = 1) and the total elastic strain energy is dissipated. These results highlight the role of the fracture threshold (ψ c ) introduced in the fracture energy description. Two time scales can be observed for the elastic strain energy dissipation rate. Damage occurs on a very short time scale, leading to the dissipation of the total elastic energy (ψ el + → 0), while viscoplasticity occurs on a longer time scale, leading to non-zero elastic strain energy at equilibrium (ψ el + → ψ ∞ ). A dichotomous search procedure is used to numerically determine the critical shear [γ xy /A] c and tensile [ εyy /A] c ratios required to reach the fracture threshold. Starting from a large interval in which the ratio is expected, we progressively reduce it (depending on whether the failure of the element is reached or not) 22 Fig. 8. Transition between brittle and ductile response as a function of Norton's exponent (n). In the shear case (a) the numerical solution is compared with the analytical solution given by Eq. ( 47). In the tensile case (b), only the numerical solution is shown because there is no explicit analytical solution.

until the size of the interval is less than 1% of its central value.

To change the ratio, we vary the applied strain rate and keep Norton's constant (A) unchanged. ratios as a function of Norton's exponent (n). In the shear case, the analytical solution (red continuous line)

given by Eq. ( 47) is consistent with its numerical counterpart. The critical shear and tensile ratios increase exponentially with Norton's exponent. For a ratio less than the critical value, viscoplasticity dominates the material's behavior and fracture is not observed. By increasing the ratio, the response becomes more and more brittle and fracture appears when the critical ratio is reached.

This first case provided a deeper understanding of the interaction between fracture and viscoplasticity in the proposed model.

Mode I: Tensile opening

The mode I tensile opening test is used to analyze the effect of viscoplasticity on the fracture behavior. is finely meshed with a maximum element size of h = l c /10 to improve the precision of the simulations and reduce the effect of spatial mesh on the results (Molnár et al., 2020a). η = 50% is sufficient to ensure good convergence of the simulation in this case. for a 0 = l c and n = {1, 3}. By increasing the ratio [ εyy /A], σ 0,c yy increases but U t decreases. This result is consistent with experimental observations [START_REF] Verleysen | Quasi-static and high strain rate fracture behaviour of Ti6Al4V[END_REF]. For a high [ εyy /A] ratio, σ 0,c yy converges to the critical load (σ el,0,c yy ) and U t to the energy density until failure (U el t ) in the case of an elastic brittle material. These results show that when the viscoplastic strain rate ( εvp (A)) is comparable to the applied strain rate ( εyy ), DR is observed because the material undergoes viscoplastic strain and absorbs significant energy before fracture. However, if the viscoplastic strain rate is small compared to the applied strain rate, a BR is observed. The material does not have enough time to dissipate the energy through viscoplastic activity.

The correlation between σ 0,c yy and U t is shown in Fig. 11(c). For a given Norton's exponent (n), σ 0,c yy increases as U t decreases. Interestingly, it can be seen that by acting simultaneously on Norton's exponent (n) and on Norton's constant (A), both σ 0,c yy and U t can be improved.

The critical load as a function of the normalized initial crack length is shown in Fig. 12. The elastic (red circle symbol) and viscoplastic (black star symbol) cases are represented. In the elastic case (A = 0), for a 0 � l c , σ 0,c yy varies linearly with a 0 in log-log scale. This result is consistent with the original Griffith's solution. When a 0 is smaller than l c , σ 0,c yy converges to a plateau. This result is similar to that of Molnár et al. (2020a) and is consistent with the size effect found experimentally [START_REF] Bažant | Size Effect in Blunt Fracture: Concrete, Rock, Metal[END_REF]. This mean that the length scale (l c ) introduced in the fracture phase-field approach does not spoil the original Griffith's description.

In the viscoplastic cases, by increasing the ratio [ εyy /A], the viscoplastic solution converges to its elastic counterpart. By reducing the ratio [ εyy /A], viscoplasticity dominates the material's behavior and the effect of the initial crack tends to disappear. The viscoplastic activity at the crack tip mitigates the effect of the initial crack.

Effect of residual stresses on the interaction between fracture and phase transformation

During the growth of an oxide layer, which is the phase formed during the transformation of an alloy, two major phenomena can influence its fracture behavior. These are the induced residual stresses and the evolution of the effective properties of the alloy.

Critical load-

Normalized initial crack length-

� yy A =10 2 [ ] � yy A =10 3 [ ] � yy A =3.10 3 [ ] � yy A =27.10 3 [ ] � yy A [ ] elastic (A=0) viscoplastic
a 0 <l c a 0 >l c Fig. 12. Critical load obtained by phase-field simulations as function of the initial crack length (a 0 ) normalized by l c ; Norton's exponent was set to n = 3 for viscoplastic cases.

Effect of residual stresses on the fracture behavior: homogeneous solution

During the growth of an oxide layer in an alloy, residual stresses can be induced by the growth strain of the oxide layer. This growth strain may be due to the incompatibility of the molar volumes between the newly formed phase and the parent phase [START_REF] Panicaud | On the growth strain origin and stress evolution prediction during oxidation of metals[END_REF]. The growth strain is modeled by an additional eigenstrain (ε * ). The effect of this eigenstrain on the fracture behavior can be understood by a single homogeneous elastic 2D plane strain element. A square plate with dimensions of 1 × 1 is subjected to a simple vertical tensile load. Young's modulus was set to E = 210 with ν = 0.3, g c = 2 × 10 -3 and l c = 2 × 10 -2 . The energy increment constraint was set to η = 1%. A planar eigenstrain is embedded within the element ε * = ε * δ i j , where δ i j is the 2 nd order identity tensor.

Fig. 13(a) and Fig. 13(b) show respectively the evolution of the normalized tensile elastic strain energy (ψ el + ) and the corresponding stress as a function of the applied strain (ε yy ). The reference case (ε * = 0) and the elastic case (ε * = 0.01) are portrayed. For the latter, it can be seen that the material remains in compression as long as ε yy < ε * and the corresponding tensile elastic strain energy which controls the fracture remains at zero. When ε yy > ε * , the material goes into tension and tensile elastic strain energy starts to increase. ε * does not change the critical load compared to the reference case, but induces an offset on the critical strain.

Using the fracture criterion and the stress-strain relation, the critical load (σ c,homo yy ) for this homogeneous case can be obtained by the following equation:

(49) σ c,homo yy = � 2λµ λ + 2µ + 2µ � � � � � � ψ c 2λµ 2 (λ + 2µ) 2 + µ , where λ = Eν (1 + ν)(1 -2ν) , µ = E 2(1 + ν)
are Lame's parameters.

As highlighted in Fig. 13(b), compared to the reference case, the elastic case induces an improvement of U t corresponding to the energy density in compression that does not contribute to fracture.

Heterogeneous case with planar interface

This section investigates the influence of the growth of a misfitting oxide layer (ε * ) on its fracture behavior. We proceed in two steps. (i) first, an oxide layer is generated in a matrix under stress-free conditions.

Despite such conditions, due to the oxide layer growth strain, the local stress in the oxide layer and matrix are non zero during this growth step. (ii) Second, a macroscopic mechanical load is applied at high rate until failure of the previously formed oxide layer (loading step). The effect of internal stresses on growth kinetics has already been studied by [START_REF] Ammar | Modelling inheritance of plastic deformation during migration of phase boundaries using a phase field method[END_REF][START_REF] De Rancourt | Homogenization of viscoplastic constitutive laws within a phase field approach[END_REF]. In this example the mechanical contribution to the growth kinetics is neglected. This is done by choosing the parameters k i large enough.

The illustration of the mechanical problem combined with the chemical problem is depicted in Fig. 14. An equilibrium hyperbolic tangent-like shape is used to set the initial structure:

(50) φ(x, t = 0) = 1 2 � 1 -tanh � θ(x -f o ) δ �� ,
θ is given by Eq. ( 40.2) and δ refer to the interface thickness, f o = 2 × δ is the initial fraction of the oxide layer (φ = 1). The matrix is initially in chemical equilibrium = c m = 0) and a linear gradient is prescribed within the oxide layer. The oxide layer growth is therefore due to the supersaturation prescribed by the Diriclet condition (c = 1.6, note that c is the normalized concentration, therefore its value can exceed 1) on the left side. During the loading step, the mechanical load is applied using a multi-point constraint (MPC) approach as shown in Fig. 14. The average stress corresponding to the force applied to the top face is given by:

(51)

�σ yy � = 1 L � σ yy dx.
L = 1 is the system size, �σ yy � is applied with a sufficiently high rate compared to the viscoplastic strain rate and the oxide layer growth kinetics. This avoids growth of the oxide layer and viscoplastic activity during the loading step.

The material parameters used for the computation are summarized in Table 3. g c of the matrix is set 500 times higher than that of the oxide layer.

The system is meshed with 1000 linear quadrangular elements in the x direction and 1 element in the y direction. Ten elements in the interface are sufficient to resolve the interface [START_REF] De Rancourt | Homogenization of viscoplastic constitutive laws within a phase field approach[END_REF].

Regarding the time step control, the energy constraint of η = 1% is needed during the loading step to ensure good accuracy of the simulations. Dimensionless material parameters used for the simulation of the growth of a misfitting rigid oxide layer on a soft matrix. The subscripts o, m refer to the oxide layer (φ = 1) and the matrix (φ = 0), respectively.

Elasticity:

E m 58 E o 210 ν m = ν o = ν 0.3 ε * 0.01 Viscoplasticity: n m = n o 3 Fracture: g m = 500 × g o 1. l c 2 × 10 -2 Chemical: ω 2 × 10 1 δ 0.01 M 10 6 D m = D o 1 c m 0. c o 1. k m = k o 10 5
layer/matrix interface is indicated by the vertical dotted line. In the elastic case the stress remains homogeneous in both matrix and oxide layer while a gradient with high compressive stresses at the interface is observed in the oxide layer in the viscoplastic case. This stress gradient obtained by simulation is consistent with experimental observations [START_REF] Godlewski | Stress Distribution Measured by Raman Spectroscopy in Zirconia Films Formed by Oxidation of Zr-Based Alloys[END_REF][START_REF] Ly | Understanding crack formation at the metal/oxide interface during corrosion of zircaloy-4 using a simple mechanical model[END_REF] and is due to the interplay between stress relaxation effects and oxide layer growth. The oxide layer that has been present for a longer period of time has had more time than the newly formed oxide layer to relax the stresses.

In the reference case there is obviously no stress generated during the growth step. The position of the As expected for diffusion-controlled processes the kinetics follow a square-root law ( f = √ Kt). K = 1.02 corresponds to the best fit of the data by least square methods. The kinetics is the same in the reference, elastic and viscoplastic cases because as mentioned above, the material parameters are chosen in such a way that the internal stresses do not have any impact on the growth kinetics.

ε x y u y (y=0) =0 L=1 0.1 u x (x=1) =0 u x (x=0) =u 0x =1 =0 c 0 =1.6 c p =1 c m =0 (
The evolution of the residual stress (σ r yy ) in the oxide layer (extracted at the defined control point shown in Fig. 14) as a function of the duration of the growth step is illustrated in Fig. 16(b) for the reference, elastic and viscoplastic cases. In the elastic case, during the growth of the oxide layer, the residual compressive stress in the oxide layer is progressively relaxed and tends towards zero when the matrix is fully transformed ( f = 1). In the viscoplastic case, the stress relaxation is accelerated by the viscoplastic activity. There is obviously no residual stress in the reference case. Fig. 16(c) compares the critical load (�σ c yy �) as a function of the duration of the growth step for the three studied cases. The critical load presented here is related to the oxide layer failure. It can be seen that the critical load increases in all cases as the oxide layer grows. The growth of the oxide layer stiffens the alloy and leads to an improvement of the critical load. When f = 1 the homogeneous solution (see Eq. ( 49)) is reached.

For a given oxide layer fraction, the critical stress is maximal in the elastic case and minimal in the reference case. This shows a beneficial effect of the compressive stresses induced in the oxide layer during its growth. In the viscoplastic case, the critical load curve shifts progressively from the elastic case to the reference case. This is due to the progressive relaxation of the oxide layer residual stress (σ r yy ) during the growth of the oxide layer. Compared to the reference case, the improvement of the critical load in the elastic or viscoplastic case corresponds to the load required to compensate the compressive stress in the oxide layer.

As shown in Fig. 16(d) the observations made for the critical load are also valid for U t .

Similar simulations were carried out by applying a monotonic mechanical load during the growth step.

The load is applied with a low rate so that the oxide layer has sufficient time to grow before its failure. The prescribed stress varies linearly from �σ yy � = 0 at t = 0 to �σ yy � = 16 at t = 1. Fig. 17 shows the stress-strain curve at the control point in the elastic and viscoplastic cases. In the elastic case the stiffening induced by the growth of the oxide layer improves the critical stress. In the viscoplastic case the stiffening is observed at the beginning of the curve after which the softening effect due to the viscoplastic activity progressively takes the lead. The critical strain is thus improved at the expense of the critical stress. The time dependence of the oxide layer growth is given by a square-root law ( f = √ Kt)

while that of the viscoplastic strain is given by a power law ( εvp = Qt n n ≥ 1). This is why the softening effect of viscoplastic activity, necessarily becomes dominant over the stiffening effect of oxide layer growth.

Heterogeneous case with non planar interface : single inclusion

In this section, the interplay between fracture and morphological evolutions during the diffussion-controlled growth of a single inclusion is investigated. As in the previous section, we proceed in two steps. First, the inclusion evolves under stress-free conditions until its fraction reaches f = 0.1 . Then, a macroscopic mechanical load is applied until failure.

The problem is depicted in Fig. 18. Periodic boundary conditions (PBC) are imposed on the edges (∂V) in order to treat the problem only on a Representative Elementary Volume (REV). PBC are applied by constraining all degrees of freedom to take the same value at two opposite points of ∂V. During the loading step, a vertical mean strain is imposed by constraining the displacement difference (Δu) along the y-direction between two opposite points respectively located on the bottom and top faces. The applied PBC are summarized by: (52)

                                         u y (x, y = 1, t) -u y (x, y = 0, t) = Δu, u x (x, y = 1, t) -u x (x, y = 0, t) = 0, u x (x = 0, y, t) -u x (x = 1, y, t) = 0, d(x, y = 1, t) -d(x, y = 0, t) = 0, d(x = 0, y, t) -d(x = 1, y, t) = 0.
The mean strain is applied with a relatively high rate compared to the phase transformation kinetics. This avoids a morphological change of the inclusion during the loading step.

A single circular inclusion which an intial radius of r 0 = 2 × 10 -2 is placed in the center of the matrix.

Initially, the inclusion is saturated (c = c i ) and evolves under the matrix supersaturation (c = c m + 0.5).

The fracture and chemical parameters used for the computation are summarized in Table 3. For sake of brievity, the inclusion and the matrix are supposed to have the same g c . The other parameters are specified for each studied case. A regular mesh of 250 × 250 linear quadrangular elements is used and the energy constraint is set to η = 1%. The viscoplastic activity and residual stresses induced by the inclusion growth strain have no influence on the crack topology in this example.

However an effect of viscoplasticity and residual stresses is observed when examining �σ c yy � and U t as illustrated in Table 4. It can be seen that the values of �σ c yy � and U t in the viscoplastic case are always between the reference and the elastic cases. Stress relaxation induced by viscoplastic activity introduces an intermediate fracture behavior between the reference and elastic behavior.

Heterogeneous case with non planar interface: several inclusions

Here, the model of the previous example is analyzed by considering several initial rigid inclusions that evolve in a soft matrix. This problem, known as the Ostwald ripening phenomenon, is a process that can Table 3 Dimensionless material parameters used for the simulation of the growth of an inclusion. The subscripts i, m respectively refers to the inclusion (φ = 1) and to the matrix (φ = 0). Elastic:

ν m = ν i 0.3 fracture: g m = g i 2 × 10 -3 l c 8 × 10 -3 chemical: ω 5 × 10 -2 δ 0.04 M 10 6 D m = D i = D 1 c m 1.33 c i 1.75 k m = k i 100 E i E m intiation intiation E m E i intiation intiation (a) E i =210, E m =58 (a) E i =58, E m =210
Fig. 19.

Fracture topology during the growth of a single inclusion for the reference case (ε * = 0, A i = A m = 0). In (a) the inclusion is stiffer than the matrix. In (b) the inclusion is softer than the matrix. The arrows indicate the propagation direction from crack initiation. be observed in heterogeneous solids. During such a process, small inclusions dispersed in a matrix evolve to minimize the total free energy of the system [START_REF] Alloyeau | Ostwald Ripening in Nanoalloys: When Thermodynamics Drives a Size-Dependent Particle Composition[END_REF][START_REF] Kabalnov | Ostwald Ripening and Related Phenomena[END_REF]. The evolution takes place mainly by the coalescence of neighboring inclusions or by the growth of large inclusions at the expense of the smallest inclusions as already illustrated by [START_REF] Ammar | Modelling inheritance of plastic deformation during migration of phase boundaries using a phase field method[END_REF]. The aim of this section is to study the effect of such a process on the crack topology.

A two-step analysis is performed. First, the inclusions evolve under stress-free conditions (Ostwald ripening or growth step). Then, a vertical mean strain is applied until failure (loading step). Fig. 20 gives an illustration of the problem.

Several circular inclusions with an initial radius of r 0 = 2 • 10 -2 are randomly distributed in the matrix.

Initially, the inclusions are saturated (c = c i ) and evolve under the matrix supersaturation (c = c m + 0.5).

The material parameters used for the computation are summarized in Table 5. In this example, the chemical Table 4 Energy density until failure (U t ) and critical load (�σ c yy �) after the growth of a single inclusion for the reference case (ε * = 0, A m = A i = 0), elastic case (ε * = 0.05, A m = A i = 0) and viscoplastic case (ε * = 0.05, A m = A i = 100). In all cases, the growth step is stopped when the inclusion fraction reaches f = 0.1. Then, the loading is applied until failure. properties are chosen so that the internal stress has a significant effect on the growth morphology.

Table 5

Dimensionless material parameters used for the simulation of the Ostwald ripening process. The subscripts i, m respectively refer to the inclusion (φ = 1) and the matrix (φ = 0).

elasticity: during the growth process. Interestingly, we observe that despite the initial defect, cracking can be initiated between two particles when they are very close to each other.

E m 58 E i 210 ν m = ν i = ν 0.3 fracture: g m = g i 2 × 10 -3 l c 8 × 10 -3 chemical: ω 5 × 10 -4 δ 0.04 M 10 6 D m = D i 1 c m 1.33 c i 1.75 k m = k i 1
This last example demonstrates the robustness of the proposed multiphase-field model which makes it possible to solve complex shape evolution and crack topology using classical numerical methods.

Conclusion

This work presented a multiphase-field model to describe the fracture behavior in viscoplastic non-static binary alloys. For this purpose, the fracture phase-field and phase-field for phase transformation approaches were succesfully coupled. Several benchmark examples were shown and results based on the critical loading, energy density until failure and crack topology were discussed.

First, the effect of the loading rate on the fracture behavior of a homogeneous viscoplastic material was studied. At low strain rates ductile fracture response is observed. There was sufficient time for significant energy dissipation by developing viscoplastic deformation, whereas at high strain rates brittle fracture response is observed. The material does not have time to dissipate energy by viscoplastic strain. Plastic activity at the crack tip was found to mitigate the effect of pre-existing defects on mechanical strength.

Second, the effect on the fracture behavior of a rigid oxide layer growing on a soft matrix was also studied. It was shown that the stiffening induced by the growth of an oxide layer improves the macroscopic mechanical strength of the alloy formed. At low strain rates, when viscoplasticity is taken into account, the softening effect induced by viscoplasticity ends up dominating over the stiffening effect of oxide layer growth over time. It was also shown that viscoplasticity introduces an intermediate fracture behavior between the reference behavior (without residual stresses induced by the phase transformation) and the elastic behavior (with residual stresses taking into account only the elasticity).

The effect of phase morphology and relative stiffness was also studied. With respect to the loading direction, when a particle is more rigid than the matrix, cracking occurs at the poles of the particle. On the other hand, when the particle is less rigid, the initiation occurs at the equator. In case of several rigid inclusions growing in a soft matrix, it has been observed that the area between neighboring inclusions aligned along the loading direction are preferential zones for crack initiation. The crack bypasses the rigid inclusions encountered on its path. Depending on the size of the encountered inclusion, the crack may bifurcate in the case of small inclusions or bypass without bifurcation in the case of large inclusions.

This work presents only simulated results. The question of the experimental validation of the presented model is an important topic that has not been addressed in this paper. In our forthcoming work, this model will be used to simulate the thermomechanical behavior and oxidation of a fuel rod cladding under pressure (A.10)

B =                                       
                                       
In Eq. (A.2) the tangents matrices are calculated as follows:

(A.11) [K uu ] i j = � V � B T ik � ∂ 2 ψ el ∂ε 2 � kl B l j � dV (A.12) � K uφ � i j = � V � B T ik � ∂ 2 ψ el ∂φ∂ε � k N j � dV (A.13) � K e φu � i j = - � V � N i � ∂ 2 ψ ∂ε∂φ � k B k j � dV (A.14) � K φc � i j = - � V � ∂ 2 ψ ∂c∂φ N i N j � dV (A.15) � K φφ � i j = - � V �� 1 MΔt + ∂ 2 ψ ∂φ 2 � N i N j + αB T ik B k j � dV (A.16) [K dd ] i j = � V � N i N j � g c l c + 2H � + g c l c B T ik B k j � dV (A.17) [K cc ] i j = � V � 1 Δt N i N j + L(φ)B T ik �� ∂ 3 ψ ∂c 3 � (∇c) k N j + � ∂ 3 ψ ∂c∂φ∂c � (∇φ) k N j + ∂ 2 ψ ∂c 2 B k j

��

Note that the global stiffness matrix is not symmetric. The problem is solved with a non symmetric solver. where C is the stiffness matrix.

Appendix C.2. Elastic strain energy derivatives

The derivatives of the elastic strain (ε el ), required for the calculation of the derivatives of the Helmholtz free energy are given below:

(B.7) ∂ε el ∂φ = h � (φ) � ε * b + ε vp b -ε * a -ε vp a � (B.8) ∂ 2 ε el ∂φ 2 = h �� (φ) � ε * b + ε vp b -ε * a -ε vp a �

Fig. 1 .

 1 Fig. 1. (a) Shematic illustration of the idealized physical problem and its corresponding regularization (b). l c and δ respectively control the width of the damaged and the interface regions.

dFig. 2 .

 2 Fig. 2. Overview of the modelled couplings.

Fig. 3 .

 3 Fig.3. Illustration of the staggered scheme for solving the fracture and displacement problems in elastoviscoplastic heterogeneous solids.

Fig. 4

 4 Fig. 4 illustrates this Helmholtz free energy density in Eq. (30) without the elastic part. There are two local minima at c = c a and c = c b which correspond to the equilibrium concentrations of each phase.

Fig. 4

 4 Fig. 4(b) displays the 2D projection onto the (ψc) plane of the 3D representation in Fig. 4(a). The black curve represents the chemical free energy density of phase b (φ = 0) and the red curve corresponds to the free energy density of phase a (φ = 1).

Fig. 4 .

 4 Fig. 4. Illustration of the chemical free energy density.

Fig. 5 .

 5 Fig. 5. Flowchart of the staggered algorithm used to solve the multi-physics problem in Abaqus.

Fig. 6 .

 6 Fig. 6. Stress-strain curves for Brittle Response (BR) and Ductile Response (DR).

Fig. 7 .

 7 Fig.7. Normalized tensile elastic strain energy as a function of the ratio [γ xy /A] c . Norton's exponent was set to n = 1 for comparison with the analytical solution given by Eq. (48).

Fig. 8

 8 Fig. 8(a) and Fig. 8(b) respectively show the evolution of the threshold shear [γ xy /A] c and tensile [ εyy /A] c

Fig. 9 Fig. 9 .Fig. 10 .

 9910 Fig. 9(a) illustrates the problem. The material is homogeneous and only fracture and viscoplasticity are considered. Symmetry boundary conditions make it possible to reduce the problem to a quarter of the geometry Fig. 9(b). The initial crack surface is defined by a Dirichlet boundary condition d = 1 and is leftfree to move. To avoid any significant size effect, the initial crack length is limited to a 0 < 10 -2 . Young's modulus was is set to E = 210 with ν = 0.3, g c = 2 × 10 -3 and l c = 2 × 10 -4 . The region near the initial crack

Fig. 10

 10 Fig. 10 shows the stress-strain curve until fracture for a ratio of [ εyy /A] = 3 × 10 3 and [ εyy /A] = 4 × 10 4 . Norton's exponent and the initial crack length were respectively set to n = 3 and a 0 = l c . For [ εyy /A] = 3 × 10 3 , a Brittle Response (BR) characterized by a high critical load (σ 0,c yy ) and a low energy density until failure (U t ) is observed, whereas for [ εyy /A] = 4 × 10 4 a Ductile Response (DR) characterized by a low σ 0,c yy

Fig. 11 .

 11 Fig. 11. Critical load (a) and energy density until failure (b) as function of the ratio [ εyy /A]. The critical load (σ el,0,c yy ) and the energy density until failure (U el t ) in the case of an elastic brittle material are also shown. The direction of the increase in viscoplastic strain (ε vp ) at the crack tip is illustrated. (c) Correlation between the critical load and the energy density until failure. For all simulations, the initial crack length was set to a 0 = l c .

Fig. 13 .

 13 Fig. 13. (a) Normalized tensile elastic strain energy and (b) stress as a function of the applied tensile strain ε yy The reference case (ε * = 0) and the elastic case (ε * = 0.01) are shown.

Fig. 15

 15 Fig. 15 depicts the stress profile during the growth of the oxide layer for the reference (ε * = A o = A m = 0), elastic (ε * = 0.01, A o = A m = 0) and viscoplastic (ε * = 0.01, A o = A m = 1) cases. The oxide

FFig. 14 .Fig. 15 .

 1415 Fig. 14. Growth of a rigid misfitting oxide layer in a soft matrix. A multi-point constraint (MPC) is applied to the top and left surfaces to enforce them to remain flat. Simulations are carried out in two steps. (i) Firstly, the oxide layer grows without external applied stress (F = 0). (ii) Secondly, a macroscopic mechanical load is applied (F � = 0). The initial concentration is highlighted in red.

Fig. 16 .Fig. 17 .

 1617 Fig. 16. (a) Parabolic oxide layer growth kinetics. The parabolic constant K = 1.02 corresponds to the best fit of the data by least square methods. (b) The residual stress (σ r yy ) at the control point after the oxide layer growth step for the reference (ε * = 0, A o = A m = 0), elastic (ε * = 0.01, A o = A m = 0) and viscoplastic (ε * = 0.01, A o = A m = 1) cases. (c) and (d) respectively show the critical load (�σ c yy �) and U t at oxide layer failure.

  Fig.18. A unit cell of periodic binary material used to study the fracture topology during the growth of a single inclusion. We procede in two steps. Firstly, the saturated inclusion (c(t = 0) = c i ) grows under matrix supersaturation (c(t = 0) = c m + 0.5) without external applied stress (Δu = 0) until its fraction reaches f = 0.1. Secondly, a high rate macroscopic vertical mean strain is applied (Δu = u 0 ) until failure.

Fig. 21 Fig. 21 .

 2121 Fig.21shows the crack topologies obtained at different stages of the Ostwald ripening process in the reference case (ε * = 0, A i = A m = 0). It can be seen that cracks always initiate in the matrix between the nearest inclusions aligned in the loading direction (Fig.21(a)-(c)). The cracks then propagate by bypassing the inclusions encountered along their propagation path (Fig.21(d)-(i)). When the inclusion encountered is small, a bifurcation is observed (Fig.21(d) and (g)), but if the inclusion is large, the crack bypasses

Fig. 22 .

 22 Fig. 22. Crack topology after Ostwald ripening process for the reference case (ε * = 0, A m = A i = 0), elastic case (ε * = 0.01, A m = A i = 0) and viscoplastic case (ε * = 0.01, A m = A i = 1). The duration of the growth step was set to t = 0.5. The loading was applied in vertical direction. The arrows indicate the direction of the crack propagation.

Fig. 22

 22 Fig.22compares the crack topologies obtained between the reference case (ε * = 0, A m = A i = 0), elastic case (ε * = 0.01, A m = A i = 0) and viscoplastic case (ε * = 0.01, A m = A i = 1) after the same duration of the Ostwald ripening process (t = 0.5). A defect (d = 1) has been added to the central nodes of the left and right side of the system. The final morphologies and thus the crack topologies are significantly different in the three cases which shows the importance of taking into account the internal stresses generated

  ∂φ : C : ε el + h � (φ) 2 ε el : (C a -C b ) : ε el + Wg � (φ) + h � (φ)(ψ ch ael ∂φ 2 : C : ε el + 2h � (φ) ∂ε el ∂φ : (C a -C b ) : ε el + ∂ε el ∂φ : C : ∂ε el ∂φ + h �� (φ) 2 ε el : (C a -C b ) : ε el + Wg �� (φ) + h �� (φ)(ψ ch aφ ch b ), (B.4) ∂ 2 ψ ∂c∂φ = h � (φ)k a (cc a ) -h � (φ)k b (cc b ), (B.5) ∂ 3 ψ ∂c∂φ 2 = h �� (φ)k a (cc a ) -h �� (φ)k b (cc b ), (B.6) ∂ 3 ψ ∂c 2 ∂φ = h � (φ)(k a -k b ),

Acknowledgements

The authors thank Antoine Ambard from EDF for the interesting discussions. This work was financially supported by EDF.

at high temperature. In this context, improvements of the model taking into account the microstructure of the oxide film may be necessary. Improvements are also needed to better describe the mechanisms associated with ductile fracture and take into account large deformations (Miehe et al., 2016a). Furthermore, in order to simulate the oxidation under stress phenomenon generally observed at low strain rates, the mass transport and probably the phase transformation equation in the damaged region must be adjusted.

Appendix A. Second law of thermodynamics

Based on the theory of generalized stresses [START_REF] Gurtin | Generalized ginzburg-landau and cahn-hilliard equations based on a microforce balance[END_REF], we note by {σ, ξ, π, ω} the microforce system which are respectively the energetic duals of {∇ s u, ∇φ, φ, d}. σ is the classical cauchy stress tensor, ξ, π are respectively the vector and scalar micro-forces associated with the phase change (Ammar et al., 2009b) and ω the driving force associated with fracture (Miehe et al., 2010b). The total internal energy density is the linear combination of the generalized stresses by their respective dual variables:

the last term of the right-hand sides represents the fracture energy density. The local form of the second law of thermodynamics in the isothermal case is given by (Ammar et al., 2009b) :

where s is the entropy density, T is the temperature, µ is the chemical potential and J is the chemical flux.

The Helmotz free energy (ψ), internal energy (e) and entropy (s) densities are related by the following Legendre transform :

The local mass conservation equation is given by [START_REF] Ammar | Modelling inheritance of plastic deformation during migration of phase boundaries using a phase field method[END_REF]:

Combining Eq. (A.1) to Eq. (A.4) leads to the Clausius-Duhem inequality :

The time derivation of the Helmotz free energy density is given by :

Replacing Eq. (A.6) in Eq. (A.5) gives :

Assuming that the Cauchy tensor (σ), the chemical potential (µ) and the generalized stress vector (ξ) are respectively independent of εel , ċ and ∇ φ gives the following state laws (de Rancourt et al., 2016):

The dissipation rate inequality is then reduced to :

where π dis and ω dis are defined by : (A.10.1)

As done in Ammar et al. (2009b), a convex dissipation potential function is defined to ensure the dissipation rate inequality :

where M is the interface mobility and L(φ) is the homogenized diffusion coefficient (Onsager coefficient).

The complementary evolution laws are then derived from the dissipation potential :

(A.12.1)

The dissipation rate inequality is then further reduced to :

(A.13) D = Ω(σ vm , ∇µ, π dis , φ) + g c γ ≥ 0 This inequality ensures that the second law of thermodynamics is respected.

Appendix B. Finite element implementation

The strong forms of the four coupled problems (displacement,diffusion,phase transformation,fracture) to be solved are summarized as follows:

A Newton-Raphson method is used to solve the multiphysic problem iteratively:

, the corresponding residue vectors read as:

In Eq. (A.4), J is the external chemical flux. The matrix of shape functions and its derivatives used for the discretization of the scalar quantities N, B and the vector quantities N, B are given by: