Topology optimization of particle-matrix – Global Homework Experts

Received: 29 October 2017 Revised: 26 February 2018 Accepted: 29 March 2018

DOI: 10.1002/nme.5818
R E S E A R C H A R T I C L E
Topology optimization of particle-matrix composites for
optimal fracture resistance taking into account interfacial
damage
Daicong Da1,2 Julien Yvonnet2 Liang Xia3 Guangyao Li1
1State Key Laboratory of Advanced Design
and Manufacturing for Vehicle Body,
Hunan University, Changsha, China
2Laboratoire Modélisation et Simulation
Multi Echelle (MSME), Université
Paris-Est, CNRS UMR 8208,
Marne-la-Vallée, France
3State Key Laboratory of Digital
Manufacturing Equipment and
Technology, Huazhong University of
Science and Technology, Wuhan, China
Correspondence
J. Yvonnet, Laboratoire Modélisation et
Simulation Multi Echelle (MSME),
Université Paris-Est, CNRS UMR 8208,
5 bd Descartes, 77454 Marne-la-Vallée,
France.
Email: [email protected]
G.Y. Li, State Key Laboratory of Advanced
Design and Manufacturing for Vehicle
Body, Hunan University, Changsha
410082, China.
E-mail: [email protected]
Funding information
National Natural Science Foundation of
China, Grant/Award Number: (61232014);
Institut Universitaire de France (IUF);
National Natural Science Foundation of
China, Grant/Award Number: (51705165,
51790171, 5171101743)
Summary
This paper presents a topology optimization framework for optimizing the fracture resistance of two-phase composites considering interfacial damage interacting with crack propagation through a redistribution of the inclusions phase. A
phase field method for fracture capable of describing interactions between bulk
brittle fracture and interfacial damage is adopted within a diffuse approximation
of discontinuities. This formulation avoids the burden of remeshing problem
during crack propagation and is well adapted to topology optimization purpose.
Efficient design sensitivity analysis is performed by using the adjoint method,
and the optimization problem is solved by an extended bidirectional evolutionary structural optimization method. The sensitivity formulation accounts
for the whole fracturing process involving crack nucleation, propagation, and
interaction, either from the interfaces and then through the solid phases, or
the opposite. The spatial distribution of material phases are optimally designed
using the extended bidirectional evolutionary structural optimization method to
improve the fractural resistance. We demonstrate through several examples that
the fracture resistance of the composite can be significantly increased at constant
volume fraction of inclusions by the topology optimization process.
KEYWORDS
BESO, cracks, interfacial damage, phase field method, topology optimization
1 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
604 Copyright © 2018 John Wiley & Sons, Ltd. wileyonlinelibrary.com/journal/nme Int J Numer Methods Eng. 2018;115:604–626.

order now

Fracture
Energy
Load
Displacement
Peak Load
Matrix Inclusions
Interfaces
Bulk
cracks
(A) (B)
FIGURE 1 A, Structure containing both bulk cracks and interface cracks, possibly occurring at the interfaces; B, Mechanical response of
the damageable structure [Colour figure can be viewed at wileyonlinelibrary.com]
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, eg, the work of Quan et al
1) 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 eg, a review in the work of Cadman et al
2 and the
references therein), and more recently, in the nonlinear context (see, eg, the work of Fritzen et al
3). 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 (eg, finite elements). Gu et al
4 used a modified greedy optimization algorithm for
composites made up of soft and stiff building blocks to improve material toughness. San and Waisman
5 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 Xia et al,
6 topology optimization for maximizing the fracture resistance of quasibrittle composites
has been introduced by combining the phase field method and a gradient-based bidirectional evolutionary structural
optimization (BESO) algorithm. However, in the aforementioned work, the crack propagation resistance was only evaluated on the basis of phase distribution. In most heterogeneous quasibrittle materials (eg, ceramic matrix composites and
cementitious materials), the interfacial damage plays a central role in the nucleation and propagation of microcracks.
710
The main objective of the present work is to extend the framework developed in the work of the aforementioned author6
for defining through topological optimization the optimal phase distribution in a quasibrittle composite with respect to
fracture resistance, taking into account crack nucleation both in the matrix and in the interfaces, as described in Figure 1.
To our best knowledge, such study is investigated here for the first time.
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 (FEM),
11,12
thick level-set method,13,14 and phase field method,1517 among the most recent popular techniques, have been introduced
to investigate this topic.
Based on the pioneer works of Francfort and Marigo,
15 the phase field method makes use of a variational principle
framework for brittle fracture
1820 and of a regularized description of the discontinuities related to the crack front.21,22
The method has been adapted to a convenient algorithmic setting by Miehe et al.23 The main advantages of the method
are (i) 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;
and (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 the work of Nguyen et al
24 to take into account both bulk brittle fracture and
interfacial damage.

Interface 0
1
(A) (B) (C)
Inclusion
Matrix
Crack
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 [Colour figure can be viewed at wileyonlinelibrary.com]
During the past decades, the research in topology optimization witnessed tremendous development since the seminal
paper by Bensøe and Kikuchi.
25 So far, various topology optimization methods have been proposed, eg, density-based
method,
2628 evolutionary procedures,29,30 and level-set methods.3133 The extended BESO method developed for the design
of elastoplastic structure in the work of Fritzen et al
3,34 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 crack nucleation, propagation, and merging of microcracks until complete failure of the specimen. By using
discrete design variables, the evolutionary structural optimization–type methods omit naturally the definition of additional pseudorelationships between fictitious materials and fracture toughness, resulting in a clear physical interpretation
and algorithmic advantages.
6,3537 While there do exist numerous researches on topology optimization including material
interface behavior (see, eg, other works
3841) and with the enforcement of stress constraints,4248 local buckling constraint49
or damage constraints,5055 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 this 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 and 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.
2 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
24 and allows simulating the initiation and propagation of a multiple-crack network in heterogeneous
microstructures. The main concepts and details are recalled in the following.
2.1 Regularized representation of discontinuous field
Let 𝛺 Rd be an open domain describing a solid with external boundary 𝜕𝛺. The solid contains internal material
interfaces between different phases, collectively denoted by
𝛤 I. During the loading, cracks may propagate within the
solid and can pass through the material interfaces as depicted in Figure 2A. The crack surfaces are denoted by
𝛤 c. In this
work, we adopt smeared representations of both cracks and material interfaces, ie, 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 nonevolving sharp crack defined on a surface 𝛤 c, a regularized (smeared) representation of the corresponding damage d(x) (see Figure 2C) can be obtained by solving the following equations (see the work of Miehe et al16 for
more details):
⎧⎪⎨⎪⎩
d(x) – 𝓁22d(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 Equation (1) is the Euler-Lagrange equation corresponding to the variational problem

d(x, t) = Arg {dinf Sd𝛤 d(d)} , 𝛤 d(d) = 𝛺

𝛾d(d)d𝛺, (2)
where
Sd = {d|d(x) = 1, x 𝛤} and 𝛤 d(d) represents the total crack length per unit area in 2D and total crack area of
per unit volume in 3D. In (2),
𝛤 d(d) is defined by

𝛾d(d(x)) = 1
2
𝓁
d(x)2 +
d(x) · ∇d(x). (3)

𝓁
2
When
𝓁 0, the aforementioned variation principle leads to the corresponding one with a sharp crack description in
the
𝛤convergence sense. The details of definition (3) can be found, eg, in the work of Miehe et al.23
In the work of Nguyen et al,24 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) – 𝓁22𝛽(x) = 0, in 𝛺
𝛽(x) = 1,
𝛽(x) · n = 0,
on 𝛤 I
on 𝜕𝛺.

(4)
Equation (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
2
l𝛽(x)2 + 𝓁2 𝛽(x) · ∇𝛽(x). (6)
For
𝓁 0, the aforementioned 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
24 using Taylor expansion of the displacement field around a point x located on the interface:
u(x)w(x) = u (x + h2 n𝐼) u (x h2 n𝐼) = h∇(u(x))nI, (7)
where
w(x) denotes the smoothed displacement jump approximation and nI 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 the work of Nguyen et al,24
a level-set technique has been proposed (see more detailed in the aforementioned paper). In the work of Verhoosel
and de Borst,
56 another definition using the possible modification of the interface by the damage related to bulk cracks
was introduced.

2.2 Energy functional
The following total energy functional formulated is introduced for the solid body related to both cracks and interfaces:
E = 𝛺Wue (𝜺e(u, 𝛽), d) d𝛺 + 𝛺[1 𝛽(x)]gc𝛾d(d)d𝛺 + 𝛺𝜓I(w)𝛾𝛽(𝛽)d𝛺, (8)
where
gc is the toughness and 𝜓I is a strain density function depending on the displacement jump across the interface
𝛤 I. According to the aforementioned equation, 𝜺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 the
work of Verhoosel and de Borst
56), in which case we recover the energy functional for a cracked body without considering
the interface behavior
E = 𝛺Wue (𝜺e(u), d) d𝛺 + 𝛺gc𝛾d(d)d𝛺. (10)
From Equation (8), the free energy
W can be identified as
W = We
u
(𝜺e(u, 𝛽), d) + [1 𝛽(x)]gc𝛾d(d) + 𝜓I(w, 𝜶)𝛾𝛽(𝛽). (11)
Following the work of Nguyen et al,
24 the elastic energy Wue is defined as
We
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± =
Di=1
𝜀i±ni ni. (15)
The decomposition (11) ensures that damage is only induced by traction, which removes the complications related to
autocontact in the crack. In (14) and (15),
x± = (x±|x|)∕2. 𝜀i and ni 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

ḋ 0, = –𝜕W
𝜕d .
A threshold function F() is assumed such that no damage occurs in the form
F() =  ≤ 0.
(16)
(17)

The principle of maximum dissipation requires that dissipation ḋ must be maximum under the constraint,
ie,
ḋ > 0 and F = 0. Therefore, it yields
F = –𝜕W
𝜕d = –
𝜕We
u
𝜕d – (1 𝛽)gc𝛿𝛾(d) = 0 (18)
with the functional derivative
16
𝛿𝛾(d) = d
𝓁
𝓁Δd. (19)
It follows that, when
ḋ > 0, then
2(1 d)𝜓e+ + (1 𝛽)gc𝛿𝛾(d) = 0, (20)

To handle loading and unloading, the strain history function adopted in the works of Miehe et al16 and Nguyen et al24

is employed here (x
, t) = max
𝜏∈[0,t]
{𝜓e+(x, t)}
(21)

and (20) is substituted by
2(1 d)+ (1 𝛽)gc𝛿𝛾(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)– (1 𝛽)gc
𝓁 (d 𝓁22d) = 0, in 𝛺

d(x) = 1,
d(x) · n = 0,
on 𝛤 c
on 𝜕𝛺.

(23)
The associated weak form is obtained as (see the work of Nguyen et al
24)
𝛺 {(2+ [1 𝛽]g𝓁c ) d𝛿d + [1 𝛽]gc𝓁d · ∇(𝛿d)} d𝛺 = 𝛺2𝛿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𝜺eue 𝜺e(𝛿u)d𝛺 + 𝛺 𝜕𝜓𝛿Iw(w) · 𝜕w𝛾𝛽(𝛽)d𝛺 = 𝛺f · 𝛿ud𝛺 + 𝜕𝛺
F
F̄ · 𝛿ud𝛤 = 𝛿Wext. (25)
In the absence of body forces, Equation (25) can be rewritten 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 nI and 𝛿w =
h∇(𝛿u)nI. Using 𝝈en = t, the aforementioned equation can be further rewritten as

𝛺𝝈e {𝜺e(𝛿u) + n ⊗ 𝛿w𝛾𝛽(𝛽) – ∇s𝛿u} d𝛺,
which is satisfied for an admissible strain field
(27)

𝜺e = ∇su nSw𝛾𝛽, (28)
where
(∇su)i j = (ui, j + uj,i)∕2 and (nSw)i𝑗 = (niw𝑗 + win𝑗). From (13), 𝜺̄ can be identified as 𝜺̄ = nSw𝛾𝛽.
With the aforementioned 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𝜺e1 + 2𝜇𝜺e. (29)
2.3.3 Model for interface damage
The general form of the traction vector t(w) in Equation (26) is given by
t(w) = [tn, tt]T, (30)
where
tn and tt denote, respectively, normal and tangential parts of the traction vector t across the interface oriented
by its normal
nI. In the presented work, a simplified nonlinear elastic cohesive model is used. We have shown that the
framework proposed in the work of Nguyen et al
24 allows avoiding the use of internal variables related to the interface

FIGURE 3 Illustration of the cohesive model for the interfaces
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, ie,
t(w) · nI = tn, the cohesive law can be written as

tn = gcI (w𝛿nn ) exp (w𝛿nn ) ,
where wn = w · nI. For this model, we obtain KI = 𝜕t(w)
(31)

𝜕w
. The relationship between 𝛿n, the toughness gI c, and the fracture
strength
tu is given by 𝛿n = gI c∕(tue), with e = ex𝑝(1) (see Figure 3).
Details about the linearization of the displacement problem and about FEM discretizations and numerical implementation are provided in the Appendix.
3 TOPOLOGY OPTIMIZATION METHOD
In this section, we present the topological optimization method based on BESO.30 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.
3.1 Model definitions
The considered domain 𝛺 is discretized into Ni finite elements, and each element i is assigned with a topology design
variable
𝜌i. Within the framework of the BESO method, the discrete design variables are defined as

𝜌i = 0 or 1, i = 1, 2, , Ni.
Following the multiple material interpolation model in the works of Huang and Xie57 and Da et al,58 we have
Ei = 𝜌iEinc + (1 𝜌i)Emat,
(32)
(33)

where Einc and Emat 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 (33), 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.
5962 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, ie,
J 1
2
nload
n=1
(fext (n) + fext (n1))TΔu(n), (34)

and will be used as our definition to fracture resistance in the following. According to the aforementioned equation, nload
is the total number of displacement increments, Δu(n) is the nth nodal displacement components, and fext n is the external
nodal force at the
nth 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 = 0
V(𝝆) = 𝜌ivi = Vreq
𝜌i = 0 or 1, i = 1, · · ·, Ni.
(35)
In the aforementioned equation,
vi is the volume of ith element; V(𝝆) and Vreq are the total and required material volumes,
respectively; and
denotes the nodal residual force
= fext fint. (36)
In (36),
fint 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
fint =
Nii=1
𝜌i𝛺
i
BT𝝈ed𝛺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, eg, the work of Bendsøe and
Sigmund
26), resulting in a clear physical interpretation.
3.2 Sensitivity analysis
In order to perform the topology optimization, the sensitivity of the objective function J with respect to topology design
variable
𝝆 must be computed. The adjoint method is employed to derive the sensitivity analysis. First, two Lagrangian
multipliers
𝜇(n) and 𝝀(n) of the same dimension as the vector of unknown u are introduced to enforce zero residual at
time
tn1 and tn for each term of the total mechanical work (34). Then, the objective function J can be rewritten in the
following form without modifying the original objective value as:

̂J = 1
2
n=1
{(fext (n) + fext (n1))Δu(n) + (𝝀(n))T(n) + (𝝁(n))T(n1)}
nload T
. (38)

Due to the asserted static equilibrium, the residuals (n) and (n1) must vanish. The objective value is thus invariant
with respect to the values of the Lagrangian multipliers
𝝀(n) and 𝜇(n) (n = 1, , nload), ie,
̂J (𝝆; {𝝀(n), 𝝁(n)}

n=1, ,nload) = J (𝝆) .
This equivalence also holds for the sensitivity with respect to changes of the topology design variable 𝜌i on element i
(39)

𝜕̂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
into essential (index E, associated with Dirichlet boundary conditions) and free (index F, remaining degrees of freedom)
entries. For a vector
w and a matrix M, we have

w [ wwEF ] and M [ MMEE FE MMEF FF ] . (41)
In the current context, the displacements uE on the Dirichlet boundary are prescribed. They are independent of the
current value of 𝝆. This implies that
𝜕Δu
𝜕𝜌i
=
𝜕
𝜕𝜌
i [ ΔΔuuEF ] = [ 𝜕u0F}∕𝜕𝜌i ]
(42)

 

holds for arbitrary increment loads, ie, for u = u(n) or u = u(n1). The components fext,F of the force vector fext vanish at
all times and the only (possibly) nonzero components are the reaction forces
fext,E
f (n)
ext = [ fext (0n),E ] . (43)
Equations (42) and (43) imply
(fext (m))T 𝜕Δ𝜕𝜌 ui(n) = 0. (44)
Hence, for arbitrary time indices
n, m = 1, , nload, we have
𝜕
𝜕𝜌
i ((fext (m))TΔu(n)) = (𝜕𝜕𝜌 fext (mi ) )TΔu(n). (45)
Therefore, the derivative of the modified design objective in (38) can be rewritten as
𝜕̂J
𝜕𝜌i

=
12
nload
n=1
⎧⎪⎨
(
𝜕f
𝜕𝜌 i )Δu(n) + (𝝀(n))T 𝜕+ (𝝁(n))T 𝜕𝜕𝜌 i⎫ ⎪ ⎬
⎪⎩
Recall the balance equation at each load time increment in (A1); the derivatives of


.
(46)
𝜕𝜌 i 𝜕𝜌 i
ext ext
(n) 𝜕f(n1) T (n)

+ (n1) (m) at the equilibrium of the mth load
increment with respect to
𝜌i can be expanded as
𝜕m
𝜕𝜌i
=
𝜕f m
ext
𝜕𝜌i
𝛺
i
BT(𝜎e)(m)d𝛺i K( tan m) 𝜕Δu(m)
𝜕𝜌i
, (47)
where
K(m)
tan = –
𝜕(m)
𝜕u(m) (48)
is the tangent stiffness matrix of the nonlinear mechanical system at the balance equation of the
mth load increment.
With the expression (47), (46) can be reformulated as
𝜕Ĵ
𝜕𝜌
i
=
12
nload
n=1
⎧⎪⎨⎪⎩
(
𝜕𝜕𝜌 fext (ni) )T (Δu(n) + 𝝀(n)) + (𝜕f𝜕𝜌 ext (ni 1) )T (Δu(n) + 𝝁(n))
(𝝀(n))T (𝛺iBT(𝝈e)(n)d𝛺i + K( tan n) 𝜕Δ𝜕𝜌 ui(n) )
(𝝁(n))T (𝛺iBT(𝝈e)(n1)d𝛺i + K( tan n1) 𝜕Δ𝜕𝜌 u(ni 1) )} .
(49)
As aforementioned, 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 (43), the first two terms in (49) can be omitted by setting
𝝀(n)
E = -Δu( En) and 𝝁( En) = -Δu( En). (50)
Accounting further for the structure of the sensitivities of
u in (42) and for the symmetry of the stiffness matrices, we have
𝜕Ĵ
𝜕𝜌
i
=
12
nload
n=1
{(𝝀(n))T𝛺iBT(𝝈e)(n)d𝛺i (𝝁(n))T𝛺iBT(𝝈e)(n1)d𝛺i
(K( tan n) ,FE𝝀( En) + K(ntan ) ,FF𝝀( Fn))T 𝜕Δ𝜕𝜌 ui( Fn)
(K( tan n,1FE ) 𝝁( En) + K(ntan 1),FF𝝁( Fn))T 𝜕Δ𝜕𝜌 u( Fni 1) } .
(51)

To avoid the evaluation of the unknown derivatives of u(n)
F and u( Fn1), ie, eliminating the last two lines of (51), the
values of
𝝀(n)
F and 𝝁( Fn) are sought as following by solving the adjoint systems with the prescribed values 𝝀( En) = -Δu( En) and
𝝁
(n)
E = -Δu( En) at the essential nodes:
𝝀(n)
F = (K(ntan ) ,FF)1

K(i)tan,FEΔu( En), (52) and
𝝁
(n)
F = (K(ntan 1),FF)1K(ntan 1),FEΔu( En).
(53)

The two relations (52) and (53), together with (50), fully determine the values of the Lagrange multipliers 𝝀(n) and 𝜇(n).
Finally, the objective function
𝜕̂J𝜕𝜌i can be computed via
𝜕Ĵ
𝜕𝜌
i
= –
12

nload
n=1
{(𝝀(n))TBT(𝝈e)(n)d𝛺i + (𝝁(n))TBT(𝝈e)(n1)d𝛺i} .
(54)
𝛺i 𝛺i

The computation of the sensitivity consists in solving two linear systems of Equations (52) and (53). Note that, because
the proportional loading is increased at a constant rate, ie,
Δu(n)
E =
Δ
t(n)
Δt(n1) Δu( En1), (55)
the solution of the second linear system (53) can also be omitted by means of the recursion formula
𝝁
(n)
F =
Δ
t(n)
Δt(n1) 𝝀( Fn1). (56)
3.3 Extended BESO method
The extended BESO method recently developed in the work of Xia et al34 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 the
following using the element sensitivity computed from (54):

𝛼i =
⎧⎪
(
𝜕𝜌 𝜕̂Ji )
(57)

⎨⎪⎩
𝜂, for 𝜌i = 1
0
, for 𝜌i = 0,
in which 𝜂 is a numerical damping coefficient (the same as the one applied in the optimality criteria method for
density-based methods
63). When 𝜂 = 1, we recover the conventional sensitivity numbers for linear elastic designs.64 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 the works of Fritzen et al3 and Xia et al.34
To avoid checkerboard patterns, sensitivity numbers are first smoothed by means of a filtering scheme63

𝛼i =

N 𝑗=i 1 wi
, (58)

N 𝑗=i 1 wi𝑗𝛼𝑗
𝑗
where wi j is a linear weight factor
wi𝑗 = max(0, rmin – Δ(i, 𝑗)), (59)
determined according to the prescribed filter radius
rmin 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.

(A)
Matrix
Inclusion
Crack

(B)
15.83 20
25
(C)
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 [Colour figure can be viewed at wileyonlinelibrary.com]
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 the work of Huang and Xie
30)
𝛼
(l)
i
(𝛼i(l) + 𝛼i(l1))
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 aforementioned 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.
4 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 Ei = 52 Gpa, Em = 10.4 Gpa, and vi = vm = 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
gc = gI c = 1 × 104 kN/mm and the interface
fracture strength is chosen as
tu = 102 Gpa.
4.1 Design of a plate with one initial crack under traction
The problem geometry of this example is depicted in Figure 4A. The dimensions of the plate are 50 × 100 mm2, 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. During the crack propagation, interfacial damage can occur and interact with the propagation of the
pre-existing matrix crack.
Figure 4B 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 the initial crack notch (up to 2 times of the length scale parameter 𝓁) is treated as a nondesignable region to avoid nonphysical
designs with the inclusion material added within the already existing crack (see a discussion in the work of Xia et al
6). The
final structural topology of the inclusion phase is shown in Figure 4C. It can be observed that the material on the right
side of the reinforcement inclusion moves up and down on the left side, and holes are generated to tailor the topology
of the inclusion phases. The fracture resistance improvement of the resultant composite structure is evaluated by comparing the initial and new designs response in Figure 5. In the present work, a staggered procedure has been employed
for solving the coupled displacement-phase field problems formulated in Section 2. Then, for one load increment, the
number of subiterations, ie, the number of times the displacement and phase field problems are solved alternatively, has
an effect on the solution. Most authors (see, eg, the works of Miehe et al
16 and Nguyen et al,24,65 among many others)
only use one subiteration 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 two, three, four, and five subiterations in the staggered
scheme on the load-displacement response of the structure. In Figure 5, we denote by “FD” and “ID” final and initial
designs, respectively, and “s-it” means subiterations. We can observe from Figure 5A that the number of subiterations 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 Figure 5B to I). We note that, for four subiterations, the load-displacement curve is roughly
initial design
J = 5.30 mJ
final design
J = 7.89 mJ
(A) (B) (C)
(D) (E) (F) (G) (H) (I)
initial design
J = 4.48 mJ
initial design
J = 3.97 mJ
initial design
J = 3.78 mJ
final design
J = 6.73 mJ
final design
J = 6.04 mJ
final design
J = 5.61 mJ
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 subiteration; B, Initial design (ID) with 2 subiterations (s-it); C, Final design (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

(A) (B) (C) (D) (E)
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
converged. With these observations in mind, we then conduct the whole optimization procedure for obtaining the final
design with only two subiterations to reduce the computational costs, as it has been seen that the number of subiterations has small influence on the final shape of the inclusion. Then, we use the obtained final design and recompute the
load-displacement curve with four subiterations to avoid underestimating the fracture energy. In the present example, we
can see from Figure 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 Figure 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.
4.2 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 Figure 7A where the inclusion phase occupies a volume fraction of 5
% of the sample. The
optimized geometry of the inclusion phase is depicted in Figure 7B. Detailed propagations of the phase field cracks of the
initially and finally designed composite structures subjected to incremental traction load are given in Figures 8 and 9,
100
(A)
Matrix
Inclusion
(B)
15.83
12.5
20
50
FIGURE 7 A plate without initial crack subjected to incremental traction loads. A, Geometry of the initial design; B, Final design [Colour
figure can be viewed at wileyonlinelibrary.com]

(A) (B) (C) (D) (E)
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
(A) (B) (C) (D) (E)
FIGURE 9 Crack propagation of the initial design composite structure with one initial crack subjected to incremental traction loads. A,
Ū = 0 mm; B, Ū = 0.035 mm; C, Ū = 0.055 mm; D, Ū = 0.075 mm; E, Ū = 0.105 mm
respectively. The cracks are first generated around the upper and lower material interfaces for the initial guess design,
while it nucleates in the middle of the inner inclusion phase and 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 Figure 10. In this example,
the structure with optimal design is 15
% more resistant to fracture than the initial composite structure.
4.3 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 Figure 11A. The dimensions of the plate are 100
× 100 mm2, 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 Figure 11A. The
final design of the inclusion phase is shown in Figure 11B. 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 Figure 12. Here, the interface cracks first initiate at the
interface (see Figure 12B). Subsequently, the left and right interface cracks propagate vertically and merge. The complete
fracture patterns and design objective values are shown in Figure 13B for the initial design, and in Figure 13C, for the

initial design
J = 8.34 mJ
final design
J = 9.52 mJ
(A) (C) (B)
FIGURE 10 Fracture resistance comparison of two composite structures without initial crack subjected to incremental traction loads
(A) (B)
Matrix
Inclusion
100
100
33.33
15
FIGURE 11 A square plate without initial crack subjected to uniaxial tension. A, Geometry of the initial design; B, Final design [Colour
figure can be viewed at wileyonlinelibrary.com]
(A) (B) (C) (D)
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 mm; D, Ū = 0.08 mm
optimal design. Both responses are compared in Figure 13A. Here, the fracture resistance of the final design structure has
been increased by 33
% as compared with the initial design.
4.4 Design of a plate with one initial crack in three-point bending
The purpose of this example is to design a plate subjected to three-point bending with one initial crack. The problem
geometry of the square plate is depicted in Figure 14A. The dimensions of the plate are 50
× 100 mm2. The domain is
uniformly discretized into 60
× 120 square-shaped bilinear elements. The load consists into a prescribed displacement at
the center of the beam on the top edge. The left bottom corner node is fixed, while the node at the right bottom corner the

initial design
J = 19.10 mJ
final design
J = 25.52 mJ
(A) (B) (C)
FIGURE 13 Fracture resistance comparison of two composite structures without initial crack subjected to uniaxial tension
(A) (B)
50
100
12.5
Matrix
Inclusion
Crack
(C)
25
20
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 [Colour figure can be viewed at wileyonlinelibrary.com]
y-displacement is fixed, and the x-displacement is free. For this case, an initial pre-existing crack is shown in Figure 14B,
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 nondesigned region. Figure 14C 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 Figure 15. The initial matrix crack
first 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 Figure 16A. The final crack patterns and the total required
fracture energies for complete failure for the initial and final designs are shown in Figure 16B and Figure 16C, respectively. In this example, the fracture resistance of the structure with final design has been increased by 76
% as compared
with the initial design.

(A) (B)
(C) (D)
FIGURE 15 Crack propagation of the initial design composite structure with one initial crack subjected to uniaxial tension. A, Ū = 0 mm;
B,
Ū = 0.11 mm; C, Ū = 0.14 mm; D, Ū = 0.16
initial design
J = 5.24 mJ
final design
J = 9.20 mJ
(A) (C)
(B)
FIGURE 16 Fracture resistance comparison of two composite structures with one initial crack subjected to three-point bending
4.5 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 Figure 17A.
On the lower and upper ends, the
y-displacements are fixed, while the x-displacements are free. On the left and right
ends, the
y-displacements 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 nondesignable to avoid
nonphysical designs. The final design is presented in Figure 17C. 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 Figure 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 and 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 Figure 19. In this
example, the resistance to fracture has been increased by 53
% in comparison with the initial design.

(A) (B) (C)
100

Matrix
Inclusion
Crack
10
10
25
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 [Colour figure can be viewed at wileyonlinelibrary.com]
(A) (B) (C) (D)
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
initial design
J = 9.77 mJ
final design
J = 14.94 mJ
(A) (C) (B)
FIGURE 19 Fracture resistance comparison of two composite structures containing multiple inclusions and one initial crack subjected to
tensile loads
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. (i) 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; (ii) initiation of cracks can be included in the
analysis, which would be not possible with techniques like extended FEM
11; (iii) including interfacial damage is simple
in the present framework, using, eg, extensions of the phase field as proposed in the work of Nguyen et al.
24 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.

5 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 crack 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 BESO 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 with 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.
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, and 5171101743).
ORCID
Julien Yvonnet http://orcid.org/0000-0002-9365-6601
Liang Xia http://orcid.org/0000-0001-6928-0854
REFERENCES
1. Quan Z, Larimore Z, Wu A, et al. Microstructural design and additive manufacturing and characterization of 3D orthogonal short carbon
fiber/acrylonitrile-butadiene-styrene preform and composite.
Compos Sci Technol. 2016;126:139-148.
2. Cadman JE, Zhou S, Chen Y, Li Q. On design of multi-functional microstructural materials.
J Mater Sci. 2013;48(1):51-66.
3. Fritzen F, Xia L, Leuschner M, Breitkopf P. Topology optimization of multiscale elastoviscoplastic structures.
Int J Numer Methods Eng.
2016;106(6):430-453.
4. Gu GX, Dimas L, Qin Z, Buehler MJ. Optimization of composite fracture properties: method, validation, and applications.
J Appl Mech.
2016;83(7):1-7.
5. San B, Waisman H. Optimization of carbon black polymer composite microstructure for rupture resistance.
J Appl Mech. 2016;84(2):1-13.
6. Xia L, Da DC, Yvonnet J. Topology optimization for maximizing the fracture resistance of quasi-brittle composites.
Comput Methods Appl
Mech Eng
. 2018;332:234-254.
7. Lamon J, Carrere N, Martin E. The influence of the interphase and associated interfaces on the deflection of matrix cracks in ceramic
matrix composites.
Composites A. 2000;31(11):1179-1190.
8. Tvergaard V. Model studies of fibre breakage and debonding in a metal reinforced by short fibres.
J Mech Phys Solids. 1993;41(8):1309-1326.
9. Nguyen TT, Yvonnet J, Bornert M, Chateau C. 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.
J Mech Phys Solids.
2016;99:320-350.
10. Narducci F, Pinho ST. Exploiting nacre-inspired crack deflection mechanisms in CFRP via micro-structural design.
Compos Sci Technol.
2017;153:178-189. https://doi.org/10.1016/j.compscitech.2017.08.023
11. Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing.
Int J Numer Methods Eng.
1999;46(1):131-150.
12. Sukumar N, Moes N, Moran B. Extended finite element method for three-dimensional crack modelling.
Int J Numer Methods Eng.
2000;48(11):1549-1570.

13. Bernard PE, Moes N, Chevaugeon N. Damage growth modeling using the thick level set (TLS) approach: efficient discretization for
quasi-static loadings.
Comput Methods Appl Mech Eng. 2012;233-236:11-27.
14. Cazes F, Moes N. Damage growth modeling using the thick level set (TLS) approach: efficient discretization for quasi-static loadings.
Int
J Numer Methods Eng
. 2015;103(2):114-143.
15. Francfort GA, Marigo JJ. Revisiting brittle fracture as an energy minimization problem.
J Mech Phys Solids. 1998;46(8):1319-1342.
16. Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: robust algorithmic implementation
based on operator splits.
Comput Methods Appl Mech Eng. 2010;199(45-48):2765-2778.
17. Kuhn C, Müller R. A continuum phase field model for fracture.
Eng Fract Mech. 2010;77(18):3625-3634.
18. Bourdin B, Francfort GA, Marigo JJ. The variational approach to fracture.
J Elast. 2008;91(1-3):5-148.
19. Pham K, Marigo JJ. The variational approach to damage. I: the foundations.
Comptes Rendus Mec. 2010;338:191-198.
20. Bourdin B, Francfort GA, Marigo JJ. Numerical experiments in revisited brittle fracture.
J Mech Phys Solids. 2000;48(4):797-826.
21. Mumford D, Shah J. Optimal approximations by piecewise smooth functions and associated variational problems.
Commun Pure Appl
Math
. 1989;42(5):577-685.
22. Ambrosio L, Tortorelli VM. Approximation of functional depending on jumps by elliptic functional via
𝛾-convergence. Commun Pure Appl
Math
. 1990;43(8):999-1036.
23. Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: variational principles and multi-field
fe implementations.
Int J Numer Methods Eng. 2010;83(10):1273-1311.
24. Nguyen TT, Yvonnet J, Zhu QZ, Bornert M, Chateau C. A phase-field method for computational modeling of interfacial damage interacting
with crack propagation in realistic microstructures obtained by microtomography.
Comput Methods Appl Mech Eng. 2016;312:567-595.
25. Bendsøe MP, Kikuchi N. Generating optimal topologies in structural design using a homogenization method.
Comput Methods Appl Mech
Eng
. 1988;71(2):197-224.
26. Bendsøe M, Sigmund O.
Topology Optimization: Theory, Methods and Applications. Berlin, Germany: Springer-Verlag Berlin Heidelberg;
2003.
27. Nguyen-Xuan H. A polytree-based adaptive polygonal finite element method for topology optimization.
Int J Numer Meth Eng.
2017;110(10):972-1000.
28. Chau K, Chau K, Ngo T, Hackl K, Nguyen-Xuan H. A polytree-based adaptive polygonal finite element method for multi-material topology
optimization.
Comput Methods Appl Mech Eng. 2018;332:712-739.
29. Xie YM, Steven GP. A simple evolutionary procedure for structural optimization.
Comput Struct. 1993;49(5):885-896.
30. Huang X, Xie YM. Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method.
Finite
Elem Anal Des
. 2007;43(14):1039-1049.
31. Sethian JA, Wiegmann A. Structural boundary design via level set and immersed interface methods.
J Comput Phys. 2000;163(2):489-528.
32. Wang MY, Wang X, Guo D. A level set method for structural topology optimization.
Comput Methods Appl Mech Eng.
2003;192(1-2):227-246.
33. Allaire G, Jouve F, Toader AM. Structural optimization using sensitivity analysis and a level-set method.
J Comput Phys.
2004;194(1):363-393.
34. Xia L, Fritzen F, Breitkopf P. Evolutionary topology optimization of elastoplastic structures.
Struct Multidiscip Optim. 2017;55(2):569-581.
35. Xia L, Xia Q, Huang X, Xie Y. Bi-directional evolutionary structural optimization on advanced structures and materials: A comprehensive
review.
Arch Comput Methods Eng. 2018;25(2):437-478. https://doi.org/10.1007/s11831-016-9203-2
36. Xia L, Breitkopf P. Concurrent topology optimization design of material and structure within FE
2 nonlinear multiscale analysis framework.
Comput Methods Appl Mech Eng. 2014;278:524-542.
37. Xia L, Breitkopf P. Multiscale structural topology optimization with an approximate constitutive model for local material microstructure.
Comput Methods Appl Mech Eng. 2015;286:147-167.
38. Vermaak N, Michailidis G, Parry G, Estevez R, Allaire G, Brechet Y. Material interface effects on the topology optimization of multi-phase
structures using a level set method.
Struct Multidiscip Optim. 2014;50(4):623-644. https://doi.org/10.1007/s00158-014-1074-2
39. Lawry M, Maute K. Level set topology optimization of problems with sliding contact interfaces.
Struct Multidiscip Optim.
2015;52(6):1107-1119.
40. Liu P, Luo Y, Kang Z. Multi-material topology optimization considering interface behaviour via XFEM and level set method.
Comput
Methods Appl Mech Eng
. 2016;308:113-133.
41. Behrou R, Lawry M, Maute K. Level set topology optimization of structural problems with interface cohesion.
Int J Numer Methods Eng.
2017;112(8):990-1016. https://doi.org/10.1002/nme.5540
42. Duysinx P, Bendsøe MP. Topology optimization of continuum structures with local stress constraints.
Int J Numer Methods Eng.
1998;43(8):1453-1478.
43. Guo X, Zhang WS, Wang MY, Wei P. Stress-related topology optimization via level set approach.
Comput Methods Appl Mech Eng.
2011;200(47-48):3439-3452.
44. Xia Q, Shi T, Liu S, Wang MY. A level set solution to the stress-based structural shape and topology optimization.
Comput Struct.
2012;90-91:55-64.
45. Bruggi M, Duysinx P. Topology optimization for minimum weight with compliance and stress constraints.
Struct Multidiscip Optim.
2012;46(3):369-384.

46. Luo Y, Wang MY, Kang Z. An enhanced aggregation method for topology optimization with local stress constraints. Comput Methods Appl
Mech Eng
. 2013;254:31-41.
47. Cai S, Zhang W, Zhu J, Gao T. Stress constrained shape and topology optimization with fixed mesh: a b-spline finite cell method combined
with level set function.
Comput Methods Appl Mech Eng. 2014;278:361-387.
48. Cai S, Zhang W. Stress constrained topology optimization with free-form design domains.
Comput Methods Appl Mech Eng.
2015;289:267-290.
49. Sun Z, Li D, Zhang WS, Shi S, Guo X. Topological optimization of biomimetic sandwich structures with hybrid core and CFRP face sheets.
Compos Sci Technol. 2017;142:79-90.
50. Challis VJ, Roberts AP, Wilkins A. Fracture resistance via topology optimisation.
Struct Multidiscip Optim. 2008;36(3):263-271.
51. Amir O, Sigmund O. Reinforcement layout design for concrete structures based on continuum damage and truss topology optimization.
Struct Multidiscip Optim. 2012;47(2):157-174.
52. Amir O. A topology optimization procedure for reinforced concrete structures.
Comput Struct. 2013;114-115:46-58.
53. Jansen M, Lombaert G, Schevenels M, Sigmund O. Topology optimization of fail-safe structures using a simplified local damage model.
Struct Multidiscip Optim. 2013;49(4):657-666.
54. James KA, Waisman H. Failure mitigation in optimal topology design using a coupled nonlinear continuum damage model.
Comput
Methods Appl Mech Eng
. 2014;268:614-631.
55. Kang Z, Liu P, Li M. Topology optimization considering fracture mechanics behaviors at specified locations.
Struct Multidiscip Optim.
2016;55(5):1847-1864.
56. Verhoosel CV, de Borst R. A phase-field model for cohesive fracture.
Int J Numer Methods Eng. 2013;96(1):43-62.
57. Huang X, Xie YM. Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials.
Comput
Mech
. 2009;43(3):393-401.
58. Da DC, Cui XY, Long K, Li GY. Concurrent topological design of composite structures and the underlying multi-phase materials.
Comput
Struct
. 2017;179(1):1-14.
59. Maute K, Schwarz S, Ramm E. Adaptive topology optimization of elastoplastic structures.
Structural Optimization. 1998;15(2):81-91.
60. Schwarz S, Maute K, Ramm E. Topology and shape optimization for elastoplastic structural response.
Comput Methods Appl Mech Eng.
2001;190(15-17):2135-2155.
61. Huang X, Xie YM. Topology optimization of nonlinear structures under displacement loading.
Eng Struct. 2008;30(7):2057-2068.
62. Xia L, Breitkopf P. Recent advances on topology optimization of multiscale nonlinear structures.
Arch Comput Methods Eng.
2017;24(2):227-249.
63. Sigmund O. A 99 line topology optimization code written in MATLAB.
Struct Multidiscip Optim. 2001;21(2):120-127.
64. Huang X, Xie YM.
Topology Optimization of Continuum Structures: Methods and Applications. Chichester, UK: John Wiley & Sons; 2010.
65. Nguyen TT, Yvonnet J, Zhu QZ, Bornert M, Chateau C. A phase field method to simulate crack nucleation and propagation in strongly
heterogeneous materials from direct imaging of their microstructure.
Eng Fract Mech. 2015;139:18-39.
How to cite this article: Da D, Yvonnet J, Xia L, Li G. Topology optimization of particle-matrix composites for
optimal fracture resistance taking into account interfacial damage.
Int J Numer Methods Eng. 2018;115:604–626.
https://doi.org/10.1002/nme.5818
APPENDIX
A.1 Linearization of the displacement problem
Even though the phase field problem is linear in the staggered scheme, ie, 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 (13) and the interface cohesive model in (31). 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
= 𝛺𝝈e 𝜺e(𝛿u)d𝛺 + 𝛺𝛾𝛽(x)t(w, 𝛼) · 𝛿wd𝛺 𝛺f · 𝛿ud𝛺 𝜕𝛺
F
F̄ · ud𝛤 = 0, (A1)
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(uk, d) = –(uk, d) = 0, (A2)

where uk is the displacement solution from the kth iteration. The displacements at the current iteration are given by
uk+1 = uk + Δu. (A3)
From (A2), we obtain
DΔu(uk) = 𝛺 𝜕𝜕𝝈𝜺ee 𝜺e𝜺) ∶ 𝜺e(𝛿𝜺) + 𝛺 𝜕t𝜕(ww) ∶ Δw 𝛿wd𝛺 (A4)
with
Δw(x) = h∇Δu(x) ∇𝜙(x)
||𝜙(x)||
(A5)
and
𝜕[𝝈e]
𝜕[𝜺e]
=
C(u, d) = [(1 d)2 + k] {𝜆R+[1]T[1] + 2𝜇P+} + {𝜆R[1]T[1] + 2𝜇P} , (A6)
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.
A.2 Finite element method 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 the work of Nguyen et al
65). Using FEM, the phase field and phase
field gradient in one element are approximated by
d(x) = Nd(x)de, d(x) = Bd(x)de, (A7)
where
de are nodal phase field values in one element and Nd(x) and Bd(x) are matrices of shape functions and of shape
functions derivatives associated to phase field variable, respectively. Introducing the aforementioned FEM discretization
into the weak form (24), the following linear discrete system of equations can be obtained:
Kdd̃ = Fd, (A8)
where
Kd = 𝛺 {(g𝓁c (1 𝛽) + 2) NT d Nd + (1 𝛽)gc𝓁BT d Bd} d𝛺 (A9)
and
Fd = 𝛺2NT d(un)d𝛺. (A10)
A.3 Finite element method discretization of the displacement problem
Similarly, the displacement field and incremental displacement field can be expressed using the FEM approximations

u = Nue, Δu = NΔue, (A11)
where N denotes the matrix of shape functions associated to displacement variables and ue and Δue are nodal displace
ment components and nodal incremental displacement components in one element. Furthermore, we have
[𝜺](Δu) = BuΔue
(A12)

where Bu is a matrix of shape function derivatives. From (7), the diffuse jump approximation vector and its incremental
counterparts can be discretized as
w = hNB̃ uue, Δw = hNB̃ uΔue, (A13)
where
N = [ n0 0 1 n2 n0 0 1 n2 ] , (A14)

and n1 and n2 are the x– component and y– component of the normal vector. The smoothed jump strain at the interfaces
is defined by

[𝜺̄] =
⎡⎢
⎢⎣
𝜺̄11 = 𝛾𝛽(x)
⎡⎢⎢⎣
𝜔1n1
𝜔2n2
. (A15)
𝜺̄22
2
𝜺̄12 1
2
(𝜔1n2 + 𝜔2n1)

⎤⎥⎥⎦
⎤⎥⎥⎦
Then, it yields
[𝜺̄u)] = h𝛾𝛽(x)MB̃ uΔue (A16)
with

M =
⎡⎢⎢
n2
1
n1n2 0
0
0 0 n1n2 n2 2
.
(A17)


1
2
n1n2
1

2
n2
2
1

2
n2
1
1

2
n1n2
⎤⎥⎥⎦
With the aforementioned FEM discretization, the linear tangent problem reduces to the following linear system of

algebraic equations: KtanΔũ = –R (
ũ k) , (A18)

where
Ktan = 𝛺 [BT u h𝛾𝛽(x)B̃ T uMT] C(x) [Bu h𝛾𝛽(x)B̃ uM] d𝛺 + 𝛺h2𝛾𝛽(x)B̃ T uNTKINB̃ ud𝛺, (A19)
and
R = 𝛺 [BT u h𝛾𝛽(x)B̃ T uMT] C(x) [Bu h𝛾𝛽(x)B̃ uM] (ue)kd𝛺
+𝛺h𝛾𝛽(x)B̃ T uNTt (wk) d𝛺 + 𝛺fNTd𝛺 + 𝛺FN ̄ Td𝛤 . (A20)
Overall numerical algorithm for crack propagation The overall algorithm is described as follows.
1. Set the initial displacement field
u0(x), the phase field d0(x), and the strain-history function 0.
2. Compute the phase field
𝛽(x).
3.
For all loading increments: (at each time tn+1), given dn, un, and n(x).
(a) Compute the history function
(tn+1) according to (21).
(b) Compute the crack phase field
dn+1(x) by solving linear problem (A8).
(c) Compute
un+1(x) :
(i) Initialize
uk = unWhile||Δuk+1|| > 𝜖, 𝜖 ≪ 1:
(ii) Compute
Δue k+1 by (A18).
(iii) Update
uk+1 = uk + Δue k+1.
(iv)
(.)n+1 (.)n and go to (a).
End
End