Structural Topology Optimization typically features continuum-based descriptions of the investi- gated systems. In Part 1 we have proposed a Topology Optimization method for discrete systems and tested it on quasi-static 2D problems of energy minimization, assuming linear elastic material. However, discrete de- scriptions become particularly convenient in the fail- ure and post-failure regimes, where discontinuous pro- cesses take place, such as fracture, fragmentation, and collapse. Here we take a rst step towards failure prob- Enrico Masoero School of Engineering, Newcastle University, Newcastle upon Tyne NE1 7RU, United Kingdom; E-mail: enrico.masoero@newcastle.ac.uk Connor O'Shaughnessy School of Engineering, Newcastle University, Newcastle upon Tyne NE1 7RU, United Kingdom; Peter D. Gosling School of Engineering, Newcastle University, Newcastle upon Tyne NE1 7RU, United Kingdom; Bernardino M. Chiaia Department of Structural, Building and Geotechnical Engi- neering - DISEG, Politecnico di Torino, 10129, Turin, Italylems, testing Discrete Element Topology Optimization for systems with nonlinear material responses. The in- corporation of material nonlinearity does not require any change to the optimisation method, only using ap- propriately rich interaction potentials between the dis- crete elements. Three simple problems are analysed, to show how various combinations of material nonlinear- ity in tension and compression can impact the optimum geometries. We also quantify the strength loss when a structure is optimized assuming a certain material be- havior, but then the material behaves di erently in the actual structure. For the systems considered here, as- suming weakest material during optimization produces the most robust structures against incorrect assump- tions on material behavior. Such incorrect assumptions, instead, are shown to have minor impact on the service- ability of the optimized structures. Keywords Topology Optimization Discrete Element MethodMaterial Nonlinearity 2 Enrico Masoeroet al. 1 Introduction 1 Structural Topology Optimisation (TO) is a family of 2 methods to distribute mass within a design domain 3 and maximize the utilization of material under a set 4 of imposed loads and constraints [1, 2]. First ideas of 5 structural of TO date back to the early 20thcentury, 6 with the analytical work of Michell [3]. The advent of 7 modern computers rapidly expanded the capabilities 8 and scope of TO; homogenisation-based TO came rst 9 [4, 5, 6] and was then mostly replaced by the Solid 10 Isotropic Material with Penalisation (SIMP) approach 11 [4, 7]. A crucial part of any TO algorithm is the calcu- 12 lation of objective functions depending on the distribu- 13 tion of material during the optimisation process. This 14 step requires structural analyses, for which continuum- 15 based methods such as the the Finite Element Method 16 (FEM) are the norm. In Part 1, we have coupled SIMP- 17 based TO with Discrete Element (DE) analyses NOTE 18 TO EDITOR: please cite or otherwise link Part 1 here. 19 The resulting Discrete Element Topology Optimisation 20 (DETO) o ers a pathway to consider processes that 21 challenge continuum-based descriptions, such as granu- 22 lar behaviors [8], fracture [9, 10], fragmentation [11, 12], 23 and collapse [13, 14, 15]. DETO may therefore impact 24 scienti c communities whose favor for discrete analyses 25 has been precluding access to topology optimisation, 26 for example in soil mechanics or in nanoscale materials 27 modelling. 28 An advantage of DETO is that, in principle, it could 29 be used to target performance indicators that involve 30 discontinuous structural behaviors, such as resistance to 31fracture or collapse. In such near-failure or even post- 32 failure regimes, the stress-strain behavior of many en- 33 gineering materials is nonlinear. Material nonlinearity 34 has been included in continuum-based TO since the mid 35 1980's [16, 17], mostly focusing on elastoplasticity: see 36 review in Ref. [18]. From a structural design perspec- 37 tive, an important nding has been that structures opti- 38 mized for serviceability performance, vizassuming lin- 39 ear elasticity, may be signi cantly sub-optimal towards 40 failure, when the material behaves nonlinearly [19, 20]. 41 For discrete systems, however, DETO has only been 42 applied to linear elasticity thus far; geometric nonlin- 43 earity from large displacements is included too in Part 1 44 NOTE TO EDITOR: please cite or otherwise link Part 45 1 here as it is naturally captured in DE analyses, where 46 particle interactions are always computed in there de- 47 formed state. However, material nonlinearity in DETO 48 is still to be addressed. 49 This manuscript presents a rst application of DETO 50 to problems with material nonlinearity. Section 2 de- 51 scribes the DETO method for a general problem of 52 interaction energy minimization under imposed forces, 53 and then particularises it for the 2D structures in this 54 manuscript. We show how nonlinear material behaviors 55 can be included seamlessly, without any change to the 56 methodological framework. We then discuss how, un- 57 der imposed displacements instead of forces, the same 58 optimisation scheme maximises the structural ductil- 59 ity, which is a typical objective function when consid- 60 ering material nonlinearity. The last part of Section 2 61 presents four interaction potentials to be used through- 62 out the manuscript, describing linear elastic materi- 63 Topology Optimization using the Discrete Element Method. Part 2: Material Nonlinearity 3 als as well as nonlinear materials undergoing strain- 64 hardening (i.e.becoming less sti ) or strain-sti ening 65 (i.e. becoming sti er) under tension, compression, or 66 both. Section 3 considers three optimization problems 67 including material nonlinearty. Two of these problems 68 are known in the literature on continuum-based TO 69 [19], whereas the third one is proposed here to further 70 appreciate the e ect of material nonlinearity. Particu- 71 lar attention is paid to the scenario in which a structure 72 is optimized assuming a certain material behavior, but 73 then the material turns out behaving di erently from 74 what was planned. Our results con rm the applicabil- 75 ity of DETO to systems with nonlinear material behav- 76 iors, and prompt a discussion on service performance, 77 strength loss, and robustness against incorrect material 78 assumption, all of which are important for structural 79 design. 80 2 Methodology 81 The DETO method has been presented in Part 1 and 82 compared with classical SIMP-based TO using Finite 83 Element analyses NOTE TO EDITOR: please cite or 84 otherwise link Part 1 here. A MATLAB implementa- 85 tion is available on GitHub [21], but the repository also 86 includes a C++ version which is more ecient and rec- 87 ommended for replicating the simulations in the present 88 manuscript. Here we rst present DETO for a problem 89 of interaction energy minimization in a system of parti- 90 cles under imposed forces and interacting viaa generic 91 potential. The problem is then particularised into one 92 of strain energy minimisation in a 2D structure made of 93 closely packed, pairwise-interacting disks. We then dis- 94cuss how the same optimisation algorithm maximises 95 structural ductility when displacements are imposed in- 96 stead of forces. Finally we present several interaction 97 potentials to describe material behaviors that will later 98 be used in Section 3. 99 2.1 DETO: Energy minimization with generic interac- 100 tions 101 Consider a system of Ninteracting particles under a set of imposed external forces and constraints to motion; these latter may represent structural supports such as pins or rollers. Each particle ihas an associated vari- able i2[0;1]. Particles with = 0 interact with zero intensity with the others, e ectively representing voids. Particles with i= 1 interact with full intensity, thus representing full solid. All the per-particle iare gath- ered into a vector . An optimization problem to min- imize the total interaction energy in the system, Utot, implies the possibility to redistribute the values in under a set of constraints, for example: min :c( ) =Utot( ;r; ) (1) subject to :V( ) V0=f (2) : 0< min i1 (3) cis the objective function to minimize, here identi ed 102 withUtot. Eqs. 2 and 3 are the constraints on . Eq. 2 103 xes the target solid fraction fof a nal system where 104 all particles have either i= 0 or i= 1, viza 0-1, 105 void-solid only system. Speci cally, V( ) =PN i=1 i 106 andV0=N. Eq. 3 bounds the values of ibetween 107 1 and a small value min; in principle one could use 108 4 Enrico Masoeroet al. min= 0 to represent void, however small but nite 109 values, such as 103, are frequently preferred to avoid 110 issues with later parts of the optimization process, in 111 particular the ltering step (see Part 1 for more discus- 112 sion on this point NOTE TO EDITOR: please cite or 113 otherwise link Part 1 here). The total interaction en- 114 ergyUtotin Eq. 1 is a generic function of the particle 115 positions rand orientations . Following the concept 116 ofpenalisation , which is central to the SIMP approach, 117 the dependence of Utoton should be set to favour 0-1 118 only structural solutions and to penalise gray solutions 119 with intermediate ivalues. 120 The optimization problem in Eqs. 1-3 is quite generic 121 but some notes on its scope and underlying assumptions 122 are due. Discrete Element analyses typically include 123 velocity-dependent dissipative terms [8]; here we do not 124 consider them because the problem refers to static equi- 125 librium conditions. Extensions to dynamic problems is 126 possible but beyond the scope of this manuscript. Eq. 1 127 also assumes that Utotis history-independent, with no 128 irreversible processes such as bond breakages. Irreversible 129 events can be included in DE analyses and they moti- 130 vate in part the development of DETO. However, min- 131 imization of Utotmay not be a meaningful problem 132 when such irreversibilities are included. One issue that 133 may arise, for example, is that Utotwould be minimum 134 when all the bonds between particles are broken; such 135 a fragmented system is unlikely to be the desired opti- 136 mum. Similar and other issues may arise depending on 137 the functional form of Utot. Therefore the problem in 138 Eqs. 1-3 may not be appropriate, or might require ad- 139 ditional constraints, depending on: the functional form 140ofUtot, the possible presence of irreversible processes, 141 and the type of analyses to be conducted (static or dy- 142 namic).The calculations in this manuscript will target 143 static equilibrium, without irreversibilities, and with 144 expressions of Utotthat ensure the applicability and 145 relevance of the problem in Eqs. 1-3. 146 A slightly more specialized expression of Utot, but still quite generic and more convenient, may be: Utot=X iX j6=iij( i; j)Uij(ri;rj) + +X iX j6=iX k6=i;jijk( i; j; k)Uijk(ri;rj;rk) +::: (4) wherei,j, andkare the indexes of the interacting parti- 147 cles. The rst term describes pairwise interactions, the 148 second three-body interactions, and the series expan- 149 sion may continue to include more multi-body terms. 150 This expansion is commonly used in materials mod- 151 elling, e.g.in atomistic simulations [22]. The spatial 152 dependence of the interactions is described by the Uij 153 andUijkfunctions; to simplify the notation they are 154 assumed to depend only on particle positions r, but 155 the treatment of orientation-dependent terms would be 156 analogous. The penalising functions depend on the 157 values of the particles involved in the interaction. Spe- 158 ci c forms of Uandwill be introduced later in the 159 manuscript; now we continue presenting DETO for the 160 still rather general form of Utotin Eq. 4. 161 To solve the optimization problem in Eqs. 1-3 we use the same approach as detailed in Part 1 NOTE TO EDITOR: please cite or otherwise link Part 1 here. At the generic optimization step, a DE analysis provides Topology Optimization using the Discrete Element Method. Part 2: Material Nonlinearity 5 the static equilibrium con guration corresponding to the current values in under the imposed forces and constraints to motion. The updating scheme for be- tween two subsequent optimization steps is: new i= old i@c @ i (5) =1 2is a parameter to improve convergence. is a parameter that rescales the predicted values of to re- spect the constraints in Eqs. 2 and 3, as well as another typical constraint on the maximum allowed change of ibetween two subsequent steps, viz new i old i (this maximum change is always set to 0.1 in this manuscript). See Part 1 for more discussion on this point NOTE TO EDITOR: please cite or otherwise link Part 1 here.@c @ i is the sensitivity of the cost function with respect to the design variables. Since c=Utothere, the de nition ofUtotin Eq. 4 gives: @c @ i=X iX j6=i@ij( i; j) @ iUij(ri;rj) + +X iX j6=iX k6=i;j@ijk( i; j; k) @ iUijk(ri;rj;rk) + +::: (6) Eq. 6 shows that spatial complexity in Utotdoes not add 162 complexity to the calculation of sensitivity. Indeed, the 163 derivatives in Eq. 6 apply only to the penalising func- 164 tions, whereas the values of UijandUijkcan be di- 165 rectly taken from the results of the DE analysis. Eqs. 5 166 and 6 come from the Optimality Criteria approach to 167 the solution of the optimisation problem, adapted from 168 Ref. [23]. This method is particularly ecient for the 169energy minimization problem here; other methods may 170 better suit other problems. 171 Often the updating scheme in Eq. 5 does not feature directly the sensitivity from Eq. 6, but rather a coarse- grained version of it: d@c @ i=nfX k=1@c @ kWk k infX k=1Wk(7) wherenfis the number of particles within a distance 172 rminfrom the center of particle i, including particle i 173 too.Wk=rminrik0 is a weight function ensuring 174 that particles closer to icontribute most to its coarse- 175 grained sensitivityddc d i. The coarse-graining process in 176 Eq. 7, known as ltering , is commonly used in Finite 177 Element based TO to avoid the checkerboarding prob- 178 lem [24, 25]. Part 1 shows that DETO does not su er 179 from checkerboarding, but also that ltering still im- 180 proves the quality of the optimum solutions NOTE TO 181 EDITOR: please cite or otherwise link Part 1 here. 182 2.2 DETO: Strain energy minimization in 2D 183 The problem in the previous section is now particular- ized for a 2D system that is relevant for all the case studies in this manuscript. The system is a rectangu- lar design domain initially lled with nely layers ofnelx particles each (actually, nelx for odd layers and nelx1 for even layers). The particles are hexagonally closed packed disks of diameter D: see Fig. 1. External forces and constraints to motion may be imposed, as shown for example in the gure. At rst, all disks have same =f. The interaction energy in Eq. 4 is assumed to 6 Enrico Masoeroet al. Fig. 1 Example of a 2D system of interacting disks, subjected in this case to two external forces and constrained by two pinned supports (the black triangles). feature only a pairwise term that depends only on the interparticle distance rij. Forwe use the penalisation scheme that we proposed in Part 1 NOTE TO EDI- TOR: please cite or otherwise link Part 1 here: ij= p i p j (8) This leads to: Utot=X iX j6=i p i p jUij(rij) (9) As discussed in Part 1, also here we take p= 2. A re- 184 striction that we will respect on the spatial function 185 Uijis that the corresponding interaction force Fij= 186 dUij drwill be monotonically increasing with rij. This 187 means that we can only consider material nonlinearity 188 where the sti ness of the material may decrease (strain- 189 hardening) or increase (strain-sti ening) with strain, 190 but always remains positive. Strain-softening nonlin- 191 earity, where the sti ness becomes negative, will not 192 be considered here as it may lead to strain localisation 193 and non-uniqueness of the solutions; this would require 194 changes to the problem formulation in Eqs. 1-3 and re- 195 lated solution method. 1962.3 Ductility maximisation under imposed displacements 197 The solution approach in Eqs. 5 and 6 has been pre- 198 sented for the energy minimisation problem in Eqs. 1-3, 199 under imposed external forces. However, when a ma- 200 terial behaves nonlinearly, imposed forces may cause 201 mechanical instability and compromise the optimisa- 202 tion process. To mitigate this issue, nonlinear structures 203 are often studied under displacement-control, vizim- 204 posing displacements to certain points and computing 205 the corresponding forces ( e.g.in Ref. [19, 20]). Under 206 such conditions, it is desirable for a structure to re- 207 quire high strain energy to reach the imposed displace- 208 ments. Therefore, the minimization problem in Eq. 1 209 should be converted into an energy maximization prob- 210 lem instead. However, below we show that the solution 211 procedure in Eqs. 5 and 6 leads indeed to energy max- 212 imization when displacements are imposed instead of 213 forces. 214 Eqs. 5 and 6 e ectively direct solid ( i.e.increase ) 215 where a structure is most strained, and remove solid 216 (i.e.reduce ) elsewhere. When penalisation schemes 217 such as in Eqs. 8 and 9 are used, larger imply stronger 218 and sti er interactions. Therefore, the solution proce- 219 dure maximizes the structural sti ness for the set of 220 imposed forces and constraints. For this reason, the 221 problem in Eqs. 1-3 is often referred to as `compli- 222 ance minimization' rather than strain energy minimiza- 223 tion. When displacements are imposed instead of forces, 224 Eqs. 5 and 6 still maximize structural sti ness. How- 225 ever, in this case sti ness maximization implies that the 226 Utotto reach the imposed displacements will be max- 227 Topology Optimization using the Discrete Element Method. Part 2: Material Nonlinearity 7 imum; this problem is also known in the literature as 228 ductility maximization [19, 20]. To immediately appre- 229 ciate the impact of sti ness maximization on Utot, con- 230 sider a linear elastic system where forces and displace- 231 ments are collected in vectors Fandu, which are linked 232 by a structural sti ness matrix Kso that F=K u. 233 The strain energy is Utot=1 2uTKu. If forces are im- 234 posed and Kis maximized, then uwill be minimum, 235 and so will Utot. By contrast, if displacements are im- 236 posed and Kis maximized, the required forces Fwill 237 be maximized too, and so will Utot. 238 Running DE analyses with imposed displacements, 239 as part of DETO, requires some dedicated arrange- 240 ments. Consider a generic optimisation step: the cur- 241 rent values in provide the penalisation functions ij242 for each pair of interacting particles. These ijare used 243 in the DE analysis that returns the static equilibrium 244 values of rand thus of Uijunder the imposed displace- 245 ments. These values of Uijare then used to compute the 246 sensitivity in Eq. 6. Under imposed forces, the DE anal- 247 ysis is simply an energy minimization process, carried 248 out using e.g.damped dynamics algorithms such as the 249 Sheppard algorithm [26]. Imposed displacements uimp250 require additional care to control the strain rate during 251 energy minimization. One possibility, adopted in this 252 manuscript, is to start each DE analysis from the under- 253 formed con guration and to progressively increase the 254 imposed displacements with a suitably small rate _uimp. 255 Energy minimization should proceed while increasing 256 the imposed displacements, until the target uimpare 257 reached. After that, energy minimization should con- 258 tinue until satisfying appropriate convergence criteria, 259such asUcurrent totUprevious tot Uprevious tot< etol whereetolis a user- 260 provided small tolerance. In this manuscript, the change 261 in energy is averaged over 10 steps before comparing 262 withetol. 263 2.4 Interaction potentials 264 We consider four types of interaction potentials Uij, 265 each representing a di erent material behavior. Table 1 266 shows the expressions of each potential, along with the 267 corresponding interaction forces (positive when repul- 268 sive). Fig. 2.a compares the Uij(rij) for the various ma- 269 terials and for a set of k0,a. andDparameters that are 270 relevant for this manuscript. As expected from strain 271 energies, all the Uijcurves are zero in the undeformed 272 staterij=D, and positive elsewhere. Fig. 2.b shows 273 theFij(rij) curves, from which the strain-hardening 274 and strain-sti ening regimes can be appreciated. The 275 Fij(rij) curves are proportional to the stress-strain be- 276 havior of the material, which can be quantitatively esti- 277 mated assuming that D= 1 mm, that the box thickness 278 in the third direction is tz=D, and that the contact 279 area between two particles is one sixth of the lateral 280 surface area of the disk, viz1 6Dtz. Under these as- 281 sumptions, the strain between particles in 103units 282 is equal to the elongation in min Fig. 2, whereas the 283 maximum stress between particles in strain-hardening 284 regimes, whenjrijDj 0, is capped tok0 a6 D2= 477 285 MPa (assuming k0= 100 kN/mm and a= 400 mm1286 as in Fig. 2). Simulation results will later con rm this 287 estimation. 288 8 Enrico Masoeroet al. Table 1 Interaction potentials for the case of linear elastic materials ( Lin), symmetric strain-hardening material in tension and in compression ( Weak ), asymmetric material hardening in tension and sti ening in compression ( Weak-T ), and asymmetric material hardening in compression and sti ening in tension ( Weak-C ). Potential name Uij(rij) Fij(rij) =dUij drij Lin1 2k0(rijD)2k0(rijD) Weakk0 a2lnfcosh[a(rijD)]g -k0 aftanh [a(rijD)]g Weak-Tk0 a1 aexp [a(rijD)] + (rijD) k0 a2 -k0 af1exp [a(rijD)]g Weak-Ck0 a1 aexp [a(rijD)](rijD) k0 a2 -k0 afexp [a(rijD)]1g 0246810 -8-6-4-2 02468(a) Energy U ij ( N mm ) Elongation r ij - D ( µm )Lin Weak 0246810 -8-6-4-2 02468(a) Energy U ij ( N mm ) Elongation r ij - D ( µm )Weak-T Weak-C -1-0.5 0 0.5 1 -6 -4 -2 0 2 4 6(b) Force F ij ( kN ) Elongation r ij - D ( µm ) Fig. 2 (a) Interaction potentials from Table 1 for some of the parameters used in this manuscript: k0= 100 kN, D= 1 mm,a= 400 mm1; (b) Corresponding force-elongation curves from Table 1, which are proportional to the stress- strain behaviors of the materials.The interactions in Fig. 2 capture material nonlin- 289 earity under strain. However, the potentials are all elas- 290 tic, with same stress-strain responses upon loading and 291 unloading. Irreversible deformations could be consid- 292 ered in principle, and indeed elastoplasticity and elasto- 293 plastic interactions are within the current capabilities of 294 continuum-based TO and of DE analyses ( e.g.[19, 13]). 295 Such irreversibilities would impact the results if the DE 296 analyses involved dynamic or cyclic loads, or if buck- 297 ling instabilities or strain localization, e.g.due to ma- 298 terial softening or fracture, led to stress relaxation in 299 some parts of the structure. In this manuscript, ma- 300 terial softening and fracture are not considered (only 301 hardening and sti ening as per Fig. 2), buckling will not 302 occur, and imposed loads or displacements will always 303 induce monotonically increasing strain everywhere in 304 the structure. Under these conditions, the reversible in- 305 teractions in Fig. 2 are as representative of large-strain 306 material behaviors as elastoplastic interactions would 307 be. Therefore in this manuscript, for the sake of brevity, 308 we will use the term \plastic ow" even if the mate- 309 rial is not strictly plastic, because its behavior under 310 quasi-statically and monotonically increasing strain is 311 the same. 312 Topology Optimization using the Discrete Element Method. Part 2: Material Nonlinearity 9 2.5 Simulation details 313 All the simulations in this manuscript are based on the 314 parameters in Table 2. The values of D,k0, andain 315 the table are those returning the interactions in Fig. 2. 316 The ltering length rminis used for the weight fac- 317 tors in Eq. 7. The domain thickness tzis only used to 318 estimate the cross section between two interacting par- 319 ticles inDtz=D2, which has been used to relate the 320 force{elongation curves in Fig. 2.b to the stress-strain 321 behavior of the materials. Table 2 Input parameters for the DETO simulations in this manuscript. In some explicitly mentioned cases, ltering will be deactivated and smaller strain rates _ uimpwill be used. D 1 mm dt 0.02s k0 100 kN/mm dmax 0.02 mm a 400 mm1m 1 mg rmin 1.1 mm etol 108 tz 1 mm min 103 _uimp 104mms1 322 The DE-based energy minimization algorithm to nd 323 con gurations at static equilibrium is based on Shep- 324 pard's damped dynamics [26]. The algorithm uses a 325 time stepdt, particle mass m, and caps the maximum 326 particle displacement to dmaxat each Euler integration 327 step; convergence is judged based on the energy tol- 328 eranceetol, as already explained in Section 2.3. The 329 DETO process is considered as converged when new i 330 old i<0:004 for all particles. 331 The stresses per particle ab;i(withab= xx, xy, or yy) can be computed from the interaction forces us- ing the virial approach [27] (see Part 1 for more de- tailsNOTE TO EDITOR: please cite or otherwise linkPart 1 here): ab;i=1 2ViX j6=i(ra;iFij;b;i+ra;jFij;b;j) (10) Here we use Vi=Vtot N, whereVtotis the total volume of the design domain and Nis the number of particles in the system. From these stresses we can then compute Von Mises deviatoric stresses per particle: dev=q 2xxxxyy+2yy+ 32xy (11) having dropped the subscript ito simplify the notation. 332 3 Results 333 This section presents three optimization problems where 334 material nonlinearity may signi cant impact the result- 335 ing topologies. Results are rst obtained for the Linand 336 Weak interactions in Table 1, which are analogous to 337 the linear elastic and elastoplastic materials considered 338 in Ref. [19]. Two of the systems are also tested using the 339 asymmetric potentials Weak-C andWeak-T in Table 1. 340 3.1 Three-support system with imposed displacement 341 from the top 342 The system is shown in Fig. 3; it is analogous to one 343 originally analysed in Ref. [19] using FEM-based TO. 344 Fig. 4.a shows the topology resulting from DETO when 345 the material is linear elastic. Most of the structure gets 346 concentrated into a central pillar, which provides the 347 shortest and sti est path to transfer the load from the 348 point A down to the central support. The bene t of 349 increasing the cross section of the pillar is limited by 350 10 Enrico Masoeroet al. H = 59 mmL = 119 mmL /10L /10L /10uimp= 0.4 mmAL /5H 20 mmL = 80 mmuimp= 0.4 mm AH/2FixedFixed(a)(b) Fig. 3 Optimization problem for a beam on three-supports with imposed displacement at point A at midspan. The value ofuimphas been ne-tuned to obtain an appreciable impact of material nonlinearity. the size of the support, to the extent that for the target 351 solid fraction used here, f= 0:3, additional sti ness is 352 gained by creating diagonal branches that reach for the 353 lateral support, despite such branches are longer than 354 the central pillar and thus contribute less eciently 355 to the overall sti ness. A similar result was obtained 356 in Ref. [19] using FEM-based TO; in that work, how- 357 ever, the optimum structure did not feature the diag- 358 onal branches. This di erence may be due to the dif- 359 ference between an FE-based description and our DE- 360 based one. Another possible explanation lies in di er- 361 ent optimization procedures, e.g.the di erent updat- 362 ing schemes for or parameters such as the maximum 363 change of iallowed between subsequent optimization 364 steps. In any case, additional simulations not presented 365 here have shown that the overall sti ness changes only 366 very slightly when the mass is all concentrated into the 367 central pillar, rather than being partly distributed to 368 the thin diagonal branches in Fig. 4.a. 369Fig. 4.b shows the optimization result for the sym- 370 metrically nonlinear material. The limiting factor for 371 Utotin this case is that some pairs of particles may 372 reach the maximum asymptotic value of their interac- 373 tion force (see Fig. 2), thus entering into the analogous 374 of a plastic ow regime. This happens near the supports 375 and under the plate applying the imposed displacement, 376 as shown by the sharp diagonal fronts of large devia- 377 toric stress in Fig. 4.b. The magnitude of such stresses 378 is consistent with the quantitative estimation in Sec- 379 tion 2. Because of these mechanisms, the response of 380 the system is controlled by the thinnest cross section 381 across which the load is transferred. Therefore, if the 382 thick central pillar in Fig. 4.a was retained, all its mass 383 in excess to its smallest cross section, vizthe size of 384 the support below it, would not contribute to the max- 385 imumUtot. Therefore, when the material nonlinearity is 386 considered in the optimization process, the excess mass 387 is removed from the central pillar and used to thicken 388 the diagonal branches, exploiting as much additional 389 area from the lateral supports as possible. The result is 390 analogous to that obtained in Ref. [19]. 391 Fig. 4.c shows the evolution of the objective func- 392 tion,Utot, during the optimization process. As expected, 393 the weaker nonlinear material ends up with signi cantly 394 lowerUtot. The snapshots within the gure show how 395 the systems in Figs. 4.a and b appear after 8 optimiza- 396 tion steps only. Both systems then feature thick diago- 397 nal branches, but with the key di erence that the Weak 398 system is already clearly utilizing the branches (light 399 color meaning intense von Mises stresses in them), whereas 400 theLinsystem is not utilizing them signi cantly (dark 401 Topology Optimization using the Discrete Element Method. Part 2: Material Nonlinearity 11 Fig. 4 Optimization results for the three-support system in Fig. 3 with target solid fraction f= 0:3, assuming (a) linear elastic and (b) symmetrically strain-hardening material, as per Linand Weak expressions in Table 1. The colors represent the intensity of local deviatoric von Mises stress; (c) Evolution of the objective function Utot,vizthe total strain energy of the systems during the optimization. The base case with inputs in Table 2 is compared with cases with no ltering and with smaller _uimp; (d) Evolution of force{displacement curves during the optimization. color meaning little stress). As a result, at this step 402 during the optimization mass tends to move away from 403 the branches in the Lincase, whereas it tends to move 404 towards the branches in the Weak case. 405 Fig. 4.d shows the force{displacement curves for the 406 LinandWeak materials, evolving during the optimiza- 407 tion process. Clearly the nal solutions are much stronger 408 than the initial ones, where all particles had i=f= 409 0:3. At small displacements the two systems feature 410 similar sti ness, whereas the nonlinearity caused by the 411 material in the Weak system becomes evident at larger 412 uimp. 413 For both types of material, Fig. 4.c and d compare 414 results for three di erent cases: the base case with in- 415 put data in Table 2, the base case but without lter- 416ing, and the base case but with a smaller loading rate 417 _uimp= 105(instead of 104mms1in the base 418 case). For the linear material all cases give identical re- 419 sult. For the nonlinear material, instead, the case with- 420 out ltering reaches a less optimal solution with lower 421 Utot, whereas the other two cases returns the same evo- 422 lution ofUtot. A close scrutiny of the force{displacement 423 curves for the Weak system indicates that the curves 424 for to the base case are the highest, suggesting a more 425 optimum outcome. However, when reaching the target 426 uimp, the base systems continues to minimize its strain 427 energy which causes a drop of force while uimpremains 428 xed at 0.4 mm. The case with lower loading rate fea- 429 tures a lower curve but with no further relaxation at 430 uimp= 0:4 mm. As a result, both the base case and 431 12 Enrico Masoeroet al. the slower one attain the same nal value of force, and 432 thus ofUtot, atuimp= 0:4 mm; this explains why their 433 Utotare identical in Fig. 4.c. By contrast, the force{ 434 displacement curves for the case without ltering are 435 intermediate between the base and slower cases but, 436 whenuimp= 0:4 mm is reached, a signi cant further 437 relaxation sees the force dropping below those of the 438 other cases (see the thin vertical lines at uimp= 0:4 in 439 Fig. 4.d). This explains why Utotin Fig. 4.c is smaller in 440 this case than for the others. In terms of geometry evo- 441 lution, what limits the un ltered case is that the system 442 rapidly converges to a con guration with all i0 or 1, 443 getting e ectively stuck into a local energy minimum. 444 This hinders a full transfer of mass towards the lat- 445 eral branches, hence a full exploitation of the supports. 446 A similar e ect of un ltered systems getting stuck into 447 sub-optimal local minima was already observed and dis- 448 cussed in Part 1 for systems with linear material NOTE 449 TO EDITOR: please cite or otherwise link Part 1 here. 450 Fig. 5 shows results that are particularly relevant 451 for structural design. The Lin from Weak series ex- 452 plores how the structure in Fig. 4.b, optimized for a 453 nonlinear material ( vizfor best performance approach- 454 ing failure), behaves in the linear elastic range. The 455 results show that the sti ness of the structure is lower 456 than that in the Linstructure, which was originally op- 457 timised assuming a linear elastic material. The loss in 458 sti ness is 11%, from a gradient of 53.5 kN/mm in the 459 Lincase to 47.75 kN/mm in the Lin from Weak case. 460 In the same gure, the Weak from Lin series explores 461 how the structure in Fig. 4.a, optimized for a linear 462 material ( vizfor maximum sti ness in service condi- 463 0 5 10 15 20 25 0 0.1 0.2 0.3 0.4Ftot (kN) uimp (mm)Lin Weak 0 5 10 15 20 25 0 0.1 0.2 0.3 0.4Ftot (kN) uimp (mm)Lin from Weak Weak from LinFig. 5 Force{displacement curves for the con guration in Fig. 4.b assuming linear elastic material ( Lin from Weak ), and for the con guration in Fig. 4.a assuming nonlinear ma- terial ( Weak from Lin ). The curves are compared with the base cases for linear and nonlinear materials already shown in Fig. 4.d (solid curves). All curves here were obtained using loading rate _ uimp= 105mms1. tions) behaves when approaching failure. The results 464 show that the maximum force and the strain energy 465 atuimp= 0:4 mm are both substantially smaller than 466 in the Weak structure, which was originally optimized 467 assuming nonlinear material. The maximum force and 468 strain energy go from 7.83 kN and 2.34 kN mm for the 469 Weak case, to 6.12 kN and 2.11 kN mm for the Weak 470 from Lin case, decreasing by 22% and 10% respectively. 471 A 15% loss in maximum force was obtained in Ref. [19] 472 for a system with same geometry, but using FEM-based 473 TO and elastoplastic material. An 11% loss of sti ness 474 in service conditions is likely to be less problematic than 475 a 22% loss of strength approaching failure. Therefore, 476 the designer should use TO with linear elastic materials 477 carefully and favour optimization using realistic mate- 478 rial behaviors when addressing strength and structural 479 failure. 480 Topology Optimization using the Discrete Element Method. Part 2: Material Nonlinearity 13 3.2 Three-support system with mid support settlement 481 The system in Fig. 6 has very similar geometry as the 482 previous one in Fig. 3. The di erences are that the dis- 483 placement is imposed at the mid support instead of 484 above the beam, and that the lateral supports are only 485 half as wide as before. This problem was also originally 486 addressed in Ref. [19], there using FEM-based TO with 487 elastoplastic material. Fig. 6 Optimization problem for a beam on three supports with imposed settlement of the mid support. 488 Fig. 7 shows the optimum geometries obtained from 489 DETO. In all cases, the resisting mechanism is akin to 490 that in the seminal work of Michell [3], where the cen- 491 tral ties connect the settling plate to the compressed 492 arch above, which transfers the load to the stable lat- 493 eral supports. The linear elastic Lincase produces a 494 structure that is very similar to that in Ref. [19], de- 495 spite the already mentioned methodological di erences. 496 A material that is strain-hardening both in tension and 497 in compression produces the Weak structure in Fig. 7.b, 498 with a attening of the arch at is its top and with fewer 499 thicker ties linking the settling mid support with the 500 compressed arch. Another important detail is that the 501Weak structure concentrates more mass near the later 502 supports, which are instead not fully utilized in the Lin 503 case. An analogous tendency to fully exploit the sup- 504 ports has been already discussed in the previous section, 505 and was also observed in Ref. [19] for this case study. 506 The structures in Fig. 4, in the previous section, 507 were fully under compression when loaded, hence con- 508 sidering asymmetric materials in tension and compres- 509 sion was not useful then. Here instead, Fig. 7 shows how 510 asymmetric material behaviors lead to di erent opti- 511 mum structures. In particular, Fig. 7.c shows that a ma- 512 terial that is weak, vizstrain-hardening, in compression 513 and strong, vizstrain-sti ening, in tension produces a 514 structure with thin central ties under high stress, and 515 a thicker compressed arch that fully utilizes the lateral 516 supports. By contrast, in Fig. 7.d, a material that is 517 weak in tension and strong in compression creates thick 518 central ties and a shallower and thinner compressed 519 arch which utilizes the lateral supports only in part. 520 In this latter case, the limiting factor is the size of the 521 settling central support, which controls the maximum 522 cross section in tension and thus the maximum force 523 that the structure can carry. 524 Fig. 8 shows the force-displacement curves for the 525 four systems in Fig. 7. As expected, all curves start 526 with the same gradient in the initial linear regime. The 527 Weak-T system displays an initially increasing gradi- 528 ent, due to strain{sti ening in the compressed arch. 529 At displacements over 0.5 mm, however, the strain- 530 hardening behavior of the central ties takes over and 531 plastic ow caps the maximum force. The Weak sys- 532 tem features the smallest strength, but the Weak-C is 533 14 Enrico Masoeroet al. Fig. 7 Optimum geometries for the problem in Fig. 6, with solid fraction f= 0:25 and for di erent material behaviors as per Table 1: (a) linear elastic, (b) symmetrically strain-hardening, (c) strain-hardening in compression and strain-sti ening in tension, and (d) strain-hardening in tension and sti ening in compression. The colors represent the intensity of the deviatoric von Mises stress. 0 0.2 0.4 0.6 0.8 1 1 10 100Utot (J) Optimization step 0 0.2 0.4 0.6 0.8 1 1 10 100Utot (J) Optimization step 0 2 4 6 8 0 0.2 0.4 0.6 0.8Ftot (kN) uimp (mm)Lin Weak 0 2 4 6 8 0 0.2 0.4 0.6 0.8Ftot (kN) uimp (mm)Weak-T Weak-C Fig. 8 Force-displacement curves for the structures in Fig. 7, each with their respective material behavior. only marginally better, as opposed to the signi cantly 534 stronger Weak-T system. This happens because the 535 strength-controlling element in the Weak and Weak- 536 Csystem is the compressed arch. The Weak-C system 537 can transfer a bit more mass from the central ties into 538the arch, but eventually the minimum cross sectional 539 area of the arch is limited by the size of the lateral 540 supports, which both the Weak and Weak-C systems 541 utilize in full or almost. In the Weak-T system, instead, 542 strength is controlled by the central ties and therefore 543 the system has more freedom to move mass away from 544 the compressed arch and alter the overall geometry to 545 maximize its strain energy. 546 The di erent optimum solutions in Fig. 7 raise the 547 question of how much an incorrect assumption of ma- 548 terial behavior in the TO process may a ect the struc- 549 tural performance. As an example, consider a structure 550 where the elements under compression are con ned us- 551 ing FRP to induce strain-sti ening in a material that 552 would otherwise be symmetrically strain-hardening. In 553 our model, this means turning a Weak system into a 554 Weak-T one. If optimized assuming Weak-T behavior, 555 Topology Optimization using the Discrete Element Method. Part 2: Material Nonlinearity 15 the geometry in Fig. 7.d would be obtained. However, if 556 the FRP system failed in the actual structure, the ma- 557 terial behavior would go back to Weak , for which the 558 optimum geometry would be that in Fig. 7.b instead. 559 This raises two questions: how much strength loss may 560 be caused by an incorrect assumption of material be- 561 havior? Which of the four material behaviors considered 562 here would produce the most robust structure, in case 563 the material ends up behaving di erently? The results 564 in Fig. 9 address these questions. 565 Each sub gure in Fig. 9 shows how one of the op- 566 timized structures in Fig. 7 would behave for any of 567 the four material types in this manuscript. A rst take 568 is that all four structures, irrespective of the material 569 assumption underlying them, feature a similar force- 570 displacement curve when the material behaves linearly 571 (compare the black solid Lin curves across the four 572 sub gures in Fig. 9). This means that, for the struc- 573 tural system considered here, sti ness is not sensitive 574 to the geometric details and the risk of losing service 575 performance due to an incorrect material assumption 576 is low. A second take is that, for all the material be- 577 haviors considered here, the structure that has been 578 optimized assuming the correct type of material is the 579 one featuring highest strength. For example, consider 580 theWeak-T curves in all the sub gures in Fig. 9: the 581 one reaching the highest force is that in Fig. 9.d, where 582 the structure was indeed optimized assuming a Weak- 583 Tmaterial. The same applies to the other three mate- 584 rial types, con rming and extending the result in Fig. 5 585 in the previous section. 586 0246800.20.40.60.8Ftot(kN)LinWeak0246800.20.40.60.8Ftot(kN)Weak-TWeak-C 0246800.20.40.60.8Ftot(kN)0246800.20.40.60.8Ftot(kN) 0246800.20.40.60.8Ftot(kN)0246800.20.40.60.8Ftot(kN) 0246800.20.40.60.8Ftot(kN)uimp(mm)0246800.20.40.60.8Ftot(kN)uimp(mm)(a) orig. = Lin(b) orig. = Weak(c) orig. = Weak-C(d) orig. = Weak-T Fig. 9 Force-displacement curves for di erent material be- haviors and for structures originally optimized assuming the following material types: (a) linear elastic Lin, (b) symmet- rically strain-hardening Weak , (c) strain-hardening in com- pression and sti ening in tension Weak-C , and (d) strain- hardening in tension and sti ening in compression Weak-T . The snapshots of the optimized structures are identical to those in Fig. 7. 16 Enrico Masoeroet al. To address the question on strength loss from un- 587 expected material behavior, consider Fig. 9.d. Going 588 back to our example with the FRP, a structure opti- 589 mized assuming Weak-T material should resist a force 590 ofca.5 kN, if the material behaves as predicted. How- 591 ever, if the FRP system fails and the material ends up 592 behaving as Weak , the maximum force drops to 1 kN, 593 with an 80% strength loss that would likely entail col- 594 lapse. An analogous loss of strength would occur for 595 structures optimized assuming LinorWeak-C mate- 596 rials, in Fig. 9.a and c, albeit less pronounced in the 597 latter case due to the aforementioned, similar resist- 598 ing mechanisms in the Weak and Weak-C cases. The 599 only case not involving strength loss is that of a struc- 600 ture optimized assuming Weak material, in Fig. 9.b. 601 At rst sight, this may be simply reduced to a \design 602 for the worst-case scenario" message. However, design- 603 ing for the worst case is a way to de ne suitably large 604 cross sections for the various structural elements. Here 605 the problem is di erent, as optimization with xed f 606 implies that any increase in cross section at one place 607 requires a reduction of cross section elsewhere. Under 608 this constraint, it is a nontrivial nding that the geom- 609 etry optimized assuming Weak material gives the most 610 robust structure with respect to other possible material 611 behaviors. 612 3.3 Doubly xed beam 613 In the previous section, the load was transferred to the 614 lateral supports viaa serial arrangement of ties working 615 in tension, followed by the arch working in compression. 616 In this section, a problem is devised to obtain elements 617in tension working in parallel with elements in com- 618 pression. The system in Fig. 10 is proposed to this end; 619 it is similar, but not identical, to the system that was 620 used in Part 1 to highlight the impact of geometric non- 621 linearity NOTE TO EDITOR: please cite or otherwise 622 link Part 1 here. H = 59 mmL = 119 mmL /10L /10L /10uimp= 0.4 mmAL /5H = 20 mmL = 80 mmuimp= 0.4 mm AH /2FixedFixed(a)(b) H = 59 mmL = 119 mm L /20L /10L /20uimp= 0.4 mmAH =39 mmL = 135 mmuimp= 0.6 mm AL/15FixedH/5H/5Fixed Fig. 10 Optimization problem for a beam partly xed on both ends and with imposed settlement at the mid point. 623 Fig. 11 shows the optimum geometries obtained with 624 di erent assumptions on material behavior. Despite the 625 symmetry of the system in Fig. 10, the structure op- 626 timized assuming Linmaterial is asymmetric with re- 627 spect to the horizontal axis: see Fig. 11.a. The asymme- 628 try stems from geometric nonlinearity, which generates 629 additional tensile stresses and thus favors concentra- 630 tion of material in the lower half of the structure (see 631 Part 1 for more discussion on this point NOTE TO 632 EDITOR: please cite or otherwise link Part 1 here). 633 Fig. 11.b shows the optimum structure for a symmet- 634 rically strain-hardening material, Weak . The material 635 nonlinarity enhances the asymmetry caused by the ge- 636 ometric nonlinearity, while mass is more concentrated 637 in the main compressed arch and lower deck in tension, 638 removing some of the diagonal struts that were present 639 in Fig. 11.a. Fig. 11.c shows the optimum structure for 640 a material that is weak in compression only, Weak- 641 C. This case features further concentrates mass in 642 Topology Optimization using the Discrete Element Method. Part 2: Material Nonlinearity 17 Fig. 11 Optimum geometries for the problem in Fig. 10, with solid fraction f= 0:4 and for di erent material behaviors as per Table 1: (a) linear elastic, (b) symmetrically strain-hardening, (c) strain-hardening in compression and strain-sti ening in compression, and (d) strain-hardening in tension and sti ening in compression. The colors represent the intensity of the deviatoric von Mises stress. the lower deck, which is now fully exploited in tension, 643 whereas the weaker compressed arch is signi cantly re- 644 duced in size. An almost specular geometry, except for 645 a slight asymmetry due to geometric nonlinariy, is ob- 646 tained for the Weak-T material, as shown in Fig. 11.d. 647 Fig. 12 shows the force-displacement curves for the 648 four structures in Fig. 11. The results agree conceptu- 00.20.40.60.81 110100Utot(J)Optimizationstep00.20.40.60.81 110100Utot(J)Optimizationstep 02468 00.30.60.91.21.5Ftot(kN)uimp(mm)LinWeak 02468 00.30.60.91.21.5Ftot(kN)uimp(mm)Weak-TWeak-C Fig. 12 Preliminary force displacment curves for the doubly xed beam case. 649ally with those in the previous section, with all mate- 650 rials providing similar sti ness at small deformations, 651 and with signi cant di erences emerging at larger uimp. 652 As expected, the Weak material results in the lowest 653 strength. The Weak-C and Weak-T materials lead to 654 very similar force-displacement curves, which well re- 655 ect their almost specular geometries, combined with 656 their specular material behaviors (see Fig. 2 in Sec- 657 tion 2). Both structures with Weak-C andWeak-T ma- 658 terials overshoot the Lincurve atuimp<0:9 mm; this 659 is due to the strain-sti ening behavior of the Weak- 660 CandWeak-T materials respectively in tension and in 661 compression, which is eventually overtaken by strain- 662 hardening in compression and tension. 663 Fig. 13 explores how robust the structure in Fig. 11 664 are with respect to wrong assumptions of material be- 665 havior. The results in Fig. 13 corroborate those in Fig. 9 666 in the previous section. Namely, all structures feature 667 a similar sti ness, meaning comparable performance in 668 service conditions. By contrast, strength is sensitive to 669 material behavior and geometry. Out of the structures 670 18 Enrico Masoeroet al. 0246800.30.60.91.21.5Ftot(kN)LinWeak0246800.30.60.91.21.5Ftot(kN)Weak-TWeak-C 0246800.30.60.91.21.5Ftot(kN)0246800.30.60.91.21.5Ftot(kN) 0246800.30.60.91.21.5Ftot(kN)0246800.30.60.91.21.5Ftot(kN) 0246800.30.60.91.21.5Ftot(kN)uimp(mm)0246800.30.60.91.21.5Ftot(kN)uimp(mm) (a) orig. = Lin(b) orig. = Weak(c) orig. = Weak-C(d) orig. = Weak-T Fig. 13 Force-displacement curves for di erent material be- haviors and for structures originally optimized assuming the following material types: (a) linear elastic Lin, (b) symmet- rically strain-hardening Weak , (c) strain-hardening in com- pression and sti ening in tension Weak-C , and (d) strain- hardening in tension and sti ening in compression Weak-T . The snapshots of the optimized structures are identical to those in Fig. 7.considered here, only the structure assuming Weak ma- 671 terial preserves a similar strength if the material ends 672 up behaving di erently: see Fig. 13.b. Instead, struc- 673 tures optimized assuming Weak-C orWeak-T materi- 674 als, in Figs. 13.c and 13.d, would end up with as little 675 as half their design strength if the material turns out 676 to feature a di erent type of nonlinearity. This means 677 that, also for the parallel tension-compression system 678 considered here, assuming the weakest material behav- 679 ior for the optimization leads to the structure that is 680 most robust against other unexpected material behav- 681 iors approaching failure. 682 4 Conclusion 683 The manuscript has presented a rst application of DE- 684 based TO to materials with nonlinear stress-strain be- 685 havior. In the DEM, the material behavior is captured 686 by the interaction potentials between particles. The op- 687 timization problem under imposed displacement has 688 been formulated for complex interaction potentials, and 689 then particularized for potentials that model linear, strain- 690 hardening, and strain-sti ening materials under mono- 691 tonically increasing strain. 692 Four di erent material behaviors have been anal- 693 ysed, exploring their impact on three structural prob- 694 lems. General ndings are that: 695 {Structural optimization with linear or strain-sti ening 696 materials may lead to a partial utilization of the sup- 697 ports, whereas full utilization is to be sought when 698 materials are expected to enter a strain-hardening 699 regime near the supports themselves. This con rms 700 Topology Optimization using the Discrete Element Method. Part 2: Material Nonlinearity 19 previous literature results from FEM-based TO with 701 elastoplastic materials [19]. 702 {Assuming the correct material behavior for the opti- 703 mization will produce the strongest structure, if the 704 material eventually behaves as expected. Therefore, 705 structural otpimization with the frequently used as- 706 sumption of linear material behavior, is likely to 707 generate signi cantly sub-optimum structures to- 708 ward failure, when material nonlinearity takes over. 709 On the other hand, structures optimized assuming 710 nonlinear material have been shown to produce sim- 711 ilar sti ness in service conditions as structures that 712 have been optimized assuming linear material. There- 713 fore the safest approach in design is to optimize 714 structures for nonlinar material behaviours and large 715 deformations. 716 {If instead the material ends up behaving di erently 717 from what is assumed during TO, the structure may 718 end up being signi cantly weaker than predicted in 719 design, causing signi cant risk of failure. Out of the 720 material behaviors considered here, for both struc- 721 tures with serial and parallel tension-compression 722 load paths, assuming the weakest material for TO 723 has produced the most robust structural geometry 724 against unexpected material behaviors. 725 Overall, this manuscript has shown how the newly 726 proposed DETO method can incorporate material non- 727 linearity, and has provided some design-relevant insights 728 from several case studies. So far, material nonlinar- 729 ity has only been considered viareversible interactions 730 between particles. Our results are therefore represen- 731 tative of structures that approach failure in a quasi- 732static regime of monotonically increasing strain. De- 733 spite this limitation, the inclusion of material nonlin- 734 earity is an important contribution towards future in- 735 corporation of irreversible behaviors, including for ex- 736 ample inter-particle collisions and bond failure. This 737 outlines a path where DETO will eventually enable the 738 optimization of systems undergoing discontinuous and 739 post-failure processes, which are challenging the current 740 optimization methods based on the FEM or on other 741 continuum descriptions. 742 Acknowledgements The authors would like to acknowledge 743 Newcastle University Rocket High Performance Computing 744 service where simulations were run. 745 Con ict of interest 746 The authors declare that they have no con ict of inter- 747 est. 748