Daicong Da

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

Liang Xia

Guangyao Y Li
email: gyli@hnu.edu.cn

Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage

Keywords: Cracks, Phase field method, Interfacial damage, Topology optimization, BESO

Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage

Daicong Da, Julien Yvonnet, Liang Xia,

INTRODUCTION

Mechanical and physical properties of complex heterogeneous materials are determined on one hand by the composition of their constituents, but can on the other hand be drastically modified, at a constant volume fraction of heterogeneities, by their geometrical shape and by the presence of interfaces. For example, it is well known that in the context of linear properties, modifying the geometrical topology of phases at a constant volume fraction can induce a change in the maximum values of effective moduli, but also induces anisotropy. In that context, topological optimization of microstructures can help designing materials with higher effective properties while maintaining the volume fraction of constituents, or to obtain new properties which are not naturally available (metamaterials). Recently, the development of 3D printing techniques and other additive manufacturing processes (see e.g. [START_REF] Quan | Microstructural design and additive manufacturing and characterization of 3d orthogonal short carbon fiber/acrylonitrile-butadiene-styrene preform and composite[END_REF]) have made possible to manufacture directly the designed materials from a numerical file, opening routes for totally new designs. Topological optimization of microstructures has been extensively studied in the linear context (see e.g. a review in [START_REF] Cadman | On design of multi-functional microstructural materials[END_REF] and the reference therein), and more recently in the nonlinear context (see e.g. [START_REF] Fritzen | Topology optimization of multiscale elastoviscoplastic structures[END_REF]). One major challenge is to use topological optimization to improve the fracture resistance of heterogeneous materials, taking into account the heterogeneities and their interfaces in the material. Optimization design of composite materials accounting for fracture resistance remains relatively unexplored so far, mainly due to the lack of robust numerical methods for fracture propagation in presence of complex heterogeneous media and interfaces, until recently. In addition, these numerical simulation models should be formulated in a context compatible with topological optimization (e.g. finite elements). In [START_REF] Gu | Optimization of composite fracture properties: Method, validation, and applications[END_REF], Gu et al. used a modified greedy optimization algorithm for composites made up of soft and stiff building blocks to improve material toughness. San and Waisman [START_REF] San | Optimization of carbon black polymer composite microstructure for rupture resistance[END_REF] explored the optimal location of carbon black particle to maximize the rupture resistance of polymer composites by a genetic algorithm. In a recent work by the authors [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF], topology optimization for maximizing the fracture resistance of quasi-brittle composites has been introduced by combining the phase field method and a gradient-based (BESO) algorithm. However, in the mentioned work, the crack propagation resistance was only evaluated on the basis of phase distribution. In most heterogeneous quasi-brittle materials (e.g. ceramic matrix composites, cementitious materials), the interfacial damage plays a central role in the nucleation and propagation of microcracks [START_REF] Lamon | The influence of the interphase and associated interfaces on the deflection of matrix cracks in ceramic matrix composites[END_REF][START_REF] Tvergaard | Model studies of fibre breakage and debonding in a metal reinforced by short fibres[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] Narducci | Exploiting nacre-inspired crack deflection mechanisms in cfrp via micro-structural design[END_REF]. The main objective of the present work is to extend the framework developed in [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF] for defining through topological optimization the optimal phase distribution in a quasi-brittle composite with respect to fracture resistance, taking into account crack nucleation both in the matrix and in the interfaces, as described in Fig. 1. To our best knowledge, such study is investigated here for the first time.

Fracture Energy

Simulating interfacial damage and its interaction with matrix crack for complex heterogeneous materials is a highly challenging issue for meshing algorithms. Many numerical methods such as eXtended Finite Element Method (XFEM) [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF][START_REF] Sukumar | Extended finite element method for three-dimensional crack modelling[END_REF], Thick Level-set method (TLS) [START_REF] Bernard | Damage growth modeling using the thick level set (TLS) approach: Efficient discretization for quasi-static loadings[END_REF][START_REF] Cazes | Damage growth modeling using the thick level set (TLS) approach: Efficient discretization for quasi-static loadings[END_REF], and Phase Field Method (PFM) [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] Kuhn | A continuum phase field model for fracture[END_REF], among the most recent popular techniques, have been introduced to investigate this topic.

Based on the pioneer works of Marigo and Francfort [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], the phase field method makes use of a variational principle framework for brittle fracture [START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Pham | The variational approach to damage: I. the foundations[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] and of a regularized description of the discontinuities related to the crack front [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF][START_REF] Ambrosio | Approximation of functional depending on jumps by elliptic functional via γconvergence[END_REF]. The method has been adapted to a convenient algorithmic setting by Miehe et al. in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF]. The main advantages of the method are: (a) the crack paths are mesh-independent; (b) initiation and propagation of multiple, complex cracks patterns can be easily handled; (c) the method is convergent with respect to the mesh size and not sensitive to the mesh regularity; (d) it is stable due to its inherent gradient-based formulation. In the present work, we use the extension of the phase field to interfacial damage as 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] to take into account both bulk brittle fracture and interfacial damage.

During the past decades, the research in topology optimization witnessed tremendous development since the seminal paper by Bensoe and Kikuchi [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF]. So far, various topology optimization methods have been proposed, e.g., density-based method [START_REF]Topology Optimization: Theory, Methods and Applications[END_REF][START_REF] Nguyen-Xuan | A polytree-based adaptive polygonal finite element method for topology optimization[END_REF][START_REF] Chau | A polytree-based adaptive polygonal finite element method for multi-material topology optimization[END_REF], evolutionary procedures [START_REF] Xie | A simple evolutionary procedure for structural optimization[END_REF][START_REF] Huang | Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method[END_REF], level-set methods [START_REF] Sethian | Structural boundary design via level set and immersed interface methods[END_REF][START_REF] Wang | A level set method for structural topology optimization[END_REF][START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF]. The extended bi-directional evolutionary structural optimization (BESO) method developed for the design of elastoplastic structure in [START_REF] Fritzen | Topology optimization of multiscale elastoviscoplastic structures[END_REF][START_REF] Xia | Evolutionary topology optimization of elastoplastic structures[END_REF] is adopted in this work to carry out the topology optimization. A computationally efficient adjoint sensitivity formulation accounting for the complete fracturing process is derived, involving cracks nucleation, propagation and merging of microcracks until complete failure of the specimen. By using discrete design variables, the ESO-type methods omit naturally the definition of additional pseudo-relationships between fictitious materials and fracture toughness, resulting in a clear physical interpretation and algorithmic advantages [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF][START_REF] Xia | Bi-directional evolutionary structural optimization on advanced structures and materials: A comprehensive review[END_REF][START_REF] Xia | Concurrent topology optimization design of material and structure within FE 2 nonlinear multiscale analysis framework[END_REF][START_REF] Xia | Multiscale structural topology optimization with an approximate constitutive model for local material microstructure[END_REF]. While there do exist numerous researches on topology optimization including material interface behavior (see e.g. [START_REF] Vermaak | Material interface effects on the topology optimziation of multi-phase structures using a level set method[END_REF][START_REF] Lawry | Level set topology optimization of problems with soliding contact interfaces[END_REF][START_REF] Liu | Multi-material topology optimization considering interface behaviour via xfem and level set method[END_REF][START_REF] Behrou | Level set topology optimization of structural problems with interface cohesion[END_REF]), and with the enforcement of stress constraints [START_REF] Duysinx | Topology optimization of continuum structures with local stress constraints[END_REF][START_REF] Guo | Stress-related topology optimization via level set approach[END_REF][START_REF] Xia | A level set solution to the stress-based structural shape and topology optimization[END_REF][START_REF] Bruggi | Topology optimization for minimum weight with compliance and stress constraints[END_REF][START_REF] Luo | An enhanced aggregation method for topology optimization with local stress constraints[END_REF][START_REF] Cai | Stress constrained shape and topology optimization with fixed mesh: A b-spline finite cell method combined with level set function[END_REF][START_REF] Cai | Stress constrained topology optimization with free-form design domains[END_REF], local buckling constraint [START_REF] Sun | Topological optimization of biomimetic sandwich structures with hybrid core and cfrp face sheets[END_REF] or damage constraints [START_REF] Challis | Fracture resistance via topology optimisation[END_REF][START_REF] Amir | Reinforcement layout design for concrete structures based on continuum damage and truss topology optimization[END_REF][START_REF] Amir | A topology optimization procedure for reinforced concrete structures[END_REF][START_REF] Jansen | Topology optimization of fail-safe structures using a simplified local damage model[END_REF][START_REF] James | Failure mitigation in optimal topology design using a coupled nonlinear continuum damage model[END_REF][START_REF] Kang | Topology optimization considering fracture mechanics behaviors at specified locations[END_REF], topology optimization for maximizing the fracture resistance, taking into account interactions between interfacial damage and bulk brittle fracture for complete fracturing process is to our best knowledge explored for the first time in the present work.

The layout of the paper is organized as follows. Section 2 gives the detailed phase field framework incorporating bulk fracture and cohesive interface. Numerical details and FEM discretization details are presented. In Section 3, we formulate the topology optimization model, containing detailed sensitivity derivation as well as the updating scheme with a damping on the sensitivity numbers. In Section 4, several numerical benchmark tests are presented to demonstrate the potential of the proposed method. Finally, conclusions are drawn in Section 5.

PHASE FIELD MODELING OF BULK CRACK AND COHESIVE INTERFACES

In this section we described the numerical method which is used by the topological optimization algorithm to obtain the fracture energy of the sample. The technique is based on the phase field model for fracture extended to interfacial damage by Nguyen et al. [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], and allows simulating the initiation and propagation of a multiple cracks network in heterogeneous microstructures. The main concepts and details are recalled in the following.

Regularized representation of discontinuous field

Let Ω ∈ R d be an open domain describing a solid with external boundary ∂Ω. The solid contains internal material interfaces between different phases, collectively denote by Γ I . During the loading, cracks may propagate within the solid and can pass through the material interfaces as depicted in Fig. 2 (a). The crack surfaces are denoted by Γ c . In this work, we adopt smeared representations of both cracks and material interfaces: , the cracks are approximated by an evolving phase field d(x, t).

Interfaces between different material phases are described by a fixed scalar phase field β(x). The material interfaces do not evolve during the loading. The regularized parameters describing the actual widths of the smeared cracks and material interfaces are respectively denoted by d and β .

In the following, the same regularization length = β = d is adopted for cracks and material interfaces for the sake of simplicity. Given a non-evolving sharp crack defined on a surface Γ c , a regularized (smeared) representation of the corresponding damage d(x) (see Fig. 2 (c)) can be obtained by solving the following equations (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) -2 ∇ 2 d(x) = 0, in Ω d(x) = 1, on Γ c ∇d(x, t) • n = 0, on ∂Ω (1)
where ∇ 2 (.) is the Laplacian, is the regularization parameter, and n is the unit outward normal to the external boundary ∂Ω. It can be shown that Eq. ( 1) is the Euler-Lagrange equation corresponding to the variational problem:

d(x, t) = Arg{ inf d∈S d Γ d (d)}, Γ d (d) = Ω γ d (d) dΩ, (2)
where

S d = {d | d(x) = 1, ∀x ∈ Γ }
γ d (d(x)) = 1 2 d(x) 2 + 2 ∇d(x) • ∇d(x). (3)
When → 0, the above variation principle leads to the corresponding one with a sharp crack description in the Γconvergence sense. The details of definition (3) can be found e.g. in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF].

In Nguyen et al. [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], an interface phase field has been introduced to describe in the same manner the discontinuities related to the damage of interfaces and obtained by solving the problem:

       β(x) -2 ∇ 2 β(x) = 0, in Ω β(x) = 1, on Γ I ∇β(x) • n = 0, on ∂Ω. (4)
Eq. ( 4) corresponds to the Euler-Lagrange equation associated with the variational problem

β(x, t) = Arg{ inf β∈S β Γ β (β)}, Γ β (β) = Ω γ β (β) dΩ (5)
where

S β = {β | β(x) = 1, ∀x ∈ Γ I }, Γ β
represents the total interface length, and γ β is defined by

γ β (β) = 1 2l β(x) 2 + 2 ∇β(x) • ∇β(x). (6)
For → 0 the above variation principle leads to a sharp interface description. This function will be used as an indicator to particularize the damage model to the interfaces or to the bulk in the formulation described in following sections.

In addition, it is necessary to introduce an approximation for the displacement jump at the interfaces, to associate a damage model specific to the interfaces and different from the bulk. For this purpose, the following approximation has been proposed [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] using Taylor expansion of the displacement field around a point x located on the interface:

[[u(x)]] w(x) = u(x + h 2 n I ) -u(x - h 2 n I ) = h∇(u(x))n I , (7)
where w(x) denotes the smoothed displacement jump approximation and n I is an approximation of the normal to the interface Γ I at a point x. Several techniques are possible to define this normal.

For example, 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 level-set technique has been proposed (see more detailed in the mentioned paper). In [START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF], another definition using the possible modification of the interface by the damage related to bulk cracks was introduced.

Energy functional

The following total energy functional formulated is introduced for the solid body related to both cracks and interfaces:

E = Ω W e u (ε e (u, β), d)dΩ + Ω [1 -β(x)]g c γ d (d)dΩ + Ω ψ I (w)γ β (β)dΩ, (8)
where g c is the toughness and ψ I is a strain density function depending on the displacement jump across the interface Γ I . Above, ε e is the bulk part of the infinitesimal strain tensor ε which satisfies the following relationship:

ε = ε e + ε, (9)
where ε is the strain part induced by the smoothed jump at the interfaces such that ε → 0 away from the interfaces (see [START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF]), in which case we recover the energy functional for a cracked body without considering the interface behavior:

E = Ω W e u (ε e (u), d)dΩ + Ω g c γ d (d)dΩ. ( 10
)
From Eq. ( 8), the free energy W can be identified as

W = W e u (ε e (u, β), d) + [1 -β(x)]g c γ d (d) + ψ I (w, α)γ β (β). (11)
Following [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], the elastic energy W e u is defined as

W e u = ψ + e (ε e )[g(d) + k] + ψ - e (ε e ) (12)
where

ε e = ε e+ + ε e-, (13)
and

ψ ± e (ε) = λ tr[ε e ] 2 ± /2 + µtr[ε e± ] 2 , (14)
and

ε e± = D i=1 ε i ± n i ⊗ n i . ( 15
)
The decomposition [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF] ensures that damage is only induced by traction, which removes the complications related to auto-contact in the crack. In ( 14)-( 15), x ± = (x ± |x|)/2. ε i and n i are respectively the eigenvalues and eigenvectors of ε e . λ and µ are initial Lamé coefficients and k is a very small value to maintain the well-posedness of the system. The degradation function g(d) is assumed to have a simple form g(d) = (1 -d) 2 . A reduced Clausius-Duhem inequality form related to the evolution of the damage parameter d can be written as

A ḋ ≥ 0, A = - ∂W ∂d . (16)
A threshold function F (A) is assumed such that no damage occurs in the form:

F (A) = A ≤ 0. ( 17
)
The principle of maximum dissipation requires that dissipation A ḋ must be maximum under the constraint, i.e. ḋ > 0, F = 0. Therefore, it yields

F = - ∂W ∂d = - ∂W e u ∂d -(1 -β)g c δγ(d) = 0 (18)
with the functional derivative [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] δγ

(d) = d -∆d. (19)
It follows that when ḋ > 0 then

-2(1 -d)ψ + e + (1 -β)g c δγ(d) = 0, (20)
To handle loading and unloading, the strain history function adopted in [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 for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF] is employed here:

H(x, t) = max τ ∈[0,t]
{ψ + e (x, t)} [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] and ( 20) is substituted by

-2(1 -d)H + (1 -β)g c δγ(d) = 0. ( 22
)
2.3. Displacement and phase field problems 2.3.1. Weak form of the phase field problem Using ( 19), the evaluation of the crack field d(x, t) can be determined by solving the following phase field problem

       2(1 -d)H -(1 -β) gc (d -2 ∇ 2 d) = 0, in Ω d(x) = 1, on Γ c ∇d(x) • n = 0, on ∂Ω. (23)
The associated weak form is obtained as (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]):

Ω (2H + [1 -β] g c )dδd + [1 -β]g c ∇d • ∇(δd) dΩ = Ω 2HδddΩ. ( 24
)
2.3.2. Weak form of the displacement problem Using the variational principle for minimizing the total energy E with respect to the displacement u, the weak form associated with the displacement problem can be formulated as

Ω ∂W e u ∂ε e : ε e (δu)dΩ + Ω ∂ψ I (w) δw • ∂wγ β (β)dΩ = Ω f • δudΩ + ∂Ω F F • δudΓ = δW ext .
(25) In the absence of body forces, Eq.( 25) can be re-written as

Ω σ e : ε e (δu)dΩ + t(w) • δwγ β (β)dΩ - Ω σ e : ∇ s δudΩ = 0, (26)
where σ e = ∂we ∂ε e is the Cauchy stress, and t(w) is the traction vector acting on the interface Γ I oriented by n I and δw = h∇(δu)n I . Using σ e n = t, the above equation can be further re-written as which is satisfied for an admissible strain field:

Ω σ e : {ε e (δu) + n ⊗ δwγ β (β) -∇ s δu} dΩ (27)
ε e = ∇ s u -n⊗ S wγ β (28)
where

(∇ s u) ij = (u i,j + u j,i )/2 and (n⊗ S w) ij = (n i w j + w i n j ). From (13), ε can be identified as ε = n⊗ S wγ β .
With the above description of strain energy function, the Cauchy stress now reads:

σ e = ∂ψ + e ∂ε e [g(d) + k] + ∂ψ - e ∂ε e = [(1 -d) 2 + k] λ trε e + 1 + 2µε e+ + λ trε e -1 + 2µε e-. (29)

Model for interface damage

The general form of the traction vector t(w) in Eq. ( 26) is given by t

(w) = [t n , t t ] T ( 30
)
where t n and t t denote respectively normal and tangential parts of the traction vector t across the interface oriented by its normal n I . In the presented work, a simplified nonlinear elastic cohesive model is used. We have shown that 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] allows avoiding the use of internal variables related to the interface damage to handle loading and unloading by exploiting the damage phase field itself as the history variable. In addition, by only taking into account the normal traction, i.e. t(w) • n I = t n , the cohesive law can be written as

t n = g I c ( w n δ n )exp(- w n δ n ) ( 31
)
where w n = w • n I . For this model, we obtain K I = ∂t(w) ∂w . The relationship between δ n , the toughness g I c and th fracture strength t u is given by δ n = g I c /(t u e), with e = exp(1) (see Fig. 3). Details about the linearization of the displacement problem and about FEM discretizations and numerical implementation are provided in Appendix 6.

TOPOLOGY OPTIMIZATION METHOD

In this section, we present the topological optimization method based on the bi-directional evolutionary structural optimization (BESO) [START_REF] Huang | Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method[END_REF]. In the proposed procedure, the geometry of inclusions is optimized so as to maximize the fracture resistance of the sample. This procedure involves evaluating the sensitivity of the whole fracturing process (initiation of multiple cracks, propagation and complete failure of the sample) with respect to changes in the geometry.

Model definitions

The considered domain Ω is discretized into N i finite elements and each element i is assigned with a topology design variable ρ i . Within the framework of the bi-directional evolutionary structural optimization (BESO) method, the discrete design variables are defined as

ρ i = 0 or 1, i = 1, 2, . . . , N i . ( 32
)
Following the multiple material interpolation model in [START_REF] Huang | Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials[END_REF][START_REF] Da | Concurrent topological design of composite structures and the underlying multiphase materials[END_REF] we have

E i = ρ i E inc + (1 -ρ i )E mat (33)
where E inc and E mat are the Young's moduli of the inclusion and matrix phases, respectively. The Poisson's ratios of the two material phases are assumed identical. The design variables can thus be interpreted as an indicator for the generalized material density. From [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF], the density value of zero corresponds to the matrix phase, whereas one corresponds to the inclusion phase.

In nonlinear topology optimization, it is conventional to adopt displacement-controlled loading due to stability considerations [START_REF] Maute | Adaptive topology optimization of elastoplastic structures[END_REF][START_REF] Schwarz | Topology and shape optimization for elastoplastic structural response[END_REF][START_REF] Huang | Topology optimization of nonlinear structures under displacement loading[END_REF][START_REF] Xia | Recent advances on topology optimization of multiscale nonlinear structures[END_REF]. For a prescribed displacement load, the objective of structural fracture resistance maximization is equivalent to maximization of the mechanical work during the fracturing process. In practice, the total mechanical work J is calculated by using numerical integration, i.e.

J ≈ 1 2 n load n=1 f (n) ext + f (n-1) ext T ∆u (n) , (34)
and will be used as our definition to fracture resistance in the following. Above, n load is the total number of displacement increments, ∆u (n) is the n-th nodal displacement components and f n ext is the external nodal force at the n-th load increment.

During the optimization design, the material volume fractions of matrix and of inclusion phases are prescribed. Then the topology optimization problem subjected to balance equation and inclusion volume constraint can be formulated as

max ρ : J(ρ, u, d, β) subjected to : R = 0 : V (ρ) = ρ i v i = V req : ρ i = 0 or 1, i = 1, . . . , N i . (35)
In the above, v i is the volume of i-th element, V (ρ) and V req are the total and required material volumes, respectively and R denotes the nodal residual force:

R = f ext -f int . (36)
In [START_REF] Xia | Concurrent topology optimization design of material and structure within FE 2 nonlinear multiscale analysis framework[END_REF], f int is defined in each element as the internal force vector given in terms of the associated topology design variable ρ i and the Cauchy stress as:

f int = Ni i=1 ρ i Ωi B T σ e dΩ i . (37)
It should be mentioned that in this framework the discrete topology design variables are adopted to indicate the existence of the associated material phase (matrix/inclusion) of the element i. Then there is no need to define the physical meaning of the intermediate densities as is the case in the continuous density-based models (see e.g., [START_REF]Topology Optimization: Theory, Methods and Applications[END_REF]) resulting in a clear physical interpretation.

Sensitivity analysis

In order to perform the topology optimization, the sensitivity of the objective function J with respect to topology design variables ρ must be computed. The adjoint method is employed to derive the sensitivity analysis. Firstly, two Lagrangian multipliers µ (n) , λ (n) of the same dimension as the vector of unknowns u are introduced to enforce zero residual R at time t n-1 and t n for each term of the total mechanical work [START_REF] Xia | Evolutionary topology optimization of elastoplastic structures[END_REF]. Then the objective function J can be rewritten in the following form without modifying the original objective value as

J = 1 2 n load n=1 f (n) ext + f (n-1) ext T ∆u (n) + λ (n) T R (n) + µ (n) T R (n-1) . ( 38
)
Due to the asserted static equilibrium the residuals R (n) and R (n-1) must vanish. The objective value is thus invariant with respect to the values of the Lagrangian multipliers λ (n) and µ (n) (n = 1, . . . , n load ), i.e.

J ρ; λ (n) , µ (n) n=1,...,n load = J (ρ) . ( 39
)
This equivalence also holds for the sensitivity with respect to changes of the topology design variable

ρ i on element i ∂ J ∂ρ i = ∂J ∂ρ i . ( 40
)
In the following the derivative ∂ J/∂ρ i is computed with properly determined values of λ (n) and µ (n) leading to simplifications in the derivation. To formally describe these derivations, we introduce a partitioning of all degrees of freedom (DOF) into essential (index E; associated with Dirichlet boundary conditions) and free (index F; remaining DOF) entries. For a vector w and a matrix M we have

w ∼ w E w F and M ∼ M EE M EF M FE M FF . (41)
In the current context, the displacements u E on the Dirichlet boundary are prescribed. They are independent of the current value of ρ. This implies that

∂∆u ∂ρ i = ∂ ∂ρ i ∆u E ∆u F = 0 ∂{∆u F }/∂ρ i (42)
holds for arbitrary increment loads, i.e. for u = u (n) or u = u (n-1) . The components f ext,F of the force vector f ext vanish at all times and the only (possibly) non-zero components are the reaction forces

f ext,E f (n) ext = f (n) ext,E 0 . ( 43
)
Eqs. ( 42) and ( 43) imply

f (m) ext T ∂∆u (n) ∂ρ i = 0. ( 44
)
Hence for arbitrary time indices n, m = 1, . . . , n load , we have

∂ ∂ρ i f (m) ext T ∆u (n) = ∂f (m) ext ∂ρ i T ∆u (n) . (45)
Therefore, the derivative of the modified design objective in [START_REF] Vermaak | Material interface effects on the topology optimziation of multi-phase structures using a level set method[END_REF] can be rewritten as:

∂ J ∂ρ i = 1 2 n load n=1    ∂f (n) ext ∂ρ i + ∂f (n-1) ext ∂ρ i T ∆u (n) + λ (n) T ∂R (n) ∂ρ i + µ (n) T ∂R (n-1) ∂ρ i    . (46)
Recall the balance equation at each load time increment in [START_REF] Huang | Topology optimization of nonlinear structures under displacement loading[END_REF], the derivatives of R (m) at the equilibrium of the m-th load increment with respect to ρ i can be expanded as

∂R m ∂ρ i = ∂f m ext ∂ρ i - Ωi B T (σ e ) (m) dΩ i -K (m) tan ∂∆u (m) ∂ρ i (47)
where

K (m) tan = - ∂R (m) ∂u (m) (48)
is the tangent stiffness matrix of the nonlinear mechanical system at the balance equation of the m-th load increment. With the expression ( 47), ( 46) can be reformulated as

∂ J ∂ρ i = 1 2 n load n=1 ∂f (n) ext ∂ρ i T ∆u (n) + λ (n) + ∂f (n-1) ext ∂ρ i T ∆u (n) + µ (n) -λ (n) T Ωi B T (σ e ) (n) dΩ i + K (n) tan ∂∆u (n) ∂ρ i -µ (n) T Ωi B T (σ e ) (n-1) dΩ i + K (n-1) tan ∂∆u (n-1) ∂ρ i . ( 49
)
As mentioned before, the aim is to find proper values of the Lagrangian multipliers λ (n) and µ (n) such that the sensitivities can be explicitly and efficiently computed. From [START_REF] Guo | Stress-related topology optimization via level set approach[END_REF], the first two terms in [START_REF] Sun | Topological optimization of biomimetic sandwich structures with hybrid core and cfrp face sheets[END_REF] can be omitted by setting

λ (n) E = -∆u (n) E and µ (n) E = -∆u (n) E . ( 50
)
Accounting further for the structure of the sensitivities of u in [START_REF] Duysinx | Topology optimization of continuum structures with local stress constraints[END_REF] and for the symmetry of the stiffness matrices we have

∂ J ∂ρ i = 1 2 n load n=1 -λ (n) T Ωi B T (σ e ) (n) dΩ i -µ (n) T Ωi B T (σ e ) (n-1) dΩ i -K (n) tan,FE λ (n) E + K (n) tan,FF λ (n) F T ∂∆u (n) F ∂ρ i -K (n-1) tan,FE µ (n) E + K (n-1) tan,FF µ (n) F T ∂∆u (n-1) F ∂ρ i . ( 51
)
To avoid the evaluation of the unknown derivatives of u

F and u

(n-1) F

, i.e. eliminating the last two lines of ( 51), the values of λ

(n) E = -∆u (n) E and µ (n) E = -∆u (n) E at the essential nodes: λ (n) F = K (n) tan,FF -1 K (i) tan,FE ∆u (n) E , (52)
and µ

(n) F = K (n-1) tan,FF -1 K (n-1) tan,FE ∆u (n) E . ( 53
)
The two relations [START_REF] Amir | A topology optimization procedure for reinforced concrete structures[END_REF] and [START_REF] Jansen | Topology optimization of fail-safe structures using a simplified local damage model[END_REF] together with [START_REF] Challis | Fracture resistance via topology optimisation[END_REF] fully determine the values of the Lagrange multipliers λ (n) and µ (n) . Finally, the objective function ∂ J/∂ρ i can be computed via

∂ J ∂ρ i = - 1 2 n load n=1 λ (n) T Ωi B T (σ e ) (n) dΩ i + µ (n) T Ωi B T (σ e ) (n-1) dΩ i . ( 54
)
The computation of the sensitivity consists in solving two linear systems of equations ( 52) and [START_REF] Jansen | Topology optimization of fail-safe structures using a simplified local damage model[END_REF]. Note that because the proportional loading is increased at a constant rate, i.e.

∆u (n) E = ∆t (n) ∆t (n-1) ∆u (n-1) E , ( 55
)
the solution of the second linear system (53) can also be omitted by means of the recursion formula

µ (n) F = ∆t (n) ∆t (n-1) λ (n-1) F . ( 56
)

Extended BESO method

The extended BESO method recently developed in [START_REF] Xia | Evolutionary topology optimization of elastoplastic structures[END_REF] using an additional damping treatment on sensitivity numbers is adopted in this work. It has been shown that this treatment can improve the robustness and effectiveness of the method, especially in dealing with nonlinear designs in presence of dissipative effects. In this method, the sensitivity numbers associated with the relative ranking of the element sensitivities are chosen to determine material phase exchange. When uniform meshes are used, the sensitivity numbers for the considered objective are defined as following using the element sensitivity computed from ( 54)

α i = ∂ J ∂ρi η , for ρ i = 1 0 , for ρ i = 0. ( 57
)
in which η is a numerical damping coefficient (the same as the one applied in the Optimality Criteria method for density-based methods [START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF]). When η = 1, we recover the conventional sensitivity numbers for linear elastic designs [START_REF] Huang | Topology Optimization of Continuum Structures: Methods and Applications[END_REF]. In the presence of dissipative effects, the sensitivity numbers vary by several orders of magnitude resulting in instabilities of the topology evolution process, especially when removing certain structural branches. For this reason, the sensitivity numbers are damped in this work with "η = 0.5" as suggested in [START_REF] Fritzen | Topology optimization of multiscale elastoviscoplastic structures[END_REF][START_REF] Xia | Evolutionary topology optimization of elastoplastic structures[END_REF].

To avoid checkerboard patterns, sensitivity numbers are firstly smoothed by means of a filtering scheme [START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF]

α i = Ni j=1 w ij α j Ni j=1 w ij , ( 58
)
where w ij is a linear weight factor

w ij = max(0, r min -∆(i, j)), (59)
determined according to the prescribed filter radius r min and the element center-to-center distance ∆(i, j) between elements e and j. We note that the filter (58) is also responsible for material exchange from the matrix phase (ρ i = 0) to the inclusion phase (ρ i = 1) by attributing filtered sensitivity number values to design variables that are associated to the matrix phase.

Due to the discrete nature of the BESO material model, the current sensitivity numbers are needed to be averaged with their historical information to improve the design convergence (see [START_REF] Huang | Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method[END_REF])

α (l) i ← (α (l) i + α (l-1) i ) 2 . ( 60
)
In this work, the material volume fraction of reinforced inclusion phases is kept constant during the optimization process. Interfaces between different material phases or even the topology of the reinforced inclusion materials will be tailored through the redistribution of the quantitative inclusion phases. The above extended BESO method is adopted to solve the optimization problem by using the modified sensitivity numbers, which account for the whole fracturing process, involving crack nucleation, propagation and interaction until complete failure of the considered heterogeneous materials, so as to improve the fracture resistance of the specimens.

NUMERICAL EXAMPLES

In this section, several numerical examples are presented to demonstrate the potential of the proposed topology optimization framework. In all tests, regular meshes using quadrilateral bilinear elements are adopted, and plane strain condition is assumed. The same finite element discretization is adopted for both displacement and crack phase fields. The regularization parameter describing the width of smeared crack and interface is chosen as two times the finite element size = 2 e . The material parameters of each phase are taken as:

E i = 52 Gpa, E m = 10.4 Gpa, v i = v m = 0.3,
where the indices i and m correspond to the matrix and inclusion materials, respectively. These parameters are those of a mortar composed of a cement paste (matrix) and sand (inclusion). The toughness is g c = g I c = 1 × 10 -4 kN/mm and the interface fracture strength is chosen as t u = 10 -2 Gpa.

Design of a plate with one initial crack under traction

The problem geometry of this example is depicted in Fig. 4(a). The dimensions of the plate are 50 × 100 mm 2 , and the domain is uniformly discretized into 60 × 120 square shaped bilinear elements. The boundary conditions are as follows: on the lower end, the vertical displacements are fixed while the horizontal displacements are free, and the left bottom corner node is fixed in both directions. On the upper end, the horizontal displacements are free, while the vertical displacements are prescribed with monotonic displacement increments Ū = 0.005 mm during the simulation. The incremental loading process continues until the reaction forces is below a prescribed value indicating that the structure is completely broken. It is worth noted that during the crack propagation, interfacial damage can occur and interact with the propagation of the pre-existing matrix crack.

Fig. 4(b) is the initial guess design and consists of a single square inclusion occupying a volume fraction of 10% of the sample. In all next examples, the material volume fraction of the inclusion phase is maintained constant during the optimization process. By the extended BESO method, the inclusion phase will be redistributed based on their sensitivity numbers so as to improve the fracture resistance of the considered structure. The pre-existing crack notch is simulated by prescribing Dirichlet conditions on the crack phase field with d = 1 along the crack. The surrounding area of procedure has been employed for solving the coupled displacement-phase field problems formulated in section 2. Then, for one load increment, the number of sub-iterations, i.e. the number of times the displacement and phase field problems are solved alternatively has an effect on the solution. It is worth noting that most authors (see e.g. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], and also our own previous works, [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], among many others) only use one sub-iteration by assuming that the load increments are small enough, but thus requires a preliminary convergence study. In the following, we have studied the effects of using 2, 3, 4 and 5 sub-iterations in the staggered scheme on the load-displacement response of the structure. In Fig. 5, we denote by "FD" and "ID" Final and Initial Designs, respectively, and "s-it" means sub-iterations. We can observe from Fig. 5 (a) that the number of sub-iterations in the whole optimization procedure does have an effect on the final load-displacement curve but almost no effects on the final optimized shape (see Fig. 5 (b)-(i)). We note that for 4 sub-iterations, the load-displacement curve is roughly converged. With these observations in mind, we then conduct the whole optimization procedure for obtaining the final design with only 2 sub-iterations to reduce the computational costs, as it has been seen that the number of sub-iterations as small influence on the final shape of the inclusion. Then, we use the obtained final design and re-compute the loaddisplacement curve with 4 sub-iterations to avoid underestimating the fracture energy. In the present example, we can see from Fig. 5 that the total required fracture energies for complete failure is 3.97 mJ for the initial design and 6.04 mJ for the final design, which means that the final structure is 52% more resistant to fracture than the initial guess design. Detailed propagation of the phase field crack of the final design composite structure with one pre-existing crack notch subjected to incremental traction loads is given in Fig. 6. The initial crack propagates into the inner supporting structure and is blocked by the reinforced inclusion phase during the previous incremental loads. Then, two interface cracks nucleate and propagate along the upper and lower material interfaces. Inclusion materials which are redistributed to up and down on the left side try to prevent the vertical propagation of matrix crack. Finally, the interface and matrix cracks are intersected and propagated horizontally until the structure is fully broken.

Design of a plate without initial cracks for traction loads

The problem setting of this example is the same as in the last section except that there is no initial crack. The geometry of the plate is depicted in Fig. 7 (a), where the inclusion phase occupies a volume fraction of 5% of the sample. The optimized geometry of the inclusion phase is depicted in Fig. 7 (b). Detailed propagations of the phase field cracks of the initially and finally designed composite structures subjected to incremental traction load are given in Fig. 8 and Fig. 9, respectively. The cracks are firstly generated around the upper and lower material interfaces for the initial guess design, while it nucleates in the middle of the inner inclusion phase as well as in the interfaces for the new design. Then, similar interface cracks around the upper and lower material interfaces are generated. Finally, cracks propagate horizontally until the structure is fully broken in both cases. The final crack patterns, the total required fracture energies for complete failure, and the load-displacement curves for the initial and optimized designs are given in Fig. 10.

In this example, the structure with optimal design is 15% more resistant to fracture than the initial composite structure.

Design of a square plate without initial cracks in tensile loading

This example aims to design a square composite plate without initial crack subjected to uniaxial tension. The problem geometry of the square plate is depicted in Fig. 11 (a). The dimensions of the plate are 100 × 100 mm 2 , and the domain is uniformly discretized into 120 × 120 square bilinear elements. The boundary conditions are as follows: on the upper and lower end, the vertical displacements are fixed while the horizontal displacements are free. On the left and right end, the horizontal displacements are prescribed to an increasing uniform value of Ū = 0.005 mm during the simulation. The incremental loading process continues until the reaction force is below a prescribed value indicating that the structure is fully broken. The initial guess involves an inclusion representing 5% of the volume fraction as shown in Fig. 11 (a). The final design of the inclusion phase is shown in Fig. 11 (b). The width of the final design is larger than the initial guess to resist the x-directional tension, and the inclusion phase topology is changed. Detailed propagation of the phase field crack in the composite structures with optimal design is given in Fig. 12. Here, the interface cracks firstly initiate at the interface (see Fig. The left bottom corner node is fixed, while node at the right bottom corner the y-displacement is fixed and the x-displacement is free. For this case, an initial pre-existing crack is shown in Fig. 14 (b), and the initial guess design with the inclusion phase occupying 10% volume fraction of the domain area. The computation is performed with monotonic displacement increments of Ū = 0.01 mm until the reaction forces is below a prescribed criterion value. The displacements are prescribed along the y-direction while the displacement along x is free. The pre-existing crack is simulated by prescribing Dirichlet conditions with d = 1 along the crack. The surrounding area of the initial crack notch is treated as a non-designed region. Fig. 14 (c) shows the final design of inclusion topology. Detailed propagation of the phase field crack of the finally designed composite structure with one pre-existing crack notch subjected to three-point bending with monotonic displacement increments is shown in Fig. 15. The initial matrix crack firstly propagates vertically and is blocked by the reinforced inclusion phase materials. During the next steps, it tries to spread along the horizontal direction but is blocked by the inclusion material redistributed by using the proposed topological optimization method. Eventually, the matrix crack propagates along the loading direction until the structure is fully broken. The two load-displacement curves are compared in Fig. 16 (a). The final crack patterns as well as the total required fracture energies for complete failure for the initial and final designs are shown in Fig. 16 (b) and Fig. 16 (c), respectively. In this example, the fracture resistance of the structure with final design has been increased by 76% as compared to the initial design.

Design of a plate containing multiple inclusions

Finally, a plate containing 20 periodically distributed square inclusions is considered. The plate is modeled as a square domain whose length is 100 mm. The length of the inclusions is computed such that the volume fraction of inclusion phase is equal to 20%. The boundary conditions are the same as described in Section 4.3 and is shown in Fig. 17 y-displacements are fixed, while the x-displacements are free. On the left and right ends, the ydisplacements are free, while x-displacements are prescribed, with a uniform increasing value of Ū during the simulation. The computation is performed with monotonic displacement increments Ū = 0.005 mm until the structure is complete broken. Here again, the surrounding region of the initial crack is assumed to be non-designable to avoid nonphysical designs. The final design is presented in Fig. 17 (c). Detailed propagation of the phase field crack of the final design of the considered composite structure containing multiple inclusions and one initial crack subjected to tensile loads is shown in Fig. 18. It can be observed that the initial crack propagates vertically and new cracks around intermediate interfaces of the structure are generated. The inclusion materials are redistributed to prevent further propagation of critical cracks. In the subsequent loading steps, the pre-existing crack meets the reinforced inclusion materials and is blocked at the beginning. Then, both matrix and interface cracks propagate toward the vertical direction as well as along the material interfaces until the structure is fully broken. The comparisons of final crack paths, the total required fracture energies for complete failure and load-displacement curves between initial and optimal designs are provided in Fig. 19. In this example, the resistance to fracture has been increased by 53% in comparison with the initial design.

It is worth noting that in such optimization process, the crack pattern can involve multiple cracks, whose trajectories and number fully depends on the geometry of the inclusion phases. Then, it is very difficult to guess a priori the optimal geometry of inclusions to increase the fracture resistance. Finally, the main advantages of the phase field method in topological optimization process are summarized as follows: (a) it is easy to couple the phase field method with a topological optimization algorithm as a fixed mesh can be employed for the fracture simulation; (b) initiation of cracks can be included in the analysis, which would be not possible with techniques like XFEM [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF]; (c) including interfacial damage is simple in the present framework, using e.g. extensions of the phase field as 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]. Using cohesive elements in that case would make the analysis much more complex when defining the sensitivity with respect to the material density in the topological algorithm. Then, the present framework seems to be very promising to design new composite materials with enhanced fracture resistance.

CONCLUSION

In this work, a topology optimization framework for improving fracture resistance of two-phase composites through a redistribution of reinforced inclusion phases, considering interactions between bulk brittle fracture and interfacial damage, has been developed. One particular merit of the adopted phase field method is the regularized description of discontinuous fields, which avoids the burden of remeshing during the crack propagation which is fully adapted to topology optimization. The sole and unique phase field is employed to describe both bulk brittle fracture and interface cracking, and thus allows interaction between the two types of cracks. Before performing the topology optimization, a computationally efficient adjoint sensitivity formulation is derived to account for the whole fracturing process, involving cracks nucleation, propagation and interaction, either from the interfaces and then through the solid phases, or the composite. Based on the determined sensitivity numbers, the constant amount of reinforced inclusion materials are redistributed by using the extended bi-directional evolutionary structural optimization method. Several benchmark tests have been presented to demonstrate the potential of the proposed design framework. It has been shown that significant improvement of the fracture resistance of the considered composite structures can be achieved for final designs accounting for full failure when compared to the initial guess designs. To our best knowledge, the topology optimization for fracture resistance taking into account interactions between interfacial damage and bulk brittle fracture for complete fracturing process has been done for the first time in this work. The presented method then provides a very promising design tool to improve the fracture resistance of heterogeneous materials where both interfacial damage and matrix crack propagation occur, and could constitute a basis for improving other physical properties involving interfacial damages.

APPENDIX

Linearization of the displacement problem

Even though the phase field problem is linear in the staggered scheme, i.e. for a fixed value to u, it should be mentioned that for a fixed crack phase field value d, the mechanical problem ( 25) is nonlinear since the computation of eigenvalues of ε e in [START_REF] Bernard | Damage growth modeling using the thick level set (TLS) approach: Efficient discretization for quasi-static loadings[END_REF] and the interface cohesive model in [START_REF] Sethian | Structural boundary design via level set and immersed interface methods[END_REF]. A linear procedure to solve this nonlinear problem by the Newton method is introduced in the following. From ( 25) and ( 29), the balance equation can be rewritten as

R = Ω σ e : ε e (δu)dΩ + Ω γ β (x)t(w, α) • δwdΩ - Ω f • δudΩ - ∂Ω F F • udΓ = 0, (61)
where ε e (δu) = ∇ s δu -n⊗ S wγ β . In a standard Newton method, the displacements are updated for each loading by solving the following tangent equation:

D ∆u R(u k , d) = -R(u k , d) = 0, (62)
where u k is the displacement solution from the k-th iteration. The displacements at the current iteration are given by u k+1 = u k + ∆u.

with ∆w(x) = h∇∆u(x) ∇φ(x) ||∇φ(x)|| [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] and

∂[σ e ] ∂[ε e ] = C(u, d) = [(1 -d) 2 + k] λR + [1] T [1] + 2µP + + λR -[1] T [1] + 2µP - (66)
where [σ e ] and [ε e ] are the vector forms for the second-order tensors σ e and ε e , respectively and C is the matrix form corresponding to the fourth-order tensor C.

FEM discretization of the phase field problem

A staggered solution procedure is adopted in this work, where the phase field and the mechanical problems are solved alternatively. At each increment, given the displacement field from the mechanical problem, the phase field problem is linear using a shifted algorithm (see more details in [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF]). Using FEM, the phase field and phase field gradient in one element are approximated by:

where de are nodal phase field values in one element, N d (x) and B d (x) are matrices of shape functions and of shape functions derivatives associated to phase field variable, respectively. Introducing the above FEM discretization into the weak form [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], the following linear discrete system of equations can be obtained:

K d d = F d (68)
where (70)

K d = Ω ( gc (

FEM discretization of the displacement problem

Similarly, the displacement field and incremental displacement field can be expressed using the FEM approximations where N denotes the matrix of shape functions associated to displacement variables, u e and ∆u e are nodal displacement components and nodal incremental displacement components in one element. Furthermore, we have

[ε](∆u) = Bu∆u e (72)
where Bu is a matrix of shape function derivatives. From [START_REF] Lamon | The influence of the interphase and associated interfaces on the deflection of matrix cracks in ceramic matrix composites[END_REF], the diffuse jump approximation vector and its incremental counterparts can be discretized as

and n 1 and n 2 are the x-and y-components of the normal vector. The smoothed jump strain at the interfaces is defined by

[ ε] =   ε11 ε22 √ 2 ε12   = γ β (x)    ω 1 n 1 ω 2 n 2 1 √ 2 (ω 1 n 2 + ω 2 n 1 )    . (75)
Then it yields

[ ε(∆u)] = hγ β (x)M Bu∆u e (76)
with

M =    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    . ( 77
)
With the above FEM discretization, the linear tangent problem reduces to the following linear system of algebraeic equations (80)

K tan ∆ũ = -R(ũ k ), (78) where
6.3.1. Overall numerical algorithm for crack propagation The overall algorithm is described as follows.

1. Set the initial displacement field u 0 (x), the phase field d 0 (x), and the strain-history function H 0 .

2. Compute the phase field β(x).

3. For all loading increments: (at each time t n+1 ), given dn, un, and Hn(x).

(a) Compute the history function H(t n+1 ) according to [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF]. ii. Compute ∆u e k+1 by (78). iii. Update u k+1 = u k + ∆u e k+1 . iv. (.) n+1 → (.)n and go to (a).

End End

Figure 1 .

 1 Figure 1. (a) structure containing both bulk cracks and interface cracks, possibly occurring at the interfaces, (b) mechanical response of the damageable structure.

Figure 2 .

 2 Figure 2. Illustration of the regularized representation of cracks and interfaces: (a) a solid containing interfaces and cracks; (b) regularized representation of the interfaces; (c) regularized representation of the cracks.

Figure 3 .

 3 Figure 3. Illustration of the cohesive model for the interfaces.

F

  are sought as following by solving the adjoint systems with the prescribed values λ

Figure 4 .

 4 Figure 4. A plate with one pre-existing crack notch subjected to incremental traction loads: (a) problem geometry; (b) initial guess design; (c) final design.

Figure 5 .

 5 Figure 5. Fracture resistance comparison of two composite structures with one initial crack subjected to incremental traction loads: (a) load-displacement curves with different numbers of sub-iteration; (b) ID with 2 s-it; (c) FD with 2 s-it; (d) ID with 3 s-it; (e) FD with 3 s-it; (f) ID with 4 s-it; (g) FD with 4 s-it; (h) ID with 5 s-it; (i) FD with 5 s-it.

Figure 6 .Figure 7 .

 67 Figure 6. Crack propagation of the final design composite structure with one initial crack subjected to incremental traction loads: (a) Ū = 0 mm; (b) Ū = 0.05 mm; (c) Ū = 0.06 mm; (d) Ū = 0.08 mm; (e) Ū = 0.095 mm.

Figure 8 .Figure 9 .Figure 10 .

 8910 Figure 8. Crack propagation of the initial design composite structure with one initial crack subjected to incremental traction loads: (a) Ū = 0 mm; (b) Ū = 0.025 mm; (c) Ū = 0.065 mm; (d) Ū = 0.085 mm; (e) Ū = 0.105 mm.

Figure 11 .

 11 Figure 11. A square plate without initial crack subjected to uniaxial tension: (a) geometry of the initial design; (b) final design.

Figure 12 .Figure 13 .

 1213 Figure 12. Crack propagation of the initial design composite structure with one initial crack subjected to uniaxial tension: (a) Ū = 0 mm; (b) Ū = 0.025 mm; (c) Ū = 0.04 (d) Ū = 0.08 mm.

Figure 14 .Figure 15 .

 1415 Figure 14. A plate with one initial crack subjected to three-point bending: (a) geometry and boundary condition; (b) geometry of the initial design and crack; (c) final design.

  12 (b)). Subsequently, the left and right interface cracks propagate vertically and merge. The complete fracture patterns and design objective values are shown in Fig.13 (b) for the initial design and in Fig.13 (c) for the optimal design. Both responses are compared in Fig.13 (a). Here, the fracture resistance of the final design structure has been increased by 33% as compared to the initial design.

Figure 16 .

 16 Figure 16. Fracture resistance comparison of two composite structures with one initial crack subjected to three-point bending.

  (a): on the lower and upper ends, the

Figure 17 .

 17 Figure 17. A plate containing multiple inclusions and one initial crack for tensile loads: (a) geometry and boundary condition; (b) geometry of the initial design and crack; (c) final design.

Figure 18 .Figure 19 .

 1819 Figure 18. Crack propagation of the final design composite structure with one initial crack subjected to uniaxial tension: (a) Ū = 0.01 mm; (b) Ū = 0.02 mm; (c) Ū = 0.05 mm; (d) Ū = 0.09 mm.

  ), we obtain D ∆u R(u k ) = Ω ∂σ e ∂ε e : ε e (∆ε) : ε e (δε) + Ω ∂t(w) ∂w : ∆w : δwdΩ

  d(x) = N d (x)d e , ∇d(x) = B d (x)d e ,

  1 -β) + 2H)N T d N d + (1 -β)gc B T d B d dΩ

u=

  Nu e , ∆u = N∆u e (71)

w

  = hN Buu e , ∆w = hN Bu∆u e

  K tan = Ω [B T u -hγ β (x) BT u M T ]C(x)[Bu -hγ β (x) BuM]dΩ + Ω h 2 γ β (x) BT u N T K I N BudΩ,(79)andR = Ω [B T u -hγ β (x) BT u M T ]C(x)[Bu -hγ β (x) BuM](u e ) k dΩ + Ω hγ β (x) BT u N T t(w k )dΩ + Ω fN T dΩ + Ω FN T dΓ.

( b )

 b Compute the crack phase field d n+1 (x) by solving linear problem (68). (c) Compute u n+1 (x) : i. Initialize u k = un While ||∆u k+1 || > , 1:

  , and Γ d (d) represents the total crack length per unit area in 2D and total crack area of per unit volume in 3D. In[START_REF] Cadman | On design of multi-functional microstructural materials[END_REF], Γ d (d) is defined by

ACKNOWLEDGEMENTS

This work was supported by State Key Program of National Natural Science Foundation of China (61232014) and Institut Universitaire de France (IUF). L. Xia thanks the financial support of the National Natural Science Foundation of China (51705165, 51790171, 5171101743 ).