Pengfei Li
email: lipengfeihhu@163.com

Julien Yvonnet

Yi Wu

Improved fracture resistance of 3D-printed elastoplastic structures with respect to their topology and orientation of deposited layers

Keywords: Topology optimization, Phase field method, Anisotropic fracture, Elastoplasticity, 3D printing

3D-printed structures may be characterized by anisotropic fracture behavior because of their layered nature. Depending on the orientation of the sample during the layer deposition, a completely different mechanical response can be obtained, ranging from quasi-brittle to elastoplastic, and with large variations in the maximum stress to failure. In this study, an optimization framework is proposed for 3D-printed samples to maximize their resistance to fracture with respect to both the orientation of the deposited layers during the process and the topology of the sample. To achieve this, a phase-field anisotropic elastoplastic fracture model is combined with a Bidirectional Evolutionary Structural Optimization topology optimization. The model makes it possible to predict the response of the structure until failure with respect to the orientation of the deposited layers in the 3D-printing process and then optimize this orientation to maximize the mechanical response. A large increase in fracture resistance can be obtained by optimizing the orientation, and a significant increase in fracture resistance can be achieved using the present nonlinear anisotropic topology optimization compared with the use of linear topology optimization.

Introduction

Currently, three-dimensional (3D)-printing technologies are revolutionizing the design and manufacturing of products in many engineering fields, including aerospace [START_REF] Seabra | Selective laser melting (slm) and topology optimization for lighter aerospace componentes[END_REF][START_REF] Guanghui | An aerospace bracket designed by thermo-elastic topology optimization and manufactured by additive manufacturing[END_REF], biotechnology [START_REF] Kumar | Topology optimization and 3D printing of large deformation compliant mechanisms for straining biological tissues[END_REF], civil engineering [START_REF] Vantyghem | 3D printing of a post-tensioned concrete girder designed by topology optimization[END_REF][START_REF] Suiker | Mechanical performance of wall structures in 3D printing processes: Theory, design tools and experiments[END_REF], and automotive engineering [START_REF] Chuang | Topology optimization with additive manufacturing consideration for vehicle load path development[END_REF][START_REF] Reddy | Application of topology optimization and design for additive manufacturing guidelines on an automotive component[END_REF]. The advantages of these techniques include the possibility of considering new shapes that cannot be obtained directly with classical manufacturing processes. With the rapid development of topology optimization (TO), new designs can be used to optimize the ratio of weight to mechanical properties. However, 3D-printing processes may induce additional defects in the manufactured samples. For example, in Fused Deposition Modeling from thermoplastic materials [START_REF] Masood | Intelligent rapid prototyping with fused deposition modelling[END_REF] or Selective Laser Sintering [START_REF] Lü | Selective laser sintering[END_REF] from polymer/metal powder, the parts are built layer by layer (see Fig. 1). The presence of pores or heterogeneities between the layers may induce a strongly anisotropic behavior, and, more specifically, a preferential orientation of cracks during the damage of the samples. As a result, the apparent resistance of the sample may be completely different depending on the orientation of the deposited layers and the external load, as shown experimentally e.g. in [START_REF] Zou | Isotropic and anisotropic elasticity and yielding of 3D printed material[END_REF][START_REF] Zhu | Backbone cup-a structure design competition based on topology optimization and 3d printing[END_REF][START_REF] Huynh | Property-map of epoxy-treated and as-printed polymeric additively manufactured materials[END_REF][START_REF] Somireddy | Anisotropic material behavior of 3d printed composite structures-material extrusion additive manufacturing[END_REF][START_REF] Nurizada | Induced anisotropy in the fracturing behavior of 3d printed parts analyzed by the size effect method[END_REF]. Therefore, it is important to consider the orientation of the layers in the design of 3D-printed workpieces. However, this point has rarely been considered in the literature.

The use of TO [START_REF] Bendsøe | Optimal shape design as a material distribution problem[END_REF][START_REF] Rozvany | Generalized shape optimization without homogenization[END_REF][START_REF] Bendsøe | Material interpolation schemes in topology optimization[END_REF][START_REF] Wang | A level set method for structural topology optimization[END_REF][START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF][START_REF] Xie | A simple evolutionary procedure for structural optimization[END_REF][START_REF] Rozvany | A critical review of established methods of structural topology optimization[END_REF][START_REF] Van Dijk | Level-set methods for structural topology optimization: a review[END_REF][START_REF] Zhu | Topology optimization in aircraft and aerospace structures design[END_REF][START_REF] Deaton | A survey of structural and multidisciplinary continuum topology optimization: post 2000[END_REF][START_REF] Sigmund | Topology optimization approaches[END_REF][START_REF] Djourachkovitch | Multiscale topology optimization of 3D structures: A microarchitectured materials database assisted strategy[END_REF][START_REF] Valm | Topology optimization of imperfect frame structures with improved manufacturability[END_REF] for increasing the fracture resistance of structures is an emergent and exciting topic in computational mechanics. In the first series of studies, linear fracture mechanics objectives in TO, such as the elastic energy released by the crack in tension [START_REF] Challis | Fracture resistance via topology optimization[END_REF] or stress intensity factors [START_REF] Kang | Topology optimization considering fracture mechanics behaviors at specified locations[END_REF], were used. In another group of studies, the incremental damage response of the structure during a full load, from initiation to damaged/cracked structures, was considered using various crack modeling techniques, such as damage, enhanced gradient damage, or extended finite-element method techniques [START_REF] Amir | Reinforcement layout design for concrete structures based on continuum damage and truss topology optimization[END_REF][START_REF] Amir | A topology optimization procedure for reinforced concrete structures[END_REF][START_REF] Kato | Multiphase layout optimization for fiber reinforced composites considering a damage model[END_REF][START_REF] Hilchenbach | Optimization of multiphase structures considering damage[END_REF][START_REF] James | Failure mitigation in optimal topology design using a coupled nonlinear continuum damage model[END_REF][START_REF] Liu | Multi-material topology optimization considering interface behavior via XFEM and Level-Set method[END_REF][START_REF] Li | Design of fracture resistant energy absorbing structures using elastoplastic topology optimization[END_REF][START_REF] Li | Topology optimization of energy absorbing structures with maximum damage constraint[END_REF][START_REF] Li | Failure resistant topology optimization of structures using nonlocal elastoplastic-damage model[END_REF][START_REF] Wu | Comprehensive implementations of phase-field damage models in abaqus[END_REF][START_REF] Seleš | Crack propagation prediction in heterogeneous microstructure using an efficient phase-field algorithm[END_REF][START_REF] Zhang | Topological design of all-ceramic dental bridges for enhancing fracture resistance[END_REF]. In [START_REF] Suresh | Topology optimization for transversely isotropic materials with high-cycle fatigue as a constraint[END_REF], Suresh et al. developed a TO method for designing transversely isotropic materials (e.g., additive manufactured components) in which high-cycle fatigue constraints are considered. Recently, the phase field method [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The Variational Approach to Fracture[END_REF][START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF][START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF][START_REF] Nguyen | Initiation and propagation of complex 3D networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF][START_REF] Zhou | Phase field modelling of crack propagation, branching and coalescence in rocks[END_REF][START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF][START_REF] Mandal | A length scale insensitive anisotropic phase field fracture model for hyperelastic composites[END_REF][START_REF] Zhang | Phase field modeling of fracture in fiber reinforced composite laminate[END_REF] (also called the variational approach to fracture in the literature) for simulating the quasi-brittle fracture of structures and materials was combined with TO [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF][START_REF] Da | Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage[END_REF][START_REF] Da | Topology optimization for maximizing the fracture resistance of periodic quasi-brittle composites structures[END_REF][START_REF] Russ | Topology optimization for brittle fracture resistance[END_REF][START_REF] Russ | A novel topology optimization formulation for enhancing fracture resistance with a single quasibrittle material[END_REF][START_REF] Wu | Level-set topology optimization for maximizing fracture resistance of brittle materials using phase-field fracture model[END_REF][START_REF] Wu | Topology optimization for enhanced dynamic fracture resistance of structures[END_REF][START_REF] Li | A SIMP-Phase field topology optimization framework to maximize quasi-brittle fracture resistance of 2D and 3D composites[END_REF]. This is an important step for TO, with the objective of increasing the fracture resistance. The phase field can be used in regular meshes, which is advantageous in TO, and the approach can address the initiation, propagation, and interaction of complex, multiple 2D and 3D fractures in arbitrary geometric configurations.

Deposited layers

In this article, a TO framework is proposed to improve the fracture resistance of 3D-printed structures with anisotropic elastoplastic behavior, considering the orientation of the layers induced by the 3D printing process. The optimization is performed with respect to the fracture resistance of the sample, defined as the cumulative external work of the structure during the entire loading cycle (from zero load until full failure). Thus a model proposed previously [START_REF] Li | Anisotropic elastoplastic phase field fracture modeling of 3D printed materials[END_REF] is used to consider the fracture behavior in layered elastoplastic structures and is combined with the bidirectional evolutionary structural optimization (BESO) TO method.

Several numerical examples are employed to illustrate the methodology and the significance for the design of 3D-printed structures. A topology optimization with objective of fracture resistance maximization, taking into account an orientation-dependent fracture behavior due to the 3D printing process is proposed here for the first time to our best knowledge.

Phase field modeling of anisotropic elastoplastic damage in 3D printed materials

In this section, an anisotropic elastoplastic fracture model developed previously [START_REF] Li | Anisotropic elastoplastic phase field fracture modeling of 3D printed materials[END_REF] for 3D printed parts is described briefly. The proposed model is intended to describe any elastoplastic material associated to a process for which a layered structure is induced. Such structure has been mainly observed e.g. in polymer or metallic SLS. Note that Fused Deposition Modeling (FDM)/Fused Filament Fabrication (FFF) may lead to more complex local anisotropic behaviors due to the deposition path and may require extensions as compared to the present framework. The model relies on the phase field method and is based on a variational principle with a regularized description of discontinuities (see [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF] and the review in [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF] for classical quasi-brittle models). In the present model, two variables are used to describe the damage: one associated with bulk damage d and the other describing the distribution of micro interfacial damage along the layers induced by the 3D-printing process, denoted by α, which is associated with an anisotropic damage in the direction of the weak planes (see Fig. 1). Considering a solid defined in a domain Ω ⊂ R 2 , with boundary ∂Ω, the total energy of the system is defined by

W = Ω ψ e + ψ p + ψ d + ψ α dΩ - ∂Ω F F • udS, (1)
where ψ e denotes the elastic strain density energy function, ψ p denotes the plastic strain density energy function, ψ d is the bulk layer fracture density function, ψ α is the micro interfacial fracture density function , ∂Ω F denotes the Neumann boundary, and F denotes the prescribed forces. We define

ψ e = 1 2 (ε -ε p ) : C (α, d) : (ε -ε p ), (2)
where ε = 1 2 (∇u + ∇ T u) and ε p denote the total and plastic strain tensors, respectively, and the elastic strain is given by

ε e = ε -ε p . ( 3
)
Plastic incompressibility is assumed, i.e., Tr(ε p ) = 0, where Tr(.) is the trace operator. In (2), C (α, d) is the damage-dependent elasticity tensor, defined as

C (α, d) = g (d) C α (α) , (4)
where we have chosen

g (d) = (1 -d) 2 , (5)
and the matrix form associated with C α can be expressed as a function of the layer orientation angle θ (see Fig. 1) as

C α (α) = Q(θ) T C (α)Q(θ) (6)
with the transformation matrix

Q(θ) =      c 2 s 2 2cs s 2 c 2 -2cs -cs cs c 2 -s 2      , (7)
and c = cos(θ), s = sin(θ). In ( 6) C (α) is the effective elasticity matrix in the frame of the layers, given in the present model by the following quadratic function of α:

C (α) =      a 1 + a 2 (1 -α) 2 b 1 + b 2 (1 -α) 2 0 b 1 + b 2 (1 -α) 2 c 1 + c 2 (1 -α) 2 0 0 0 d 1 + d 2 (1 -α) 2      , (8)
where a 1 , a 2 , ..., d 2 are material constants that can be identified either experimentally or numerically by RVE (Representative Volume Element) calculations (see [START_REF] Li | Anisotropic elastoplastic phase field fracture modeling of 3D printed materials[END_REF]).

The second-order Cauchy stress tensor is defined as

σ = ∂ψ e ∂ε e = (1 -d) 2 C α : ε e . (9)
In (1), ψ p is here chosen as

ψ p = σ y p + 1 2 Hp 2 . ( 10
)
This describes a linear isotropic hardening law for the plasticity evolution, where σ y and H > 0 denote the yield stress and hardening modulus parameters, respectively. In [START_REF] Zou | Isotropic and anisotropic elasticity and yielding of 3D printed material[END_REF], p is an equivalent plastic strain described by the evolution equation

ṗ = 2 3 εp , (11)
where ( .) denotes the time derivative. In (1), ψ α and ψ d are defined as

ψ α = g α c γ α (α) , ψ d = g d c γ d (d) , (12)
where g α c and g d c represent the fracture toughness with respect to α and d, respectively. In 3D printed structures, interphase slip does not seem to be the dominant mechanism for failure and is then not taken into account in the present model. However, models such as cohesive tractions have have been included in previous phase field models [START_REF] Li | An extension of the phase field method to model interactions between interfacial damage and brittle fracture in elastoplastic composites[END_REF] and could be added in future extension.

The bulk layer crack surface density function for d is defined by

γ d (d) = d 2 2 d + d 2 ∇d • ∇d, (13)
where d is a length-scale regularization parameter associated with d.

To enforce a preferential direction related to the damage induced by the micro interfacial failure, an anisotropic crack surface density function is introduced [START_REF] Clayton | Phase field modeling of directional fracture in anisotropic polycrystals[END_REF][START_REF] Nguyen | Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials[END_REF]:

γ α (α) = α 2 2 α + α 2 ω α : (∇α ⊗ ∇α) , ( 14
)
where α is the length scale regularization parameter associated with α, and

ω α (θ) = 1 + ξ α (1 -n α (θ) ⊗ n α (θ)) , (15)
where ξ α 1 is a parameter used to penalize the damage along the direction normal to n α (see Fig. 1). Employing the variational framework for fracture, as introduced in [START_REF] Li | An extension of the phase field method to model interactions between interfacial damage and brittle fracture in elastoplastic composites[END_REF][START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF][START_REF] Rodriguez | A variational approach to the phase field modeling of brittle and ductile fracture[END_REF], i.e., minimizing the energy W with respect to u, d, and α with ḋ ≥ 0 and α ≥ 0, yields a set of coupled equations to determine the fields d(x) and α(x) and the displacements u(x) and ∀x ∈ Ω:

         ∇ • σ = 0 in Ω, u = ū on ∂Ω u , σn = F on ∂Ω F , ( 16
)
where n is the unitary normal vector to ∂Ω, ∂Ω u is the Dirichlet boundary, and u are the prescribed displacements on ∂Ω u .

         g d c d d -2 d ∆d = 2 (1 -d) H in Ω, ∇d • n = 0 on ∂Ω, d = 1 on Γ d . ( 17
)
The history strain energy density function is employed to prescribe the irreversibility of the bulk damage field [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] and is defined by

H = max s∈[0,t] ψ e 0 (u, s) with ψ e 0 = 1 2 ε e : C α : ε e . ( 18
)
The irreversibility condition α ≥ 0 is enforced algorithmically (see [START_REF] Li | Anisotropic elastoplastic phase field fracture modeling of 3D printed materials[END_REF]). The equations associated with the micro interfacial damage distribution α are given by

         1 2 (1 -d) 2 ε e : ∂C α ∂α : ε e + g α c α α -g α c α ω α : ∇∇α = 0 in Ω, ∇α • n = 0 on ∂Ω, α = 1 on Γ α , (19)
where (∇∇α) ij = ∂ 2 α ∂xi∂xj . Considering J 2 -plasticity, the yield criterion is given by

F p = 3 2 σ dev -(σ y + Hp) ≤ 0. (20)
Finally, εp = ṗ 3 2

σ dev σ dev with ṗ ≥ 0, (21)
where σ dev := σ -1 3 Tr (σ) 1 is the deviatoric stress tensor. Equations ( 16)-( 17)-( 19) are discretized by finite elements and solved in a staggered scheme, i.e., through a fixed-point algorithm where the mechanical problem ( 16) and the phase field problems ( 17) and ( 19) are solved separately. The mechanical problem ( 16) is nonlinear owing to the elastoplastic behavior and is solved using a Newton algorithm. The equivalent plastic strain is updated using a return-mapping algorithm [START_REF] Simo | Computational Inelasticity[END_REF]. More details about the solving procedure can be found in [START_REF] Li | Anisotropic elastoplastic phase field fracture modeling of 3D printed materials[END_REF].

Topology optimization formulation

In this section, a topological optimization is presented for maximizing the fracture resistance of elastoplastic anisotropic material described in Section 2. The method is based on the BESO method [START_REF] Huang | Evolutionary topology optimization of continuum structures: methods and applications[END_REF][START_REF] Fritzen | Topology optimization of multiscale elastoviscoplastic structures[END_REF][START_REF] Xia | Evolutionary topology optimization of elastoplastic structures[END_REF]. In the proposed procedure, the fracture resistance is defined as the cumulative external work of the structure for a complete loading, from the initiation of cracks to failure, under the constraint of material volume fraction.

Optimization problem statement

The design domain Ω is discretized into N e finite elements, and each element e is associated with a topology design variable ρ e . The vector {ρ} = ρ 1 , ρ 2 , ..., ρ Ne containing the discrete values for all elements of the mesh is defined. Following the well-known BESO formulation [START_REF] Huang | Evolutionary topology optimization of continuum structures: methods and applications[END_REF], the design variables take the values ρ e ∈ [0; 1], ρ e = 1 for solid material, whereas ρ e = 0 corresponds to voids.

For stability considerations, displacement-controlled loading is adopted. For a prescribed displacement load, the fracture resistance maximization is equivalent to the maximization of the mechanical work. The total (cumulated) mechanical work J is approximated by numerical integration using the trapezoidal rule,

J ≈ 1 2 n load n=1 f n ext + f n-1 ext T ∆u n , (22)
where n load is the total number of displacement increments, ∆u n denotes the prescribed load increment, and f n ext is the external nodal force vector at load n. During optimization, the material volume fraction is prescribed. Then, the optimization problem discretized with N e elements can be formulated as follows [START_REF] Fritzen | Topology optimization of multiscale elastoviscoplastic structures[END_REF][START_REF] Xia | Evolutionary topology optimization of elastoplastic structures[END_REF][START_REF] Huang | Topology optimization of nonlinear structures under displacement loading[END_REF]:

Maximize : J (23) subject : r = 0 (24) Ne e=1 ρ e v e /( Ne e=1 v e ) = f (25
)

ρ e = 0 or 1, e = 1, 2, ..., N e . (26)
Here, v e is the volume (area in 2D) of the e-th element, and f is the material volume fraction. Vector r denotes nodal residual forces

r = f ext - Ne e=1 ρ e Ωe B T u σdΩ (27)
where Ω e denotes the element domain, and B u is the matrix of the displacement shape function derivatives. In this model, the discrete topology design variable ρ e = 0; 1 indicates the existence of the associated solid element e, whereas the constitutive behavior is assumed to be independent of ρ e . This assumption does not require the definition of supplementing pseudo-relationships between intermediate densities and their constitutive behaviors, as in the case of models of the solid isotropic material with penalization type [START_REF] Russ | Topology optimization for brittle fracture resistance[END_REF][START_REF] Bogomolny | Conceptual design of reinforced concrete structures using topology optimization with elastoplastic material modeling[END_REF], thus avoiding the construction of C α in Eq. ( 8) for all intermediate densities and being well-suited for the present phase field model.

Sensitivity analysis

The derivation of the sensitivity requires the use of the adjoint method (e.g., [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF][START_REF] Buhl | Stiffness design of geometrically nonlinear structures using topology optimization[END_REF]). Assuming that the displacement problem has been solved, the Lagrangian is introduced:

J = 1 2 n load n=1 f n ext + f n-1 ext T ∆u n + (λ n 1 ) T r n + (λ n 2 ) T r n-1 . (28)
Here, r n and r n-1 are the residuals of ( 27) at the n-th and (n -1)-th load increments, respectively. Furthermore, λ n 1 and λ n 2 are Lagrange multipliers with the same dimensions as the displacement vector u. The final objective derivatives can be obtained as

∂J ∂ρ e = - 1 2 n load n=1 (λ n 1 ) T Ωe B T u σ n dΩ + (λ n 2 ) T Ωe B T u σ n-1 dΩ (29)
with

λ n 1,E = -∆u n E , λ n 1,F = K n tan,FF -1 K n tan,FE ∆u n E (30) λ n 2,E = -∆u n E , λ n 2,F = K n-1 tan,FF -1 K n-1 tan,FE ∆u n E ( 31
)
where we define two sets among all degrees of freedom: essential (index E; associated with Dirichlet boundary conditions) and free (index F) nodal values, e.g.. For a vector v and a matrix M we have:

v ∼   v E v F   and M ∼   M EE M EF M FE M FF   . (32)
The tangent stiffness matrix of the nonlinear mechanical system at the equilibrium of the m-th load increment is defined as

K m tan = - ∂r m
∂u m (33)

BESO updating scheme

The BESO variable updating scheme for this work follows part of the basic procedure presented in [START_REF] Xia | Evolutionary topology optimization of elastoplastic structures[END_REF][START_REF] Huang | Topology optimization of nonlinear structures under displacement loading[END_REF]. In addition, a normalization strategy introduced previously [START_REF] Zhou | A normalization strategy for beso-based structural optimization and its application to frequency response suppression[END_REF] is used to improve the TO convergence.

The target material volume fraction f is reached by removing the material from the entire design domain, i.e.,

f k = max f , (1 -c er ) f k-1 , (34)
where f k is the material volume fraction at the k-th design iteration, and c er is an evolutionary ratio that determines the percentage of material to be removed from the design of the previous iteration. The actual volume fraction obtained during the design iterations usually deviates slightly from the prescribed value.

Once the target material volume fraction f is reached, the optimization algorithm only modifies the topology but keeps the volume fraction constant. At each design iteration, the sensitivity numbers, which denote the relative ranking of the element sensitivities, are used to determine the material removal and addition. Thus, without changing the relative ranking of the element sensitivities, a normalization strategy introduced elsewhere [START_REF] Zhou | A normalization strategy for beso-based structural optimization and its application to frequency response suppression[END_REF] is employed to force the element sensitivities φ e to distribute in an interval of [0, 1], i.e.,

φ e = X e -X min X max -X min , e = 1, 2, ..., N e (35)
where X e is the element sensitivity obtained in [START_REF] Kang | Topology optimization considering fracture mechanics behaviors at specified locations[END_REF], and X max and X min denote the maximum and minimum values of X e , e = 1, 2, ..., N e , respectively. With this normalization strategy, the magnitude of the sensitivities can be efficiently unified to the same order while maintaining a history-averaging technique, thus avoiding the possible numerical instabilities of the TO process caused by the strong nonlinearity.

After normalization, the new sensitivities are smoothed to avoid mesh dependency by using a filtering scheme [START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF]

φ e = Ne j=1 w ej φ j Ne j=1 w ej , (36)
in which w ej is a linear weight factor

w ej = max (0, r min -∆ (e, j)) (37)
which is determined by the prescribed filter radius r min and the element center-to-center distance ∆ (e, j) between element e and j. Then a history-averaging technique [START_REF] Huang | Topology optimization of nonlinear structures under displacement loading[END_REF][START_REF] Huang | Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method[END_REF][START_REF] Huang | Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials[END_REF] is employed to improve the convergence of the TO process.

Numerical examples

In this section, several illustrative examples of 3D-printed elastoplastic structures subjected to cracking are presented. The proposed framework is used to design the optimal topology of the structures to resist cracking, considering the orientation of deposited layers during the process. In this work, the material coefficients in C (see Eq. ( 8)) were identified numerically by RVE calculations from an RVE corresponding to a layered RVE (see Fig. 2 (b)). The material parameters for the matrix were chosen as follows: E = 10 GPa, ν = 0.25. We denote the Young's modulus of the interphase E I with respect to the micro interfacial damage parameter α according to:

E I (α) = G(α)E, G(α) = (1 -α) 2 χ -(χ -1) (1 -α) 2 with χ ≥ 1. ( 38
)
The Poisson's ratio ν I for the interphase is assumed to be independent on α: ν I = ν. The evolution of the identified componenents of C is shown in Fig. 2. The corresponding coefficients in C are: In the present model, the defects between the layers induced by the 3D printing process are taken into account by the presence of the weak interphase (see Fig. 2 (a)-(b)) and the related interfacial damage model, which induces the anisotropic fracture behavior. The morphology and distributions of pores at the interfaces between the deposited layers can be explicitly described in the RVE (see e.g. [START_REF] Li | Anisotropic elastoplastic phase field fracture modeling of 3D printed materials[END_REF]). If experimental observations are available to provide information about the layers and porous interphase thicknesses, they can be included in the description of the RVE. If not, the coefficients in (8) can be fitted from macroscopic calculations on samples for different orientations. The parameters of the numerical model are as follows. In all examples, a uniform mesh of quadrilateral bilinear elements was employed. The characteristic length scale parameters d and α were both set to be twice the typical finite-element size d = α = 2h e . The parameter c er in [START_REF] James | Failure mitigation in optimal topology design using a coupled nonlinear continuum damage model[END_REF] was set to 0.02. The filter radius r min was set to r min = 6h e , and ξ α = 30 was used.

a 1 = 9, a 2 = 1.67, b 1 = c 1 = d 1 = 0, b 2 = 2.
4.1. Topology optimization of a L-shaped bracket: Type I As a first example, a bracket defined in an L-shaped design domain is considered, as shown in Fig. 3(a). The domain was uniformly discretized into 3600 square bilinear elements with h e = 0.5 mm.

We have shown in our previous work [START_REF] Li | Anisotropic elastoplastic phase field fracture modeling of 3D printed materials[END_REF] that even though the results depend on the grid size, the results are convergent when the grid is refined for the proposed anisotropic fracture model. In addition, a common point for any topology optimization technique is that the final topology depends on the size of the grid used during the topology optimization procedure. Analyzing the dependence on the final geometries with respect to the grid size is an interesting point, and could be analyzed in future studies.

The boundary conditions were as follows: on the upper end, the x, y-displacements were fixed. On the left side, the vertical displacements were prescribed with fixed incremental displacement loads of ∆U = 0.002 mm. The target volume fraction of the solid material in the final structure was set to f = 60%. First, the fracture process in the initial domain (without TO) was investigated. Fig. 4 shows the final interfacial damage (α), bulk crack (d), and equivalent plastic strain (p) fields for three different layer orientations. When θ = 0 • and θ = 60 • , the failure of the structure was mainly caused by the propagation of micro interfacial cracks (represented by the α distribution) along the layer orientation. However, for θ = 90 • , the activation of the bulk layer damage d indicates that the fracture of the bulk layers was the main mechanism that caused the failure of the structure. One can also observe a significant distribution of the equivalent plastic strain field in that case, which did not occur for the other orientations. It is worth noting that different fracture modes can be observed in this example with respect to the orientation angle. For example, mode I occurs for θ = 0 • , while for θ = 60 • , mode II occurs. The corresponding load-displacement curves for different layer orientations are shown in Fig. 5. The peak load increases with an increase in the layer orientation angle θ. For θ = 90 • , a long plastic stage can be observed. This first investigation illustrates that the proposed anisotropic phase field model can consider the orientation of the layers to describe the fracture process and describe the transition from quasi-brittle to elastoplastic behavior with damage according to the layer orientation. In addition, even without considering TO, the proposed model makes it possible to optimize the orientation of the sample in the 3D printing process with respect to the layer deposition to maximize the mechanical strength. Figs. 678show, respectively, the design history using the present TO framework and the corresponding final crack patterns for θ = 0 • , θ = 60 • , and θ = 90 • . The objective value fluctuates and globally decreases as the material is removed until the prescribed material volume fraction is reached. In Fig. 7, the objective value slightly improves while the material is removed from initial 100% to 80%. This indicates that, for the same fracture resistance performance, the required amount of material can be saved using this proposed TO process. Fig. 8 shows that oscillations occur in the transition from elastoplastic to quasi-brittle behaviors before convergence. Good convergence of the proposed scheme was achieved in all cases.

J=58.7 mJ

Fig. 9 shows the final interfacial damage (α), bulk crack (d), and equivalent plastic strain (p) fields for the corresponding optimized final designs, for different orientations of the layers. For all orientations, the main fracture mechanism appears to be the interfacial damage because the field d does not significantly develop. The corresponding load-displacement curves in Fig. 10 reveal that the response for the case θ = 90 • induces a larger resistance to fracture, but there are no significant differences in the quasi-brittle nature of the response. Compared to the load-displacement curve for initial design domain shown in Fig. 5, the loaddiaplacement curve for θ = 90 • optimized design does not show a significant plastic stage. Several reasons can cause the difference: (i) the optimized design is 60% volume fraction of the initial design domain, which is obtained by removing material from the initial design domain; (ii) the optimized design fails earlier because of the interfacial crack initiation and propagation which are much easier than the bulk crack initiation and propagation.

p p To evaluate the gains provided by the proposed framework, an analysis was performed. First, a classical TO of the structure was performed considering a linear behavior with the material volume fraction constraint and the same parameters. The obtained topologies are called "isotropic linear" (IL). However, the fully nonlinear framework TO described in section 3 was used to obtain the final design, called "anisotropic elastoplastic fracture" (AEF). Then, both designs were used to perform a fracture simulation using the anisotropic fracture model described in section 2. Fracture patterns for different orientations are shown in

Topology optimization of a L-shaped bracket: Type II

In the next example, the problem defined in Fig. 3(b), which depicts the geometry of the design domain and the boundary conditions, was considered. The geometry is similar to that in the previous example, but the boundary conditions induce a different distribution of strains within the structure. Again, the domain was uniformly discretized into 3600 square bilinear elements with h e = 0.5 mm. The boundary conditions were as follows: on the lower end (y = 0), the x, y-displacement was fixed. On the upper end, the horizontal displacements were prescribed with fixed incremental displacement loads ∆U = 0.002 mm. The target volume fraction in the optimized structure was set to f = 60%.

The final interfacial damage, bulk crack, and equivalent plastic strain field in the initial design domain for the three different layer orientations are shown in Fig. 13. In this case, when θ = 0 • and θ = 120 • , the failure of the structure was mainly caused by the propagation of micro interfacial cracks (represented by the α distribution) along the layer orientation. However, for θ = 90 • , the activation of the bulk layer damage d indicates that the fracture of the bulk layers was the main mechanism causing the failure of the structure in this case. The corresponding load-displacement curves for different layer orientations are shown in Fig. 14. Here again, the orientation θ = 90 • exhibits a very different response, with a flattening of the curve, which can be attributed to the high development of plasticity. In this case, the proposed model makes it possible to select the optimal direction for improved mechanical resistance as θ = 90 • . It is not always intuitive to define the optimal angle orientation by simply considering the geometry of the sample and loading orientation. Thus, such a model can greatly aid the design of 3D-printed samples. To illustrate this point better, the maximum load before failure and the cumulative external work, interpreted here as the fracture resistance, were computed as a function of the orientation angle θ for both structures (Type I bracket and Type II bracket structures), as shown in Figs. 15 and16, respectively. These figures show that there is an optimal resistance with respect to the angle, which corresponds to approximately θ = 95 • for the Type I bracket and θ = 80 • for the Type II bracket. The experimental observations in [START_REF] Zou | Isotropic and anisotropic elasticity and yielding of 3D printed material[END_REF][START_REF] Kiendl | Controlling toughness and strength of FDM 3D-printed PLA components through the raster layup[END_REF] also illustrate that the fracture resistance for 3D-printed samples significantly depends on the printing angle, and the maximal fracture resistance tends to occur around the crack-vertical printing orientation. For the Type I bracket, an increase in the external work (fracture resistance) of approximately a factor 4 can be obtained by appropriately choosing the orientation of the sample. Such optimization clearly brings large gains at no additional cost.

Figs. 17-19 show, respectively, the design evolutions and associated final crack patterns for θ = 0 • , θ = 90 • , and θ = 120 • . for the θ = 0 • and θ = 120 • cases, the objective values gently improved while the material was removed for the first 13 design iterations, and both had a sharp decrease and reached the final convergence. This resulted from material removal in the path of the crack in the sample, as two "legs" were formed on the bottom of the sample (see, for example, Fig. 17). In all cases, global convergence was achieved. Fig. 20 shows the final interfacial damage, bulk crack, and equivalent plastic strain field for the corresponding final designs. The final designs for θ = 0 • and θ = 120 • failed mainly because of micro interfacial cracking mechanisms. However, for θ = 90 • , the activation of the bulk damage d indicates that the mechanism was, in this case, the cracking of bulk layers. The corresponding load-displacement curves are shown in Fig. 21. In this case, it is obvious that the orientation angle θ = 90 • is again the optimal orientation of the sample in the 3D-printing process to maximize the mechanical resistance. As in Section 4.1, the improvements provided by the proposed analysis compared with a design obtained by TO were demonstrated under a linear assumption. Fracture patterns for different orientations are shown in Figs. [START_REF] Van Dijk | Level-set methods for structural topology optimization: a review[END_REF]. The corresponding comparisons of the load-displacement curves and design objective values for the IL and AEF designs are shown in Figs. [START_REF] Zhu | Topology optimization in aircraft and aerospace structures design[END_REF]. For θ = 0 • and θ = 120 • , a reasonable improvement in the external work ( 20% and 18%, respectively) was found. However, for the case θ = 90 • , an improvement of 51% was obtained over the simple IL analysis, fully justifying the consideration of the anisotropic elastoplastic damage model in the proposed TO framework.

J=29.0 mJ

J=31.5 mJ

Conclusion

A TO framework that considers the orientation of the deposited layers in 3D-printed structures and an elastoplastic anisotropic fracture model was proposed. The model can describe the transition from quasibrittle to elastoplastic fracture behavior according to the layer orientation. A large increase in the fracture resistance (cumulated external work and maximum stress to failure) can be obtained by optimizing the deposited layer orientation (angle in 2D) in the 3D printing process. Furthermore, the present framework makes a significant increase in the fracture resistance possible for a minimal amount of material owing to the TO process in comparison with simpler, isotropic linear TO. Such a framework should facilitate the design of 3D-printed parts with increased fracture resistance while lowering costs because of material saving, which is of major interest to industry. Experimental validations of the proposed simulations constitute an exciting perspective for this work.

Figure 1 :

 1 Figure 1: Different orientations of deposited layers during a 3D printing SLS (Selective Laser Sintering) process for different orientations of the sample during the process.

Figure 2 :

 2 Figure 2: (a) Layer structure in a 3D printed material. The layers are separated by porous interphases (white dotted lines); (b) numerical mesh of a periodic cell associated with the material; (c) evolution of the effective elastic moduli computed numerically as a function of the degradation of the interface α.

  67, c 2 = 10.67, and d 2 = 4.

Figure 3 :

 3 Figure 3: Two 3D printed structures; the angle θ indicates the orientation of the layers induced by the 3D printing process. (a): Type I L-shaped bracket; (b) Type II L-shape bracket.

Figure 4 :

 4 Figure 4: Interfacial damage (α), bulk damage (d) and equivalent plastic strain (p) in the Type I L-shaped bracket structure subjected to a load on the y-direction. The angle θ indicates the orientation of the layers induced by the 3D printing process: (a) θ = 0 • ; (b) θ = 60 • ; (c) θ = 90 • .

5 :

 5 Load-displacement curves for a Type I L-shaped bracket (non-optimized geometry). The angle θ denotes the orientation of the layers in the structure due to the 3D printing process.J=55.7 mJ

Figure 6 :

 6 Figure 6: Evolution of the design during topology optimization of the Type I L-shape bracket, for an orientation of the 3D printed layers θ = 0 • . Final crack patterns (interfacial damage field α) are depicted for different steps of the optimization procedure. The evolutions of the final external work and volume of solid are depicted in blue and red colors, respectively.

Figure 7 :

 7 Figure 7: Evolution of the design during topology optimization of the Type I L-shape bracket, for an orientation of the 3D printed layers θ = 60 • . Final crack patterns (interfacial damage field α) are depicted for different steps of the optimization procedure. The evolutions of the final external work and volume of solid are depicted in blue and red colors, respectively.

Figure 8 :

 8 Figure 8: Evolution of the design during topology optimization of the Type I L-shape bracket, for an orientation of the 3D printed layers θ = 90 • . Final crack patterns (interfacial damage field α and bulk damage d) are depicted for different steps of the optimization procedure. The evolutions of the final external work and volume of solid are depicted in blue and red colors, respectively.

Figure 9 :

 9 Figure 9: Optimized designs for the Type I L-shape structure. The final interfacial damage fields α, bulk crack d and equivalent plastic strains p are indicated for different orientation angles θ of the deposited layers during the 3D printing process: (a) θ = 0 • ; (b) θ = 60 • ; (c) θ = 90 • .

Figs. 11 . 10 :

 1110 The corresponding comparisons of the load-displacement curves and design objective values for the IL and AEF designs are shown in Figs. 12. The gains in fracture resistance with the increase in external work were evaluated for the IL and AEF designs. Increases of 18%, 15%, and 8% were obtained for the orientations θ = 0 • , θ = 60 • , and θ = 90 • , respectively. A higher gain was obtained in the case of θ = 0 • . Load-displacement curves for a Type I L-shaped bracket (optimized design). The angle θ denotes the orientation of the layers in the structure due to the 3D printing process.

Figure 11 :

 11 Figure 11: Final crack patterns (interfacial crack damage α) for optimized designs of the type I L-shape bracket structure."IL" indicates that the topology optimization has been conducted using Isotropic Linear fracture behavior. "AEF" indicates that the topology optimization has been conducted using anisotropic fracture behavior.

Figure 12 :

 12 Figure 12: Load-displacement curves for optimized designs of the type I L-shape bracket structure. The red curves correspond to the optimized design obtained using a classical linear isotropic fracture behavior (IL). The blue curves are obtained using the present fully anisotropic Fracture behavior (AEF). The values J indicate the values of the objective function. The curves are depicted for different orientation angles of deposited layers during the 3D process.

Figure 13 :

 13 Figure 13: Interfacial damage (α), bulk damage (d) and equivalent plastic strain (p) in the Type II L-shaped bracket structure subjected to a load on the y-direction. The angle θ corresponds to the orientation of the layers induced by the 3D printing process: (a) θ = 0 • ; (b) θ = 60 • ; (c) θ = 90 • .

14 :

 14 Load-displacement curves for a Type II L-shaped bracket (non-optimized geometry). The angle θ denotes the orientation of the layers in the structure due to the 3D printing process.

Figure 15 :Figure 16 :

 1516 Figure 15: Maximum load to failure in two structures (Type I and Type II L-shape structures) as a function of the layer orientation angle θ chosen during the 3D printing process.

Figure 17 :Figure 18 :

 1718 Figure 17: Evolution of the design during topology optimization of the Type II L-shape bracket, for an orientation of the 3D printed layers θ = 0 • . Final crack patterns (interfacial damage field α) are depicted for different steps of the optimization procedure. The evolution of the final external work and volume of solid are depicted in blue and red colors, respectively.

Figure 19 :

 19 Figure 19: Evolution of the design during topology optimization of the Type II L-shape bracket, for an orientation of the 3D printed layers θ = 120 • . Final crack patterns (interfacial damage field α) are depicted for different steps of the optimization procedure. The evolution of the final external work and volume of solid are depicted in blue and red colors, respectively.

Figure 20 : 21 :

 2021 Figure 20: Optimized designs for the Type II L-shape structure. The final interfacial damage fields α, bulk crack d and equivalent plastic strains p are indicated for different orientation angles θ of the deposited layers during the 3D printing process: (a) θ = 0 • ; (b) θ = 60 • ; θ = 90 • .

Figure 22 :Figure 23 :

 2223 Figure 22: Final crack patterns (interfacial crack damage α and bulk crack damage d) for optimized designs of the type II L-shape bracket structure. "IL" indicates that the topology optimization has been conducted using Isotropic Linear fracture behavior. "AEF" indicates that the topology optimization has been conducted using anisotropic fracture behavior.

Acknowledgements

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

Conflict of interest statement

On behalf of all authors, the corresponding author states that there is no conflict of interest.