Topology Optimization for Maximizing the Fracture Resistance of Periodic Quasi-Brittle Composites Structures

Topology optimization for maximizing the fracture resistance of particle-matrix composites is investigated. The methodology developed in our previous works, combining evolutionary topology optimization and phase field method to fracture embedding interfacial damage, is applied and extended to periodic composites and multiple objectives. On one hand, we constrain the periodicity of unit cells geometry and conduct their topology optimization for one given load prescribed over the whole structure. On the other hand, we consider a single unit cell whose topology is optimized with respect to the fracture energy criterion when subjected to multiple loads. Size effects are investigated. We show that significant enhancement of the fracture resistance can be achieved for the studied composite structures by the present method. In addition, a first attempt to fracture resistance enhancement of a unit cell associated with a material is investigated for multiple loads, exhibiting a complex optimized microstructure.


Introduction
Over the past few decades, topology optimization (TO) has undergone a tremendous development not only regarding academic investigations but also in various industrial applications since the seminal paper in [1]. The key merit of topology optimization over conventional size and shape optimization is that the former can provide more design freedom, consequently leading to the creation of novel and highly efficient designs. With the increasing development of contemporary manufacturing processes such as the popular additive manufacturing/3D printing techniques [2], it is nowadays possible to directly fabricate topologically optimized designs from digital files with arbitrary geometrical shapes.
Topology optimization in the context of minimizing structural compliance and maximizing the elastic mechanical properties of composite materials has been extensively studied [3][4][5][6][7][8][9][10][11]. In most cited works, only linear properties have been investigated. One new and challenging topic which has emerged recently is TO for maximizing fracture resistance of structures and materials. Fracture resistance to mechanical load is amongst the most important criteria to design a structure in engineering and the use of TO to help designing structures and materials with respect to this criterion is of critical importance. By resistance to fracture, we mean here maximizing the total mechanical work during the whole fracturing process, from initiation up to propagation of internal micro cracks and finally failure of the structure (see Figure 1). This is in contrast with, e.g., maximizing limit stress state in an elastic context, as it requires simulating the full nonlinear fracturing process within the optimization scheme. Early works attempting to maximize fracture resistance using TO can be traced back to [12] in the context of Level-set TO and linear fracture mechanics, and have been recently developed in a similar framework in [13]. Improvements of this idea by taking into account an incremental damage response during the load have been proposed, incorporating damage mechanics [14,15], nonlocal damage theory [16], or interfacial damage [17].
Another important progress was to include brittle fracture propagation within TO. First attempts have made use of fracture mechanics like XFEM [18], but such approaches require the existence of an initial crack and are limited to moderately simple crack morphologies. The phase field method [19][20][21][22][23][24][25][26][27][28], has many advantages to describe brittle fractures while maintaining the simplicity and versatility of simple finite elements, and allows for initiating cracks, handling complex and arbitrary, 3D crack networks and can be extended to more complex behaviors like plasticity [29], interfacial damage [28,30], or coupled multiphysics [31], among many others. Another very important point is that, within phase field, crack paths are not mesh-dependent, and a regular mesh can be employed, which might be advantageous in TO procedures.
Phase field and genetic algorithms have been used in [32] to find the optimal location of particles to maximize fracture mechanics indicators such as the peak force, maximum deformation at failure point, and maximum fracture energy during an incremental procedure. The first works to our best knowledge combining phase field and TO were introduced in [33][34][35], where the Bi-directional Evolutionary Structural Optimization (BESO) TO [36] was employed to optimize the fracture resistance of two-phase structures with respect to inclusion shapes, including cracks in both bulk and interfaces. In [37,38], Solid Isotropic Material with Penalization (SIMP) TO and phase field were combined to optimize the fracture energy in one-phase material structures, and a Level-Set TO-phase field approach was proposed in [39] to optimize the fracture resistance of composites. A SIMP-phase field procedure for two-phase materials with comparison with BESO has been developed in [40]. An experimental study involving the fracture of 3D-printed topology-optimized beams can be found in [41]. A comparative review of the different TO approaches can be found in [10].
This paper presents developments of the framework proposed in [34], combining BESO topology optimization and the phase field method to fracture embedding interfacial damage. The contributions include extensions to periodic composites, where the topology optimization is constrained to be periodic over the different unit cells composing the structure. Another extension to multiple loads is proposed, applied to the fracture resistance of one unit cell subjected to possible multiple loads. Finally, we investigate strength size effects in the optimized periodic composite structures.
The paper is organized as follows. We first review the phase field method to fracture in Section 2. The topology optimization framework and the introduced extensions to maximize the fracture resistance, combining the incremental simulations and BESO method, is presented in Section 3. Numerical examples are provided in Section 4, to show the potential of the method and to investigate strength size effects in such optimized structures.

Phase Field Modeling of Bulk Crack and Interfaces
In this section, we briefly review the phase field model for fractures interacting with interfacial damage as presented in [28]. This model will basically serve as an engine for the topology optimization procedure. Linearization and implementation aspects are briefly reviewed in Appendixes A and B. More details can be found in our previous works [28,34].
Let Ω ∈ R d be an open domain describing a solid with external boundary ∂Ω (see Figure 2a). The solid contains internal material interfaces between different phases, collectively denoted by Γ I , as illustrated in Figure 2a. The heterogeneous material is composed of two phases, the matrix and inclusions, associated respectively to the domains Ω m and Ω i , and such that Ω = Ω m ∪ Ω i . During the loading, cracks may propagate within the solid and can pass through the material interfaces. The crack surfaces are collectively denoted by Γ. In this work, we adopt smeared representations of both cracks and material interfaces: the interfaces between different material phases are described by a fixed scalar phase field β(x) (see Figure 2b) while cracks are approximated by an evolving phase field d(x, t) (see Figure 2c). 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. The macroscopic energy of a structure containing sharp cracks and interfaces is defined by (see e.g., [28]): where Ψ e is the elastic bulk strain energy density, ε is the second-order strain tensor, u is the displacement field, Ψ I is an interface strain density function, g c is the toughness, and [[u]] is the displacement jump across the interface Γ I . In [28], we have proposed a smoothed representation of both discontinuities associated with bulk and interface cracks. In this framework, the above energy functional is substituted by the following expression: where γ d is a crack density function associated with the smoothed representation of bulk cracks defined by  In (2), Ψ e refers to the strain energy density function depending on the part of the strain ε e defined by: where (∇ s u) ij = (u i,j + u j,i )/2 and (n⊗ S w) ij = (n i w j + w i n j )/2. The strain density function is defined such that the damage is associated with tensile strain only and is given in the form: where Ψ e+ and Ψ e− denote the parts of the elastic strain density function associated with tensile and compressive parts, respectively, and are expressed for an isotropic material as: Above, ε e± are the tensile and compressive parts of the strain tensor, respectively, defined in the following (see Equation (13)) and λ and µ are the Lamé's elastic parameters. Note that we have omitted the dependance of the local parameters to the local phase, i.e., the local elasticity parameters are defined as λ i , µ i and λ m , µ m , where the indices i and m refer to the inclusion and matrix phase materials, respectively. Similarly, local toughness are denoted by g i c and g m c , respectively. The field β(x) is a smooth indicator for interfaces which satisfies: where β is the regularization parameter describing the width of the regularized interfaces, with Problem (7) can be solved by classical finite elements (see [28]). Above, w(x) is a smooth representation of the displacement jump function across the interface, defined by: Using variational principle for minimizing E with respect to displacements (assuming the crack phase field d and interface β are fixed), i.e., F · udΓ with f and F being body forces and prescribed traction over the boundary ∂Ω F , we obtain the weak form for u(x) ∈ S u : (Ω) and where the Cauchy stress σ e is expressed as with k << 1 a small scalar parameter, ε e (δu) = ∇ s δu − n⊗ S wγ β and δw = h∇(δu)∇φ/ ∇φ . Above, ε e± are defined according to: where ε i and n i are the eigenvalues and eigenvectors of ε e , i.e., satisfying ε e n i = ε i n i . In (12) and (13), The cohesive traction law at the interfaces is given in 2D by: where t n and t t denote normal and tangential parts of the traction vector t across the interface Γ I oriented by its normal n I . We chose the simple model [42]: where δ n is related to the interface fracture toughness g I c and the interface fracture strength t u by: w n = w · n I , t t = 0 for the sake of simplicity, and g I c is the toughness associated with the interface. Note that the problem (11) is nonlinear. A Newton iterative scheme is then required to solve it, as described in Appendix A.
The phase field problem is obtained by minimization of energy under the constraint thatḋ ≥ 0. Here, we follow the approach of Miehe et al. [23] and extended in [28] where the weak form associated with this process is given by: and (Ω) and where the history functional is introduced to handle the Irreversibility condition: In the present work, we used a staggered procedure to solve alternatively problems (11) and (17). Numerical details and finite element discretization of the mechanical problem and phase field problem are provided in Appendix A and B, respectively. The overall algorithm is summarized 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. For all loading increments: (at each time t n+1 ), given d n , u n , and H n (x).

Topology Optimization Problem to Fracture Resistance Maximization
In this section, the topology optimization problem to fracture resistance maximization, first defined in [33,34], is extended here to periodic composite structures and multiple objectives.

Periodic Structures
The topology optimization problem is conducted with respect to a density variable ρ(x) which is associated with the inclusion phase. In other words, ρ(x) = 1 ∀x ∈ Ω i and ρ(x) = 0 ∀x ∈ Ω m (see Figure 2a). We consider periodic structures composed of a repetition of N x × N y unit cells along x− and y− directions. The dimensions of the unit cell are x 0 × y 0 . The optimization problem is then defined as follows: with where R 1 is given by (11) and R 2 is given by (17), t max denotes the maximum simulation pseudo-time associated with the final prescribed displacement at the failure step, f ext is the external force response at the load point, and f is the prescribed volume fraction of the inclusion phase andŪ is the prescribed displacement load vector. Above, S u and S d are the appropriate vector spaces associated with the unknown fields u and d defined in Section 2. The objective function is the total mechanical work during the fracturing process. Above, H = [i 0; 0 j], i = 1, 2, ..., N x , j = 1, 2, ..., Ny is a matrix which defines the periodicity with respect to the unit cell and x 0 is a point coordinate in the periodic unit cell Ω (see Figure 1). The pseudo-density ρ can be interpreted as an indicator such that the value of one corresponds to the inclusion phase, whereas zero corresponds to the matrix phase. Note that, in the present framework, quasi-static simulations are performed, and the time is here only related to the load evolution. In addition, the elasticity parameters and toughness of the phases are defined with respect to the density as: Note that evaluation of (24) requires solving the full fracture problem, from initiation to complete failure of the structure.

Multiple Objectives
In the case when the topology is optimized with respect to k = 1, 2, ..., N loads, the optimization problem is modified according to: and with t max k the maximum simulation time for the k − th load.

Discrete Topology Optimization Problem
In what follows, we only define the solution procedure for the periodic composite problem, as the idea is very similar for the multiple loading case. Problems (19)-(24) are solved using Finite Elements based on the same mesh than the one used to solve the displacement and fracture problems defined in Section 2. A regular mesh of bi-linear 4-node elements is employed.
We define the vector of discrete density values {ρ}= {ρ 1 , ρ 2 , ..., ρ N e } with N e the number of Finite Elements. The discrete form of (22) and (23) is then defined as: where n load is the number of load steps and J h in (34) is approximated by: where ∆u (n) denotes the prescribed load increment at load n and R 1 is approximated by

Sensitivity Analysis and Extended BESO Procedure
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 adjoint vectors µ (n) , λ (n) of the same dimension as the vector of unknowns u are introduced to enforce zero residual R 1 at load steps n − 1 and n for each term of the total mechanical work (39). To be solved by the adjoint method, the modified expression for the objective function is introduced as follows: Due to the asserted static equilibrium, the residuals R 1 (n) and R 1 (n−1) must vanish. The objective value is thus invariant with respect to the values of the adjoint vectors λ (n) and µ (n) (n = 1, . . . , n load ), i.e., This equivalence also holds for the sensitivity with respect to changes of density ρ e on element e In the following, the derivative ∂ J/∂ρ e 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 a and a matrix M, we have It can be shown (see e.g., [35]) that where K (m) is the tangent stiffness matrix of the nonlinear mechanical system at the balance equation of the m−th load increment.
To avoid the evaluation of the unknown derivatives of u F are sought as following by solving the adjoint systems with the prescribed values at the essential nodes: and λ (n) and µ (n) The two relations (48) and (49) together with (47) fully determine the values of the adjoint vectors λ (n) and µ (n) . Finally, the objective function ∂ J/∂ρ e can be computed via The computation of the sensitivity consists in solving two linear systems of Equations (48) and (49). The extended BESO method recently developed in [43,44] 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 efficiency of the method, especially in dealing with nonlinear designs in the presence of dissipative effects. Note that, in the present work, due to the staggered approach used to solve the problems (11) and (17), we only take into account the constraint (20) through the sensitivity analysis, assuming that (17) is verified at the previous iteration. Details about BESO implementation can be found in [34].

Numerical Examples
In all following numerical examples, bi-material composites are investigated. The material parameters of each phase are taken as E i = 52 GPa, E m = 10.4 GPa, ν i = ν m = 0.3, where E denotes the Young's modulus and ν the Poisson's ratio. The indices i and m refer to the inclusion and matrix materials, respectively. The bulk toughness g c is the same for the matrix and the inclusions and equal to g m c = g i c = g c = 0.1 N/mm. The interface fracture strength is chosen as t u = 10 MPa. In addition, we use g c = g I c = 1 × 10 −4 kN/mm. In all tests, regular meshes using quadrilateral bilinear elements are adopted. The same finite element discretization is employed for both displacement and crack phase fields. The regularization parameter d = β describing the width of bulk cracks and interface cracks is chosen as twice the finite element size d = β = 2 e , with e the characteristic finite element size.
During the topology optimization process, the volume fraction of inclusions is maintained constant. We define the gain in fracture resistance as G = 100 × (J opt − J init )/J init , where J opt is the fracture resistance of the optimized structure while J init is the initial structure design.

Periodic Composite under Three-Point Bending
We first consider the composite structure depicted in Figure 3a. The composite is made of N s = 7 × 3 periodic unit cells. The initial guess design comprises the unit cells consisting of circular inclusions in a square domain, such that the volume fraction is equal to f = 30%. The dimensions of the structure are 350 × 150 mm 2 . The domain is uniformly discretized into 350 × 150 4-node bilinear elements. The boundary conditions are defined as follows: the left bottom corner node is fixed, while in the right bottom corner node, the y-displacement is fixed and the x-displacement is free. On the upper end, the load consists of a prescribed displacementŪ at the center point. The displacement is prescribed along the y-direction while the displacement along x is free. The computation is performed with monotonic displacement increments ∆Ū = 0.01 mm until the corresponding resultant vertical force F y at the displacement application point is below a prescribed criterion value, indicating that the structure is completely broken. As mentioned in Section 3, all unit cells are constrained to keep the same geometry during the topology optimization process. Figure 3b shows the obtained final design for this composite structure and this specific load. We observe that, while the initial design was composed of disconnected phases, the optimal design is formed with connected phases along the x-direction and forms a kind of wavy layered structure.
A comparison between the crack paths obtained within the structures with initial and optimal designs are depicted in Figures 4 and 5, respectively. In both cases, micro cracks nucleate from interfaces or from the surface and merge to form a macro crack, which is more complex in the case of the optimal design, with more interfacial cracks, and forming many branches, which may explain that the structure requires more energy for complete failure.
The displacement-force response curve of both structures with initial and optimal designs are depicted in Figure 6a. The gain in fracture resistance for the optimized structure is equal to G = 124%.

Size Effects
In this section, size effects are investigated. We analyze the influence of the size on the strength of the structure when one fixed optimal design is used. More specifically, we use the optimal design obtained in the example of Section 4.1 and define four structures composed of 7 × 3, 14 × 6, 21 × 9 and 28 × 12 cells, corresponding to length L equal to 350 mm, 700 mm, 1050 mm, and 1400 mm, respectively. The associated heights H are 150 mm, 300 mm, 450 mm, and 600 mm. In all cases, the vertical displacement is prescribed as an increasing uniform value of ∆Ū = 0.01 mm during the simulation, until the reaction force goes below a prescribed value indicating that the periodic composite is fully broken. One convenient way to describe the size effects in quasi-brittle structures is to plot the nominal stress σ N with respect to the height H in log scale. The nominal stress is defined as [45] where F is the critical force to failure (estimated numerically in the load-displacement curve as the maximum load), L is the length of the beam, and b is an arbitrary thickness, taken here as b = 1 m. Results are presented in Figure 7. The size effects are clearly shown in Figure 7a, where a decrease of nominal stress with the size is noticed. This is consistent with classical experimental results for fracture in quasi brittle beams [45]. In Figure 7b, we can see the effects of keeping the same topology and increasing the size of the structure. What is very interesting is that, when using the optimized topology obtained on the smallest structure, then, at relatively low computational costs, its use on large structure still bring a large gain (40%) in fracture resistance. (b) (a) Figure 7. Size effects associated with the strength of the structure (a) nominal stress to failure with respect to the height H for initial and optimized topologies; (b) gain between initial and optimized designs with respect to the size. For each size, the same topology is kept.
As an illustration, the crack propagation in both initial and optimized structures for the structure involving 6 × 14 cells is depicted in Figures 8 and 9, respectively. We can observe the same type of patterns in both structures rather than in the smaller structure. A comparison between the reference and the optimized solution in that case is provided in Figure 6b

Design of a Periodic Composite under Non-Symmetric Three-Point Bending
We consider the same structure as in the former example, but using a non-symmetric three-point bending load. The dimensions, geometry, and boundary conditions are depicted in Figure 10a, showing the initial design, here again composed of circular inclusions as defined in the previous example. The considered composite is composed of 9 × 3 unit cells which are repeated periodically along each space directions. The dimension of the composite are 450 × 150 mm 2 , and the domain is uniformly discretized into 450 × 150 square shape bilinear elements. The optimized design is depicted in Figure 10b. We can see that a change in the position of the load drastically modifies the optimized design. The crack propagation steps are illustrated for the initial and optimized structures in Figures 11 and 12, respectively. In both cases, cracks initiate from the interface, and then merge in complex crack patterns at the failure step. The crack network is much more complex in the case of the optimized structure and induces a more ductile behavior, as shown in Figure 13. In this example, the gain in fracture resistance for the optimized design is G = 100.6%.

Multi-Objective Topology Optimization of a Material Unit Cell with Randomly Distributed Inclusions
In this part, we extend the topology optimization to multi objectives in order to find the optimal topology which maximizes the fracture resistance of a material unit cell possibly subjected to several representative loads. An initial design for the unit cell is described in Figure 14a, consisting of randomly distributed squares, such as their volume fraction is f = 12%. Three load cases are considered, whose corresponding boundary conditions are described in Figure 14a-c. Both first loads correspond to traction in the xand y-directions, while the third load corresponds to shear. Here, the optimization problem is then conducted with respect to the response obtained by these three loads to maximize the objective functions J = ∑ i J i , where J i , i = 1:N are the objective function defined by (41). The crack paths for the three corresponding loads in the RVE with initial design are shown in Figures 15-17, respectively. In all cases, the cracks first initiate from the interfaces and merge into a macro crack which crosses the RVE. We first consider only the first two loads in the optimization process, i.e., N = 2. The final design is shown in Figure 18a. We can observe a complex, heterogeneous re-distribution of the inclusion phase. The crack paths for the optimized design under loads 1 and 2 are depicted in Figures 19 and 20.  The force-displacement curves for loads 1 and 2, for initial and optimized designs, are shown in Figure 21a,b, respectively. A gain of G = 38.8% improvement of the fracture resistance is achieved as compared with the initial design. We now consider the three loads in the optimization process, i.e., N = 3. The obtained new distribution of inclusions is shown in Figure 18b. As compared to the previous optimized design, the new topology is more complex here, with a kind of multiscale phase distribution. As shown in Figure 22, the new distribution of inclusions creates a first concentration of cracks within the inclusions, before there is a merging of the different cracks. Comparisons between responses of initial and optimized designs for the three loads are depicted in Figure 23. Here, the fracture energy has been improved by 23.8% for case 1, 23.7% for case 2, and 31.2% for case 3. Compared with the design for considering only the first two loading cases in the previous example, the improvement of fracture energy has been reduced; however, we have considered three objectives herein, and the average increase of the required fracture energy is 27.1%.

Conclusions
In this paper, we have investigated the topology optimization of particle-matrix composites with the objective to maximize its fracture resistance. The framework, initially proposed in [34], combined a phase field method to fracture with interfacial damage and BESO topology optimization, taking into account the full fracture initiation and propagation process, has been extended to periodic composites and multi-objective optimization. We have shown evidence of significant fracture resistance improvements in composite structures under various loading conditions. Size effects have been investigated, showing similar trends as in classical fracture experiments in quasi-brittle structures. More importantly, we have shown that the topology optimization can be performed on a structure with a small number of inclusions at low computational costs. Then, using the optimized topology for larger structures still brings a large gain in the fracture resistance. Finally, we have investigated the effects of multiple loads such that a unit cell of material requires maximizing its fracture for the different possible loads.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Linearization and Discretization of the Mechanical Problem
Classical FEM discretization is used for discretizing the nodal displacement and nodal displacement increments as: u with n 1 and n 2 being the xand y-components of the normal vector to the interface Γ I , andB u is a matrix of shape functions derivatives such that The smoothed jump strain at the interfaces is defined by Using first-order Taylor expansion of (11) and assuming a staggered solving procedure (i.e., that d is fixed when the mechanical problem is solved), the linearized problem is given by where D v f (.) is the directional derivative, and u k is the displacement solution from the k-th iteration.
The displacements at the current iteration are updated according to Using (11)  where P ± is such that [ε ± ] = P ± [ε].
With the above FEM discretization, the linear tangent problem reduces to the following linear system of algebraic equations K tan ∆ũ = −R(ũ k ) (A16)

Appendix B. Discretization of the Phase Field Problem
Using FEM discretization, the phase field and phase field gradient in one element are approximated by: where d e is a vector of nodal phase field values in one element, and N d (x) and B d (x) are matrices of shape functions and of shape functions derivatives associated with the phase field variable, respectively. Introducing the above FEM discretization into the weak form (17), the following linear discrete system of equations can be obtained: and