Yi Wu

Julien Yvonnet
email: julien.yvonnet@univ-eiffel.fr

Pengfei Li

Zhi-Cheng He

Topology optimization for enhanced dynamic fracture resistance of structures

Keywords: Topology optimization, Phase-field method, Dynamics, Structural optimization, SIMP, Fracture

A topology optimization framework for improving the dynamic fracture resistance of structures is proposed. The phase field method for fracture is combined with Solid Isotropic Material with Penalization (SIMP) topology optimization. The topology optimization problem is defined as minimizing the fracture energy during the whole dynamic loading process, from initiation of cracks to full failure of the structure, under volume and compliance constraints. Semi analytical expressions of sensitivities in a dynamic context are provided to solve the topology optmization problem efficiently. Numerical examples involving structures subjected to impact loading are investigated.

It is shown that the present framework allows a significant reduction of the fracture energy as compared to designs obtained by static optimization.

Introduction

Since the late 1980s [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF], topology optimization (TO) has evolved into one of the most powerful numerical design methods. In a given design domain, topology optimization approaches [START_REF] Xie | A simple evolutionary procedure for structural optimization[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] Guo | Doing topology optimization explicitly and geometrically-a new moving morphable components based framework[END_REF] generate the optimal topological designs that minimize or maximize an objective function with certain design constraints. The problems solved by TO have gradually evolved from simple linear global optimization problems to localized stress problem [START_REF] Duysinx | Topology optimization of continuum structures with local stress constraints[END_REF], to dynamic problems [START_REF] Zhang | Dynamic topology optimization of piezoelectric structures with active control for reducing transient response[END_REF], nonlinear problems [START_REF] Xia | Recent advances on topology optimization of multiscale nonlinear structures[END_REF], uncertainties [START_REF] Wu | Robust concurrent topology optimization of structure and its composite material considering uncertainty with imprecise probability[END_REF], multidisciplinary integration [START_REF] Zhu | Topology optimization in aircraft and aerospace structures design[END_REF] or towards large-scale problem [START_REF] Aage | Giga-voxel computational morphogenesis for structural design[END_REF], among many others. These achievements have demonstrated the potential of TO in both academic and industrial applications.

Recently, an exciting new branch of TO has emerged to enhance the mechanical resistance to damage or cracks in structures and materials. The pioneering work on this topic might be traced back to Challis et al. [START_REF] Challis | Fracture resistance via topology optimization[END_REF], in which a level-set method was proposed to maximize the structural fracture resistance. In [START_REF] Kang | Topology optimization considering fracture mechanics behaviors at specified locations[END_REF], Kang et al. used a J-integral approach to predict crack opening at predefined locations. In [START_REF] Amir | A topology optimization procedure for reinforced concrete structures[END_REF][START_REF] Amir | Reinforcement layout design for concrete structures based on continuum damage and truss topology optimization[END_REF], non-local damage field was considered for the first time in TO to obtain the optimal mechanical resistance design of concrete structures and their reinforcement. Similarly, Kato and Ramm [START_REF] Kato | Multiphase layout optimization for fiber reinforced composites considering a damage model[END_REF] investigated fiber-reinforced composites considering a damage model, in which the layout of the multi-phase materials was optimized. James and Waisman [START_REF] James | Failure mitigation in optimal topology design using a coupled nonlinear continuum damage model[END_REF] developed a non-local damage-TO-coupled algorithm for failure reduction, in which the maximal damage was constrained, similarly to what is usually done in stress-constrained problems. In [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], Li et al. investigated TO methods involving stored energy while constraining the elastoplastic-damage. More recently, Russ and Waisman [START_REF] Russ | A novel elastoplastic topology optimization formulation for enhanced failure resistance via local ductile failure constraints and linear buckling analysis[END_REF] proposed a method for the structural resistance of both ductile failure and buckling in a new aggregated optimization objective with local ductile failure constraints. Liu et al. [START_REF] Liu | Multi-material topology optimization considering interface behavior via xfem and level set method[END_REF] investigated multi-material fracture resistance TO including cohesive models. Note that the above mentioned studies do not include a complete damage or crack evolution involving the whole loading history. In [START_REF] Zhang | Topological design of all-ceramic dental bridges for enhancing fracture resistance[END_REF], Zhang et al. firstly included a full crack propagation analysis within TO for fracture resistance designs, using the X-FEM [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] method. However, the complexity of XFEM for dealing with initiation and complex cracks configurations strongly restricts its use within TO analysis. In [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF], Xia et al. combined for the first time TO with fracture phase field analysis to maximize the resistance of composite structures. The variational phase field approach to fracture [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] 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] has unique advantages to deal with initiation, propagation of multiple, complex, 3D cracks in possibly regular meshes due to an appropriate regularization process, and is highly compatible with TO analysis. Da et al. [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] Da | Fracture resistance design through biomimicry and topology optimization[END_REF] extended this work to consider fracture resistance enhancement in composite by considering both interfacial and bulk fracture. More specifically, the approaches developed in these works considered a full fracture initiation and propagation within the structure until failure and combined the phase field method with BESO [START_REF] Huang | Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method[END_REF] TO methods. Li et al. [START_REF] Li | A simp-phase field topology optimization framework to maximize quasi-brittle fracture resistance of 2d and 3d composites[END_REF] extended Xia et al.'s work to the SIMP [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]based topology optimization framework and provided a comprehensive comparison of the BESO and SIMP methods for composite design to brittle fracture resistance. Russ and Waisman [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 quasi-brittle material[END_REF] proposed two different topology optimization frameworks for the brittle fracture resistance involv-ing one-phase structure material, combining phase field and TO. In their work, the optimization problem was defined as minimizing the total volume or volume and fracture energy while ensuring fracture energy or energy dissipation constraints. Wu et al. [START_REF] Wu | Level-set topology optimization for maximizing fracture resistance of brittle materials using phase-field fracture model[END_REF] developed a level-set method [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] based topology optimization for the brittle fracture resistance of two-phase composite materials.

Taking into account dynamics in the fracture process is of extreme importance for resistance of structures to impacts. When dynamics are involved, the cracks can interact with wave propagation and the final crack patterns depend on the energy of the impactor. In addition, dynamics can lead to more complex crack configurations such as crack branching or initiation of cracks within the solid. In [START_REF] Miller | Energy dissipation in dynamic fracture of brittle materials[END_REF], Miller et al. analyzed the relationship between energy dissipation and crack paths instabilities in dynamic fracture of brittle materials.

The phase field method has been extended to dynamic problems in Borden et al. [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF], Hofacker and Miehe [START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF], and many extensions and applications have been proposed (see e.g. [START_REF] Bleyer | Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms[END_REF][START_REF] Nguyen | Modeling dynamic fracture of solids with a phase-field regularized cohesive zone model[END_REF][START_REF] Ren | An explicit phase field method for brittle dynamic fracture[END_REF][START_REF] Geelen | A phase-field formulation for dynamic cohesive fracture[END_REF][START_REF] Tian | A hybrid adaptive finite element phase-field method for quasi-static and dynamic brittle fracture[END_REF][START_REF] Raghu | Modeling of brittle fracture in thick plates subjected to transient dynamic loads using a hybrid phase field model[END_REF]).

In the present work, TO analysis for dynamic fracture resistance of structures is investigated.

As compared to available existing works, the present framework involves the folllowing new contributions: (a) dynamic phase field simulations to fracture is combined with TO for minimizing the fracture energy of a structure; (b) semi-analytical sensitivities required in the TO analysis in a dynamic fracture analysis context are provided; (c) comparisons of the obtained fracture resistance (reduction in the fracture energy) as compared to a static TO designs are discussed.

The organization of the paper is as follows. In section 2, the dynamic fracture model involved in the analysis is reviewed. In section 3, the related time-space discretizations in the context of the Finite Element Method (FEM) are provided. In section 4, the topology optimization framework involving the dynamic fracture analysis is developed. The sensitivity analysis related to the fracture energy of the system in a dynamic context as well as the whole SIMP TO algorithm is provided.

Finally, numerical examples are proposed in section 5 to validate the sensitivity analysis, and to evaluate the fracture energy reduction of different structures, and more specifically the added value of the dynamic analysis as compared to a static TO.

Dynamic phase field fracture model

In this section, the dynamic phase field method for crack propagation in quasi-brittle solids is briefly reviewed. A structure defined in a domain Ω ⊂ R D is considered, with D the space dimension, with external boundary ∂Ω ⊂ R D-1 . In the context of the phase field method, as shown in Fig. 1, the crack surfaces collectively denoted by Γ are described by a continuous damage field d ∈ [0, 1], which takes 0 value when the material is undamaged and 1 when the material is cracked.

The portions of ∂Ω, ∂Ω u and ∂Ω t denote the Dirichlet and Neumann boundaries, respectively (see Fig. 1). We define u, u = du dt and ü = d 2 u dt 2 as the displacement, velocity and acceleration vectors, respectively.

In this context, the elastic strain energy E s is defined by

E s (u, d) = Ω ψ e (ε (u), d) dΩ (1)
where ψ e is a strain density function, whose form will be specified later, and ε = 1 2 ∇u + ∇ T u is the linearized second-order strain tensor, with ∇(•) the gradient operator. The kinetic energy of the solid is defined by:

E k ( u) = Ω 1 2 ρ u • u dΩ, ( 2
)
where ρ is the material density. In the phase field method, a non-local fracture energy is defined according to

E f (d) = Ω c 1 G c ω(d) + 2 ∇d • ∇d dΩ (3)
where c 1 is a constant, G c is the Griffith-type critical energy release rate, ω(d) is a local damage density function, and is a length regularization parameter, which defines the width of the regularized crack. Finally the work of external forces is defined by

W ext (u) = ∂Ωt t • u dS + Ω f • u dΩ (4)
where t denotes prescribed traction over the portion of the boundary ∂Ω t (see Fig. 1), and f denotes body forces. The action-integral over the time interval [t 1 , t 2 ] is defined by:

A = t 2 t 1 E s (u, d) + E f (d) -E k ( u) -W ext (u) dt. (5)
In the dynamic context, the variational principle of nonlocal damage at the core of the phase field method implies minimization of the action-integral under the constraint of irreversibility of the damage field, i.e.

ḋ ≥ 0. ( 6
)
where ḋ = d(d)/dt denotes the rate of the damage field.

In the following, the different equations of the model in the case of an assumed isotropic quasibrittle solid are specified. We follow Miehe et al. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] and express the strain density function such that damage is induced by traction only as:

ψ e = (1 -d) 2 + b ψ + e + ψ - e , (7)
where ψ + e and ψ - e denote the positive and negative components of the strain density function, respectively, which can be computed from the strain tensor as

ψ ± e = λ 2 Tr [ε] 2 ± + µ ε ± : ε ± , (8)
where λ and µ are the Lamé coefficients, which can be related to the Young's modulus E and Poisson's ratio ν by

λ = Eν (1 + ν) (1 -2ν) and µ = E 2 (1 + ν) . (9)
Note that other decompositions exist (see a comparison and discussion e.g. in [START_REF] Nguyen | Implementation of a new strain split to model unilateral contact within the phase field method[END_REF]). Above,

• ± can be expressed by a ± = 1 2 (a ± |a|) and Tr [•] denotes the trace operator. The positive and negative parts ε ± can be expressed by:

ε ± = D p=1 ε p ± Q p , Q p = v p ⊗ v p , (10)
where D denotes the space dimension and ε p and v p are the eigenvalues and eigenvectors of ε, respectively. The following definitions are used: ω(d) = d 2 and c 1 = 1 2 [START_REF] Bourdin | The variational approach to fracture[END_REF]. Other choices are possible, e.g. using ω(d) = d and c 1 = 3 8 [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]. The first choice induces damage for any loading (even though very low at the beginning) while the second choice leads to a linear elastic stage before damage. Due to its simplicity, the first choice is adopted here. A more in-depth comparison of the different available models and applications in a dynamic context can be found in [START_REF] Mandal | Evaluation of variational phase-field models for dynamic brittle fracture[END_REF].

With these models at hand, stationary variation of (5) leads to the following Euler-Lagrange equations:

     ∇ • σ + f = ρü, G c d -2 ∆d = 2 (1 -d) ψ + e , (11)
where ∇ • (•) and ∆ (•) denote the divergence and Laplacian operators, respectively, and where σ is the Cauchy stress tensor σ = ∂ψe ∂ε , which is expressed under the above assumptions by:

σ = (1 -d) 2 + b ∂ψ + e ∂ε + ∂ψ - e ∂ε = (1 -d) 2 + b λ Tr [ε] + 1 + 2µε + + λ Tr [ε] -1 + 2µε -, (12)
where 1 is the second-order identity tensor and b << 1 a small numerical parameter used to maintain stability in the case of fully broken elements. Above, ε ± and Tr (ε) ± can be related to ε through the following operators:

ε ± = P ± : ε (13) Tr [ε] ± = R ± Tr [ε] (14)
in which the components of P ± are given in closed form as [START_REF] Miehe | Algorithms for computation of stresses and elasticity moduli in terms of Seth-Hill's family of generalized strain tensors[END_REF]:

P ± ijkl := ∂ε ± ∂ε = D p H (±ε p ) (Q p ) ij (Q p ) kl + 1 2 D p D q =p φ pq (Q p ) ik (Q q ) jl + (Q p ) il (Q q ) jk (15)
with

φ pq =      ε p ± -ε q ± ε p -ε q , if ε p = ε q H (±ε p ) , if ε p = ε q (16)
where H (•) denotes the Heaviside step function. The operator R ± is expressed by

R ± = 1 2 (sign (±Tr [ε]) + 1) . (17)
To achieve non-reversible evolution of cracks, a strain history functional introduced by Miehe et al. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF][START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. part i. balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids[END_REF]] is adopted to replace ψ + e in Eq. ( 11), as

H (x, t) = max s∈(0,t] ψ + e (x, s) . (18)
Then, the equations ( 11) are substituted by

     ∇ • σ + f = ρü, G c d -2 ∆d = 2 (1 -d) H, (19)
The above equations are completed with boundary conditions as

           u = ū on ∂Ω u , σ • n = t on ∂Ω t , ∇d • n = 0 on ∂Ω. ( 20
)
where n denotes the outward unitary normal vector to ∂Ω (see Fig. 1).

Multiplying Eqs. ( 19) by two different test functions, δu for the displacement problem and δd for the phase field problem, respectively, integrating the resulting expression over the domain Ω, and using the divergence theorem together with boundary conditions yields the associated weak forms:

find d ∈ H 1 (Ω) and find u ∈ D = {v|v = u on ∂Ω u , v ∈ H 1 (Ω)} such that: Ω 2H(u) + G c d δd + G c ∇d • ∇d dΩ = Ω 2H(u)δd dΩ (21)
Ω {ρü • δu + σ(u, d) : ε (δu)} dΩ - ∂Ωt t • δu dS - Ω f • δu dΩ = 0 (22) for all δd ∈ H 1 (Ω) and δu ∈ H 1 0 (Ω) = {v|v = 0 on ∂Ω u , v ∈ H 1 (Ω)}, H 1 is the usual Sobolev
space of square-integrable derivative functions.

Numerical solving procedure

Displacement field u, phase field d and their gradients ε (u) and ∇d are approximated by classical FEM interpolation in the elements of the FEM mesh according to:

u(x) = N u (x)u e , ε (x) = B u (x)u e , d(x) = N d (x)d e , ∇d(x) = B d (x)d e , (23)
where u e and d e denote the nodal displacement and phase field in an element e, N and B denote the matrices of the shape functions and their derivatives, respectively. The indices (u and d) of N and B refer to displacement and phase field variables, respectively.

A staggered scheme is adopted. At one time step t n , the phase field problem ( 21) is solved, assuming the displacement field u given. Then, the mechanical problem ( 22) is solved assuming the phase field d given. These problems are solved alternatively before solving the problems at the next time step. Note that the mechanical problem ( 22) is nonlinear due to the separated description of the strain field in [START_REF] Zhu | Topology optimization in aircraft and aerospace structures design[END_REF]. Here, we transform this problem into a linear one by expressing the projectors P ± with with respect to the displacements know from the previous time step n -1, i.e.

P ± (ε n ) P ± (ε n-1 ), R ± (ε n ) R ± (ε n-1 ).
At time t n , the strain history functional described in Eq. ( 18) can be calculated using

H n =    ψ + e n if ψ + e n -H n-1 > 0, H n-1 otherwise, (24)
Note that H n is discontinuous, which brings difficulties to the subsequent sensitivity derivations presented in section 4. To alleviate this issue, we introduce a continuous version of the history function as:

Hn H n-1 + ψ + e n -H n-1 g ψ + e n -H n-1 ( 25
)
where g is a regularized Heaviside function, defined by

g(x) = 1 2 1 + 2 π arctan x ζ , (26)
and ζ is regularization parameter. When ζ decreases, the approximation is closer to a sharp jump (see Fig. 2). In this paper, ζ = 10 -6 is adopted.

More specifically, expressing the strain and stress tensors in vector forms in 2D , i.e.

[ε] = [ε 11 , ε 22 , 2ε 12 ], [σ] = [σ 11 , σ 22 , σ 12 ]
, the constitutive law [START_REF] Challis | Fracture resistance via topology optimization[END_REF] can be expressed at time t n as:

[σ n ] = (1 -d n ) 2 + b λR + n-1 ([ε n ] • [1]) [1] + 2µP + n-1 [ε n ] +λR - n-1 ([ε n ] • [1]) [1] + 2µP - n-1 [ε n ] , (27)
where

R ± n-1 = R ± (ε n-1
) and P ± n-1 = P ± (ε n-1 ), and P ± are the matrix forms associated with the fourth-order tensors P ± .

Introducing ( 23) and ( 27) in ( 21)-( 22) we obtain a linear system of equations in the form

K n d d n = F n d , (phase field problem), (28)
M ün + K n u u n = F n u , (displacement problem) (29)
with

M = Ω ρN u T N u dΩ and F n u = ∂Ωt N u tn dS. ( 30
)
K n u = Ω B T u (1 -d n ) 2 + b λR + n-1 [1] T [1] + 2µP + n-1 B u dΩ, + Ω B T u λR - n [1] T [1] + 2µP - n B u dΩ, (31)
K n d = Ω {( G c + 2H n )N T d N d + G c B T d B d }dΩ (32)
and

F n d = Ω 2N T d H n dΩ. ( 33
)
Then, a time-stepping I = [t 0 , t 1 , t 2 , ..., t M ] is introduced, where t n -t n-1 = ∆t is a time step, assumed to be constant. An unconditionally stable implicit Newmark scheme is used to solve [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF] according to:

un = un-1 + ∆t 2 ün-1 + ün , (34)
with

u n = u n-1 + ∆t un-1 + ∆t 2 4 ün-1 + ün , (35) ün
= 4 ∆t 2 (u n -ûn ) (36) ûn = u n-1 + ∆t un-1 + ∆t 2 4 ün-1 . (37)
Introducing ( 36) into (29), we finally obtain at one iteration k and at one time step t n the linear problem to be solved:

Ku u n = F n u , (38)
with

Ku = 4 ∆t 2 M + K n u (39)
F n u = F n u + M 4 ∆t 2 u n-1 + ∆t un-1 + ∆t 2 2 ün-1 . ( 40
)
Note that in the present work iterations are used within the staggered scheme, i.e the mechanical and phase field problems are solved alternatively during one time step t n until a convergence criterion is reached. In this work, the convergence criterion is reached when the maximum value of the nodal phase field variation between two iterations is lower than a threshold. This algorithm allows larger time steps in the present implicit Newmark's scheme. In addition, such staggered scheme is a key ingredient to simplify the sensitivity analysis developed in section 4.3.

The general algorithm is summarized in Algorithm 1, where quantities at one time step t n and at one iteration k are denoted by (.) n k .

Topology optimization formulations

Material interpolation scheme

In the present paper, the SIMP topology optimization method (Solid Isotropic Material with Penalization (SIMP) method [START_REF] Bendsøe | Material interpolation schemes in topology optimization[END_REF][START_REF] Bendsøe | Optimal shape design as a material distribution problem[END_REF]) is adopted. This method belongs to the so-called densitybased topology optimization (TO) methods (see a review and classification of TO in [START_REF] Sigmund | Topology optimization approaches[END_REF]). In Algorithm 1: Staggered dynamic Phase Field algorithm solving procedure.

Initialize: u 0 , u0 , ü0 , d 0 , H 0 = 0;

for n = 1 : M (Loop over all time steps) do Input:

u n-1 Initialize u n 0 = u n-1 k = k + 1 while Err > tol (Convergence loop) do k = k + 1; Input: u n k-1 Solve K d (u n k-1 )d n k = F d (u n k-1 ) Output: d n k Input: d n k , u n k-1 Solve K u (u n k-1 , d n k )u n k = F u (u n k-1 ) Output: u n k Compute Err = max j | [d n k ] j -d n k-1 j | u n k-1 = u n k end u n = u n k end
this framework, the geometry of the structural domain is defined by a pseudo-density parameter ϕ ∈ [0, 1]. Then, a continuous description of the material properties is defined according to:

E (ϕ) = [E min + (1 -E min ) ϕ p E ] E 0 , ρ (ϕ) = [ρ min + (1 -ρ min ) ϕ pρ ] ρ 0 , G c (ϕ) = [G c,min + (1 -G c,min ) ϕ p G ] G c,0 , (41)
where E 0 , ρ 0 and G c,0 denote the material properties of the solid for ϕ = 1, and E min , ρ min and are defined as E min = 10 -6 , ρ min = 10 -6 and G c,min = 10 -2 , respectively. Note that alternative approaches, e.g. the BESO method [START_REF] Huang | Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method[END_REF] could be used (see [START_REF] Li | A simp-phase field topology optimization framework to maximize quasi-brittle fracture resistance of 2d and 3d composites[END_REF] for a comparison between SIMP and BESO in the context of fracture resistance maximization).

G c,

Optimization problem

In this section, we define the topology optimization problem related to minimizing the dynamic fracture of a structure. Following previous works on fracture resistance maximization using topology optimization [START_REF] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[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], the phase field method described in section 3 is used to describe the full fracture process of the structure, from initiation until full failure, while here the dynamic effects are taken into account. The problem can be formulated as follows:

min ϕ∈[0,1] : G f (ϕ), s. t. : K n d d n = F n d , n = 1, 2, ..., M M ün + K n u u n = F n u , n = 1, 2, ..., M C -C 0 ≤ 0, V -χV 0 ≤ 0, (42)
where G f (ϕ) denotes the fracture energy, and C = F s • u s denotes a structural static compliance.

Here, the compliance constraint is only used to ensure connectivity of the material within the structure. Without this constraint, unrealistic topologies with disconnected parts could be obtained in this dynamic context, as cracks may occur in the middle of the structure. The compliance C is evaluated by a separated static test with an external force F s = -100 N. Note that such compliance constraint has been used by several other authors in a dynamic topology optimization context, even though in a linear vibration regime (see e.g. [START_REF] Silva | A critical analysis of using the dynamic compliance as objective function in topology optimization of one-material structures considering steady-state forced vibration problems[END_REF][START_REF] Zhou | A normalization strategy for BESObased structural optimization and its application to frequency response suppression[END_REF]). Above, u s is the static displacement response of the undamaged structure under a test static external load F s and C 0 is the prescribed compliance upper bound; V is the target structural volume (area in 2D), which can be computed as V = Ne e=1 V e ϕ e , and V 0 = Ne e=1 V e is the total volume of the design domain, and χ denotes the target volume fraction constraint. We define the fracture energy over the whole loading history as:

G f = T 0 Ω 1 d • ḋ + ∇d • ∇ ḋ G c dΩ dt. ( 43
)
Using the trapezoidal rule, G f can be numerically approximated as

G f = M n=1 1 2 (d n ) T K g ḋn + d n-1 T K g ḋn-1 ∆t, ( 44
)
where K g is defined by

K g = Ω G c 1 (N d ) T N d + (B d ) T B d dΩ. (45)
Above, K g is a matrix which depends neither on damage nor on displacement fields.

Sensitivity of fracture energy

In this section, the sensitivity of fracture energy is derived in a dynamic context. The derivations require the use of the adjoint method [START_REF] Komkov | Design sensitivity analysis of structural systems[END_REF][START_REF] Cho | Design sensitivity analysis and topology optimization of displacement-loaded non-linear structures[END_REF][START_REF] Buhl | Stiffness design of geometrically nonlinear structures using topology optimization[END_REF]. For the widely known derivations related to the compliance and volume fraction, one may refer to [START_REF] Bendsoe | Topology optimization: theory, methods, and applications[END_REF][START_REF] Andreassen | Efficient topology optimization in MATLAB using 88 lines of code[END_REF] for details.

The sensitivity of G f with respect to a change in the pseudo-density is given by

∂G f ∂ϕ e = M n=1 1 2 ∂ ∂ϕ e (d n ) T K g ḋn + d n-1 T K g ḋn-1 ∆t, (46)
and involves evaluating ∂d n ∂ϕ e and ∂ ḋn ∂ϕ e . Using the chain rule, we have

∂ ḋn ∂ϕ e = ∂ ḋn ∂d n ∂d n ∂ϕ e (47)
where ∂ ḋn ∂d n can be obtained by the Newmark scheme (34)-(36)(37) as the simple expression:

∂ ḋn ∂d n = 2 ∆t . (48)
To express ∂d n ∂ϕe the adjoint method [START_REF] Buhl | Stiffness design of geometrically nonlinear structures using topology optimization[END_REF] is employed. Introducing two vectors of Lagrange multipliers (adjoint vectors) λ n and λ n-1 , and assuming that the problems

R n d = K n d d n -F n d = 0 (49) R n-1 d = K n d d n-1 -F n-1 d = 0 (50)
have been solved, then the terms (λ n ) T R n d and (λ n-1 ) T R n-1 d can be added to the objective function without change as:

G f = Ns n=1 1 2 (d n ) T K g ḋn + d n-1 T K g ḋn-1 ∆t + (λ n ) T R n d + λ n-1 T R n-1 d . (51)
In addition, using ( 8), (ψ + e ) n can be expressed using the discrete (vector) forms of strain tensor as:

ψ + e n = 1 2 ε n λR + n-1 [1] T [1] + 2µP + n-1 ε n . (52)
In [START_REF] Nguyen | Implementation of a new strain split to model unilateral contact within the phase field method[END_REF], it is worth noting that ε obviously depends on ϕ e . However, for the sake of simplicity, we assume that the term involving ∂ε ∂ϕe has small influence as compared to the other terms and neglect it. Then, the following approximation is made:

∂ (ψ + e ) n ∂ϕ e 1 2 ε n ∂λ ∂ϕ e R + n-1 [1] T [1] + 2 ∂µ ∂ϕ e P + n-1 ε n . ( 53
)
Taking the derivation of G f with respect to the pseudo-density, using (48) and combining similar terms, the following expression is obtained, after some calculations:

∂ Ĝf ∂ϕ e = Ns n=1 1 2 (d n ) T ∂K g ∂ϕ e ḋn + d n-1 T ∂K g ∂ϕ e ḋn-1 ∆t + (λ n ) T ∂K n d ∂ϕ e d n - ∂F n d ∂ϕ e + λ n-1 T ∂K n-1 d ∂ϕ e d n-1 - ∂F n-1 d ∂ϕ e + 1 2
ḋn T K g + (d n ) T K g ∂ ḋn ∂d n ∆t + (λ n ) T K n d ∂d n ∂ϕ e + 1 2 ḋn-1 T K g + d n-1 T K g ∂ ḋn-1 ∂d n-1 ∆t + λ n-1 T K n-1 d ∂d n-1 ∂ϕ e . (54)
The terms ∂d n ∂ϕe and ∂d n-1 ∂ϕe are difficult to evaluate in practice. However, as R n d = 0 and R n-1 d = 0, the vectors λ n and λ n-1 can be chosen arbitrarily. They are then chosen to eliminate the unknown terms ∂d n ∂ϕe and ∂d n-1 ∂ϕe such that:

1 2 ḋn T K g + (d n ) T K g ∂ ḋn ∂d n ∆t + (λ n ) T K n d ∂d n ∂ϕ e = 0 (55)
and

1 2 ḋn-1 T K g + d n-1 T K g ∂ ḋn-1 ∂d n-1 ∆t + λ n-1 T K n-1 d ∂d n-1 ∂ϕ e = 0. (56)
Eqs. ( 55) and ( 56) are equal to zero if the expressions under brackets on the left-hand are equal to zero, corresponding to the following systems of equations:

2K n d λ n = -K g ḋn + K g ∂ ḋn ∂d n d n ∆t, (57)
and

2K n-1 d λ n-1 = -K g ḋn-1 + K g ∂ ḋn-1 ∂d n-1 d n-1 ∆t. ( 58
)
Solving Eqs. ( 57) and ( 58), the Lagrange multipliers λ n and λ n-1 are then available. Above, the expressions of ḋn can be computed according to [START_REF] Da | Fracture resistance design through biomimicry and topology optimization[END_REF]. The sensitivity of the fracture energy is then finally obtained as:

∂G f ∂ϕ e = Ns n=1 1 2 (d n ) T ∂K g ∂ϕ e ḋn + d n-1 T ∂K g ∂ϕ e ḋn-1 ∆t +(λ n ) T ∂K n d ∂ϕ e d n - ∂F n d ∂ϕ e + λ n-1 T ∂K n-1 d ∂ϕ e d n-1 - ∂F n-1 d ∂ϕ e (59)
in which ∂Kg ∂ϕe can be derived from Eq. ( 45), as:

∂K g ∂ϕ e = Ω ∂G c ∂ϕ e 1 (N d ) T N d + (B d ) T B d dΩ. (60)
The terms

∂K n d ∂ϕ e = Ω 2 ∂H n ∂ϕ e + ∂G c ∂ϕ e (N d ) T N d + ∂G c ∂ϕ e (B d ) T B d dΩ, ∂F n d ∂ϕ e = Ω 2 ∂H n ∂ϕ e N d dΩ, (61)
where ∂Gc ∂ϕe can be obtained from the material interpolation scheme [START_REF] Wu | Level-set topology optimization for maximizing fracture resistance of brittle materials using phase-field fracture model[END_REF], and ∂H n ∂ϕe can be derived from Eq. ( 25):

∂H n ∂ϕ e = ∂ψ + e ∂ϕ e g ψ + e -H n-1 + ψ + e -H n-1 ∂g (ψ + e -H n-1 ) ∂ϕ e (62)
with

∂g ((ψ + e ) n -H n-1
)

∂ϕ e = ζ ∂(ψ + e ) n ∂ϕe π ζ 2 + ((ψ + e ) n -H n-1 ) 2 , (63)
and the term

∂(ψ + e ) n ∂ϕe
have been approximated by Eq. ( 53). Note that above adjoint vectors λ n and λ n-1 are here path-independent, in contrast to other formulations, see e.g. [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 quasi-brittle material[END_REF]. This strong assumption has the advantage to gratefully simplify the formulation and the implementation. The influence of such simplification on the accuracy of the sensitivities will be tested in the numerical examples. To improve stability, mesh independence and to eliminate so called checkerboard issues [START_REF] Sigmund | Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima[END_REF],

Optimization techniques

Initialize design variables and filter

filtering techniques are often used in topology optimization. Following [START_REF] Schevenels | Robust topology optimization accounting for spatially varying manufacturing errors[END_REF][START_REF] Andreassen | Efficient topology optimization in MATLAB using 88 lines of code[END_REF], a filtered density variable θ e is introduced as θ e = i∈Ne e V e ϑ e i∈Ne e V e [START_REF] Andreassen | Efficient topology optimization in MATLAB using 88 lines of code[END_REF] where ϑ denotes the design variable and N e denotes the set of elements whose center-to-center distance r ei to the e-th element is lower than the filter radius r min . The corresponding weighting factor ei is defined by ei = max (0, r min -r ei ).

The projection technique proposed by Guest et al. [START_REF] Guest | Achieving minimum length scale in topology optimization using nodal design variables and projection functions[END_REF] is then adopted to minimize transition regions with pseudo-density values ϕ between zero and one, as where ϕ e is the elemental pseudo-density, and η is a parameter defined by 1 in the first iteration and is doubled after every specified time steps until it reaches a chosen maximum value, taken here as 128 by numerical tests.

ϕ e = 1 -e -ηθe + θ e e -η (65)
Using this procedure, the sensitivities of the objective functions and optimization constraints with respect to the design variable can be further derived by means of the chain rule as

∂f (ϕ e ) ∂ϑ e = ∂f ∂ϕ e ∂ϕ e ∂θ e ∂θ e ∂ϑ e . (66)
The method of moving asymptotes (MMA) proposed by Svanberg [START_REF] Svanberg | The method of moving asymptotes-a new method for structural optimization[END_REF] is adopted for seeking the optimal distribution of the design variables ϑ. Following [START_REF] Guest | Eliminating beta-continuation from heaviside projection and density filter algorithms[END_REF], the η-based modification on the asymptotes are adopted for removing spurious oscillations after doubling the projection parameter η. The convergence criterion of the topology optimization is determined by the maximal change on the design variable, which should be less than 10 -3 . Fig. (3) summarizes the flow chart of the proposed topology optimization.

Numerical examples

In this section, two typical 2D structures are considered. These have been widely studied in the topology optimization community. The first one is a 3-point bending beam, also called MBB-beam in the literature [START_REF] Olhoff | On cad-integrated structural topology and design optimization[END_REF]. For the sake of computational costs, only the right half of this axisymmetric beam is considered as shown in Fig. 4(a). The left end is simply supported in the x-direction and the lower right-end corner is simply supported in the y-direction. The second structure is a cantilever beam, which is shown in Fig. 4(b). The length and width of these two structures are the same, L = 150 mm and H = 60 mm. A velocity is prescribed on a surface of length L f = 4mm.

Fig. 5 depicts the loading velocity profile, which increases from 0 to v 0 by a time t 0 , and then v t t max t 0 Critical energy release rate G c,0 22.17 N/mm remains constant until the maximum time t max is reached. The material properties adopted here are taken from the Kalthoff-Winkler experiment [START_REF] Kalthoff | Failure mode transition at high rates of shear loading, DGM Informationsgesellschaft mbH[END_REF], and are summarized in Table 1. The same geometry, loading curve and parameters will be kept in all following examples. In the example of section 5.1, the structure is discretized into a coarse 75 × 30 four-node quadrilateral elements mesh for the sake of computational costs. In the examples of sections 5.2 and 5.3, a finer mesh with 150 × 60 four-node quadrilateral finite elements is adopted.

To evaluate the added value of the present framework, two solutions are defined:

1. A so-called "S-design" solution. This solution is obtained by static topology optimization with minimization of compliance under volume constraint with a static force chosen as F s = -100N. Then, the design is remained unchanged during the dynamic fracture simulation. ϕ ∆ϕ e v 0 ∆t t 0 t max tol 0.5 10 -6 4 mm 40m/s 1 µs 2 µs 80 µs 10 -5

The obtained design for the two problems studied in the next examples, namely the half MBB-beam and the cantilever beam are depicted in Fig. 6.

2. A so-called "DF-design" solution. In that case, the problem 42 is solved to define the topology: at each iteration of the algorithm, a static problem is firstly solved to prescribe the compliance constraint, then a full dynamic fracture simulation is performed to evaluate the fracture energy, and compute the sensitivities to update the topology.

Validation of sensitivity analysis

First the sensitivity analysis developed in section 4.3 is validated. Both half MBB-beam and cantilever beam are considered.

The central finite difference method is employed to provide a reference solution to be compared with our semi-analytical sensitivities expressions, according to: df (ϕ) dϕ e ≈ f (ϕ 1 , ..., ϕ e + ∆ϕ e , ..., ϕ Ne ) -f (ϕ 1 , ..., ϕ e -∆ϕ e , ..., ϕ Ne ) 2∆ϕ e

where ∆ϕ e is a pseudo-density perturbation parameter. The value of the numerical parameters are listed on Table 2.

error = ξ dif -ξ ana ξ ana (68)
where ξ denotes the vector of element sensitivity values. The superscripts dif and ana indicate the finite difference method and semi-analytical method, respectively, and ξ ana denotes the maximum element sensitivity obtained by the semi-analytical method. Figs. 7(e) and (f) show a comparison between elemental sensitivities associated with the fracture energy of these two structures. A good agreement between our analytical expressions of sensitivities and the reference finite difference solution is noticed. Fig. 8 depicts the sensitivity validation on a structure with random distribution of densities ϕ e ∈ [0, 1] in the elements. After generating the densities using a uniform probability of distributions, a filter is then applied. The other parameters are provided in Table 2.

We can note that even though the absolute values of sensistivities are good, the relative errors might locally be high, even though localized, associated with the approximation made in Eq. ( 53).

However, these errors remain acceptable. In addition, it will be shown in the next examples that the made approximation allows a large simplification of the whole methodology, while keeping important dynamic fracture reduction results.

Half-MBB beam

In this example, the presented methodology is applied to the Half-MBB beam (see Fig. 3. When the loading rate increases, and thus the related dynamic effects, the reduction of the fracture energy as compared to the one obtained by static analysis is even larger.

Fig. 11 depicts the iterative process plots of the above topology optimizations. Regardless of the jumps caused by the variation of the projection parameter η, a good convergence is appreciated.

All the optimization constraints are verified, except the compliance constraint for v 0 = 60 m/s, which might be too strict to be reached in this case. The competition between minimizing the fracture energy and satisfying the compliance constraint might be one possible reason for the observed oscillations. For the case v 0 = 20m/s, 647 iterations were necessary, for a total of 8.7 h on a single processor for the whole optimization process.

Next, a longer loading period is investigated, with t max = 100 µs. Fig. 12 depicts the topological designs and their final fracture patterns under different loading rates. Compared to the designs for t max = 40 µs shown in Fig. 9, the obtained designs show an obvious difference, and the final fracture patterns also change accordingly. Fig. 13 depicts the iterative processes of these topological designs. In this case, although the objective function remains oscillatory, these oscillations remain small and around a stable value. The computational time for the case v 0 = 20m/s is 23 h for 672 iterations. The computational times are here proportional to the chosen loading period. Fig. 14 compares the fracture energy evolution in a period of time [0-100 µs] obtained by the present DF-designs for t max = 100 µs and the S-designs. Once again, an important decrease of the fracture energy using the DF-design for all loading rates is appreciated. Corresponding comparisons and fracture energy reductions at t = 100 µs are indicated in Table 4.

Cantilever beam

In this section, the cantilever cantilever beam shown in Fig. 4 It is worth noting that there remain some gray elements in the DF-designs. This issue is a classical one found by several other authors in dynamic topology optimization. For example, it is discussed as a key issue in [START_REF] Silva | A critical analysis of using the dynamic compliance as objective function in topology optimization of one-material structures considering steady-state forced vibration problems[END_REF] and found in other works such as in [START_REF] Yoon | Structural topology optimization for frequency response problem using model reduction schemes[END_REF][START_REF] Zhao | Topology optimization for minimizing the maximum dynamic response in the time domain using aggregation functional method[END_REF]. As the main objective of this paper is to present the new topology optimization algorithm with fracture minimization Fig. 17 depicts the iterative processes of the topology optimizations for different loading rates.

A good convergence is obtained and all the constraints are reached. The computational time for the case v 0 = 20m/s is 12.7 h for 624 iterations.

As a final remark, we can note that in most studied examples, the cracks are rather diffuse damage zones. In the present phase field framework, the cracks width depends on the mesh density. To maintain reasonable computational costs, we used meshes which do not allow very fine descriptions of cracks. However, it has been shown in many other studies (see e.g. [START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF]) that the phase field method is fully convergent with respect to the mesh density, even in the dynamic case.

Then, finer crack descriptions can be obtained if faster computational ressources are available.

Figure 1 :

 1 Figure 1: Cracked solid: (a) representation of cracks by surfaces; (b) continuous approximation of cracks by a damage field d(x).

Figure 2 :

 2 Figure 2: Regularized Heaviside function.

  min are artificial lower values to avoid numerical singularity in zero-pseudo density zones. Above, p E , p ρ and p G are penalty parameters. These parameters are chosen so as to avoid intermediate values of ϕ. Here, we choose p E = 3, p ρ = 1 and p G = 1. The values of E min , ρ min and G c,min

Figure 3 :

 3 Figure 3: Optimization flow chart.

Figure 4 :

 4 Figure 4: Geometry and boundary conditions of: (a) Half MBB-beam; (b) Cantilever beam.

Figure 5 :

 5 Figure 5: Prescribed velocity.

Figure 6 :

 6 Figure 6: S-designs obtained for :(a) the half MBB-beam; (b) the cantilever beam.

Figs. 7

 7 Figs.7(a) and (b) depict the fracture plot of the half MBB-beam and cantilever beam at time t max = 80 µs, in which only d > 0.6 is depicted for the sake of clarity. Figs.7(c) and (d)show the normalized error map of the sensitivity values, which is defined by:

Figure 7 :

 7 Figure 7: Validation of fracture energy sensitivity analysis: fracture pattern at t max = 80 µs in (a) the half MBBbeam; (b) the cantilever beam; error of normalized sensitivity values (c) the half MBB-beam; (d) the cantilever beam; sensitivity values in elements at t max = 80 in (e) the half MBB-beam and (f) the cantilever beam.

Figure 8 :

 8 Figure 8: Validation of fracture energy sensitivity analysis on a design-variable-random-distributed structure: fracture pattern at t max = 80 µs in (a) the half MBB-beam; (b) the cantilever beam; error of normalized sensitivity values (c) the half MBB-beam; (d) the cantilever beam; sensitivity values in elements at t max = 80 in (e) the half MBB-beam and (f) the cantilever beam.

4

 4

  (a)) to minimize the fracture energy with respect to the topology of the structure. Three different values of loading rates are investigated, v 0 = 20 m/s, v 0 = 40 m/s and v 0 = 60 m/s, respectively. Two maximal loading times, t max = 40 µs and t max = 100 µs, are separately considered. The compliance

Figure 9 :

 9 Figure 9: Obtained topological designs for t max = 40µs and corresponding final fracture patterns for different loading rates: (a) DF-design: final topology; (b) DF-design: final crack pattern; (c) S-design: final crack patterns.

Figure 10 :

 10 Figure 10: Comparison of time-energy curves for S-and DF-designs for t max = 40 µs at different loading rates: (a) v 0 = 20 m/s; (b) v 0 = 40 m/s; v 0 = 60 m/s

Fig. 10

 10 Fig. 10 compares the fracture energy evolution in time for the present Dynamic Fracture DFand S-optimized designs, in a period of time [0-40µs]. An important decrease of the fracture energy using the DF-design for all loading rates is appreciated, which shows the importance of

Figure 11 :

 11 Figure 11: Iterative topology optimization process for t max = 40 µs under different lading rates: (a) v 0 = 20 m/s; (b) v 0 = 40 m/s; (c) v 0 = 60 m/s.

Figure 12 :

 12 Figure 12: Topological designs for t max = 100 µs and corresponding final fracture patterns for different loading rates: (a) DF-design: final topology; (b) DF-design: final crack pattern; (c) S-design: final crack patterns.

  Fig.15(a) shows the topological designs of the cantilever beam for different loading rates. Here, the proposed method gives different topology designs for different loading rates. Fig.15(b)shows the crack pattern at t max = 60 µs. Fig.15(c) depicts the final crack pattern of the S-design for comparison. It is worth noting that there remain some gray elements in the DF-designs. This issue

Figure 13 :

 13 Figure 13: Iteration process of the topology optimization for t max = 100 µs for different lading rates: (a) v 0 = 20 m/s; (b) v 0 = 40 m/s; (c) v 0 = 60 m/s.

Figure 14 :Figure 15 :

 1415 Figure 14: Comparison of time-energy curves for S-and DF-designs for t max = 100µs at different loading rates: (a) v 0 = 20 m/s; (b) v 0 = 40 m/s; v 0 = 60 m/s.

Figure 16 :

 16 Figure 16: Comparison of time-energy curves for S-and DF-designs at different loading rates: (a) v 0 = 20 m/s; (b) v 0 = 40 m/s; v 0 = 60 m/s.

  Fig. 16 depicts the fracture energy evolution of the DF-and S-designs in a period of time [0-60 µs]. Table 5 provides the comparison of fracture energy for different loading rates of the DF-and S-designs at the final time. Again, the DF-designs show large reductions of the fracture energy. Further investigations, including comparisons with stress-based linear topology optimization, could be conducted in future studies.

Figure 17 :

 17 Figure 17: Iteration process of the topology optimization of cantilever beam for different lading rates: (a) v 0 = 20 m/s; (b) v 0 = 40 m/s; (c) v 0 = 60 m/s.

A

  SIMP topology optimization framework for maximizing the dynamic fracture resistance has been proposed. Several contributions have been introduced. The dynamic phase field method for fracture has been combined with SIMP topology optimization. Then, a topology optimization minimizing the fracture optimization as an objective function under constraints of material volume and verification of local equilibrium equations has been originally proposed. Semi-analytical expressions of sensitivities in this context have been derived, and their accuracy using numerical finite difference approximations has been validated. The algorithm involves solving at each iteration first a static problem to evaluate the compliance and then a full dynamic fracture problem from initiation to crack propagation, during a given period of time, then taking into account the whole loading history. A staggered scheme with convergence iterations has been used to solve the dynamic phase field problem thus authorizing larger time steps. Numerical examples on structural problems subjected to impacts for different loading velocities have been investigated. The examples show that the present dynamic analysis allows reducing the fracture energy as compared to the designs obtained from static classical topology optimization analysis.

Table 1 :

 1 Material parameters.

	Material properties	Symbol	Value	Unit
	Young's modulus	E 0	1.9 × 10 5	MPa
	Poisson's ratio	ν	0.3	-
	Density	ρ 0	8 × 10 -9	ton/mm 3

Table 2 :

 2 Numerical parameters for validation of sensitivity analysis.

Table 3 :

 3 Comparison of fracture energy using DF-design for t max = 40 µs and S-design at final simulation time for different loading rates.

	Loading rate (m/s) S-design DF-design Fracture energy reduction
	20	258.23	91.51	64.5%
	40	2213.33	347.92	84.3%
	60	4133.94	627.92	84.8%

Table 4 :

 4 Comparison on the fracture energy for DF-designs for t max = 100 µs and S-design at the final time for different loading rates.

	Loading rate (m/s) S-design DF-design Fracture energy reduction
	20	2312.37	313.97	86.4%
	40	5396.64	1007.77	81.2%
	60	8127.18	2530.02	68.9%
	including the dynamics in the topology optimization analysis as compared to the designs obtained
	by simple static analysis. Corresponding comparisons and fracture energy reduction at t = 40 µs
	are presented in Table

Table 5 :

 5 Comparison on the fracture energy of DF-and S-design for the cantilever beam at the final time for different loading rates.

	Loading rate (m/s) S-design DF-design Fracture energy reduction
	20	638.56	227.20	64.4%
	40	1289.30	546.75	57.6%
	60	1788.05	578.54	67.6%

objective, fully addressing this problem is reported to later studies.

Acknowledgement

We acknowledge Prof. Svanberg from KTH, Sweden for providing the code of the MMA algorithm. Yi Wu acknowledges the support from China Scholarship Council (CSC No. 201906130024).