Pengfei Li

Julien Yvonnet
email: julien.yvonnet@univ-paris-est.fr

Christelle Combescure

An extension of the phase field method to model interactions between interfacial damage and brittle fracture in elastoplastic composites

Keywords: Fracture, Phase field method, Interfacial damage, Crack propagation, Elastoplastic microstructures

An extension of the phase field method to model interfacial damage in elastoplastic composites is proposed. In the matrix, an elastoplastic phase field is employed to model the fracture process. To introduce interfacial damage between inclusions and the matrix, a strain density function depending on the jump due to decohesion is added to the total energy. Smooth indicator functions are used to maintain the regularized character of the approximation. They weight the different terms in the energy with respect to the vicinity of interfaces. Then, the different problems (mechanical and phase field problems) are derived and an algorithmic procedure is described. Numerical examples show the capabilities of the method to handle initiation, propagation and interactions between both elastoplastic fracture and interfacial cracks in complex elastoplastic composite microstructures. It is also shown that the solutions are convergent with respect to the mesh refinement.

Introduction

The evaluation of composite resistance to cracking is of major importance in engineering. Recently, several advances in formulations, numerical methods and computer capabilities have opened the route to analyze damage of composites and heterogeneous materials directly from their microstructure, by simulating microcracking initiation, propagation and merging up to global failure. Such analyses offer new possibilities in material design, such as optimization of the constituents composition/morphologies [START_REF] San | Optimization of carbon black polymer composite microstructure for rupture resistance[END_REF][START_REF] Gu | Optimization of composite fracture properties: Method, validation, and applications[END_REF][START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF][START_REF] Da | Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage[END_REF], or investigation of microcracking mechanisms directly from realistic microstructures, such as the ones arising from experimental micro-CT images [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF][START_REF] Nguyen | Initiation and propagation of complex 3d networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF].

However, several challenges are related to microstructure-based damage analyses: (i) the development of appropriate formulations and related numerical methods to handle initiation, propagation and interactions of micro cracks in complex geometrical domains related to heterogeneous microstructures; (ii) the computational complexity, which warrants describing all heterogeneities in the whole sample; (iii) the presence of interfaces which may involve different damage mechanisms as compared to the bulk cracks.

Regarding the first point, many numerical methods have been proposed in the past decades to simulate crack propagation, including: discrete approaches (peridynamics [START_REF] Silling | Reformulation of elasticity theory for discontinuities and long-range forces[END_REF], discrete elements [START_REF] Hedjazi | Application of the discrete element method to crack propagation and crack branching in a vitreous dense biopolymer material[END_REF]), cohesive models [START_REF] Needleman | An analysis of decohesion along an imperfect interface[END_REF][START_REF] Camacho | Computational modelling of impact damage in brittle materials[END_REF][START_REF] Yang | Fully automatic modelling of cohesive discrete crack propagation in concrete beams using local arc-length methods[END_REF], Augmented Finite Elements [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF], FEM remeshing techniques, XFEM [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Sukumar | Extended finite element method for three-dimensional crack modeling[END_REF], or continuum damage models with regularization [START_REF] De Borst | Simulation of strain localization: a reappraisal of the cosserat continuum[END_REF][START_REF] Peerlings | Gradient-enhanced damage for quasi-brittle materials[END_REF][START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF]. Discrete approaches are suitable for fragile cracking but they require ad-hoc calibrations and cannot be easily extended to other behaviors like elastoplasticity. In methods based on sharp fracture description (FEM, XFEM), it is very difficult to model crack initiation and extension to complex 3D patterns may be highly cumbersome. Augmented FEM suffers from stability and mesh-dependency issues. While local damage models are well-known to be associated with mesh-dependency issues and lack of energetic convergence, nonlocal models such as [START_REF] Peerlings | Gradient-enhanced damage for quasi-brittle materials[END_REF][START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] require nontrivial extensions of FEM discretizations.

The development of the variational approach to fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF], also called phase field method in the literature, offers several advantages in this context, like: the possibility to initiate cracks from undamaged configurations; the possibility to handle arbitrary crack networks (including branching, merging, in both 2D and 3D) without specific treatment and use of classical finite elements; a variational framework allowing to include many models or mechanisms, and a mesh-independence due to an appropriate regularization process. The method requires a fine mesh along the crack path and a suitable definition for regularization parameters (see a discussion in [START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF] and recent internal length-insensitive formulations in [START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF][START_REF] Mandal | Length scale and mesh bias sensitivity of phase-field models for brittle and cohesive fracture[END_REF]). However, due to the above-mentioned advantages, the phase field method has been widely developed and applied to many problems, such as, among many others: brittle fracture [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF][START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF][START_REF] Zhang | A modification of the phase-field model for mixed mode crack propagation in rock-like materials[END_REF], composite delamination [START_REF] Roy | Phase field based peridynamics damage model for delamination of composite structures[END_REF], dynamic fracture [START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF][START_REF] Doan | Hybrid phase field simulation of dynamic crack propagation in functionally graded glass-filled epoxy[END_REF][START_REF] Liu | Abaqus implementation of monolithic and staggered schemes for quasi-static and dynamic fracture phase-field model[END_REF], hydraulic fracture [START_REF] Wilson | Phase-field modeling of hydraulic fracture[END_REF][START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. part iii. crack driving forces in hydro-poro-elasticity and hydraulic fracturing of fluid-saturated porous media[END_REF][START_REF] Xia | Phase field modeling of hydraulic fracturing with interfacial damage in highly heterogeneous fluid-saturated porous media[END_REF][START_REF] Ehlers | A phase-field approach embedded in the theory of porous media for the description of dynamic hydraulic fracturing[END_REF], topology optimization for resistance to cracking [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF][START_REF] Da | Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage[END_REF], anisotropic material fracture [START_REF] Nguyen | Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials[END_REF][START_REF] Bleyer | Phase-field modeling of anisotropic brittle fracture including several damage mechanisms[END_REF][START_REF] Zhang | Phase field modeling of fracture in fiber reinforced composite laminate[END_REF], ductile fracture [START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF][START_REF] Ambati | A phase-field model for ductile fracture at finite strains and its experimental verification[END_REF][START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. part ii. coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[END_REF][START_REF] Borden | A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects[END_REF][START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF], ductile/fragile transition [START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF][START_REF] You | Incorporation of tension-compression asymmetry into plastic damage phase-field modeling of quasi brittle geomaterials[END_REF], fracture in micro tomography image-based models of microstructures [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF][START_REF] Nguyen | Initiation and propagation of complex 3d networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF][START_REF] Nguyen | Large-scale simulations of quasibrittle microcracking in realistic highly heterogeneous microstructures obtained from micro ct imaging[END_REF] and more recently adapted in machine learning strategies in [START_REF] Goswami | Transfer learning enhanced physics informed neural network for phase-field modeling of fracture[END_REF].

Another issue related to microstructural damage is the presence of many interfaces, which may constitute additional weak parts and can strongly impact the damage mechanisms. Several works have been devoted to describe interactions between bulk and interfacial damage: Paggi et al. [START_REF] Paggi | Revisiting the problem of a crack impinging on an interface: a modeling framework for the interaction between the phase field approach for brittle fracture and the interface cohesive zone model[END_REF] proposed a novel combined approach in which bulk brittle fracture is treated by phase field model and interface debonding is modelled by classical cohesive elements. Msekh et al. [START_REF] Msekh | Predictions of j integral and tensile strength of clay/epoxy nanocomposites material using phase field model[END_REF][START_REF] Msekh | Fracture properties prediction of clay/epoxy nanocomposites with interphase zones using a phase field model[END_REF] applied the phase field approach to model progressive failure in polymer-matrix composites. Hansen-Dörr et al. [START_REF] Hansen-Dörr | Phase-field modelling of interface failure in brittle materials[END_REF] proposed a phase field approach for interface failure between two dissimilar materials. In [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF][START_REF] Nguyen | Phase field modeling of interfacial damage in heterogeneous media with stiff and soft interphases[END_REF], a simple framework to extend phase field to interfacial damage has been proposed. A regularized description of interface displacement jumps has been used as well as appropriate weights indicator functions to associate different damage behavior to the bulk and to the interfaces in the energy.

In this work, we follow the framework developed in [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF][START_REF] Nguyen | Phase field modeling of interfacial damage in heterogeneous media with stiff and soft interphases[END_REF] to interfacial damage in composites whose phases are elastoplastic. For this purpose, an elastoplastic phase field is employed to model the fracture process as proposed in [START_REF] Duda | A phase-field/gradient damage model for brittle fracture in elastic-plastic solids[END_REF]. To introduce interfacial damage between inclusions and the matrix, a strain density function depending on the jump due to decohesion is added to the total energy. To maintain the regularized character of the approximation, smooth indicator functions are used to weight the different terms in the energy with respect to the vicinity of interfaces. Then, the different problems (mechanical and phase field problems) are derived and an algorithmic procedure is described.

The present paper is organized as follows. In Section 2, we introduce diffuse approximations for crack and interface displacement jumps in a regularized context. In Section 3, a phase field method for brittle fracture in elastoplastic solids with interfacial damage is proposed. Section 4 provides the weak forms and FEM discretizations and the flowcharts for the overall algorithm. Finally, numerical examples are presented in Section 5.

Diffuse approximation of discontinuous fields

Let Ω ⊂ R D be an open domain with D = 2, 3, describing a heterogeneous medium which contains internal interfaces between two elasto-plastic phases. The external boundary of Ω is denoted by ∂Ω ∈ R D-1 . During the loading, cracks may propagate in the medium phases and can pass through the interfaces as depicted in Fig. 1(a), where the crack surfaces and the interfaces are collectively denoted by Γ and Γ I , respectively. In this work, we adopt the framework proposed in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Ambrosio | Approximation of functional depending on jumps by elliptic functional via tconvergence[END_REF] for a regularized representation of discontinuities extended to interfaces as in [START_REF] Xia | Phase field modeling of hydraulic fracturing with interfacial damage in highly heterogeneous fluid-saturated porous media[END_REF][START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]. In this regularized framework, the cracks are approximately represented by a scalar phase field d (x, t) (see Fig. 1(b)) and the interfaces by a fixed scalar function β (x) (see Fig. 1(c)).

Phase field approximation of bulk cracks and interfaces

For a known fixed crack surface Γ (see Fig. 1(a), the scalar crack phase field d (x, t) can be determined through solving the following boundary value problem subjected to Dirichlet boundary conditions d = 1 on the crack (see [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] for more details):

           d (x, t) -2 d ∆d (x, t) = 0 in Ω, d (x, t) = 1 on Γ, ∇d (x, t) • n = 0 on ∂Ω, (1)
where ∆ (•) and ∇ (•) are the Laplacian and gradient operator respectively, d is a length scale parameter that governs the width of the regularization zone and gives for d → 0 the exact sharp crack in Fig. 1(a), and n the outward normal on ∂Ω. It can be shown that (1) is the Euler-Lagrange equation associated with the variational problem:

d = Arg inf d∈S d Γ d (d) , Γ d (d) = Ω γ d (d) dΩ, S d = {d | d (x) = 1, ∀x ∈Γ} , (2)
where Γ d (d) represents the total length of the crack in 2D and the total crack surface area in 3D, and γ d (d) is the crack surface density function per unit volume defined by:

γ d (d) = d 2 2 d + d 2 ∇d • ∇d, (3)
where the second term in γ d (d) penalizes high values of ∇d(x) and where d varies between 0 and 1. Note that in the absence of the second right-hand term in (3), a local damage model is found, with well-known related non-convergence issues with respect to the mesh discretization.

It must be noted that, d does not represent physically the exact crack width, but a parameter which is used to regularize the discontinuities. It has been shown that this parameter can be treated as a material parameter related to the Young's modulus, the tensile strength, and the critical energy release rate of the material in [START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF][START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Zhang | Numerical evaluation of the phase-field model for brittle fracture with emphasis on the length scale[END_REF]. In our previous work [START_REF] Nguyen | Initiation and propagation of complex 3d networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF], an inverse approach was developed to identify this parameter, by combining simulations and experiments.

The scalar interface function β (x) can also be determined through solving the following boundary value problem subjected to Dirichlet boundary conditions β = 1 on the interfaces Γ I :

           β (x) -2 β ∆β (x) = 0 in Ω, β (x) = 1 on Γ I , ∇β (x) • n = 0 on ∂Ω, ( 4
)
where β is a length scale parameter which governs the width of the regularization zone of the interface and gives for β → 0 the exact sharp interfaces in Fig. 1(a). Similarly, (4) corresponds to the Euler-Lagrange equation of the variational problem:

β = Arg inf β∈S β Γ β (β) , Γ β (β) = Ω γ β (β) dΩ, S β = β | β (x) = 1, ∀x ∈Γ I , (5)
where γ β (β) is defined by:

γ β (β) = β 2 2 β + β 2 ∇β • ∇β. (6)
In the present work, the length scale parameter d and β are interpreted as material parameters. Then, one possibility to select these parameters is to perform an inverse identification from experimental data. Such identification procedure has been described in [START_REF] Nguyen | Initiation and propagation of complex 3d networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF] in the context of quasibrittle heterogeneous materials and could be extended to the present elastoplastic framework in future studies using similar ideas. Then, in the present paper, the numerical values of d and β have been chosen arbitrarily as small values in the numerical examples of section 5.

Regularized representation of the displacement jumps within interfaces

In this section, an approximation for the displacement jump [[u]] is introduced to consider the interface debonding (see [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]). Let Γ I be the interface. We define Γ I as the zero level-set of a function φ (x), such that (see [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF] for its construction): where Ω i and Ω m denote the set of inclusions and matrix, respectively. Using a Taylor expansion of the displacement field around a point x located on the interface (see Fig. 2):

           φ (x) > 0 for x ∈Ω i , φ (x) < 0 for x ∈Ω m , φ (x) = 0 for x ∈Γ I , (7)
[[u]] w (x) = u x+ h 2 n I -u x- h 2 n I = h∇ (u (x)) n I , (8)
where n I is the normal vector to Γ I at the point x defined by:

n I = ∇φ (x) ∇φ (x) (9)
and w (x) denotes the smooth displacement jump approximation. Above, h is a small scalar parameter, chosen as h = h e to minimize the estimation error (see [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]), with h e the typical element size of the finite element mesh around the interface. A detailed description for the numerical computation of φ (x) is introduced in [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF].

It should be noted that φ (x) and n I do not change throughout the simulation because the interfaces do not evolve. For a sharp interface, the singular part of the strain along the interface can be defined by:

ε I (x) = n I (x)⊗ s [[u(x)]] δ(x) ∀x ∈ Γ. ( 10
)
Then using the above regularization framework, this expression can be approximated as [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]:

ε I (x) n I (x)⊗ s w(x)γ β (x) ∀x ∈ Ω. ( 11
)

Phase field modeling of elastoplastic damage interacting with interfacial damage

In this section, we propose a phase field model for modeling interactions between interfacial damage and bulk cracking in elasto-plastic composites. The framework is developed using the variational framework for fracture as introduced in [START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF][START_REF] Rodriguez | A variational approach to the phase field modeling of brittle and ductile fracture[END_REF].

First, we introduce the total strain as:

ε = ε e + ε p + ε I , (12)
where ε I has been defined in [START_REF] Yang | Fully automatic modelling of cohesive discrete crack propagation in concrete beams using local arc-length methods[END_REF] and ε p is the plastic strain. The plastic incompressibility is assumed, i.e. Tr(ε p ) = 0, Tr(.) being the trace operator. The cumulated plastic strain is defined as:

ṗ = 2 3 εp (13)
where ( .) denotes time derivative. In the following, we introduce the deviatoric parts of the stress and of the elastic strain tensors as s = dev(σ) and e e = dev(ε e ), respectively, with dev(.) = (.) -1 3 Tr(.)1, 1 being the first-order identity operator. Finally, we define the directional derivative as:

D v f (u) = d dh f (u + hv) h=0 . ( 14
)
The damage variable d is introduced to penalize the stifness of the material.

Total energy

Let us consider a two-phase medium whose phases are elasto-plastic with possible damage. Small strains are assumed. The total energy of the medium is defined by

W (u, p, d, β) = Ω ψ e (ε (u) , p, d) dΩ + Ω ψ p (p) + ψ d (β, d) + ψ I (w (u) , β) dΩ - ∂Ω F F • udS - Ω f • udΩ, (15)
where ψ e , ψ p , ψ d and ψ I denote the elastic strain density function, the plastic and damage dissipative potentials, and a strain density function depending on the approximated displacement jump across the interfaces, respectively. Above, f and F are body forces and prescribed traction over the boundary ∂Ω F , respectively. Note that the total energy is a function of 3 state variables u, p and d, and of one field β which plays the role of a parameter, as this field does not change during the evolution of the system.

Variational framework

We follow the framework presented in [START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF] to construct the variational principle, which involves: irreversibility condition, stability condition and energy balance. In this framework, stability condition provides mechanical balance equation, damage and plastic criteria. The energy balance provides damage consistency and plastic flow rule.

Irreversibility condition

The irreversibility condition is imposed on the damage variable to disallow material regeneration. It can be expressed as

ḋ ≥ 0, 0 ≤ d ≤ 1. ( 16
)
In the present work, the irreversibility condition is implemented by using an appropriate history function [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] (see section 3.3).

First-order stability condition

The first order stability condition (see [START_REF] Mielke | Evolution of rate-independent systems[END_REF][START_REF] Mielke | A mathematical framework for generalized standard materials in the rate-independent case[END_REF][START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF]) is expressed by:

D δu W (u, p, d) + D δp W (u, p, d) + D δd W (u, p, d) ≥ 0. ( 17
)
Applied to [START_REF] De Borst | Simulation of strain localization: a reappraisal of the cosserat continuum[END_REF], it yields:

Ω σ : ε e (δu) + ∂ψ I ∂w •w (δu) dΩ + Ω - 3 2 σ : n + ∂ψ p ∂p δpdΩ + Ω ∂ψ e ∂d + ∂ψ d ∂d δddΩ - ∂Ω F F • δudS - Ω f • δudΩ ≥ 0 (18)
where

σ = ∂ψ e ∂ε e ( 19
)
and n is a unit tensor in the direction of the plastic flow. From [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF]:

ε e = ε -ε p -ε I (20)
and then:

D δu ε e = ε e (δu) = ε(δu) -ε I (δu). (21)
From [START_REF] Camacho | Computational modelling of impact damage in brittle materials[END_REF],

ε I (δu) = n I ⊗ s w(δu)γ β , w(δu) = h∇(δu)n I . ( 22
)
From this expression, the following results stem out:

• For δp = δd = 0 we obtain:

Ω σ : ε e (δu) dΩ + Ω ∂ψ I ∂w •w (δu) dΩ - ∂Ω F F • δudS - Ω f • δudΩ = 0 (23)
which is the weak form of the equilibrium equation.

• For δd = 0 and δu = 0:

Ω - 3 2 σ : n + ∂ψ p ∂p δpdΩ ≥ 0. ( 24
)
For J 2 -plasticity, this expression leads to

Ω 3 2 s - ∂ψ p ∂p δpdΩ ≤ 0, (25)
which is the weak form of the plasticity yield criterion. The local form of the plastic yield criterion can be expressed as

F p (p) = 3 2 s - ∂ψ p ∂p ≤ 0 in Ω ( 26
)
which is the classical von Mises yield criterion.

• For δp = 0 and δu = 0:

Ω ∂ψ e ∂d + ∂ψ d ∂d δddΩ ≥ 0 (27)
which is the weak form of the damage criterion. In local form, it can be expressed as:

F d (d) = - ∂ψ e ∂d + ∂ψ d ∂d ≤ 0 in Ω. ( 28
)

Energy balance

The energy balance represents the need for the total energy to remain constant as the state variables evolve. Following a procedure analogous to the treatment of the stability condition, this condition leads to

Ω -σ : ε e ( u) - ∂ψ I ∂w •w ( u) + 3 2 σ : n - ∂ψ p ∂p ṗ - ∂ψ e ∂d + ∂ψ d ∂d ḋ dΩ + ∂Ω F F • udS + Ω f • udΩ = 0. ( 29
)
The following cases are analyzed.

• For u = 0 and ḋ = 0, and using [START_REF] Mandal | Length scale and mesh bias sensitivity of phase-field models for brittle and cohesive fracture[END_REF], the plasticity consistency condition is obtained:

F p (p) ṗ = 0. ( 30
)
• For u = 0 and ṗ = 0, and using (28), the damage consistency condition is obtained:

F d (d) ḋ = 0. ( 31
)

Alternate minimization

In this section, a staggered alternate minimization algorithm is applied, which naturally stems out from the energetic principles. This procedure takes advantage of the fact that although the global energy is non-convex, it is convex with respect to u, p and d individually [START_REF] Alessi | Variational approach to fracture mechanics with plasticity[END_REF]. With the total energy (15) at hand, the alternate minimization follows.

• Minimization with respect to the displacement field:

D δu W (u, p, d) = 0 ( 32
)
leads to

Ω σ : ε e (δu) dΩ + Ω ∂ψ I ∂w •w (δu) dΩ - ∂Ω F F • δudS - Ω f • δudΩ = 0 ( 33
)
which corresponds to the weak form of the mechanical problem to be solved for u, given d.

• Minimization with respect to the equivalent plastic strain:

D δp W (u, p, d) = Ω - 3 2 σ : n + ∂ψ p ∂p δpdΩ = 0 (34)
which is the weak form of the plastic yield criterion (24) which has to be satisfied for ṗ ≥ 0.

In the present work, this condition is handled by a return-mapping algorithm (see [START_REF] Simo | Computational Inelasticity[END_REF]). In [START_REF] Rodriguez | A variational approach to the phase field modeling of brittle and ductile fracture[END_REF] a regularization term was introduced in the total energy and the above equation was verified through solving a global problem for p. Here we do not adopt this approach and treat this criterion as a local one (at Gauss integration points). Eqs. [START_REF] Wilson | Phase-field modeling of hydraulic fracture[END_REF] and [START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. part iii. crack driving forces in hydro-poro-elasticity and hydraulic fracturing of fluid-saturated porous media[END_REF] are solved together using the return-mapping algorithm (see Algorithm 1).

• Minimization with respect to the damage field:

D δd W (u, p, d) = Ω ∂ψ e ∂d + ∂ψ d ∂d δddΩ = 0. ( 35
)
which corresponds to the global problem to be solved to find the field d(x) (phase field problem), given u and p. Note that if ψ d does not include the gradient of damage ∇d, then [START_REF] Xia | Phase field modeling of hydraulic fracturing with interfacial damage in highly heterogeneous fluid-saturated porous media[END_REF] leads to a local relationship to define the damage evolution, leading to well-known issues like non-convergence with respect to the mesh size, or dependence of the crack paths to the mesh structure and orientation.

To avoid these issues, we adopt a form similar to (3) for ψ d , following [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] but adding influence of interface weighting function β as detailed in the following.

Specialization

We now specialize the energy components introduced in [START_REF] De Borst | Simulation of strain localization: a reappraisal of the cosserat continuum[END_REF]. Inspired by [START_REF] Xia | Phase field modeling of hydraulic fracturing with interfacial damage in highly heterogeneous fluid-saturated porous media[END_REF][START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF], we define the different strain density functions as follows:

ψ e (u, d) = g (d) ψ e+ 0 (ε e (u)) + ψ e- 0 (ε e (u)) , (36)
with [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF]

ψ e+ 0 (ε e ) = 1 2 κ Tr (ε e ) 2 + + µe e : e e , (37)
ψ e- 0 (ε e ) = 1 2 κ Tr (ε e ) 2 -. (38)
Above, κ and µ denote bulk modulus and shear modulus, respectively,

x ± = 1 2 (x ± |x|). g(d) = (1 -d) 2 + k, ( 39
)
where k is a small numerical parameter to prevent loss of definite posedness of the elastic tensor in case of full damage,

ψ p (p) = σ Y p + 1 2 Hp 2 , (40)
ψ d (β, d) = (1 -β) 2 g c d 2 2 d + d 2 ∇d • ∇d , (41)
In [START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF], σ Y and H are yield stress and hardening modulus, respectively. Notice that the elastoplastic brittle fracture case model can be recovered by simply setting β = 0 in [START_REF] Ambati | A phase-field model for ductile fracture at finite strains and its experimental verification[END_REF]. Following our recent work [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF], we specialize the energy of interfacial jump component as:

ψ I (u, β) = ψ I (w(u)) γ β (β) , (42)
where the traction acting on the interface oriented by n I (see Fig. 2) is defined by:

t (w) = ∂ψ I (w) ∂w . (43)
Above, ψ I is the interface strain density depending on the regularized displacement jump w.

Note that as discussed in [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF], it is not required in this formulation to include history variables for the interfaces, the diffuse damage field being used to describe the irreversibility of the interfacial damage.

For 2D problems, the traction t (w) is in the form

t (w) = t n , t t T , ( 44
)
where t n and t t denote normal and tangential parts of the traction vector t. In this paper, we neglect the effects of the tangential component in the interface model (t t = 0), and use [START_REF] Van Den Bosch | An improved description of the exponential xu and needleman cohesive zone law for mixed-mode decohesion[END_REF]: where δ n is related to the interface fracture toughness g I c and the interface fracture strength t u by:

t n = g I c w n (δ n ) 2 exp - w n δ n , (45)
δ n = g I c t u exp (1) , (46)
w n is the normal displacement jump defined by:

w n = w•n I , (47)
and w and n I are given in ( 8) and [START_REF] Needleman | An analysis of decohesion along an imperfect interface[END_REF], respectively. This relationship is illustrated in Fig. 3, where the fracture toughness g I c denotes the value of the interface energy function ψ I at full crack opening. Relation [START_REF] You | Incorporation of tension-compression asymmetry into plastic damage phase-field modeling of quasi brittle geomaterials[END_REF] can be easily obtained by expressing the extremum of the function [START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF] (∂t n (w n )/∂w n = 0), which is found at w n = δ n , and then simply replacing this solution in [START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF].

Using the aforementioned constitutive specialization results, we can re-write ( 19), ( 26) and ( 13) as:

σ = g (d) σ + eff + σ - eff , (48)
F p (p) = 3 2 s -(σ Y + Hp) ≤ 0, ( 49
) εp = ṗ 3 2 s s with ṗ ≥ 0 ( 50
)
where σ + eff and σ - eff are the effective tensile and compressive stresses defined as

σ + eff = ∂ψ e+ 0 ∂ε e = κ Tr (ε e ) + 1 + 2µe e , σ - eff = ∂ψ e- 0 ∂ε e = κ Tr (ε e ) -1. (51)
The associated Euler-Lagrange equations to [START_REF] Wilson | Phase-field modeling of hydraulic fracture[END_REF] are given by:

           ∇ • σ -γ β (β)f I + f = 0 in Ω, u = ū on ∂Ω u , σn = F on ∂Ω F . ( 52
)
where the body force term f I is expressed by f I = -h∇t(w) • n I , and is obtained from [START_REF] Wilson | Phase-field modeling of hydraulic fracture[END_REF] by applying the divergence theorem to the term Ω γ β (β)t(w)h∇(δu)dΩ, and assuming t • n = 0 over the external boundary of the domain.

Using [START_REF] Xia | Phase field modeling of hydraulic fracturing with interfacial damage in highly heterogeneous fluid-saturated porous media[END_REF] and the property:

(∆d) δd = ∇ • (∇dδd) -∇d • ∇(δd) (53)
as well as the divergence theorem and ∇d • n = 0, we obtain the weak form of the phase field problem as:

Ω -2(1 -d)ψ e+ 0 + (1 -β) 2 g c d d δd + g c (1 -β) 2 d ∇d • ∇(δd) dΩ = 0. ( 54
)
To prescribe irreversibility, we employ the technique introduced in [START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. part i. balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids[END_REF] which consists in substituting the above weak form with:

Ω -2(1 -d)H(ε e ) + 2(1 -β) 2 ψ c d δd + 2(1 -β) 2 ψ c 2 d ∇d • ∇(δd) dΩ = 0, (55)
where

H (ε e ) = max s∈[0,t] ψ e+ 0 (ε e , s) -ψ c + (56)
and ψ c is a specific fracture energy density, which can be further related to a critical fracture stress σ c by:

ψ c = 1 2E σ 2 c , ( 57
)
where E is the Young's modulus. The corresponding Euler-Lagrange equations to (55) are given by:

           (1 -β) 2 ψ c d -2 d ∆d = (1 -d) H (ε e ) ∇d • n = 0 on ∂Ω, d = 1 on Γ. ( 58
)
Irreversibility ḋ ≥ 0, 0 ≤ d ≤ 1 Mechanical balance ∇ • σ -γ β (β)f I + f = 0 in Ω u = ū on ∂Ω u , σn = F on ∂Ω F Constitutive law σ = g (d) σ + eff + σ - eff σ + eff = κ Tr (ε e ) + 1 + 2µe e , σ - eff = κ Tr (ε e ) -1 Cohesive law t (w) = g I c w n δ 2n exp -w n δ n , 0 T Damage criterion (1 -β) 2 ψ c d -2 d ∆d -(1 -d) H(ε e ) ≥ 0 H (ε e ) = max s∈[0,t] ψ e+ 0 ( ε e , s) -ψ c + Damage consistency (1 -β) 2 ψ c d -2 d ∆d -(1 -d) H(ε e ) ḋ = 0 Plastic yield criterion F p (p) = 3 2 s -(σ Y + Hp) ≤ 0, Plastic flow rule εp = ṗ 3 2 s s with ṗ ≥ 0
The different equations of the model are summarized in Table 1. The weak forms of the problems to be solved alternatively are summarized in Table 2.

Remark: In the present work, we did not consider coupling between damage and plasticity. The reason is that currently, the available models from the literature which introduce such phenomena are not able to control the evolution of the plastic strain when damage occurs. In other words, the plastic strain continues to increases when the damage occurs within the crack. So far, it seems that the only model able to deal with this issue is the the model proposed in [START_REF] Miehe | Phase field modeling of ductile fracture at finite strains: A variational gradient-extended plasticity-damage theory[END_REF], involving gradient plasticity and gradient damage and strong plasticity-damage coupling. Extending the present interfacial damage framework with the gradient plasticity model in [START_REF] Miehe | Phase field modeling of ductile fracture at finite strains: A variational gradient-extended plasticity-damage theory[END_REF] could be a perspective for this work.

Discretization and numerical implementation

In this section, we detail the weak forms and FEM discretizations for displacement and damage problems, and finally provide the different algorithms.

Weak forms the linearized mechanical problem

Using ( 42) and ( 43), we can re-write the associated weak form for the displacement problem (33) as:

Ω σ : ε e (δu) dΩ + Ω γ β t•w (δu) dΩ = Ω f •δudΩ + ∂Ω F F•δudS, (59)
where t and σ are given in [START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF] and [START_REF] Goswami | Transfer learning enhanced physics informed neural network for phase-field modeling of fracture[END_REF], respectively. We can rewrite the balance equation ( 59) as Mechanical problem: given d, solve for u:

R = Ω σ : ε e (δu) dΩ + Ω γ β t•w (δu) dΩ - Ω f •δudΩ - ∂Ω F F•δudS = 0. ( 60
)
Ω σ : ε e (δu) dΩ + Ω γ β (β)t(w)•w (δu) dΩ -∂Ω F F • δudS -Ω f • δudΩ = 0 F p (p) = 3 2 s -(σ Y + Hp) ≤ 0 εp = ṗ 3 2 s
s with ṗ ≥ 0 Phase Field problem: given u, solve for d:

Ω H(u) + (1 -β) 2 ψ c dδd + (1 -β) 2 ψ c 2 d ∇d • ∇(δd) dΩ = Ω H(u)δddΩ.
In a standard Newton method, the displacements are updated for each loading increment by solving the tangent problem:

D ∆u R u (k) , d = -R u (k) , d = 0, ( 61
)
where u (k) is the displacement solution known from the previous iteration. The displacement corrections are obtained as

u (k+1) = u (k) + ∆u. (62)
In [START_REF] Mielke | A mathematical framework for generalized standard materials in the rate-independent case[END_REF],

D ∆u R u (k) = Ω ∂σ ∂ε e : ε e (∆ε) : ε e (δε) dΩ + Ω γ β ∂t (w) ∂w : ∆w : δwdΩ, (63)
with

∂σ ∂ε e = C s (u, d) . (64)
Using [START_REF] Zhang | Numerical evaluation of the phase-field model for brittle fracture with emphasis on the length scale[END_REF], we can re-write the associated weak form for the damage problem [START_REF] Xia | Phase field modeling of hydraulic fracturing with interfacial damage in highly heterogeneous fluid-saturated porous media[END_REF] as:

Ω H+ (1 -β) 2 ψ c dδd + (1 -β) 2 ψ c 2 d ∇d•∇ (δd) dΩ = Ω Hδd dΩ. (65)

Discretization of the displacement problem

In this work, for the sake of clarity, only 2D FEM discretization is detailed. The vector form of second-order tensors are introduced as

[ε] = ε 11 , ε 22 , √ 2ε 12 T , [σ] = σ 11 , σ 22 , √ 2σ 12
T , as well as the FEM approximations u = N u u e , δu = N u δu e , and ∆u = N u ∆ u e where u e , δu e , ∆u e and N u are nodal displacement components in one element, nodal trial function components, nodal incremental displacement components and a matrix of displacement shape function, respectively. Then we obtain:

[ε] (u) = B u u e , [ε e ] (∆u) = B w ∆u e , [ε e ] (δu) = B w δu e , (66)
and

w (u) = hNBu e , ∆w (u) = hNB∆u e , δw = hNBδu e , ( 67
)
where B u is a matrix of displacement shape function derivatives, and

N =   n 1 n 2 0 0 0 0 n 1 n 2   , (68)
B = ∂u 1 ∂x 1 ∂u 1 ∂x 2 ∂u 2 ∂x 1 ∂u 2 ∂x 2 T =   ∂ ∂x 1 ∂ ∂x 2 0 0 0 0 ∂ ∂x 1 ∂ ∂x 2   T N u , (69)
where n 1 and n 2 are the x-and y-components of the normal vector n I in ( 9) and B w is a modified spatial strain-displacement matrix defined by:

B w = B u -hγ β MB (70)
in which M is a matrix expressed by [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]:

M =       n 1 0 0 n 2 1 √ 2 n 2 1 √ 2 n 1       N =       n 2 1 n 1 n 2 0 0 0 0 n 1 n 2 n 2 2 1 √ 2 n 1 n 2 1 √ 2 n 2 2 1 √ 2 n 2 1 1 √ 2 n 1 n 2       . ( 71
)
After discretization, the linear system (61) with the displacement corrections (62) reduces to a standard Newton-type iteration:

K tan ∆u = -R u (k) , u (k+1) = u (k) + ∆u, (72)
where u (k) is the displacement field known from the previous (k-th) Newton-Raphson iteration,

K tan = Ω B T w C s B w dΩ + Ω γ β (hNB) T C I (hNB) dΩ, (73)
and

R u (k) = Ω B T w σ (k) dΩ + Ω γ β (hNB) T t w (k) dΩ - Ω N T u f dΩ - ∂Ω F N T u FdS, (74)
where

C I = ∂t (w) ∂w =   ∂t n ∂w n 0 0 0   , (75)
and where C s is the matrix form corresponding to the fourth-order elastoplastic consistent tangent operator C s in [START_REF] Simo | Computational Inelasticity[END_REF], which is determined by the classical elastic predictor and plastic corrector (return-mapping) algorithm outlined in [START_REF] De Souza Neto | Computational methods for plasticity: theory and applications[END_REF]. It should be noted that an accurate evaluation of the operator C s is crucial to guarantee the convergence of the Newton-Raphson iterative solution scheme. The analytical form of C s for a J 2 -plasticity yield function can be found in [START_REF] Ambati | Phase-field modeling of brittle and ductile fracture in shells with isogeometric nurbsbased solid-shell elements[END_REF]. The iterative update ( 72) is performed until convergence is achieved in the sense ∆u / u (k+1) -u (0) ≤tol.

Discretization of the phase field problem

We solve alternatively the damage problem and then the mechanical problem within a staggered procedure [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. The damage and damage gradient are approximated in one element by

d = N d d e , ∇d = B d d e , (76)
where N d and B d are matrices of damage shape function and of damage shape function derivatives, respectively, and d e denote nodal damage in one element. The discretization of damage problem (65) results into the following discrete system of equations:

K d d = F d ( 77
)
in which

K d = Ω H+ (1 -β) 2 ψ c N T d N d + (1 -β) 2 ψ c 2 d B T d B d dΩ (78)
and

F d = Ω N T d H dΩ, (79)
where H is given in (56).

Numerical implementation

In the present work, a staggered scheme is employed following [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], where at each load increment the displacement problem is solved for fixed damage field which is known from the previous time step. The damage problem is then solved with the new displacement field. The overall algorithm is illustrated in Algorithm 1.

Algorithm 1: Overall algorithm

Initialize u 0 , d 0 , ε p 0 , and α 0 with assumption of not plasticised and undamaged state. Compute the level-set function φ and interface damage β. Loop over load increments n for i = 1, . . . , n do Displacement problem Newton-Raphson iterative solution scheme k = 1, err = 1, tol = 10 -5 , u

i = u i-1 while err > tol do ε ep = B w u (k-1) i , w = hNBu (k-1) i (0)
Compute t (w) and C I from ( 44) and ( 75)

Return-mapping algorithm Given ε ep , ε p i-1 , α i-1 , d i-1 , Compute C s , ε p,(k) i , α (k) i from Algorithm 2
Compute K tan and R u (k-1) i from ( 73) and ( 74) Compute ∆u from (72) Update u

(k) i = u (k-1) i + ∆u u i = u (k) i , ε p i = ε p,(k) i , α i = α (k) i err = ∆u / u i -u i-1 , k = k +
i = i + 1 end
The flowchart for return-mapping algorithm on one Gauss point is provided in Algorithm 2.

Algorithm 2: Return-mapping algorithm

Input:

ε ep , ε p i-1 , α i-1 , d i-1 Output: C s , ε p,(k) i , α (k) i Elastic prediction ε e trial = ε ep -ε p i-1
Compute σ trial and F p σ trial , α i-1 , d i-1 from ( 48) and ( 49)

if F p ≤ 0 then ε p,(k) i = ε p i-1 , α (k) i = α i-1 else Plastic correction Compute ∆α from F p (σ, α i-1 + ∆α, d i-1 ) = 0 ∆ε p = ∆α 3 2 s trial s trial Update the variables ε p,(k) i = ε p i-1 + ∆ε p α (k) i = α i-1 + ∆α end ε e = ε ep -ε p,(k) i
Compute σ from (48) Compute C s with the analytical form in [START_REF] Ambati | Phase-field modeling of brittle and ductile fracture in shells with isogeometric nurbsbased solid-shell elements[END_REF]

Numerical examples

In this section, all numerical computations are performed within the finite element framework and assuming plain strain conditions. Both damage and displacement fields are discretized with 4-node bilinear elements. The material properties are shown in Table 3. Displacement controlled conditions are always assumed. We adopt the staggered solution strategy presented in 4.4. For the sake of simplicity, we recall that both length scale parameters d and β are assumed to be equal, i.e d = β = . (

) (b1) (c1) (a2) (b2) (c2) 0.6GPa c   0.8GPa c   U 0 15mm = . a1
U 0 17mm = .

U 0 197mm = .

U 0 15mm = .

U 0 17mm = .

U 0 197mm = .

U 0 2mm = .

U 0 24mm = .

U 0 316mm = .

U 0 2mm = .

U 0 24mm = .

U 0 316mm = .

Figure 5: Semicircular notched specimen: a1-c1 crack phase field and a2-c2 equivalent plastic strain field at three different prescribed displacements (see Fig. 7).

Semicircular notched specimen

In this example, we fist validate the convergence of the elastoplastic phase field with respect to mesh refinement, in absence of interfaces. Then, we analyze the influence of the critical fracture stress σ c . For this purpose, a specimen with a semicircular notch, as described in Fig. 4(a), is considered. The boundary conditions are as follows: on the lower end (y = 0), the y-displacement are fixed, while the x-displacement are free and the node (x = 0, y = 0) is fixed. On the upper end, the x-displacement are free, while the y-displacement are prescribed to an increasing value of U with ∆U = 0.001 mm during the simulation. The material parameters are those of Material I in Table 3 and = 0.5 mm. The spatial discretization of the model comprises 8953 4-node quadrilateral elements, with refinement in the central region where the crack is expected to propagate (see Fig. 4(b)). Fig. 5 shows the evolution of the crack phase field and the equivalent plastic strain field at three different prescribed displacements. As can be observed from Fig. 5 (a1)-(c1), the crack propagates horizontally towards the right-end boundary. The equivalent plastic strain as shown in Fig. 5 (a2)-( c2) is maximum at the notch root and localization branches form near the notch at an angle of about 45 • . The same simulation has been repeated on two other finite element meshes of 6834 and 11, 325 elements (coarse and fine mesh) to study the convergence with respect to the mesh size. Results are provided in Fig. 6, demonstrating the convergence of the method as the mesh is refined. In the following, the medium mesh has been used to limit the computational costs. coarse mesh medium mesh fine mesh

(b1) (c1) (a2) (b2) (c2) 0.6 GPa c   0.8 GPa c   1.0 GPa c   1.2 GPa c   (a1) (b1) (c1) (a2) (b2) (c2) (a1)
Figure 8: Semicircular notched specimen: effect of the critical fracture stress σc on the fracture process; a1-c1: crack phase field; a2-c2: equivalent plastic strain field at three different prescribed displacements (see Fig. 7).

Fig. 7 and Fig. 8 show the effect of the critical fracture stress σ c on results in terms of loaddisplacement curve, crack path and the equivalent plastic strain field. As can be observed, with the increase of the critical fracture stress σ c , the prescribed displacement corresponding to the point (maximum load) of initiation of the fracture rises, thus leading to an increase of the equivalent plastic strain field. This observation has been reported in many other works [START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. part ii. coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[END_REF][START_REF] Borden | A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects[END_REF][START_REF] Badnava | A phase field model for rate-dependent ductile fracture[END_REF]. Here we define the critical fracture stress σ c as a material parameter which controls the initiation of the fracture.

Fatigue cracking

In this next example, we include a cohesive interface and validate the convergence of the model with respect to the mesh size. Additionally we investigate the capability of the framework to handle fatigue cracking under cyclic loading. A square domain is considered. The domain contains a cohesive interface, whose geometry is depicted in Fig. 9(a). In order to provide a good balance between simulation accuracy and computational costs, a finely refined mesh is used in the region close to the cohesive interface, with finite element size: h e = 0.05 mm (see Fig. 9(b)). The material parameters are those of Material I in Table 2, in addition to = 0.1 mm, h = h e = 0.05 mm in [START_REF] Hedjazi | Application of the discrete element method to crack propagation and crack branching in a vitreous dense biopolymer material[END_REF], fracture strength and toughness t u = 10 MPa and g I c = 0.1 N/mm, respectively. In a first case, we validate the convergence with respect to the mesh refinement. For this purpose, two other finite element meshes have been used: a coarse mesh with 3500 elements and a fine mesh with 10, 028 elements. A displacement U (∆U = 0.001 mm) whose evolution is described in Fig. 10(a) is prescribed on the upper end, as depicted in Fig. 9(a). Results for 3 meshes are provided in Fig. 10(b), which show the convergence of the method with mesh refinement. In the following case, the medium mesh has been used to limit the computational costs.

In this next case, the evolution of U (∆U = 0.001 mm) is described in Fig. 11(a), involving multiple cycles. The corresponding load-displacement curve is depicted in Fig. 11(b). As shown in this figure, the irreversible strains upon unloading are well produced by this model. The evolution of the equivalent plastic strain in the element near node A (which is located in the center of the domain) and damage on node A are depicted in Fig. 12(a). The computed displacement jump along y in the element near node A is shown in Fig. 12(b). The normal traction force in the element near node A, with respect to the computed displacement jump and the prescribed displacement, are depicted in Fig. 12(c) and (d), respectively.

Traction test of a microstructure with a single fiber

In this example, we investigate a microstructure involving one fiber in an elastoplastic matrix with damageable interface. The importance of interfacial damage and the influence of its parameters on the results will be studied. The geometry and boundary conditions are depicted in Fig. 13(a). The finite element model (medium mesh with 24,688 elements) is shown in Fig. 13(b). External loading is applied by displacement control through a series of load increments with a fixed step value ∆U = 5 × 10 -4 mm. The material parameters for the matrix and fiber are those of Material I and Material II in Table 2, in addition to = 0.2 mm, h = 0.1 mm, g I c = 0.1 N/mm and t u = 0.3 MPa (model I). An illustration of the interface indicator β (x) and associated level-set φ (x) for the present microstructure is depicted in Fig. 14. In this composite example, we also validate the mesh convergence with two other finite element meshes: a coarse mesh with 13, 268 elements and a fine mesh with 34, 802 elements. The load-displacement curves are shown in Fig. 15(a). The evolution of the crack phase field is presented in Fig. 16. From these results it is observed that the crack nucleates from the interface and then kinks into the matrix when reaching a certain point. This leads to the final predicted semi-debonding angle of 65.10 • which is similar to the quasi brittle cases [START_REF] Nguyen | Discontinuous galerkin/extrinsic cohesive zone modeling: Implementation caveats and applications in computational fracture mechanics[END_REF][START_REF] Zhang | Modelling progressive failure in multi-phase materials using a phase field method[END_REF]. The difference compared with similar cases but involving a quasi brittle matrix [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF][START_REF] Zhang | Modelling progressive failure in multi-phase materials using a phase field method[END_REF] is that the initial interface debonding position in this model is about 45 • direction in the fiber surface (see Fig. 16(a)).

To evaluate the capability of the method to describe interfacial damage and investigate the effect of the interfacial cohesive model, another three simulations are performed: in the first one, called "model II": g I c = 0.1 N/mm and t u = 1 MPa, in the second one, called "model III": g I c = 0.2 N/mm and t u = 1 MPa and in the third one, called "no interface damage "model, only takes into account damage of the bulk (basic phase field method). The comparison of load-displacement curves for all four simulations is depicted in Fig. 15(b). The crack phase fields are shown in Fig. 17. For the three cohesive models, it can be observed that the cracks are similar which nucleate from the interface and then propagate into the matrix, while for the no interface damage model, the cracks nucleate and propagate within the matrix around the boundary after a very long plastic hardening process. We can then observe that the different interfacial damage models do not have a significant influence on the response, but when removed, the response of the sample is drastically changed. This shows the crucial importance of incorporating such interfacial damage model to properly describe the overall behavior of damaged microstructures.

M1

U 0 094mm = .

M2

U 0 08mm = .

M2

U 0 104mm = .

M3

U 0 08mm = .

M3

U 0 127mm = .

M4

U 0 098mm = .

M4

U 0 12mm = .

M5

U 0 094mm = .

M5

U 0 134mm = .

M6

U 0 078mm = .

M6

U 0 102mm = .

M7

U 0 097mm = .

M7

U 0 125mm = .

M8

U 0 113mm = .

M8

U 0 16mm = .

M9

U 0 103mm = .

M9

U 0 139mm = .

Traction test of a microstructure with randomly distributed fibers

A microstructure with elastoplastic matrix and damageable interfaces, containing 9 randomly distributed circular fibers (radius r = 1 mm) is considered, whose geometry and boundary conditions are illustrated in Fig. 18(a). External loading is applied by displacement control through a serious of load increments with a fixed step value ∆U = 2 × 10 -4 mm. The material parameters for the matrix and fiber are those of Material I and Material II in Table 2, in addition to = 0.2 mm and h = 0.1 mm. The cohesive model parameters are g I c = 0.1 N/mm and t u = 0.3 MPa. The distribution of interface phase field is shown in Fig. 18(b). It is observed that the diffusive interface concentrates in the nearby area of the interface, and quickly attenuates. The microcracking initiation and final distribution for 9 realizations of microstructures is depicted in Fig. 19. The corresponding load-displacement curves are depicted in Fig. 20. In the different figures of Fig. 19, we can note that in each case, the cracks nucleate at the interface between the matrix and the fibers, and then kink into the matrix. For most cases, a crack path is created between the fibers passing through the interfaces and then leading to the rupture of the specimen. This example illustrates well the capability and robustness of the technique to handle brittle crack propagation from interfaces to the bulk in elasto-plastic composites. In Fig. 20, it can be observed that all curves match with each others in the elastic and plastic stages. However, a large dispersion of the individual results in the damaged stage is observed. This shows a strong sensitivity of the overall brittle response to the local distribution of fibers in microstructure as well as the capability of the proposed method to capture these effect.

Traction test of a sample whose complex microstructure is obtained from microtomography

In this example, we investigate the capabilities of the method to simulate microcrack propagation in complex microstructures such as obtained by experimental imaging techniques, like X-ray microtomography. The geometry and mechanical boundary conditions are illustrated in Fig. 21 N/mm and t u = 0.3 MPa. External loading is applied by displacement control through 1800 load increments with a fixed step value ∆U = 1 × 10 -4 mm. The distribution of interface phase field is shown in Fig. 21(b). Comparing with the single fiber system and the randomly distributed fibers system, it is much more complicated due to material heterogeneities. The microcracking evolution in the domain is presented in Fig. 22. We can observe that with the increase of the external loading, cracks are initiated in the form of interface debonding and then migrate into the matrix in the form of matrix cracks. Subsequently, these interface cracks and matrix cracks are interconnected and then lead to the final failure of the microstructure. The crack paths are very complex and show the potential of the method to describe microcracking with interfacial damage in very complex, heterogeneous microstructures. The corresponding load-displacement curve is depicted in Fig. 23. It can be seen that the load-displacement curve shows several abrupt stress drops. These stress drops result mainly from the initiation of the microcracking in the complex, heterogeneous microstructure.

Conclusion

In this work an extension of the phase field to interfacial damage in elastoplastic composites has been proposed. In the matrix, the elastoplastic phase field extension as proposed in [START_REF] Duda | A phase-field/gradient damage model for brittle fracture in elastic-plastic solids[END_REF] has been used. To introduce interfacial damage, the energetic formulation has been modified by adding: (i) a strain density depending of the displacement jump related to matrix/inclusions decohesion; (ii) a modified description of the total energy involving a regularized approximation of the singular strain part along the interfaces and (iii) the use of a regularized description of interfaces through diffuse weighting functions which are introduced in the energetic formulation to differentiate the bulk and interfacial damage mechanisms. In that manner, different damage mechanisms can be associated with interfaces as compared to the matrix cracking mechanisms. Application of the variational principle provides the different equations to be solved in a straightforward manner. As a result, the technique allows simulating initiation, propagation and interactions between both fracture and interfacial cracks in elastoplastic composites without special treatment and using standard finite elements. Several numerical examples involving complex microstructures (composites, concrete micro tomography images) have shown the capability of the method to handle complex micro cracks interactions for arbitrary complexity of the microstructures, and convergent solutions with respect to the mesh size. One future perspective of this work would be the experimental identification of the internal length parameters d and β , i.e. following an inverse approach as described in the quasi-brittle context in [START_REF] Nguyen | Initiation and propagation of complex 3d networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF]. The present model could be extended to plasticity gradient models such as in [START_REF] Miehe | Phase field modeling of ductile fracture at finite strains: A variational gradient-extended plasticity-damage theory[END_REF] to involve full damage and plasticity coupling.

Figure 1 :

 1 Figure 1: Diffused approximation of cracks and interfaces: (a) a medium containing sharp cracks and interface; (b) diffused approximation of cracks; (c) diffuse approximation of interfaces.

Figure 2 :

 2 Figure 2: Approximation of the displacement jump across an interface Γ I at a point x, with n I the normal unit vector to Γ I .

Figure 3 :

 3 Figure 3: Cohesive model for the interfaces: evolution of the normal traction component t n with respect to the approximated jump normal displacement component w n .

Figure 4 :

 4 Figure 4: Semicircular notched specimen: (a) geometry and boundary conditions; (b) finite element model.

Figure 6 :

 6 Figure 6: Load-displacement curve of a semicircular notched specimen showing convergence with respect to the finite element mesh size (6834, 8953 and 11, 325 elements).

Figure 7 :

 7 Figure 7: Load-displacement curve of a semicircular notched specimen: sensitivity with respect to the critical fracture stress σc.

Figure 9 :

 9 Figure 9: Definition of the fatigue crack test: (a) geometry and boundary conditions; (b) finite element model (medium mesh).

Figure 10 :

 10 Figure 10: Fatigue crack test: (a) evolution of the load; (b) effect of the mesh size (3500, 6800 and 10, 028 elements) on the load-displacement curve.

Figure 11 :

 11 Figure 11: Fatigue crack test: (a) evolution of the load; (b) corresponding load-displacement curve.

Figure 12 :

 12 Figure 12: Fatigue crack test: (a) equivalent plastic strain and damage with respect to the prescribed displacement; (b) displacement jump along y with respect to the load number; (c) cohesive traction with respect to the displacement jump; (d) cohesive traction with respect to the prescribed displacement.

Figure 13 : 14 :

 1314 Figure 13: Traction test of a microstructure with a single fiber: (a) geometry and boundary conditions; (b) finite element model.

Figure 15 :

 15 Figure 15: Traction test of a microstructure with a single fiber: (a) load-displacement curves for three finite element meshes (13, 268, 24, 688 and 34, 802 elements); (b) comparison of load-displacement curves for different interfacial damage models.

16 : 17 :Figure 18 :

 161718 Figure 18: Traction test of a microstructure with randomly distributed fibers: (a) geometry of the domain and boundary conditions; (b) smeared interfaces represented by the interface phase field β (x).

Figure 19 :Figure 20 :

 1920 Figure 19: Evolution of crack networks in random elastoplastic microstructures subjected to traction. For each realization, the crack phase field is depicted at early and final stages of propagation. The microstructures are subjected to displacement traction U .

  (a). The material parameters for the matrix and inclusion are those of Material I and Material II in Table2, in addition to = 0.1 mm and h = 0.1 mm. The cohesive model parameters are g I c = 0.1

Figure 21 :

 21 Figure 21: Traction test of a sample whose geometry is obtained from microtomography: (a) geometry of the microstructure and boundary conditions; (b) smeared interface represented by the interface phase field β (x).

Figure 22 :Figure 23 :

 2223 Figure 22: Traction test of an elastoplastic sample whose microstructure is obtained from microtomography. Evolution of the crack phase field for traction displacements U = 0.09 mm, U = 0.11 mm, U = 0.13 mm and U = 0.18 mm.

Table 1 :

 1 Governing equations of the elastoplastic phase field model with interfacial damage.

Table 2 :

 2 Problems to be solved in the staggered procedure

  1 end Damage problem Compute H with u i and ε p i from (56) Compute K d and F d from (78) and (79) Compute the damage field d i from (77)

Table 3 :

 3 Material parameters used in the numerical simulations, from[START_REF] Mediavilla | Discrete crack modelling of ductile fracture driven by non-local softening plasticity[END_REF]

	Name	Symbol Material I Material II
	Shear modulus	κ	27.28 GPa 70.3 GPa
	Bulk modulus	µ	71.66 GPa 136.5 GPa
	Yield stress	σ Y	0.345 GPa 0.443 GPa
	Hardening modulus	H	0.25 GPa	0.3 GPa
	Critical fracture stress	σ c	1 GPa	2 GPa

Acknowledgements

This work has benefited from the project "Multiscale Modeling and Experimental Investigation of damage in composite components obtained by additive manufacturing process" (MMELED), ANR-16-CE08-0044-04.