Structural Topology optimization is attract- ing increasing attention as a complement to additive manufacturing techniques. The optimization algorithms usually employ continuum-based Finite Element anal- yses, but some important materials and processes are better described by discrete models, for example gran- ular materials, powder-based 3D printing, or structural collapse. To address these systems, we adapt the es- tablished framework of SIMP Topology optimization to address a system modelled with the Discrete Element Method. We consider a typical problem of strain en- ergy minimization, for which we de ne objective func- Connor O'Shaughnessy School of Engineering, Newcastle University, Newcastle upon Tyne NE1 7RU, United Kingdom; Enrico Masoero School of Engineering, Newcastle University, Newcastle upon Tyne NE1 7RU, United Kingdom; E-mail: enrico.masoero@newcastle.ac.uk Peter D. Gosling School of Engineering, Newcastle University, Newcastle upon Tyne NE1 7RU, United Kingdomtion and related sensitivity for the Discrete Element framework. The method is validated for simply sup- ported beams discretized as interacting particles, whose predicted optimum solutions match those from a clas- sical continuum-based algorithm. A parametric study then highlights the e ects of mesh dependence and l- tering. An advantage of the Discrete Element Method is that geometric nonlinearity is captured without ad- ditional complexity; this is illustrated when changing the beam supports from rollers to hinges, which indeed generates di erent optimum structures. The proposed Discrete Element Topology Optimization method en- ables future incorporation of nonlinear interactions, as well discontinuous processes such as during fracture or collapse. Keywords Topology Optimization Discrete Element MethodGranular material Geometric nonlinearity 2 Connor O'Shaughnessy et al. 1 Introduction 1 In structural mechanics, Topology Optimization (TO) 2 refers to a family of computational methods to nd 3 structural solutions that optimize a set of target per- 4 formance indicators under a set of constraints [1, 2]. 5 Common examples are to minimize compliance to cer- 6 tain loads using a prescribed amount of material, or vice 7 versa minimize mass while obtaining a target compli- 8 ance. TO is valuable in the design of more sustainable 9 structures, where target performance may include min- 10 imal material waste or maximum robustness towards 11 uncertainties, linked for example to malicious attacks 12 or climate change. The scope of TO is being further ex- 13 tended by the fast-paced development of additive manu- 14 facturing, which is streamlining the production of com- 15 plex structures that were too dicult and expensive to 16 fabricate in the past [3]. 17 The idea of optimising the use of structural mate- 18 rial originates with the seminal work of Michell, who 19 put forward the concept of fully stressed structures [4]. 20 First formulations of TO were based on homogenization 21 methods [5, 6, 7]. In the 1990s these were largely sup- 22 planted by the Solid Isotropic Material with Penalisa- 23 tion (SIMP) method, which proved more ecient espe- 24 cially for complex TO problems requiring discretization 25 and numerical solution [5, 8]. In SIMP-based TO, the 26 objective function to optimize (e.g. minimum compli- 27 ance to certain loads) is a function of one design variable 28 only: a density eld assigning a value between 0 and 1 29 to each point in the design domain, or to each element 30 in a discretized system. In the objective function, the 31value of the design variable is raised to a power p, which 32 penalises intermediate values and pushes the optimal 33 solution towards a 0-1, void-solid only con guration. 34 SIMP-based TO problems can be solved using various 35 methods; a particularly ecient and very popular one 36 is the Optimality Criteria method. [8, 9, 10, 11, 12]. In 37 recent years TO has expanded to applications in various 38 areas, such as compliant mechanism design [13], natu- 39 ral convection problems [14], and bre reinforced ma- 40 terial design [15], also including new design objectives 41 such as structural vibration [16] and damage tolerance 42 [17, 18, 19]. New approaches have been developed too, 43 notably based on the Level Set [20, 21] and on the Phase 44 Field methods [22]. 45 When solving a discretized TO problem, the objec- 46 tive function is computed multiple times, typically us- 47 ing the the Finite Element Method (FEM)1. For ex- 48 ample, when minimizing structural compliance, FEM 49 analyses are repeatedly used to compute displacements. 50 The FEM is extremely ecient for linear elastic analy- 51 ses. Including geometric and mechanical nonlinearities 52 is more demanding, but use of nonlinear FEM analyses 53 is common too [25, 26, 27]. FEM analyses become more 54 problematic as a structure approaches failure, and in- 55 deed applications of topology otpimisation to problems 56 involving fracture are only very recent [28]. The re- 57 liance on FEM analyses has hindered application of TO 58 methods to inherently discrete problems, such as frag- 59 mentation, structural collapse, and granular mechan- 60 ics across length scales (atoms, nanoparticles, colloids, 61 1An exception is the work of Gong et al.[23, 24] who have applied an element-free Galerkin method based on moving particles to solve TO problems. Title Suppressed Due to Excessive Length 3 and macroscopic grains). These problems are typically 62 analysed with discrete simulations [29, 30]. Therefore, 63 arguably, a TO scheme that uses for example the Dis- 64 crete Element Method (DEM) [31] would be desirable. 65 The DEM is a numerical approach describing the 66 mechanical behavior of assemblies of discrete, inter- 67 acting particles. The method was rst formulated by 68 Cundall [32] and used to describe granular media such 69 as sands, soils, and powders [33]. The DEM can also 70 be applied as an approximation of continua, in par- 71 ticular when describing processes that involve fracture 72 [34, 35], fragmentation [36, 37], or structural collapse 73 [38, 39, 40]. Rigid body motion, impacts, and geomet- 74 ric and mechanical nonlinearity [41] are naturally cap- 75 tured by explicitly integrating the equations of motion, 76 or by nding static equilibrium con gurations using 77 energy minimization techniques. The incorporation of 78 these features into a TO procedure would thus open 79 up applications to a new and wide range of quasi-static 80 and dynamic, discrete problems. 81 This paper couples DEM simulations with classical 82 SIMP-based topology optimization, leading to a new 83 Discrete Element Topology optimization (DETO) method. 84 Section 2 provides a brief background on SIMP-based 85 TO, solved using the OC method, and on the DEM. 86 Section 3 presents the optimization problem for a 2D 87 system of disks interacting via harmonic potentials, equiv- 88 alent to linear elastic springs. A penalization scheme is 89 proposed, whereby the interactions between neighbor- 90 ing particles are modulated by a particle-speci c design 91 variable . This variable ranges between 0 and 1, akin 92 to the density in continuum-based problems, but we 93do not call it \density" because its only role is to pe- 94 nalise mechanical interactions, not the inertial behav- 95 ior of the particles. We then formulate the objective 96 function, which here is simply the strain energy of the 97 system, and the updating scheme for . A numerical 98 implementation of the DETO method is provided. Sec- 99 tion 4 validates the new method by comparing results 100 on simple beams with established ones from the liter- 101 ature. Then a parametric study highlights the e ects 102 of mass fraction, lter length, and number of particles 103 per unit volume, which is akin to mesh neness in FEM. 104 The parametric study highlights the impact of topologi- 105 cal complexity on the optimized values of the objective 106 function, con rms that the solutions from DETO are 107 mesh-independent, and shows that the checkerboard 108 problem that may arise in FEM-based optimization is 109 absent in DETO. Finally, Section 4 shows that DETO 110 naturally accounts for geometric nonlinearity, which has 111 a marked e ect on the optimal topologies. Part II of the 112 present manuscript series (NOTE TO EDITOR: please 113 cite or otherwise link the Part 2 paper here) incorpo- 114 rates material nonlinearity in DETO, generalizing it to 115 problems of ductility maximization under large defor- 116 mations, and presenting results for several case studies 117 that provide design-relevant insights. All these results 118 indicate the viability of the proposed method, as well 119 as its potential for further development and application 120 to more complex, discrete structural problems. 121 2 Background on SIMP-TO and DEM 122 This section presents some basic concepts in SIMP- 123 based Topology optimization and in the Discrete El- 124 4 Connor O'Shaughnessy et al. ement Method. The aim is to enable a better appre- 125 ciation of the innovative aspects introduced later in 126 Section 3. Readers who are already familiar with the 127 techniques can proceed directly to Section 3. 128 2.1 Basics of SIMP-based Topology Optimization 129 SIMP-based TO is a well established technique that 130 nowadays encompasses sophisticated algorithms and ap- 131 plications, and can include multiple design variables, 132 constraints, and objectives [42, 43]. This section will 133 not cover such richness. By contrast, we will present 134 SIMP-based TO only for one simple but very in uen- 135 tial example by Sigmund [44]. Thanks to its simplicity, 136 this example became the entry point into TO for many 137 researchers. We hope for a fraction of a similar success 138 when presenting the same example in the context of our 139 new Discrete Element TO method, later in Section 3. 140 A numerical TO problem starts from de ning the 141 boundary conditions, external loads, supports, and spa- 142 tial domain within which to de ne the geometric detail 143 of the structure. The domain is discretized into individ- 144 ual elements and, in the SIMP method, each element 145 is associated with one continuously distributed design 146 variable e2[0; 1]. The design variable represents the 147 density of material at that point, between void ( e= 0) 148 and fully solid ( e= 1). A uniform density eld is usu- 149 ally chosen as a starting point and an objective function 150 c( ) is de ned, which speci es the performance indica- 151 tor to be optimized. is the vector collecting the eof 152 all the individual elements. 153A typical optimization problem is written as: min :c( ) =NX e=1uT ekeue (1) subject to :V( ) V0=f (2) : 0 min e1 (3) In Eq. 1 the objective function is twice the total, linear- 154 elastic strain energy of the system; minimizing it also 155 minimizes the compliance to a set of external loads that 156 induced equilibrium displacements ueto the various el- 157 ements. The uevectors depend on , but the notation 158 in Eq. 1 omits it for better readability. Eqs. 2 and 3 159 are constraints that the optimal solution must satisfy. 160 The rst constraint xes the target volume of solid, 161 asV( ) =P e e;V0isVwhen the whole domain is 162 solid with = 1 everywhere; f2(0;1) is a constant. 163 The second constraint sets the bounds for ebetween a 164 minimum value minand the fully solid e= 1. In prin- 165 ciple one could use min= 0, but later in this section 166 we will show that a small but nonzero min, e.g. 103, 167 is needed when applying ltering. 168 In Eq. 1, keis the element sti ness matrix. The distinguishing feature of SIMP is that kedepends on eas a power law: ke= p ek0 (4) where k0is a constant, base sti ness matrix. This pe- 169 nalisation scheme plays a central role in the solution of 170 the optimization problem; this will be discussed later in 171 this section. The dependence of keon e, also omitted 172 for clarity in Eq. 1, causes the aforementioned depen- 173 Title Suppressed Due to Excessive Length 5 dence of ueon e. At a generic step in the optimization 174 process, the structure features a certain vector and 175 therefore each element has a corresponding ke; a Fi- 176 nite Element analysis provides the uecorresponding to 177 the imposed external loads for the current distribution 178 ofke, and all this determines the current value of the 179 objective function c. 180 The optimization problem in Eqs. 1-3 can be solved using various methods. The most popular one is the Optimality Criteria method, which provides the updat- ing scheme , viz. the expressions and algorithms to up- date at the generic optimization step while respecting the imposed constraints (or tending towards a solution that respects them). Using the OC method and impos- ing the constraints at each optimization step, Sigmund obtained the following updating scheme [44]: new e= old e @c @ e (5) @c @ eis an entry of the gradient of cwith respect to , 181 which is called sensitivity . =1 2is a numerical damp- 182 ing coecient to improve convergence. is a parame- 183 ter that changes at every step of the optimization and 184 rescales the sensitivity so that newrespects the con- 185 straint on total volume in Eq. 2. Additional care must 186 be taken to also guarantee that newfalls between min 187 and 1, as per constraint in Eq. 3. This can be achieved 188 by capping the values of epredicted by Eq. 5, but this 189 would a ect the rst constraint. Therefore, an iterative 190 algorithm is usually needed to nd a value of that re- 191 spects both the imposed constraints; for example, the 192 bi-sectioning algorithm implemented by Sigmund [44]. 193The analytical expression of the sensitivity is ob- tained by combining the de nition of cis Eq. 1 with the expression of the penalised kein Eq. 4: @c @ e=p p1 euT ek0ue (6) The combined e ect of Eqs. 5 and 6 is to push material 194 away from under-utilised areas and move the design 195 towards a solid-void only solution. The penalisation 196 scheme, combined with the presented udating scheme, 197 are only e ective if p>1; however, a value of p>3 is 198 usually preferred for reasons linked to the fabrication of 199 actual structures using only void or fully solid parts: see 200 Ref. [45] for discussion on this point. Higher values of 201 penforce stricter solid-void only solutions and also im- 202 prove the speed of convergence, but reduce the ability 203 to escape local minima of c( ) and therefore increase 204 the probability of nding sub-optimal solutions. 205 2.2 Basics of the Discrete Element Method 206 The DEM describes a system of mechanically interact- 207 ing particles. The particles can be subjected to external 208 forces, and exchange interaction forces that depend on 209 their relative positions rand sometimes orientations. 210 Velocity-dependent dissipative forces are often included 211 too, but here we will not consider them for simplicity, 212 and because in this manuscript we will only address 213 performance at static equilibrium. 214 As a simple example of interaction, consider a har- monic potential representing pairs of particles connected by linear springs, as depicted in Fig. 1. The interaction potential energy Uijand the intensity of the repulsive 6 Connor O'Shaughnessy et al. Fig. 1 Schematic of the interaction force between two parti- cles emerging from a harmonic interaction potential. interaction force Fijbetween the two particles iandj are: Uij=1 2kij(rijr0)2(7) Fij=dUij drij=kij(rijr0) (8) kijis the sti ness of the connecting spring, rijis the 215 inter-particle distance, and r0is the equilibrium dis- 216 tance. 217 The displacements of the particles under external 218 and interaction forces can be solved either dynamically, 219 by explicitly integrating Newton's second law of motion 220 for all particles, or by minimizing the total interaction 221 energy of the system Utot=P i;jUij; the latter ap- 222 proach leads to a nal con guration at static equilib- 223 rium, similar to what one would obtain with a linear- 224 elastic, static FEM solution. Various quantities can be 225 computed during a dynamic solution or at nal static 226 equilibrium, including the total strain energy that de- 227 nedcin the optimization problem in Eq. 1, and that 228 in the example here coincides simply with Utot. 2293 Methodology 230 This section starts with introducing the new combi- 231 nation of Topology optimization with DEM-based so- 232 lutions of structural behavior. A simple implementa- 233 tion of the new method in MATLAB to perform opti- 234 mization by strain energy minimization in 2D problems 235 with circular disks hexagonally packed in rectangular 236 domains, is then presented. The code is available on 237 GitHub [46]. 238 3.1 Discrete Element Topology Optimization 239 Here we introduce some fundamental changes to the SIMP-TO method in Section 2.1 to apply it to a system described using the DEM. The rst change is straight- forward: the design variable e, which previously was speci ed for each nite element, here becomes a per- particle quantity i. For the task of minimizing the strain energy of the system, considering the expression ofUtotpresented in Section 2.2, the optimization prob- lem becomes: min :c( ) =1 2NX i=1NX j>ikij(rijr0)2(9) subject to :V( ) V0=f (10) : 0 min i1 (11) The di erence with respect to the FEM-based problem 240 in Eq. 1 lies in the expression of the cost function, which 241 now features a sum over all pairs of particles, instead 242 of over individual elements. The prefactor1 2is incon- 243 sequential, but we keep it here to fully identify cwith 244 the total interaction energy of the system Utot. The in- 245 Title Suppressed Due to Excessive Length 7 terparticle distance rijand the interaction sti ness kij 246 are now scalar quantities pertaining to pairs of parti- 247 cles, whereas the FEM framework featured a vector of 248 nodal displacements and a sti ness matrix pertaining 249 to individual elements. The constraints in Eqs. 10 and 250 11 are the same as for the FEM-based problem. 251 A key change in the DEM framework concerns the penalisation scheme. In the FEM-based approach, since each nite element contributes individually to c, the e of each element penalises only the sti ness matrix of the element itself, as shown in Eq. 4. In the DEM con- text, however, kijis associated with pairs of particles rather than individual ones. Therefore, we propose the following penalisation scheme: kij= p i p jk0 (12) wherek0is a constant base sti ness and iand j 252 are the design variables of the two interacting particles. 253 The penalisation exponent pplays an analogous role as 254 discussed in Section 2.1. Numerical experiments indi- 255 cated thatp= 2 provides a good compromise between 256 optimality of the solution, solid-void only result, and 257 numerical performance, for the examples that will be 258 considered later in this manuscript. 259 The penalisation scheme in Eq. 12 imposes that 260 the interaction sti ness kijis a function of iand j. 261 Therefore, under a given set of external forces, also the 262 interparticle distance rijat equilibrium (or at a generic 263 step during a simulated dynamic response) will depend 264 on . In Eq. 9 the dependence of kijandrijon have 265 been omitted for clarity of notation. 266The solution of the optimization problem in Eqs. 9- 11 can be obtained with the same updating scheme as in Eq. 5. The only di erence is that@c @ emust be replaced with a sensitivity@c @ ithat di ers from Eq. 6 and that now reads as: @c @ i=X j6=ip p1 i p jk0(rijr0)2(13) Eqs. 9 to 13 complete the formulation of DEM-based TO for the speci c case of strain energy minimization using a harmonic pairwise interaction potential. A full generalisation of the method is beyond the scope of this manuscript, but one immediate and important ex- tension concerns more complex interaction potentials. The harmonic potential considered here represents lin- ear springs connecting the particles; this is analogous to a linear elastic constitutive law in the FEM. How- ever, one can replace the potential in Eq. 7 with a more complex form such as: Uij:::=k( i; j;:::)g(ri;rj;:::) (14) wherekis now a generic function of the design vari- ables of the interacting particles, and gis a function of the position vectors rof the interacting particles. The ellipsis indicate that the interactions can involve more than pairs of particle, including three-body terms, four- body, etc. The objective function then simply becomes c=Utot=X i;jUij (15) having omitted multi-body terms and ellipsis only for clarity of notation. The resulting sensitivity would then 8 Connor O'Shaughnessy et al. be: @c @ i=X j;:::6=i@k @ ig(ri;rj;:::) (16) Mechanical nonlinearities embedded into g(ri;rj;:::) can 267 be seamlessly considered in this DEM framework. For 268 example, a potential that is widely used in nanoscale 269 materials modeling is the Lennard-Jones potential, which 270 is akin to a nonlinear spring with a softening regime 271 after a peak force in tension. For the Lennard Jones 272 potential,kcan be expressed simply as in Eq. 12 and 273 g= 4  rij12   rij6 , whereis a constant. Writ- 274 ing Eqs. 15 and 16 for the Lennard-Jones potential 275 would thus be straightforward. The details on how to 276 incoporate material nonlinearity in DETO are given in 277 the companion Part 2 manuscript (NOTE TO EDI- 278 TOR: please cite or otherwise link the Part 2 paper 279 here). 280 The numerical eciency of solving a TO problem 281 is typically governed by the solver to compute the - 282 dependentcand sensitivities at each step of the opti- 283 mization. For the energy minimization problem in this 284 manuscript, eciency is thus controlled by the FEM 285 or DEM solvers that provide the con guration of the 286 system under load: ufor FEM, rfor DEM. The e- 287 ciency of the DEM solver depends on the type of anal- 288 ysis that is performed. For the speci c case of harmonic 289 potentials and for suciently small external loads, the 290 DEM analysis could also be expressed as a linear al- 291 gebra inverse problem, as for the FEM in linear elas- 292 ticity; in such a case, the performance of DEM and 293 FEM would be similar. However, a strength of the DEM 294 is the simplicity to consider geometric and mechanical 295nonlinearities, also including fracture and rigid motion, 296 both in quasi-static and dynamic regimes. To address 297 this complexity, DEM problems are typically solved us- 298 ing energy minimization or explicit integration, which 299 are more time-consuming than linear-elastic Finite Ele- 300 ment analyses. However, the numerical eciency of the 301 DEM is again comparable with that of the FEM when 302 the latter includes geometric and material nonlineari- 303 ties, or even fracture [47]. 304 3.2 Numerical implementation of DETO 305 This section describes a simple implementation of the 306 proposed DETO method, to solve the optimization prob- 307 lem in Eqs. 9-11 for a simple 2D system of monodisperse 308 circular interacting via a harmonic potential. The ini- 309 tial con guration features nelxnely disks2arranged 310 in a close-packed hexagonal lattice that lls a rectangu- 311 lar domain of size ( Dnelx)p 3 2Dnely , whereD 312 is the disk diameter: see Fig. 2. Initially all disks have 313 i=f2(0;1). The parameters nelx,nely,D, andf, 314 are chosen and provided by the user. Various external 315 forces and constraints can be applied to each particle: 316 these will be speci ed in the next section for each of 317 the examples studied there. 318 Each particle interacts only with its immediately ad- 319 jacent neighbors in the hexagonal lattice. The harmonic 320 potential is the same as in Eq. 7, with equilibrium dis- 321 tancer0=D, sti nesskijpenalised as per Eq. 12, and 322 base constant sti ness k0chosen by the user. The har- 323 2Actually, the number of elements per row alternate be- tween nelx andnelx1 to respect horizontal symmetry, so the exact number of disks is nelx1 2 nelywhen nelyis even, and nelx1 2 (nely1) +nelx when nely is odd. Title Suppressed Due to Excessive Length 9 nelxnely Fig. 2 Hexagonal close packing of DE disks with linear elastic interaction potentials between immediate neighbors monic bonds are modelled as unbreakable and the parti- 324 cles are not allowed to create new bonds with other par- 325 ticles that initially were not among their rst neighbors. 326 This restricts the scope of the DEM, which usually deals 327 with particles that move widely across the system, cre- 328 ating new bonds or colliding with particles that initially 329 might have been far away. In this manuscript, however, 330 we will validate and test the new DETO method for the 331 classical example of a simple beam under point load, for 332 which we expect relatively small deformations. In such 333 applications, the particles will indeed interact only with 334 their initial rst neighbors. 335 Fig. 3 shows the ow chart for the program. First 336 the system geometry is generated from the inputs as ex- 337 plained above, adding also the required external forces 338 and constraints (e.g. pinned or roller supports). The ini- 339 tial neighbor list is recorded and stays the same during 340 the whole simulation, for the reason discussed above. 341 This is much faster than a general case in which the 342 neighbor list must be updated dynamically during the 343 DEM simulation. 344 Fig. 3 Flow chart of the DETO algorithm implementation. The optimization loop begins by computing the in- 345 teraction sti ness kijfor each pair of neighboring parti- 346 cles, following Eq. 12. The DEM minimization module 347 computes the particle positions at static equilibrium 348 using a damped dynamics algorithm by Sheppard et 349 al. [48]. The DEM solution is considered as converged 350 when the change in total strain energy between two suc- 351 cessive steps is suciently small:(Ucurrent totUprevious tot ) Uprevious tot 352 etol. The values of etolused in this manuscript will be 353 in the 1010108range. The energy minimization al- 354 gorithm requires two parameters: a time step dtand a 355 maximum particle displacement allowed at the generic 356 stepdmax. These should be ne tuned depending on the 357 system that the user wants to analyse. The algorithm 358 10 Connor O'Shaughnessy et al. also uses the masses of the particles, here all set to the 359 same value m. 360 After the DEM module converges, the DETO pro- gram computes the objective function c(i.e. the total interaction energy Utot) and the sensitivity to update as per Eq. 13. The stress tensor for each particle is also computed; this is not strictly needed for the optimiza- tion process, but knowing the stress eld inside the sys- tem will support the interpretation of the results. The per-particles stress tensor is based on the virial stress expression [49]: ab;i=1 2ViX j6=i(ra;iFij;b;i+ra;jFij;b;j) (17) ab;iis theab(xx, xy, or yy) stress component at par- ticlei,ra;iandra;jare thea-component (x or y) of the positions of particles iandj,Fij;b;i andFij;b;j are the b-component of the force on particle idue to the inter- action with particle jand vice versa. Viis the averag- ing volume, here taken equal to the tributary volume of particle i, viz.Vtot N, whereVtotis the total volume of the rectangular domain (assuming unit thickness in the third dimension) and Nis the number of particles in the system. In particular, in this manuscript we will compute and plot the hydrostatic and Von Mises devi- atoric components of the per-particle stress tensors: hyd=xx+yy 3(18) dev=q 2xxxxyy+2yy+ 32xy (19) When the sensitivity is computed, the program can proceed with updating . However, TO algorithms of- ten add an intermediate step of ltering. The idea isthat the updating scheme in Eq. 5 does not use directly the sensitivity from Eqs. 6 or 13. Instead, it uses a new coarse-grained sensitivityc@c @ that, for the generic ele- menteor particlei, depends also on the sensitivities of neighboring elements or particles: d@c @ i=nfX k=1@c @ kWk k infX k=1Wk(20) In Eq. 20 we used the notation for the DEM-based al- 361 gorithm, thus referring to particle i; however, the exact 362 same expression also applies to FEM-based algorithms, 363 just replacing subscript iwithe. In the equation, the 364 sums are over the nfparticles (or elements) at a dis- 365 tancerik< rminfrom particle i, including particle i 366 too3. The ltering length rminis chosen and provided 367 by the user. Wk= max 1rik rmin;0 is a factor that 368 linearly reduces the weight of neighboring particle k 369 with its distance from the centre of particle i: its value 370 is 1 for particle iand becomes zero for particles with 371 rik> rmin. The presence of iat the denominator in 372 Eq. 20 is the reason why one should enforce min>0. 373 In this manuscript we will always set min= 103, as 374 is customary in the literature [44]. 375 The reasons to include ltering are both practical 376 and numerical. The practical reason is that rminim- 377 poses a minimum size of solid and void regions in the 378 nal structure; this provides some control over the com- 379 plexity of the optimal structure, which may help with 380 3We use counter kinstead of jfor the neighboring particles to clarify that the neighbor list for ltering is not in general the same as for the interactions; for example, if rmin>2D, also second nearest neighbors in the lattice will be included in the ltering even if they do not contribute to the interaction energy. Title Suppressed Due to Excessive Length 11 fabricability. The numerical reason is that optimisa- 381 tion processes not including ltering often converge too 382 rapidly to solid-void solutions getting e ectively stuck 383 into sub-optimum local minima. Some ltering (viz. a 384 smallrminD) usually removes these local minima 385 and leads to a better solution, although one must be 386 careful as a larger rminmay also smoothen the global 387 minimum and thus a ect the optimality of the solution. 388 These e ects will be shown and discussed later in Sec- 389 tion 4. When FEM analyses are used, another bene t 390 of ltering is to remove the checkerboarding problem 391 [50, 51]: when a ne FE mesh is used, individual neigh- 392 boring elements in the optimal solution typically create 393 an alternating pattern of void and solid. The problem 394 arises from a locking e ect in certain types of nite el- 395 ements [50, 51]. We will see that DETO does not su er 396 from checkerboarding. 397 The last step in the optimization loop is to update 398 following Eq. 5, but using the ltered sensitivities 399 instead of the original ones. The optimization loop is 400 repeated until new i old i4103for every par- 401 ticle. In the nal solution, some particles will feature 402 i= minand others, especially at solid-void inter- 403 faces, might be \gray", viz. feature a ithat is inter- 404 mediate between 0 and 1. The MATLAB implentation 405 of DETO that accompanies this manscript includes an 406 optional post-processsing module to reduce the solution 407 to a solid-void only system, where all particles have ei- 408 ther i= 0 or 1, while respecting the constraint on 409 the total solid volume fraction f. In the simulations 410 for this manuscript, however, we did not perform any 411 post-processing to present straight optimization results 412without further manipulations. All the outputs from 413 the simulation are recorded in an XYZ con guration 414 le following the evolution of the structural topology, 415 and in a text le collecting values such as candV( ) 416 after each optimization step. 417 4 Results 418 This section starts by comparing results from DETO 419 with others from established optimization methods. The 420 comparison is done using a base system representing 421 a simple beam under point load. We then present a 422 parametric study where the base system is altered one 423 feature at the time. The simulations reveal how the op- 424 timal solutions from DETO are impacted by the target 425 solid volume fraction f, the ltering length rmin, and 426 the number of particles per unit volume (akin to mesh 427 resolution in the FEM). Finally, the supports of two 428 simple beams are arranged in such a way to show how 429 DETO captures the e ect of geometric nonlinearity on 430 the optimal solutions. All of the simulations snapshots 431 in this section have been obtained using the open source 432 software OVITO [52]. 433 4.1 Validation 434 New topology optimization methods are typically tested 435 on simple structures with known optimal geometries. 436 Two such structures are the simply-supported and pin- 437 supported beam systems shown in Fig. 4. The gure 438 shows results from DETO side by side with optimal ge- 439 ometries from established methods. The input param- 440 12 Connor O'Shaughnessy et al. eters for the DETO simulations are shown in Table 1. 441 The intensity of the external force is 1 kN in both cases. 442 Table 1 Input parameters for the DETO simulations on the base system in Fig. 4. nelx 135 k0 100 kN/mm nely 45 dt 0.01s f 0.6 dmax 0.01 mm D 1 mm m 1 mg rmin 1.5 mm tz 1 mm To improve the physical interpretation of the re- sults, consider that k0EA r0, whereEis the Young modulus of the material, ADtzis the contact area between two neighboring disks ( tzis the thickness of the structure in the third direction), and r0=Dis the length at rest of the spring connecting neighboring par- ticles. Rearranging the expression of k0we can estimate an equivalent Young modulus: E=k0r0 A=k0 tz(21) The values in Table 1 return E= 100 GPa, thus one can 443 consider the simulated structures as made of sintered 444 metallic powder. 445 For the simply supported beam in Fig. 4.a the op- 446 timal geometry from DETO is qualitatively similar to 447 that from the FEM code in Ref. [44], when similar in- 448 puts are provided. For the pinned structure in Fig. 4.d, 449 the optimal layout predicted by DETO is analogous to 450 the theoretical solution of fully stressed structures by 451 Michell [4]. In this case the similarity is less striking be- 452 cause Michell's solution features one-dimensional mem- 453 bers with the additional constraint that all members 454 must have equal cross section. Overall, Fig. 4 shows that 455the optimal solutions obtained using DETO are com- 456 parable to those coming from other more established 457 methods in the literature. 458 Fig. 5 shows the evolution of the objective function 459 c(viz.Utot) and of the corresponding geometry during 460 the optimization process. Signi cant changes in both c 461 and geometry take place during the rst 50 optimiza- 462 tion steps. Between steps 50 and 500, the geometry has 463 practically converged and cremains nearly constant. 464 Fig. 6 shows the distribution of hydrostatic and devia- 465 toric stress in the optimized structure from Fig. 4.a. The 466 deviatoric stress distribution is clearly bimodal, with 467 the lower deck of the structure and the lateral diagonal 468 elements carrying the tensile stresses, and with the up- 469 per arch and central diagonal struts being under com- 470 pression. The distribution of deviatoric von Mises stress 471 is instead centered around a value of 75 MPa quite uni- 472 formly distributed over the structure, with stress con- 473 centrations at midspan (top and bottom points) where 474 the axial stresses due to the bending moment are high- 475 est, and near the supports, where the shear stress gets 476 concentrated before transferring to the pointwise sup- 477 ports. 478 4.2 Parametric study: volume fraction 479 The base system in Fig. 4.c featured a nal volume frac- 480 tion of solid f= 0:6, as per Table 1. Here the system 481 is kept the same except f, for which 4 additional val- 482 ues are explored between 0.5 and 0.7. Fig. 7 shows the 483 impact offon the nal 0-1 optimized structure. As ex- 484 pected, small fvalues force the system to create fewer 485 elements and lead to less optimal solutions, with higher 486 Title Suppressed Due to Excessive Length 13 Fig. 4 Optimal structures for the beam problem under point load. (a) Simply supported system with (b) result from the FEM-based TO code in Ref. [44], with inputs: nelx = 67, nely = 39, volfrac = 0:6,penal = 3 and rmin = 2 (see Ref. [44] for details on the meaning of those inputs), and (c) result from DETO. (d) Pin-supported system with results from (e) Michell's analysis in Ref. [4] and (f) DETO. Fig. 5 Simply supported beam case in Fig. 4.a: evolution of geometry and objective function cduring the rst 500 opti- mization steps. ccompared to more topologically rich solutions at high 487 fas can be seen in Fig. 8. However, these results are 488 not sucient to determine how much the lower cvalues 489 at higherfcome from topological complexity rather 490 than just having used more material. The next section 491 will add insight to this point. 492 Fig. 6 Hydrostatic and Von Mises deviatoric stress distribu- tions in the optimized structure from Fig. 4.c. 14 Connor O'Shaughnessy et al. Fig. 7 E ect of target solid volume fraction fon optimized solutions to the simply supported beam problem in Fig. 4.a. Fig. 8 E ect of volume fraction on objective function 4.3 Parametric study: ltering length 493 Fig. 9 shows how the ltering length rminimpacts the 494 optimized geometries. The gure also shows a result for 495 the un ltered case. Predictably, the optimized topolo- 496 gies become simpler at higher rminvalues, which force 497 the solid to concentrate into fewer, thicker structural 498 elements. Reducing topological complexity by ltering, 499 however, constrains the optimization problem; as a re- 500 sult,cis expected to increase as the solutions become 501 Fig. 9 E ect of ltering length on optimized geometries. For generality, the values of rminare given in units of particle diameters D. The label o indicates the un ltered case. less optimal at larger rminvalues. This is con rmed in 502 Fig. 10, which shows cgrowing from 0.111 J to 0.117 J 503 asrminis increased from 1 :1Dto 3D. This complements 504 the discussion of Fig. 7 in the previous section, showing 505 indeed that more optimal solutions can be obtained by 506 increasing topological complexity while keeping f xed. 507 As the optimal topologies get simpler with ltering, 508 the structural geometries with f= 0:6 forrmin2D 509 in Fig. 9 end up resembling those in Fig. 7 for smaller 510 f= 0:5. Comparing these two examples con rms the 511 expected trend that similar geometries with smaller f 512 lead to higher values of c, as more porous structures 513 are less sti : cf. c= 0:132 J for the structure with 514 f= 0:5 in Fig. 7 with c= 0:114 J for the structure 515 withrmin= 2Din Fig. 9. 516 The top image in Fig. 9 shows a case without lter- 517 ing. Expectedly, the resulting topology is most complex 518 Title Suppressed Due to Excessive Length 15 Fig. 10 E ect of lter length on objective function compared to the other cases with ltering on. However, 519 less intuitively, the resulting cis higher than in most 520 ltered cases: see Fig. 10. This con rms the e ect that 521 we already mentioned in Section 3, whereby lack of l- 522 tering causes fast convergence to a local minimum of c. 523 Filtering tends to smoothen out and remove local min- 524 ima, thus leading to more optimal solutions. The top 525 structure in Fig. 9, obtained without ltering, shows 526 a few very thin elements but no checkerboard e ect, 527 viz. an alternating pattern of individual particles with 528 i= 0 and 1. This is because the locking problem lead- 529 ing to the checkerboard e ect is speci c to FEM-based 530 analyses [50]. 531 4.4 Parametric study: mesh resolution 532 In FE analyses, the size of the elements discretizing 533 the continuum is in principle arbitrary. Therefore, when 534 performing FEM-based TO one must monitor the im- 535 pact of mesh resolution on the results, as in some cases 536 the problem might display nonuniqueness and even nonex- 537 istence of the solution [51]. By contrast the DEM, in 538 its basic formulation, does not feature a mesh at all, as 539 particles represent physically distinct units. However, in 540practice, the particles in DE analyses are often coarse 541 grained representations of richer underlying microstruc- 542 tures; for example one particle might summarise a col- 543 lection of smaller grains. In other cases, like the simple 544 beams in this manuscript, the particles actually dis- 545 cretize a continuum to explore its behavior under sce- 546 narios that the FEM might struggle with, e.g. collisions, 547 fracture, fragmentation, or collapse. Therefore, also in 548 the DEM there can be some arbitrariness in deciding 549 the number and size of particles, which thus becomes 550 analogous to deciding the mesh resolution in FE anal- 551 yses. 552 To mimic the role of mesh resolution in FE analy- 553 ses, here we consider a rectangular design domain with 554 xed dimensions and we use di erent numbers of par- 555 ticlesnelx andnely initially lling the domain. When 556 solving problems with greater nelx andnely than our 557 base case, we reduced accordingly the particle diame- 558 tersDto always ll the same domain. When chang- 559 ing particle sizes in DE analyses, one should be care- 560 ful that the intensity of the interaction may depends 561 onD, as opposed to FE analyses where the constitu- 562 tive parameters describing the material are intrinsically 563 mesh-independent, e.g. the Young modulus E. Specif- 564 ically for our system, however, Eq. 21 shows that k0 565 andEare simply linked by the thickness of the simu- 566 lation domain in the z direction, tz; since we keep the 567 latter always constant and equal to 1 mm, we do not 568 need to change k0when changing D. This is not always 569 the case; for example, in a 3D simulation with spherical 570 discrete elements, k0EDand therefore k0would be 571 proportional to D. 572 16 Connor O'Shaughnessy et al. Fig. 11 presents optimization results for structures 573 with a range of mesh resolutions around our base case. 574 In all cases, a ltering length rmin= 1:5 mm is applied, 575 as in our base case. The resulting geometries are gener- 576 ally insensitive to the mesh resolution, except for small 577 changes such as an additional horizontal element ap- 578 pearing at low resolution. Fig. 12 shows that the qual- 579 ity of the solution is very similar for all the structures 580 in Fig. 11, viz. they all feature a similar value of cat 581 the end of the optimization process. This result also 582 proves that, for the 2D examples in this manuscript, it 583 is indeed correct to consider k0as independent of D. 584 Fig. 11 E ect of mesh resolution on optimum geometries and associate objective function cvalues: ltering included with length rmin= 1:5 mm in all cases. The results in Fig. 11 do not tell whether mesh inde- 585 pendence is intrinsic to the DETO method or is a result 586 of ltering. Indeed, mesh-independent ltering, with 587 xedrminirrespective of the mesh resolution, is known 588 to enforce mesh independence also in mesh-sensitive 589 FEM-based TO [50]. To clarify this point, Fig. 13 shows 590 optimization results for structures with various mesh 591 Fig. 12 E ect of mesh resolution on the evolution of cduring the optimization process, for the systems whose nal geome- tries are shown in Fig. 11. resolutions and without ltering. Also in this case the 592 solutions are essentially mesh-independent. Higher res- 593 olutions (viz. greater nelx andnely) lead to smoother 594 geometries and also remove some of the thin elements 595 appearing at lower resolutions. This indicates that mesh 596 re nement helps the system avoid or escape some of the 597 local minima, although in none of the cases in Fig. 13 598 the system is actually able to reach the more optimum 599 solutions obtained when ltering is applied. However, 600 the similarity between all the solutions in Fig. 13 shows 601 that the DETO method is not mesh-sensitive. 602 4.5 Geometric nonlinearity 603 When structures are subjected to high loads and un- 604 dergo large displacements, geometric nonlinearity should 605 impact the optimal topologies [53]. In FEM-based TO, 606 accounting for geometric nonlinearity requires additional 607 complexity in the formulation of the problem [25, 27]. 608 By contrast, an advantage of DETO is that geometric 609 nonlinearities are captured without any change to the 610 theoretical framework, because interactions are always 611 Title Suppressed Due to Excessive Length 17 Fig. 13 E ect of discretization on optimum geometry when no ltering is applied. computed with reference to the system in its deformed 612 con guration. 613 Fig. 14 highlights the potential impact of geometric 614 nonlinearity by considering two beam systems that are 615 identical to our base case study in Figs. 4.a,c, except 616 that: (i) the supports are applied to the central axis 617 instead of the bottom corners, and (ii) a larger point 618 load of 10kN is applied to the center of the beams in- 619 stead of above or below them; this larger load has been 620 chosen to induce larger displacements and thus better 621 appreciate the e ect of geometric non-linearity. Under 622 these new conditions, simulations assuming small dis- 623 placements should return identical solutions for both 624 the simply-supported and the pinned systems, because 625 the length of the central axis remains unchanged af- 626 ter small displacements. Indeed the results from linear 627 elastic FEM-based TO in Fig. 14.b,f are identical. By 628 contrast, considering large displacements should lead 629 to di erent optimum geometries, because the length of 630the central axis would remain unchanged in the simply- 631 supported case, where the rollers can move inward, but 632 would increase in the pinned scenario due to the nite 633 de ection. The impact of geometric nonlinearity is in- 634 deed captured by our DETO solutions in Figs. 14.c,g, 635 which feature very di erent geometries for the two sys- 636 tems. 637 The optimum solution for the simply-supported beam 638 in Fig. 14.c is similar to the linear-elastic solution in 639 Fig. 14.b. The reason is that the inward motion of the 640 rollers allows the structure to behave in pure bending 641 also when geometric nonlinearities are included. The 642 distribution of hydrostatic stresses in Fig. 14.d shows 643 indeed a symmetric distribution of elements working in 644 tension and in compression. The qualitative di erence 645 in Fig. 14.g stems from the catenary action induced by 646 the pinned supports. During the de ection, the central 647 axis of the pinned system is stretched and this generates 648 a tensile stress along the beam. During the optimiza- 649 tion process, this additional tensile stress drives mate- 650 rial away from the compressed regions and towards the 651 parts under tension. As a result, Fig. 14.h displays a 652 thickening of the lower deck, which carries most of the 653 catenary force, whereas the upper arch in compression 654 becomes smaller and migrates towards the centre of the 655 beam. 656 5 Conclusion 657 We have formulated SIMP-based TO for a problem of 658 strain energy minimization in a 2D system of circular 659 disks, modelled and analysed using the Discrete Ele- 660 18 Connor O'Shaughnessy et al. Fig. 14 Solutions of topology optimization problems highlighting the impact of geometric nonlinearity. (a) Simply-supported and (e) pinned systems, with supports and forces applied to the central axis of the beam; (b,f) solutions from linear elastic FEM-based TO using the code in