Thanh-Tung Nguyen

Julien Yvonnet

Danièle Waldmann

Qi-Chang He

Phase field modeling of interfacial damage in heterogeneous media with stiff and soft interphases

Keywords: Fracture, Phase field model, Interfacial cracking, Coherent interface, Spring-layer interface

A new interfacial cracking model in the phase field framework is proposed. The developed method is able to capture the effects of both stiff and soft interphases on the fracture behavior of heterogeneous materials. A dimensional-reduced model based on a rigorous asymptotic analysis is adapted to derive the zero thickness imperfect interface models from an original configuration containing thin interphase.

Then, the energetic approach is used to describe the material degradation both on the interface and in bulk within the context of the phase field method for fracture. This technique allows to effectively model the competition between the interface and bulk cracking. Furthermore, a unilateral contact condition is also proposed to ensure the physical crack propagation patterns in the case of spring imperfect interface. The complex cracking phenomena on interfaces such as initiation, delamination, coalescence, deflection are successfully predicted by the present method. The numerical implementation using a staggered algorithm provides an extremely robust approach. The performance of the proposed framework is demonstrated through numerical examples involving complex microcracking of both stiff and soft interfaces in complex microstructures.

Introduction

Design of materials and structure often requires the use of multi-components at the material level and/or at the structural level. In these conditions, interfaces or/and interphases are most often the location for crack initiation and thus play a major role in the physical/mechanical properties of these materials or structures [START_REF] Chawla | Composite materials: science and engineering[END_REF]. Consequently, an accurate modeling and simulation framework for failure of interfaces and interaction with bulk fracture is central to predicting the strength resistance of the Email address: thanh-tung.nguyen@uni.lu (Thanh-Tung Nguyen ) final products. Study of interface failure has been the topic of an intense research in the last decades, on both theoretical and numerical aspects.

Numerous analytical works on interfacial failure based on linear fracture mechanics have been carried out, including determination of near crack tip fields or intensity factors in bilayer materials and associated solutions [START_REF] Hutchinson | Mixed mode cracking in layered materials[END_REF][START_REF] He | The mechanics of crack growth in layered materials[END_REF], or determination of interface energy [START_REF] Evans | The fracture energy of bimaterial interfaces[END_REF]. Analytical models of bulk cracks interacting with bi-material interfaces have been studied in e.g. [START_REF] He | Crack deflection at an interface between dissimilar elastic materials[END_REF][START_REF] Leguillon | Interface debonding ahead of a primary crack[END_REF][START_REF] Martin | A revisited criterion for crack deflection at an interface in a brittle bimaterial[END_REF].

The idea of cohesive zone model has been introduced by Dugdale [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and Barenblatt [9] to address the stress singularity at a crack tip. With this concept, the stress ahead of the crack tip is bounded, and a traction-separation law is adopted to describe the failure process. Xu and Needleman [START_REF] Xu | Void nucleation by inclusion debonding in a crystal matrix[END_REF] used such model to investigate numerically the nonlinear decohesion of an matrix/inclusion interface and its interaction with the elastoplastic behavior of the matrix. Jin et al. [START_REF] Jin | Cohesive fracture modeling of elastic-plastic crack growth in functionally graded materials[END_REF] introduced cohesive models for functionally graded materials. Different other cohesive models have been proposed (see e.g. [START_REF] Camanho | Numerical simulation of mixed-mode progressive delamination in composite materials[END_REF]).

In [START_REF] Allix | Interlaminar interface modelling for the prediction of delamination[END_REF][START_REF] Allix | Modeling and simulation of crack propagation in mixed-modes interlaminar fracture specimens[END_REF], Allix et al. introduced a surface damage model with several internal variables to model the anisotropic interfacial damage in laminated composites. Recently, Esmaeili et al. [START_REF] Esmaeili | Coherent energetic interfaces accounting for in-plane degradation[END_REF] have proposed a damage model for treating coherent imperfect interfaces, i.e. not involving displacement jump at the interface but rather a jump of traction across interfaces and which can model thin interphase, as well as for non-coherent general imperfect interfaces in [START_REF] Esmaeili | Non-coherent energetic interfaces accounting for degradation[END_REF].

Numerically, crack propagation along interfaces has been widely investigated using cohesive elements [START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF][START_REF] Camacho | Computational modelling of impact damage in brittle materials[END_REF]. In the situation where the crack path is known (e.g. along the interfaces), cohesive elements provide a simple framework by debounding meshed interfaces where nodes are doubled (see [START_REF] Elices | The cohesive zone model: advantages, limitations and challenges[END_REF] for a critical review on the application of the relevant methods). Other techniques like XFEM [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] have been used to study local fields and stress intensity factors in bi materials in linear situations [START_REF] Nagashima | Stress intensity factor analysis of interface cracks using X-FEM[END_REF][START_REF] Liu | XFEM for direct evaluation of mixed mode SIFs in homogeneous and bi-materials[END_REF][START_REF] Sukumar | Partition of unity enrichment for bimaterial interface cracks[END_REF][START_REF] Huynh | The extended finite element method for fracture in composite materials[END_REF], where the displacement jump can be accurately reproduced through the enrichment of the displacement field with discontinuous functions, or by directly incorporating a cohesive model independently of the mesh within a nonlinear XFEM context [START_REF] Wang | Progressive delamination analysis of composite materials using xfem and a discrete damage zone model[END_REF]. However, XFEM methods can hardly handle triple junctions or complex interphases, such as those in polycristals. When interface cracks interact with bulk cracks, one appealing numerical scheme is the use of smeared description of discontinuities like in the context of the variational approach to fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF] (also widely called phase field method for fracture [START_REF] Kuhn | A continuum phase field model for 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] Karma | Phase-field model of mode III dynamic fracture[END_REF][START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF]). In this context, cracks can be initiated, propagate and merge without dependence on the mesh. In [START_REF] Verhoosel | A phase-field model for cohesive fracture[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], extensions of the phase field to interfacial damage have been proposed. In [START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF] an additional field has been introduced to model the displacement jump at the interface. 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], a different model has been proposed, which does not require an additional field and where both bulk and interface cracks are described by smeared discontinuities in the context of the phase field method, allowing complex interactions between bulk and interface cracks. A similar approach has been suggested in a recent work of Hansen-Dörr et al. [START_REF] Hansen-Dörr | Phase-field modelling of interface failure in brittle materials[END_REF] for simulating failure of perfect interfaces by distributing the adhesive interface over a finite width (smeared interface) and then assigning an interface fracture toughness. Numerical studies using the above mentioned numerical methods have been applied to studying delamination [START_REF] Turon | An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[END_REF], particle debonding in materials [START_REF] Li | Prediction of fracture toughness of ceramic composites as function of microstructure: I. numerical simulations[END_REF] interactions between interface failure and bulk cracks (also sometimes referred to as crack deflexion problems) e.g. in [START_REF] Parmigiani | The roles of toughness and cohesive strength on crack deflection at interfaces[END_REF][START_REF] Ponnusami | Cohesive-zone modelling of crack nucleation and propagation in particulate composites[END_REF], and in [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] where the phase field and cohesive elements are employed to model separately interfacial damage and bulk cracking.

In this work, we develop a phase field framework for dealing with damage of a larger class of imperfect interfaces, i.e. exhibiting either displacement jump or traction jump across an interface.

These different models are obtained by asymptotic analysis and can replace an interphase with finite thickness made of a third material by an imperfect interface of zero thickness with appropriate jump conditions. The main numerical advantage is to avoid meshing the interphase. The jumps conditions depend on the ratio of the stiffness of the interphase to those of the connected materials and lead to two well-known models in the literature: coherent imperfect interface for very stiff interphase, and linear spring-layer imperfect interface for very soft interphase [START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF][START_REF] Yvonnet | An XFEM/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites[END_REF][START_REF] Gu | Interfacial discontinuity relations for coupled multifield phenomena and their application to the modeling of thin interphases as imperfect interfaces[END_REF][START_REF] Zhu | Three-dimensional numerical modelling by XFEM of spring-layer imperfect curved interfaces with applications to linearly elastic composite materials[END_REF][START_REF] Javili | Thermomechanics of solids with lower-dimensional energetics: on the importance of surface, interface, and curve structures at the nanoscale. a unifying review[END_REF]. In this paper, we extend the framework proposed 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] by including these models within the phase field framework, so as to allow interaction between the damage in the thin interphase and the bulk cracks in an efficient and robust manner.

This paper is organized as follows. In the next section, a brief introduction to imperfect interface models is provided. The detailed description of the proposed framework is given in Section 3, where a variational approach in the context of the phase field is used to formulate the problem. In Section 4, numerical implementation details are provided. Finally, several numerical applications are presented in Section 5, demonstrating the performance of the proposed approach.

A physically-based isotropic elastic imperfect interface model

Geometric preliminaries

Let us consider a practical three-phase configuration composed of a matrix Ω (m) , an inclusion Ω (i) , and the interphase Ω (0) between them, as schematically sketched in Fig. 1(a). The matrix/interphase interface S (m) and the interphase/inclusion interface S (i) are assumed to be perfect. It has been demonstrated in many contributions [42-45, 47, 48] that the three-phase configuration mentioned above can be replaced by a two-phase configuration, where the interphase Ω (0) of thickness h is substituted by an imperfect interface of null thickness located at the middle surface Γ I of Ω (0) . Meanwhile, the material forming Ω (m) and Ω (i) separately extend to Γ I , see Fig. 1(b). The primary advantage of this approach is to avoid modeling a very thin interface, so as to provide an essential simplification to modeling such a heterogeneous medium numerically. A general method in differential geometry to describe a smooth interface Γ I (with arbitrary shape) is to define it as the zero level-set of a function φ(x) [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF], such that

         φ(x) > 0 for x ∈ Ω (i) φ(x) < 0 for x ∈ Ω/Ω (i) φ(x) = 0 for x ∈ Γ I (1)
Then, the unit normal vector filed on Γ I can be determined by

n Γ I (x) = ∇φ(x) ∇φ(x) , (2)
in which ∇(.) is the gradient operator and . denotes the Euclidian norm.

Next, the two orthogonal complementary projection tensors are introduced

N(x) = n(x) ⊗ n(x), T(x) = I -n(x) ⊗ n(x), (3)
where I is the second order unit tensor.

Geometrically, N(x) and T(x) describe the projections along the normal direction and tangential direction to Γ I , respectively. The introduction of N(x), and T(x) allows the following decomposition of any differentiable vector function w and any differentiable second-order tensor function W into the normal part (denoted by the n indice) and the tangential/surface part (denoted by the s indice), as follows

   w = w n + w s with w n = Nw, w s = Tw, ∇w = ∇ n w + ∇ s w with ∇ n w = ∇wN, ∇ s w = ∇wT, (4)
and

         W = W n + W s with W n = NWN, W s = TWT, div (W) = div n (W) + div s (W) with div n (W) = ∇ (W) : N,
and div s (W) = ∇ (W) : T.

(5)

Coherent imperfect interface model and Linear spring-layer imperfect interface model

In this subsection, we will provide a summary of the constitutive relations for the imperfect interface model. Since the two-phase and three-phase configurations described in Fig. 1, hold the same mechanical phenomenon, the jumps of both the displacement and traction vectors across the imperfect interface Γ I in the two-phase configuration must be, to within a admissible error, equal to those across the interphase Ω (0) in the three-phase configuration. The expressions for the discontinuity of mechanical fields can be derived by applying the asymptotic expansion technique to the physically realistic three-phase configuration [START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF].

This work is dedicated to studying the situation, where the interphase is either very stiff (denoted by Γ I,α ) or very compliant (denoted by Γ I,β ), so that the two well-known particular imperfect interface models are applied: coherent imperfect interface model (CI) and linear spring-layer imperfect interface model (SI). The case of the general imperfect interface, e.g., in [START_REF] Gu | Interfacial discontinuity relations for coupled multifield phenomena and their application to the modeling of thin interphases as imperfect interfaces[END_REF][START_REF] Javili | General imperfect interfaces[END_REF], is not considered here. In addition, we assume that the material forming each of the different phases is individually homogeneous, isotropic, and governed by the Hooke's law

σ (k) = C (k) : ε (k) in Ω (k) , with k = i, m, 0, (6)
where σ (k) , ε (k) are the second-order stress and strain tensors in the phases i, m, 0, respectively.

With [[.]] describes the discontinuity of the mechanical fields, the constitutive relations for the two models (CI) and (SI) are reported according to [START_REF] Yvonnet | An XFEM/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites[END_REF][START_REF] Gu | Interfacial discontinuity relations for coupled multifield phenomena and their application to the modeling of thin interphases as imperfect interfaces[END_REF][START_REF] Zhu | Three-dimensional numerical modelling by XFEM of spring-layer imperfect curved interfaces with applications to linearly elastic composite materials[END_REF] as follows Very stiff interphase: Coherent imperfect interface model (CI)

(i) The displacement jump [[u]] across Γ I,α is null [[u]] = 0 on Γ I,α (7)
(ii) The traction vector across Γ I,α is discontinuous

[[t]] = σ (m) -σ (i) n on Γ I,α (8)
(iii) The strain tensor is discontinuous but the tangent or surface strains ε s are continuous

[[ε s ]] = [[TεT]] = 0 on Γ I,α (9)
Very compliant interphase: Linear spring-layer imperfect interface model (SI)

(i) The displacement jump [[u]] across Γ I,β is non-null [[u]] = u (m) -u (i) on Γ I,β (10)
(ii) The traction vector across Γ I,β is continuous

t = σ (m) n = σ (i) n on Γ I,β (11)
(iii) The traction vector t acting on Γ I,β is proportional to the vector of displacement jump across Γ I,β , with C (si) s being the tangential interface stiffness, i.e.

t = C (si) s [[u]] on Γ I,β (12)

Regularized framework of smeared crack and smeared interface

In the studies of Yvonnet et al. [START_REF] Yvonnet | An XFEM/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites[END_REF] and Zhu et al. [START_REF] Zhu | Three-dimensional numerical modelling by XFEM of spring-layer imperfect curved interfaces with applications to linearly elastic composite materials[END_REF], the XFEM/level set approach has been successfully applied to study the surface/interface effects in micro and nanocomposite materials.

However, using such method to modeling interfacial cracking has some shortcomings when simulating crack initiation, and describing complex crack coalescence/interaction (especially in 3D). To overcome these drawbacks, the computational approach of the phase field model based on the smeared description of the interface/crack is here adopted [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 | Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials[END_REF], see Fig. 1(c).

The basic idea of the phase field model in the regularized framework is to approximate the sharp crack by a smeared surface Γ d , defined by a scalar parameter d(x), x ∈ Ω, taking a unit value on the crack surface and vanishing away from it [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] Γ Above, γ d (d, ∇d) is the crack density function per unit volume, which has been chosen following Pham et al [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF], as follows

d (d) = Ω γ d (d, ∇d) dΩ. ( 13
γ d (d, ∇d) = 3 8 d d + 3 d 8 ∇d • ∇d, (14)
where d is a regularization parameter describing the actual thickness of the smeared crack.

Then, the phase field d(x) can be determined by solving the variational problem

d(x) = Arg inf d∈S d Γ d (d) , with S d = {d|d(x) = 1 ∀x ∈ Γ d } . (15)
The smeared interfaces are defined in the same manner with the introduction of a field α(x) for stiff interface, and a field β(x) for compliant interface (also taking a unit value on the interface and vanishing away from it). Following our previous study [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], these two fields can be obtained by solving the variational problem defined by

α(x) = Arg inf α∈Sα Γ α (α) and β(x) = Arg inf β∈S β Γ β (β) , (16)
where S α = {α|α(x) = 1 ∀x ∈ Γ I,α }, and

S β = {β|β(x) = 1 ∀x ∈ Γ I,β }.
In [START_REF] Allix | Modeling and simulation of crack propagation in mixed-modes interlaminar fracture specimens[END_REF] the total interface length Γ α,β = Ω γ α,β dΩ is defined as a volume integration of the quadratic interface functional density γ (.) , with

γ (.) = 1 2 (.) (.) 2 + (.) 2 ∇(.) • ∇(.). (17)
Here, α,β can be understood as the internal lengths of the smeared interfaces. The illustration of the smeared description of a cracked heterogeneous solid containing both stiff and compliant interfaces is provided in Fig. 2.

Energetic approach

In this section, the formulation of the phase field model mixing the bulk energy and interface energy is constructed. A new failure model accounting for the unilateral contact condition in both bulk and interface cracking is then proposed. Finally, the phase field crack evolution law is derived in a thermodynamically consistent way.

Definition of the total energy

Total energy of the cracked heterogeneous solid (as described in section 2) in a standard framework of sharp discontinuity description is given by

E = Ω/Γ d ,Γ I ψ e dΩ + Γ d ψ d dΓ + Γ I,α ψ ci dΓ + Γ I,β ψ si dΓ, (18)
where ψ e is the elastic energy density; ψ d is the surface (fracture) energy density; ψ ci and ψ si describe the interfacial energy densities of stiff interface, and compliant interface, respectively.

In the regularized description: (i) the fracture energy is regularized by the crack density function γ d (d, ∇d), and the strain energy is replaced by the energy of a damageable material ψ e (ε e , d); (ii) the energies of stiff interface and soft interface are regularized by γ α (α, ∇α) and γ β (β, ∇β), respectively.

Furthermore, in the case of (SI) model, (iii) we substitute [[u]] by a smeared displacement jump v(x)

(see e.g. [START_REF] Verhoosel | A phase-field model for cohesive fracture[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][START_REF] Nguyen | Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials[END_REF] for more details). Then, the infinitesimal strain tensor ε in this framework can be decomposed into a part related to the bulk ε e and a part induced by the smoothed jump at the interfaces ε, implying ε = ε e + ε. The aforementioned energy functional can be now rewritten as

E = Ω ψ e (ε e , d) h α,β dΩ + Ω g c γ d (d, ∇d) dΩ + Ω ψ ci (ε s , d)γ α (α, ∇α) dΩ + Ω ψ si (v, d)γ β (β, ∇β) dΩ, (19)
in which, h α,β = (1 -α)(1 -β) is introduced to distinguish the elastic energy of bulk phases from the one interfaces; g c here denotes the fracture resistance for both bulk material, and stiff/compliant interfaces; the tangential strain ε s is expressed by

1 ε s = Tε e T = M t ε e , with M t = T⊗T. ( 20
)
The definition of the model in Eq. 19 ensures that for α(x), β(x) → 0 (away from the interface) then γ α,β → 0, h α,β → 1 and ε e → ε, we recover the regularized energy functional for brittle fracture without interfacial effects [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[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] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF]. For the case of α(x), β(x) → 1 (the interfacial zones) the coherent/spring interfacial damage model is then obtained.

For the later purpose, we identify

ψ = ψ e (ε e , d) h α,β + ψ ci (ε s , d)γ α + ψ si (v, d)γ β + ψ d (d, ∇d), (21)
as the total strain energy density.

Unilateral contact formulation for the interfacial/bulk cracking

In the following, we will construct the failure model for both bulk and interfacial crackings.

The unilateral contact condition is adapted to distinguish the different fracture behaviors in tension/compression for the bulk cracking and in normal/tangential direction for the interfacial cracking.

Bulk cracking

The model proposed by Miehe et al [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] with the assumption that damage induced by traction only is here used for the bulk cracking. Herein, the density of elastic energy for the isotropic materials can be written as

ψ e (ε e , d) = g(d) ψ e+ (ε e ) + ψ e-(ε e ) , (22)
where g(d) = (1 -d) 2 + (with 1) is the degradation function, satisfying g(0) ≈ 1, g(1) ≈ 0 and g (1) = 0, see [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]).

The positive and negative parts of the strain energy ψ e± (ε e ) in ( 22) are expressed as

ψ e± (ε e ) = λ 2 tr ε e ± 2 + µ tr ε e± 2 , ( 23
)
with ε e+ and ε e-are, respectively the extensive and compressive modes of the elastic strain tensor ε e = ε e+ + ε e-, and its derivatives with respect to the elastic strain defines two projection tensor

P ± (ε e ) = ∂ ε e ε e± (ε e ) . (24)
Basically, P ± (ε e ) can be determined following the approach proposed by Miehe [START_REF] Miehe | Comparison of two algorithms for the computation of fourth-order isotropic tensor functions[END_REF]. The curious readers are also referred to the ref [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF] for more practical information of the numerical calculation.

Interfacial cracking: coherent imperfect interface model

The coherent imperfect interface is based on the hypothesis that the interface is endowed with surface (tangential) stresses/strains which are related by the two-dimensional (2D) Hooke law [START_REF] Yvonnet | An XFEM/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites[END_REF]. In this contribution, we assume that the interfacial cracking is created by the whole tangential strains, implying

ψ ci = g(d)ψ ci e = 1 2 ε s : C ci s (d) : ε s , (25)
where the tangential interface stiffness C ci s (d) is formulated by

C ci s (d) = g(d) λ ci s T ⊗ T + 2µ ci s T⊗T . ( 26
)
The Lamé's constants characterizing the interface, λ ci s , and µ ci s in Eq. ( 26) are given following [START_REF] Yvonnet | An XFEM/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites[END_REF]

λ ci s = 2µ (0) λ (0) + 2µ (0) hλ (0) , µ ci s = hµ (0) , (27)
with λ (0) , µ (0) are Lamé's parameters of the investigated interphase, and h denotes the interphase thickness.

Interfacial cracking: spring imperfect interface model

The spring imperfect interface model represents the displacement jump across the smeared interface.

Hence, the required energy for the interfacial decohesion is here described as a function of this displacement jump. By introducing a decomposition of normal displacement jump into negative part and positive part

v n = v + n + v - n , with v + n = v n + , and v - n = v n - 2 .
Then, the interfacial cracking is assumed to be created by the tangential displacement jump v t and the positive part of normal displacement jump v + n . This leads to a decomposition expressed by

v = v + n + v - n + v t , (28)
In the proposed framework, ψ si is postulated as

ψ si = g(d) ψ si n+ + ψ si t + ψ si n-= 1 2 g(d) v + n • C si s • v + n + v t • C si s • v t + 1 2 v - n • C si s • v - n , (29)
with the tangential interface stiffness C si s is defined following [START_REF] Zhu | Three-dimensional numerical modelling by XFEM of spring-layer imperfect curved interfaces with applications to linearly elastic composite materials[END_REF] as

C si s = λ si s n ⊗ n + µ si s (I -n ⊗ n) . (30)
The Lamé's constants characterizing the interface, λ si s , and µ si s , can be determined through the following formulation

λ si s = 1 h λ (0) + 2µ (0) , µ si s = 1 h µ (0) , (31)
where λ (0) , µ (0) are Lamé's parameters of the investigated interphase, and h is the interphase thickness.

Furthermore, by introducing projection tensors Q ± for the decomposition of the normal displacement jump, with v ± n = Q ± v n , we assume the following definitions

Q + =      sign(v i ) 0 0 0 sign(v j ) 0 0 0 sign(v k )      and Q -= I -Q + , (32)
where sign(.) = 1 if (.) ≥ 0, and, sign(.) = 0 if (.) < 0.

It implies the general expression for the failure model of the spring imperfect interface as

ψ si (v, d) = 1 2 v• C si s (d)• v, (33)
in which, the formulation of the tangential interface stiffness accounting for damage C si s (d) is given by (see Appendix A)

C si s (d) = g(d) N T Q + T C si s Q + N + T T C si s T + N T Q -T C si s Q -N . ( 34
)
Remark The physical choice of which part will enable the interface degradation for the tangential displacement jump is still being discussed. In general, its negative part can be related to either the compression of interphase or the shear of interphase. The first situation does not allow crack growth, while the second one does. For the sake of simplicity, the present model assumes that the whole tangential displacement jump will contribute to the crack creation. This assumption is linked to the fact that the tangential displacement jump provided by interphase in compression mode is less common than by interphase in shear mode. A new formulation able to distinguish these two situations has to be pursued in future studies.

Phase field problem

Next, we will formulate a crack phase field evolution law that can guarantee the irreversibility of the process. Assuming isothermal process and without the external mircoforces, the Clausius-Duhem inequality takes the reduced form

A d ḋ ≥ 0, (35)
where A d = -δψ δd = -∂ψ ∂d + ∇. ∂ψ ∂∇d is the derivative of ψ with respect to the phase field d. At this stage, a threshold function F (A d ) such that

F (A d ) ≤ 0, ( 36
)
is introduced. Assuming the principle of maximum dissipation requires the dissipation A d ḋ to be maximum under the above constraint Eq. [START_REF] Hansen-Dörr | Phase-field modelling of interface failure in brittle materials[END_REF]. Using the method of Lagrange multipliers, which defines the following Lagrangian

L = -A d ḋ + λF (A d ), (37)
yields the Kuhn-Tucker equations

∂L ∂A d = 0, λ ≥ 0, F ≤ 0, λF = 0. ( 38
)
The first equality in Eq. [START_REF] Li | Prediction of fracture toughness of ceramic composites as function of microstructure: I. numerical simulations[END_REF] gives

ḋ = λ ∂F (A d ) ∂A d . (39)
Without loss of generality, the threshold function F (A d ) is assumed to take the form

F (A d ) = A d .
From Eq. ( 39) and using the second inequality in Eq. ( 38), we obtain

ḋ = λ ≥ 0. ( 40
)
For ḋ > 0, and from Eqs. [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] Ponnusami | Cohesive-zone modelling of crack nucleation and propagation in particulate composites[END_REF] and the third equality in Eq. [START_REF] Li | Prediction of fracture toughness of ceramic composites as function of microstructure: I. numerical simulations[END_REF], which give F = 0, implying

F = - δψ δd = - ∂ψ e ∂d h α,β - ∂ψ ci ∂d γ α - ∂ψ si ∂d γ β -g c δγ d (d, ∇d) = 0, (41)
with

δγ(d, ∇d) = 3 8 - 3 4 ∆d. ( 42
)
From Eqs. ( 21), ( 22), ( 25), [START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF], and (41) we obtain

2(1 -d) ψ e+ h α,β + ψ ci e γ α + ψ si n+ + ψ si t γ β -g c δγ d (d, ∇d) = 0. ( 43
)
As discussed in the last section, the critical fracture energy g c is here used to denote the fracture resistance of bulk material g b c , the fracture resistance of stiff interface g ci c , the fracture resistance of soft interface in the normal direction g si,n c and in the tangential direction g si,t c . The corresponding internal lengths for different fracture modes can be chosen as b , ci , si,n , si,t , respectively. However, for sake of simplicity, we adopt a single internal length = b = ci = si,n = si,t . The local crack driving force H(x, t) is then introduced following Miehe et al. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] to describe the dependence on history, and to make loading-unloading possible, and reads as

H(x, t) = max τ ∈[0,t] g b c ψ e+ (x, τ ) h α,β (x) + g si,n c ψ si n+ (x, τ ) + g si,t c ψ si t (x, τ ) γ β (x) + g ci c ψ ci e (x, τ ) γ α (x) . (44)
Basically, H(x, t) contains the maximum reference energy, or a measure for the maximum tensile strain or interfacial opening/sliding obtained in the loading history, which provides a very transparent representation of the balance equation that governs the diffusive crack topology. For more detailed information, curious readers may refer to, e.g., [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF].

From the Eqs. ( 42), [START_REF] Yvonnet | An XFEM/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites[END_REF], and (44), the evolution law for phase field is finally expressed as

2(1 -d)H - 3 8 - 3 2 4 ∆d = 0 in Ω, (45)
with a homogeneous Neumann condition ∇d(x) • n = 0 on ∂Ω.

The irreversible character of the crack evolution law described in Eq. ( 45) can be verified as follows.

As H(x, t) ≥ 0, the Eq. ( 43) implies

δγ d (d, ∇d) ≥ 0. ( 46
)
We can thus check that due to Eq. ( 40),( 46) the variation of crack length

Γl = Ω δγ d (d, ∇d) ḋ dΩ ≥ 0, (47)
satisfying the criteria of irreversible evolution of cracks.

The associated weak form for the phase field problem is obtained by taking the variation in the phase field δd as

Ω 2Hdδd + 3 2 4 ∇d∇δd dΩ = Ω 2H - 3 8 + δd dΩ. ( 48
)

Displacement problem

The weak form of mechanical problem is obtained by taking an arbitrary variation in the displacement field δu and in the auxiliary displacement jump field δv, while the phase field is kept constant.

In the absence of body forces, it is expressed as

Ω h α,β σ (ε e , d) : ε e (δu) + γ α δε s : C ci s (d) : ε s + γ β δv • C si s (d) • v + ζ ∂v ∂x n • ∂δv ∂x n dΩ = ∂Ωt t • δu dS + ∂Γα γ α F • Tδu dl, ( 49
)
where t is the prescribed traction on ∂Ω t ; F = σ s n is an applied force on ∂Γ α for the case of open interfaces;

x n = (x -x Γ β ) • n Γ β and x Γ β = argmin y∈Γ ( y -x ), the positive constant ζ is introduced
to ensure that the auxiliary displacement jump field is constant in the normal direction [START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF].

The Cauchy stress σ(ε e , d) in the bulk material is defined as

σ(ε e , d) = ∂ψ e ∂ε e = g(d) λ tr ε e + 1 + 2µε e+ + λ tr ε e -1 + 2µε e-, (50)
herein the expression of the elastic strain ε e is given by (see [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF])

ε e = ∇ s u -n Γ β ⊗ s vγ β , (51)
with (∇ s u) ij = 1 2 (u i,j + u j,i ) and n Γ β ⊗ s v ij = 1 2 (n i v j + v i n j ). In the case of the crack curvature is assumed to be small, it yields the BVP

                               ∇ • σ(ε e , d) = 0 ∀x ∈ Ω, ∇ s • σ s = -σ (m) + σ (i) n Γα ∀x ∈ Γ α , γ β C si s • v -n Γ β • σ(ε e , d) = ζ ∂ 2 v ∂(x n ) 2 ∀x ∈ Γ β , ∂v ∂x n = 0 ∀x ∈ ∂Γ β , u = u ∀x ∈ ∂Ω u , n t • σ = t ∀x ∈ ∂Ω t , (52)

Discretization and numerical implementation

In the present work, a staggered procedure is employed, i.e., we alternatively solve the phase field problem and the mechanical problem.

FEM discretization of the phase field problem

The phase field and phase field gradient are approximated in one element by

d = N d d e and ∇d = B d d e , (53)
where N d and B d are vector of shape functions and matrix of shape functions derivatives, respectively;

Introducing this FEM discretization in [START_REF] Hashin | Thin interphase/imperfect interface in conduction[END_REF] results in a linear system of equations

K d d = F d , ( 55
)
where d is a column vector consisting of the unknown nodal values of d.

K d = Ω 2HN T d N d + 3 2 4 B T d B d dΩ (56)
and

F d = Ω 2H - 3 8 + N T d dΩ. ( 57
)

FEM discretization of the displacement problem

We introduce the following classical FEM discretization

u = N u u e , δu = N u δu e v = N v v e , δv = N v δv e , (58)
where u e , δu e ,v e , δv e denote nodal displacement components, nodal trial function displacement components, nodal displacement jump components, and nodal trial function jump components, respectively.

Similarly, the strains are defined by

[ε e ] (u) = B u u e , [ε e ] (δu) = B u δu e , [ε] (v) = B v v e , [ε] (δv) = B v δv e , (59)
and

∂v x n = G v v e . ( 60
)
Using the above discretization yields the following linear system of equations

K   u v   = F, (61)
where

K =   K uu K uv K vu K vv   , F =   F u F v   , (62)
with

K uu = Ω h α,β B T u C(d)B u dΩ + Ω γ α B T u M T C ci s (d)MB u dΩ, K uv = Ω -γ β B T u C(d)B v dΩ, K vu = Ω -γ β B T v C(d)B u dΩ, K vv = Ω γ 2 β B T v C(d)B v + γ β N T v C si s (d)N v + ζG T v G v dΩ,
and

F u = ∂Ωt N T u t dS + ∂Γα γ α N T u F • T dl, F v = 0,
where C ci s (d), C si s (d) are the matrices corresponding to the fourth-order tensor of tangential interface stiffness C ci s (d) in Eq. ( 26), and the second-order tensor of tangential interface stiffness C ci s (d) in Eq. [START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF], respectively.

By introducing the indicator R ± , with

R + (ε e ) = sign (tr ε e ) , and R -(ε e ) = 1 -R + (ε e ) , (63)
where the operator sign(.) is defined as the same as in Eq. [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF].

The matrix form C(d) of the elastic tensor in Voigt notation is expressed as follows

C(d) = g(d) λR + (ε e ) [1] T [1] + 2µP + + λR -(ε e ) [1] T [1] + 2µP -, (64)
in which P ± are the matrix form of the fourth-order tensor P ± defined in Eq. [START_REF] Huynh | The extended finite element method for fracture in composite materials[END_REF].

In general, to avoid the non-linearity due to the spectral decomposition, P ± and R ± are evaluated based on results of the previous loading step within the context of the incremental scheme, see, e.g., [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 | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF] for the practical details.

Numerical examples

Benchmark test of the SI model

In this first example, the proposed SI model will be investigated. We consider a benchmark problem consisting in a square domain which contains two phases and an initial crack, as depicted in Fig. 3.

The size of the square domain is

L × B = 1 × 1 mm 2 .

Initial crack

In te rf ac e The interphase is very compliant as compared to the matrix and inclusion phases. Hence the linear spring-layer imperfect interface model is here adopted. Then, the null thickness interface is regularized by the smeared interface γ β (β, ∇β). The material properties for bulk phases and interphase are given in the Table 1. The tangential Lamé's constants for the equivalent interface are determined by using the relation in Eq. ( 31), we obtain λ s = 87.5 GPa/mm, µ s = 25 GPa/mm for the case of an interphase with a thickness h = 1 µm. The fracture resistance of the interface is assumed to be different in the tangential and normal directions, see Table 1. Basically, the critical fracture energy in the proposed model can be characterized by using an inverse analysis, i.e determining g eq c so that the null-thickness interface model can hold the same fracture behaviors as the thin interphase one. Our analysis gives g eq c ≈ g c . Hence, in all numerical examples, we directly use the fracture resistance of the interphase for the null-thickness interface.

The computation is performed with monotonic displacement increment ∆u y . The displacements are prescribed along the y-direction for upper edge (y = 1) mm while the displacement along x is free.

Along the lower edge (y = 0) mm, the displacements along y are fixed to zero, while the displacements along x are free. Plane strain condition is assumed. In this first analysis, the influence of the load increment ∆u y in the numerical simulation on the mechanical response is investigated. The material parameters given in Table 1 are used. The internal length scale for the smeared crack and for the smeared interface are chosen to be the same = β = 0.0075 mm. The structure is meshed with triangular elements with h max e = 0.03 mm and h min e = 0.0035 mm for the critical region of expected crack path.

Three load increments are analyzed, with ∆u y = 5 × 10 -5 mm, ∆u y = 2 × 10 -5 mm, and ∆u y = 10 -5 mm. The overall stress-displacement curve is plotted in Fig. 4. The convergence of the mechanical response is obtained. The variation is bellow 1.5 % when the increment goes from 2 × 10 -5 to 10 -5 mm. This confirms the stability of the numerical algorithm based on the staggered scheme provided sufficiently small loading steps are used.

Mesh size effects

Next, the effects of mesh size on the numerical solution of the proposed model will be clarified.

The material properties along with the internal length , β are taken to be the same as in the previous example. The adaptive load increment is used, in which the monotonic displacement increments of ∆u y = 3 × 10 -5 mm have been prescribed as long as d < 0.9 in all nodes and ∆u y = 10 -5 mm, as soon as d ≥ 0.9 in one node. varies between h e = 2 µm and h e = 6 µm. Fig. 5 presents the obtained crack paths for different meshes at the chosen loading of U y = 6 × 10 -3 mm, wherein we also provide the zoom of mesh into regions close to the interface for each case. The observed fracture phenomena are the same. The initial crack propagates reaching the interface region. Then the interfacial delamination is captured.

The convergence of the mechanical behavior with respect to the mesh size is demonstrated in Fig. 6.

The solutions become mesh-independent when the conditions , β ≥ h e are fulfilled.

Influences of interfacial properties

In the following, the influences of the interfacial fracture resistance on the cracking behavior will be considered. The material properties λ (i,m,0) , µ (i,m,0) , g

(i,m) c
, internal length , β , as well as the adaptive load condition are taken to be the same as in the previous example. The fixed adaptive mesh is used with h max e = 0.03 mm and h min e = 0.0035 mm for the region of expected crack path.

Let us define the fracture energy mismatch ratio between matrix phase and interface f i by

g si,n c = f i × g (m) c and g si,t c = 1.05 × f i × g (m) c ( 65
)
We have performed three simulations corresponding to three values of f i , i.e, f i = 0.035, f i = 0.35 and f i = 0.85. Different fracture phenomena are observed in Fig. 7. A switching from the interfacial cracking mode (interfacial delamination) to bulk cracking mode (crack penetrates through the second phase) is clearly observed by increasing the strength of interface. The obtained result is in good agreement with the observation in the literature works [START_REF] He | Crack deflection at an interface between dissimilar elastic materials[END_REF][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]. The mechanical behavior of the structure is strongly influenced by the crack penetration/deflection mode at the interface, as depicted in Fig. 8.

More interestingly, the interface with lower fracture strength f i = 0.035 exhibits a better post-cracking resistance. Due to the interfacial delamination, the second phase is still uncracked at the end of the simulation, inducing the material resistance as observed in Fig. 8. In the contrary, the interfaces with higher fracture resistance (f i = 0.35, f i = 0.85), lead to the dominant bulk fracture mode. Hence the structure is totally collapsed for the loading U y = 8 × 10 -3 mm. This interesting phenomenon will be clarified in the following numerical experiment.

Influences of bulk material properties

The competition between bulk cracking and interfacial cracking is also affected by the material mismatch ratio of the two phases. In order to get in-depth insight into this behavior, we perform several simulations by keeping the same configuration of mesh, loading, and internal length as the last example. The interfacial fracture resistance f i = 0.85 is chosen. The material parameters of the first phase λ (m) , µ (m) , g

(m) c
are kept constant as given in Table 1, while the ones of the second phase

λ (i) , µ (mi) , g (i)
c are taken to vary as follows

λ (i) = f m × λ (m) and µ (i) = f m × µ (m) and g (i) c = f m × g (m) c . ( 66
) (a) f m = 0.25 (b) f m = 2.5 (c) f m = 5 (d) f m = 7.5
Figure 9: Crack paths obtained by changing the material mismatch between phase m and phase i at the loading U y = 7.9 µm. In the case,where phase i is softer than the phase m, the crack penetration with a dominant bulk cracking mode is captured. Contrarily, when the phase i is stiffer than the phase m, it will prevent the crack penetration from entering the second phase, resulting in the main interfacial cracking mode.

Four values of the material mismatch are chosen: f m = 0.25, f m = 2.5, f m = 5.0, and f m = 7.5.

The results of crack paths and corresponding mechanical response are plotted in Fig. 9, and Fig. 10. In general, the increase of the strength of the second phase: (i) prevents the crack penetration phenomena, (ii) forces the interfacial cracking mode, and (iii) improves the resistance of the structure.

Influences of length scales

The internal length which controls the thickness of the smeared crack and of the smeared interface, is also an important parameter of the regularized framework. Herein, the smeared crack width is often considered as a material parameter, and can be identified from the tough relation with the mechanical properties of the considered material [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF]. Hence, the individual effects of each length scale are not considered here for the sake of brevity. We are interested here influences of the ratio mismatch between β and on the principle of the regularized framework.

The ratio of the interface fracture resistance f i = 0.85 and the material mismatch ratio f m = 2.5 are selected. The smeared crack width = 0.01 mm is kept constant, while the one of the smeared interface β is taken to vary as follows

β = f , (67)
where f denotes the the ratio between the smeared interface width and crack width.

We perform five simulations corresponding to five different values of f i.e. f = [0.3, 0.5, 1.0, 1.5, 2.0].

The obtained crack patterns are provided in Fig. 11. In the case of β the interfacial effects are vanishing. It is consistent with the conception of the regularized scheme. The crack deflection behavior is quasi stable β ≥ , and we also archive the convergence of the mechanical response as depicted in Fig. 12. This means the internal length of the smeared interface should be chosen greater than one of the smeared crack to avoid the breakdown of the regularized scheme. However, it seems that the choice of β also slightly alters the mechanical response of the material. A detailed investigation on the relation of this parameter with the material parameters should be considered in future works.

Benchmark test of the CI model

In this example, the proposed CI model will be assessed. We consider a benchmark problem consisting in a square domain which contains a circular inclusion and an initial crack, as depicted in interface is very stiff as compared to the matrix and inclusion phases. Hence the coherent imperfect interface model is here adopted. The sharp interface is then regularized by the smeared interface γ α (α, ∇α). The material properties of bulk phases and interphase are given in the Table 2.

The tangential Lamé's constants for the interface can be determined by using Eq. ( 27), resulting in λ s = 13.714 GPa.mm, µ s = 16 GPa.mm for an interphase thickness of h = 2 µm. The internal length scale for the smeared crack is chosen = 0.0075 mm, while the internal length scale for the smeared interface is taken as α = 0.008 mm. The computation is performed with monotonic displacement increment ∆u y . The displacements are prescribed along the y-direction for upper edge (y = 1) mm while the displacements along x are here free. Along the lower edge (y = 0) mm, the displacements along y are fixed to zero, while the displacements along x are free. Plane strain condition is assumed.

The influences of the load increment and the mesh size on the numerical solution of the mechanical response are assessed. First, for the consideration of increment load effects, we perform three simulations on a fixed adaptive mesh of the element size h e ≈ 2 µm, and the displacement increment is taken to vary with ∆u y = 3 × 10 -5 mm, ∆u y = 8 × 10 -6 mm, and ∆u y = 4 × 10 -6 mm. Second, for the investigation of mesh size effects, three meshes are used to conduct the simulations, which are h e ≈ 2 µm, h e ≈ 3.5 µm, and h e ≈ 6 µm (the increment load is fixed ∆u y = 4 × 10 -6 mm). The obtained overall stress -displacement curves are depicted in Fig. 14. The convergence of the numerical simulations is clearly observed with to the load increment and the mesh size. Fig. 15 shows the evolution of the crack at several prescribed loads. Note that, a post-processing procedure is applied to highlight the effects of the heterogeneity and interface on the fracture problem.

The geometry of inclusion is represented and marked in green color and used as a background, which is then overlaid with the phase field. Moreover, the cracked area (phase field d > 0.98) is removed from the plot to better describe the fracture morphologies. The obtained results indicate that for the case of very stiff interface, the fracture will preferably take place at the matrix phase. Although the stress concentration can take place at the interface, it still does not allow for crack initiation at this region due to its very high fracture resistance. More interestingly, this phenomenon of crack deflection due to the strong interface has been experimentally captured in the literature studies, e.g., Sinien et al. [START_REF] Sinien | Microdamage and interfacial adhesion in glass bead-filled high-density polyethylene[END_REF] have demonstrated that when the interfacial bonding is very strong, the particles cannot be separated from the matrix, thus the crack has to climb over the particles. The work of Yang et al. [START_REF] Yang | Singlecrystal sic nanowires with a thin carbon coating for stronger and tougher ceramic composites[END_REF] shows that the too strong nanowire/matrix bonding strength avoids the debonding and pullout of the nanowires, or in the research of Zhu et al. [START_REF] Zhu | Microstructure and mechanical properties of cf/sic composites reinforced with boron nitride nanowires[END_REF] the authors also report that the stiff interface between BNNWs and the matrix leads to the main bulk cracking. The effects of stiff interface are totally in opposite to the last case of soft interface. A detailed comparison for the fracture behavior between these two situations will be provided in the next example. In order to highlight the effects of interface behavior on the fracture phenomena, the inclusion is considered to be five times stiffer than matrix phase. This choice ensures that the cracking will preferably take place at either the matrix phase or interfacial zone. Both situations of very stiff The tangential Lamé's constants for the interface can be determined by using Eqs. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], and ( 27):

with the interphase thickness h = 1 µm, we obtain λ s si = 105 GPa/mm, µ s si = 30 GPa/mm and λ s ci = 6.857 GPa.mm, µ s ci = 8 GPa.mm. The displacements are prescribed along the y-direction for upper edge (y = B) while the displacements along x are free. On the lower edge (y = 0), the displacements along y are fixed to zero, while the displacements along x are free. The monotonic displacement increments of ∆u y = 2 × 10 -5 mm have been prescribed via 1000 time steps. Plane strain condition is assumed. The obtained results of crack propagation are shown in Fig. 17 for (SI1) the first case of soft interface, in Fig. 18 for (SI2) the second case of soft interface, and in Fig. 19 for very stiff interface. In general, the following main fracture phenomena are observed:

• (SI1) Crack nucleation and initiation take place at the interface. Then the major interfacial cracking mode is captured.

• (SI2) Crack initiation always takes place at the interface. The interfacial cracking still present but less than the case (SI1).

• (CI) Cracks initiate at the matrix phase (zone close to the interface). The bulk cracking is the main fracture mode. The interfacial properties have clearly a strong impact on the fracture behavior of the heterogeneous material. Again, a switching from the interfacial cracking mode to bulk cracking mode is captured by increasing the interface strength. More specifically, the very compliant interface has created the weak zones, that induce the crack to initiate and propagate along the interface (a weaker path), rather than through the bulk material. On the contrary, the very strong interface can prevent the interfacial cracking problem. The reason for the crack to choose a convenient fracture mode is dependent on the ratios of the crack driving force (history function in Eq. ( 44)) to the fracture resistance of two kinds of failure modes [START_REF] Leguillon | The strengthening effect caused by an elastic contrast? part I: the bimaterial case[END_REF][START_REF] Nguyen | Role of interfacial transition zone in phase field modeling of fracture in layered heterogeneous structures[END_REF]. In the situation of the very stiff interface, its fracture resistance is very high. Thus, the cracks preferably propagate in bulk rather than in the first case of propagation along the interface. This will strongly alter the mechanical response of the material as depicted in Fig. 20.

The presence of very stiff interphase significantly enhances the material stiffness. More interestingly, it however exhibits a poorly post-cracking behavior. This phenomenon is experimentally observed in the work of Krishnan and Xu [START_REF] Krishnan | Experimental studies on the interaction among cracks, notches and interfaces of bonded polymers[END_REF]. The reason behind that can be related to the complex stress state at the crack tip in the case of the interfacial delamination. Regarding the loading direction, the stress acting ahead of the interfacial crack tip can be in the compressive state, which will increase the fracture resistance of the interface. Therefore, this interfacial crack needs very high energy (driving force) to exceed the fracture toughness, and propagates again along the interface. As a result, the applied load increases to a higher value in order to provide enough external work for increasing the crack driving force.

The very complex behavior of the interfacial cracking is reproduced by the present model. The competition and interaction between bulk cracking and interfacial cracking are successively simulated. With this numerical experiment, we expect to capture the complex behavior of interfacial cracking in 3D. Thus, only very soft interface will be considered. The material properties of the matrix, inclusion and interphase are given in Table 5. The tangential Lamé's constants for the interface can be determined by using Eq. ( 31), reaching λ s si = 3.5 GPa/mm, µ s si = 1 GPa/mm for the interphase thickness of h = 1 µm. The fracture resistance is taken g (m) c = 50 N/m for the matrix phase and g (i) c = 300 N/m for the inclusion. The resistance of fracture in the tangential direction and normal direction is assumed to be the same. The internal length scale for the smeared crack is chosen as = 0.01 mm, while the internal length scale for the smeared interface is taken as β = 0.015 mm. The incremental displacements are prescribed along the y-direction for upper surface (y = 1.2 mm) while the displacements along x, z are here free. Along the lower surface (y = 0 mm), the displacements along y are fixed to zero, while the displacements along x, z are free. The adaptive load increment is used, in which the monotonic displacement increments of ∆u y = 6 × 10 -5 mm have been prescribed

as long as d < 0.9 in all nodes and ∆u y = 1.75 × 10 -5 mm, as soon as d ≥ 0.9 in one node. • The (L-IFR) situation presents the crack initiation at the interface only, while the (H-IFR) situation shows the crack onset at both interface and pre-notched. It means a competition between bulk cracking and interfacial cracking, which is well consistent with experimental observation in literature, e.g. the debonding of glass bead from high-density polyethylene matrix [START_REF] Sinien | Microdamage and interfacial adhesion in glass bead-filled high-density polyethylene[END_REF].

• The resistance of the (H-IFR) case is significantly higher than (L-IFR) one, in particular its In addition, the semi-elliptical crack formation is captured for both cases. It seems that the propagation tendency agrees with the theory of semi-elliptic surface crack growth [START_REF] Newman | An empirical stress-intensity factor equation for the surface crack[END_REF][START_REF] Smith | A finite element model for the shape development of irregular planar cracks[END_REF][START_REF] Kikuchi | Surface crack growth simulation under mixed mode cyclic loading condition[END_REF]. This is clearly demonstrated in the evolution of pre-cracked zone. The crack first propagates to obtain the semi-elliptical form, then propagates in both directions (depth and surface) in order to keep this kind of morphology. More interestingly, a very complex morphology of the interfacial cracking is predicted by the present model as illustrated in Fig. 24. The physical interpretation of this phenomenon can be related to the mixed mode fracture (see the unilateral contact condition), which induces the heterogeneous crack surfaces. This example has demonstrated the robustness of the present model, which is possibly applied to study the interfacial cracking in the full 3D model. More importantly, a very complicated fracture behavior is reproduced by using this framework, including the complex interfacial decohesion, interaction of bulk/interfacial cracking. Although the prediction accuracy should be validated by the experiment in the future work, the proposed method is very promising to provide a new computational tool to help material design.

Figure 1 :

 1 Figure 1: Description of the proposed modeling methodology: (a) three-phase configuration; (b) two-phase configuration with sharp interface; and (c) two-phase configuration with smeared interface and smeared crack.

Figure 2 :

 2 Figure 2: Illustration of the regularized scheme: (a) sharp description of interfaces and crack; (b) smeared description of interfaces and crack.

d

  e are the nodal values of d. The same discretization is employed for the test function δd = N d δd e and ∇δd = B d δd e .

Figure 3 :

 3 Figure 3: Illustration of the Benchmark problem for the SI model: geometry and boundary conditions.

Figure 4 :

 4 Figure 4: Comparison of the mechanical response for different load increments. The convergence of the solution with respect to the increment is obtained.

  Mesh size h e = 6 μm (b) Mesh size h e = 3.5μm (c) Mesh size h e = 2 μm

Figure 5 :Figure 6 :

 56 Figure 5: Comparison of crack paths obtained from different mesh size. Upper row: zoom of mesh around the interface; Lower row: crack paths. The same fracture phenomena are observed for all three mesh sizes.

Figure 7 :Figure 8 :

 78 Figure 7: Crack paths obtained by changing the interfacial fracture resistance at the loading U y = 8.0 µm.A transition from crack deflection to crack penetration is observed when the interfacial strength is increased. This is consistent with previous studies[START_REF] He | Crack deflection at an interface between dissimilar elastic materials[END_REF][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].

Figure 10 :

 10 Figure 10: Comparison of the mechanical response for different material mismatch ratios.

Figure 11 :Figure 12 :

 1112 Figure 11: Crack paths obtained by varying the ratio mismatch between the smeared interface width β and the smeared crack width .

Fig. 13 .Figure 13 :

 1313 Fig. 13. The size of the square domain is L × B = 1 × 1 mm 2 .

Figure 14 :

 14 Figure 14: Assessment of the coherent interface model: convergence of the mechanical response with respect to the load increment and the mesh size.

Figure 15 : 5 . 3 .

 1553 Figure 15: Fracture behavior of the structure containing a very stiff interface. The results of crack propagation for the case of the mesh size h e = 2 µm, and the load increment ∆u y = 4 × 10 -6 mm.

  interface (CI) and very soft interface (SI) are considered. Moreover, in the case of the soft interface, we also perform the simulation on two sub-cases: one is devoted to very weakly interface strength (SI1) with g si,n c = 3.75 × 10 -5 kN/mm and g si,t c = 3.18 × 10 -5 kN/mm, and the other one (SI2) is these values about seven times stronger, i.e. g si,n c = 2.50 × 10 -4 kN/mm and g si,t c = 2.13 × 10 -4 kN/mm.

Figure 16 :

 16 Figure 16: 2D microstructure simulation.

Inclusion 3 × 10 - 3 Figure 17 :Figure 18 :

 31031718 Figure 17: Complex crack paths obtained from the soft interface situation (SI1) with g si,n c = 0.625 N/m,

Figure 19 :

 19 Figure 19: Fracture phenomena obtained from the stiff interface situation (CI). The major bulk cracking is captured.

Figure 20 : 5 . 4 .

 2054 Figure 20: Comparison of mechanical response for the different interfacial properties. The stiff interface strongly enhances the stiffness of the material, but exhibit a poorly post-cracking behavior.

Figure 21 :

 21 Figure 21: Full 3D model of interfacial cracking: geometry and boundary conditions.

- 3 Figure 22 :

 322 Figures 22 and 23 depict the crack propagation for the case of low interfacial fracture resistance (L-IFR) and for the case of high interfacial fracture resistance (H-IFR), respectively. The following phenomena are noted:

  (a) U y = 5.11 µm (b) U y = 5.51 µm (a) U y = 6.00 µm (b) U y = 8.08 µm

Figure 23 :

 23 Figure 23: High interfacial fracture resistance (H-IFR): complex 3D crack morphology, both bulk and interfacial cracking are simultaneously observed.

Figure 24 :

 24 Figure 24: Complex fracture morphologies predicted by using the present model.

Table 1 :

 1 Material properties of the spring interface model

	Parameter Matrix Inclusion	Interphase	Unit
	λ	12	36	3.75 × 10 -2	GPa
	µ	8	24	2.5 × 10 -2	GPa
	g c	50	150	g si,n

c = 17.5 g si,t c = 18.375 N/m 5.1.1. Incremental loading effects

Table 2 :

 2 Material properties of the coherent interface model.

	Parameter Matrix Inclusion Interphase Unit
	λ	12	36	12 × 10 3	GPa
	µ	8	24	8 × 10 3	GPa
	g c	50	150	500	N/m

Table 3 :

 3 Material properties of the spring/coherent interface model.

	Parameter Matrix Inclusion Soft interphase Stiff interphase Unit
	λ	18	60	4.5 × 10 -2	6 × 10 3	GPa
	µ	12	32	3 × 10 -2	4 × 10 3	GPa

Table 4 :

 4 Fracture resistance of the spring/coherent interface model.

	Phases	Fracture resistance [kN/mm]
	Matrix	5 × 10 -4

Table 5 :

 5 Material properties of the coherent interface model.

	Parameter Matrix Inclusion Interface Unit

The tensor product operation A⊗B is defined by (A⊗B) ijkl = (A ik B jl + A il B jk ) /2 for any two second-order tensors A and B.

The positive/negative operators of a vector is defined as follows[v i , v j , v k ] ± = [ v i ±, v j ±, v k ±].

Acknowledgments

This research is supported by the project SeRaMCo within the programme INTERREG NWE.

Appendix A. Derivation of a stiffness tensor for the interfacial cracking model

Summary of the computational time

For all cases, a workstation with 18 cores, 112 Go Ram and 2.4 GHz processor was used (see https://hpc.uni.lu/systems/iris/). The present code has been implemented in Matlab.

Conclusion

A robust interfacial cracking model in the phase field framework is proposed in this work. The new method is based on the mixing of bulk energy with the energy of null-thickness interface which is derived from the asymptotic analysis of the very thin interphase. The new model is able to study the complex interfaces with various properties. More specifically, the effects of both stiff and soft interphases have been effectively accounted for, providing a powerful computational tool to accurately simulate the mechanical behavior of composite materials and structures. Furthermore, the numerical implementation using one-pass staggered algorithm results in two linear problems for the successive update of the phase field and the displacement field within a typical time step. That constructs one of the most efficient models to study interfacial cracking phenomena.

The proposed model is applied to investigate the interfacial cracking in composite materials and to determine the influences of thin interphase properties on the global behavior of advanced material and structures. The complex cracking phenomena on interfaces such as initiation, delamination, coalescence, deflection are successfully predicted by the present method. Especially, an extremely complicated interfacial-decohesion (heterogeneous crack surface) is observed in the three-dimensional simulation for the first time to the best of our knowledge. The strong impacts of the interface on the mechanical performance of materials are also captured. The method then constitutes a very promising modeling and simulation tool to accurately predict the behavior of a broad class of composite materials and could form a basis for many extensions of interfacial damage models involving other physical phenomena.

Appendix A.1. Bulk cracking model

With the use of the unilateral contact condition, the density of elastic energy for the isotropic materials is expressed as

By replacing tr ε e ij ± = R ± tr ε e ij and ε e± ij = P ± ijkl ε e kl , Eq. (A.1) can be rewritten as

From Eqs. (A.1), and (A.2), the elastic tensor is defined by

. Interfacial cracking model

In the case of imperfect spring interface, with the introduction of the projection tensor Q ± , the following relations are obtained for the displacement jump across interfaces

The density of interfacial energy can now be reformulated as follows

This implies the expression for the tangential interface stiffness accounting damage, C si s (d), as .6)