Topology optimization using the discrete element method. Part 2: Material nonlinearity

Structural Topology Optimization typically features continuum-based descriptions of the investigated systems. In Part 1 we have proposed a Topology Optimization method for discrete systems and tested it on quasi-static 2D problems of stiffness maximization, assuming linear elastic material. However, discrete descriptions become particularly convenient in the failure and post-failure regimes, where discontinuous processes take place, such as fracture, fragmentation, and collapse. Here we take a first step towards failure problems, testing Discrete Element Topology Optimization for systems with nonlinear material responses. The incorporation of material nonlinearity does not require any change to the optimization method, only using appropriately rich interaction potentials between the discrete elements. Three simple problems are analysed, to show how various combinations of material nonlinearity in tension and compression can impact the optimum geometries. We also quantify the strength loss when a structure is optimized assuming a certain material behavior, but then the material behaves differently in the actual structure. For the systems considered here, assuming weakest material during optimization produces the most robust structures against incorrect assumptions on material behavior. Such incorrect assumptions, instead, are shown to have minor impact on the serviceability of the optimized structures.


Introduction
Structural topology optimization (TO) is a family of methods to distribute mass within a design domain and maximize the utilization of material under a set of imposed loads and constraints [1,2]. First ideas of structural of TO date back to the early twentieth century, with the analytical work of Michell [3]. The advent of modern computers rapidly expanded the capabilities and scope of TO; homogenisation-based TO came first [4][5][6] and was then mostly replaced by the Solid Isotropic Material with Penalisation (SIMP) approach [4,7]. A crucial part of any TO algorithm Abstract Structural Topology Optimization typically features continuum-based descriptions of the investigated systems. In Part 1 we have proposed a Topology Optimization method for discrete systems and tested it on quasi-static 2D problems of stiffness maximization, assuming linear elastic material. However, discrete descriptions become particularly convenient in the failure and post-failure regimes, where discontinuous processes take place, such as fracture, fragmentation, and collapse. Here we take a first step towards failure problems, testing Discrete Element Topology Optimization for systems with nonlinear material responses. The incorporation of material nonlinearity does not require any change to the optimization method, only using appropriately is the calculation of objective functions depending on the distribution of material during the optimization process. This step requires structural analyses, for which continuum-based methods such as the the Finite Element Method (FEM) are the norm. In Part 1, we have coupled SIMP-based TO with Discrete Element (DE) analyses [8]. The resulting Discrete Element Topology Optimization (DETO) offers a pathway to consider processes that challenge continuum-based descriptions, such as granular behaviors [9], fracture [10,11], fragmentation [12,13], and collapse [14][15][16]. DETO may therefore impact scientific communities whose favor for discrete analyses has been precluding access to topology optimization, for example in soil mechanics or in nanoscale materials modelling.
An advantage of DETO is that, in principle, it could be used to target performance indicators that involve discontinuous structural behaviors, such as resistance to fracture or collapse. In such near-failure or even post-failure regimes, the stress-strain behavior of many engineering materials is nonlinear. Material nonlinearity has been included in continuum-based TO since the mid 1980's [17,18], mostly focusing on elastoplasticity: see review in Ref. [19]. From a structural design perspective, an important finding has been that structures optimized for serviceability performance, viz assuming linear elasticity, may be significantly sub-optimal towards failure, when the material behaves nonlinearly [20,21]. For discrete systems, however, DETO has only been applied to linear elasticity thus far; geometric nonlinearity from large displacements is included too in Part 1 [8] as it is naturally captured in DE analyses, where particle interactions are always computed in the deformed state. However, material nonlinearity in DETO is still to be addressed.
This manuscript presents a first application of DETO to problems with material nonlinearity. Section 2 describes the DETO method for a general problem of interaction energy maximization under imposed forces, and then particularises it for the 2D structures in this manuscript. We show how nonlinear material behaviors can be included seamlessly, without any change to the methodological framework. We then discuss how, under imposed displacements instead of forces, the same optimization scheme maximises the structural ductility, which is a typical objective function when considering material nonlinearity. The last part of Sect. 2 presents four interaction potentials to be used throughout the manuscript, describing linear elastic materials as well as nonlinear materials undergoing strain-hardening (i.e. becoming less stiff) or strain-stiffening (i.e. becoming stiffer) under tension, compression, or both. Section 3 considers three optimization problems including material nonlinearty. Two of these problems are known in the literature on continuum-based TO [20], whereas the third one is proposed here to further appreciate the effect of material nonlinearity. Particular attention is paid to the scenario in which a structure is optimized assuming a certain material behavior, but then the material turns out behaving differently from what was planned. Our results confirm the applicability of DETO to systems with nonlinear material behaviors, and prompt a discussion on service performance, strength loss, and robustness against incorrect material assumption, all of which are important for structural design.

Methodology
The DETO method has been presented in Part 1 and compared with classical SIMP-based TO using Finite Element analyses [8]. A MATLAB implementation is available on GitHub [22], but the repository also includes a C++ version which is more efficient and recommended for replicating the simulations in the present manuscript. Here we first present DETO for a problem of interaction energy maximization, which leads to maximum ductility in a system of particles under imposed displacements and interacting via a generic potential. The problem is then particularised into one of interaction energy maximization in a 2D structure made of closely packed, pairwise-interacting disks. Finally we present several interaction potentials to describe material behaviors that will later be used in Sect. 3.

DETO: Energy maximization with generic interactions
Consider a system of N interacting 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 i has an associated variable i ∈ [0, 1] . Particles with = 0 interact with zero intensity with the others, effectively representing voids. Particles with i = 1 interact with full intensity, thus representing full solid. All the per-particle i are gathered into a vector . An optimization problem to maximize the total interaction energy in the system, U tot , implies the possibility to redistribute the values in under a set of constraints, for example: c is the objective function to maximize, here identified with U tot . Equations 2 and 3 are the constraints on . Equation 2 fixes the target solid fraction f of a final system where all particles have either i = 0 bounds the values of i between 1 and a small value min ; in principle one could use min = 0 to represent void, however small but finite values, such as 10 −3 , are frequently preferred to avoid issues with later parts of the optimization process, in particular the filtering step (see Part 1 for more discussion on this point [8]). The total interaction energy U tot in Eq. 1 is a generic function of the particle positions and orientations . Following the concept of penalisation, which is central to the SIMP approach, the dependence of U tot on should be set to favour 0-1 only structural solutions and to penalise gray solutions with intermediate i values.
The optimization problem in Eqs. 1-3 is quite generic but some notes on its scope and underlying assumptions are due. Discrete Element analyses typically include velocity-dependent dissipative terms [9]; here we do not consider them because the problem refers to static equilibrium conditions. Extensions to dynamic problems is possible but beyond the scope of this manuscript. Equation 1 also assumes that U tot is history-independent, with no irreversible processes such as bond breakages. Irreversible events can be included in DE analyses and they motivate in part the development of DETO. However, maximization of U tot may not be a meaningful problem when such irreversibilities are included and other quantities, such as local stress or number of broken bonds for example, may be targeted instead. Similar and other issues may arise depending on the functional form of U tot . Therefore the problem in Eqs. 1-3 may not be appropriate, or might require additional constraints, depending on: the functional form of U tot , the possible presence of irreversible processes, and the type of analyses to be conducted (static or dynamic). The calculations in this manuscript will target static equilibrium, without irreversibilities, and with expressions of U tot that ensure the applicability and relevance of the problem in Eqs. 1-3.
A slightly more specialized expression of U tot , but still quite generic and more convenient, may be: where i, j, and k are the indexes of the interacting particles. The first term describes pairwise interactions, the second three-body interactions, and the series expansion may continue to include more multibody terms. The terms may describe nearest neighbour interactions, such as contact forces or harmonic bonds, as well as nonlocal ones, such as long-range electrostatic interactions. This expansion is commonly used in materials modelling, e.g. in atomistic simulations [23]. The spatial dependence of the interactions is described by the U ij and U ijk functions; to simplify the notation they are assumed to depend only on particle positions , but the treatment of orientation-dependent terms would be analogous. The penalising functions depend on the values of the particles involved in the interaction. Specific forms of U and will be introduced later in the manuscript; now we continue presenting DETO for the still rather general form of U tot in Eq. 4.
To solve the optimization problem in Eqs. 1-3 we use the same approach as detailed in Part 1 [8]. At the generic optimization step, a DE analysis provides the static equilibrium configuration corresponding to the current values in under the imposed forces and constraints to motion. The updating scheme for between two subsequent optimization steps is: 2 is a parameter to improve convergence. is a parameter that rescales the predicted values of to respect the constraints in Eqs. 2 and 3, as well as another typical constraint on the maximum allowed change of i between two subsequent steps, viz (this maximum change is always set to 0.1 in this manuscript). See Part 1 for more discussion on this point [8]. dc d i is the sensitivity of the cost function with respect to the design variables. Since c = U tot here, the sensitivity would be the total derivative of U tot with respect to the design variables . In the Supplementary Information of Part 1 we showed how such total derivatives can be computed using a numerical perturbation method, which starts from the equilibrated system at the generic optimization step and requires an additional equilibration for each i after altering the value of this latter by a small quantity . In Part 1 we showed also how, for the structural problems considered there, the perturbation method returned full sensitivities that were very similar to approximate sensitivities obtained considering only the partial derivative of c with respect to . This is equivalent to assuming that the energy change caused by a variation of i is governed by the change in stiffness induced by the while keeping the deformed configuration fixed to its pre-perturbation state. In the present manuscript we adopt a similar approximate version of sensitivities, comparing them with full sensitivities from numerical perturbation in Sect. A of the Supplementary Information. Because the relationship between stiffness and is determined by the penalization scheme, the penalising functions are the only ones to be derived when computing approximate sensitivities, hence the expression of sensitivity that we employ here, using the definition of U tot in Eq. 4, is: Equation 6 shows that spatial complexity in U tot does not add complexity to the calculation of sensitivity.
Indeed, the derivatives in Eq. 6 apply only to the penalising functions , whereas the values of U ij and U ijk can be directly taken from the results of the DE analysis. Equations 5 and 6 come from the Optimality Criteria approach to the solution of the optimization problem, adapted from Ref. [24]. This method is particularly efficient for the energy maximization problem here; other methods may better suit other problems.
Often the updating scheme in Eq. 5 does not feature directly the sensitivity from Eq. 6, but rather a coarse-grained version of it: where n f is the number of particles within a distance r min from the center of particle i, including particle i too. W k = r min − r ik ≥ 0 is a weight function ensuring that particles closer to i contribute most to its coarsegrained sensitivity dc d i . The coarse-graining process in Eq. 7, known as filtering, is commonly used in Finite Element based TO to avoid the checkerboarding problem [25,26]. Part 1 shows that DETO does not suffer from checkerboarding, but also that filtering still improves the quality of the optimum solutions [8].

DETO: Interaction energy maximization in 2D
The problem in the previous section is now particularized for a 2D system that is relevant for all the case studies in this manuscript. The system is a rectangular design domain initially filled with nely layers of nelx particles each (actually, nelx for odd layers and nelx − 1 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 figure.
At first, all disks have same = f . The interaction energy in Eq. 4 is assumed to feature only a pairwise term that depends only on the interparticle distance r ij . For we use the penalisation scheme that we proposed in Part 1 [8]: This leads to: As discussed in Part 1, also here we take p = 2 . A restriction that we will respect on the spatial function U ij is that the corresponding interaction force F ij = − dU ij dr will be monotonically increasing with r ij . This means that we can only consider material nonlinearity where the stiffness of the material may decrease (strain-hardening) or increase (strain-stiffening) with strain, but always remains positive. Strainsoftening nonlinearity, where the stiffness becomes negative, will not be considered here as it may lead to strain localisation and non-uniqueness of the solutions; this would require changes to the problem formulation in Eqs. 1-3 and related solution method.

Ductility maximisation under imposed displacements
The solution approach in Eqs. 5 and 6 has been presented for the energy maximization problem in Eqs. 1-3, which leads to a maximization of structural ductility under imposed displacements [20,21]. Imposing displacements rather then forces is good practice when dealing with nonlinear materials, because for these imposed forces may cause mechanical instability and compromise the optimization process.
Running DE analyses with imposed displacements, as part of DETO, requires some dedicated arrangements.
Consider a generic optimization step: the current values in provide the penalisation functions ij for each pair of interacting particles. These ij are used in the DE analysis that returns the static equilibrium values of and thus of U ij under the imposed displacements. These values of U ij are then used to compute the sensitivity in Eq. 6. Under imposed forces, the DE analysis is simply an energy minimization process, carried out using e.g. damped dynamics algorithms such as the Sheppard algorithm [27]. Imposed displacements imp require additional care to control the strain rate during energy minimization. One possibility, adopted in this manuscript, is to start each DE analysis from the underformed configuration and to progressively increase the imposed displacements with a suitably small rate ̇ imp . Energy minimization should proceed while increasing the imposed displacements, until the target imp are reached. After that, energy minimization should continue until satisfying appropriate convergence criteria, such as small tolerance. This approach for imposing displacements is shown by the two inner loops in the scheme in: Fig. 2. The other parts of the scheme in the figure describe the overall optimization problem, which is the same as for the case of imposed forces. In this manuscript, the change in energy is averaged over 10 steps before comparing with etol. If suitably small, the value of ̇ imp does not impact the optimization process, because sensitivities are computed with respect to equilibrium configurations, at which point the pseudo-dynamic effect of imposing deformation is complete. In this manuscript, we will show one example of checking the suitability of the chosen ̇ imp . Other approaches to impose a desired displacement are possible, for example adding imp as constraints to the DE energy minimization problem that provides the equilibrium configuration at the generic TO step. The resulting DE energy minimization problem can then be solved with any of the methods used in constrained optimization [28]. Choosing a fast method to perform DE energy minimization can significantly change the numerical performance of DE-based optimization as compared to FE-based optimization, because DE energy minimization is the most computationally expensive part of the algorithm. In this manuscript we did not look for the best method to perform DE energy minimization; we settled instead with a method that was sufficiently fast for problems that we wanted to consider.

Interaction potentials
We consider four types of interaction potentials U ij , each representing a different material behavior. This manuscript uses only first-neighbour interactions, but inclusion of longer-range nonlocal interactions would not require any change to the methodology presented here. Table 1 shows the expressions of each potential, along with the corresponding interaction forces (positive when repulsive).
One can immediately notice how the interaction energy and force do not diverge in the r ij → 0 + limit; potentials that are commonly used in microstructural simulations, such as the Lennard-Jones potential, feature instead diverging energy and force in such limit. However, the interactions proposed here are meant for macroscopic systems experiencing strain levels limited to few percent. For such systems, typical interactions used in Discrete Element simulations do not diverge in the r ij → 0 + , e.g. Hertz contact forces or Hookean bonds. Figure 3a compares the U ij (r ij ) for the various materials and for a set of k 0 , a. and D parameters that are relevant for this manuscript. As expected from strain energies, all the U ij curves are zero in the undeformed state r ij = D , and positive elsewhere. Figure 3.b shows the F ij (r ij ) curves, from which the strain-hardening and strain-stiffening regimes can be appreciated. The F ij (r ij ) curves are proportional to the stress-strain behavior of the material, which can be quantitatively estimated assuming that D = 1 mm, that the box thickness in the third direction is t z = D , and that the contact area between two particles is one sixth of the lateral surface area of the disk, viz 1 6 D t z . Under these assumptions, the strain between particles in 10 −3 units is equal to the elongation in m in Fig. 3, whereas the maximum stress between particles in strain-hardening regimes, when |r ij − D| ≫ 0 , is capped to k 0 a 6 D 2 = 477 MPa (assuming k 0 = 100 kN/mm and a = 400 mm −1 as in Fig. 3); simulation results will later confirm this estimation. Figure 3 also shows how, for the materials proposed here, nonlinearity become important at approximately 0.1-0.2% strain: this is representative of various metals at the macroscale, for example steel.
The interactions in Fig. 3 capture material nonlinearity under strain. However, the potentials are all elastic, with same stress-strain responses upon loading and unloading. Irreversible deformations could  [14,20]). Such irreversibilities would impact the results if the DE analyses involved dynamic or cyclic loads, or if buckling instabilities or strain localization, e.g. due to material softening or fracture, led to stress relaxation in some parts of the structure. In this manuscript, material softening and fracture are not considered (only hardening and stiffening as per Fig. 3), buckling will not occur, and imposed loads or displacements will always induce monotonically increasing strain everywhere in the structure. Under these conditions, the reversible interactions in Fig. 3 are as representative of large-strain material behaviors as elastoplastic interactions would be. Therefore in this manuscript, for the sake of brevity, we will use the term "plastic flow" even if the material is not strictly plastic, because its behavior under quasi-statically and monotonically increasing strain is the same.

Simulation details
All the simulations in this manuscript are based on the parameters in Table 2. The values of D, k 0 , and a in the table are those returning the interactions in Fig. 3. The filtering length r min is used for the weight factors in Eq. 7. The domain thickness t z is only used to estimate the cross section between two interacting particles in D t z = D 2 , which has been used to relate the force-elongation curves in Fig. 3.b to the stress-strain behavior of the materials.
The DE-based energy minimization algorithm to find configurations at static equilibrium is based on Sheppard's damped dynamics [27]. The algorithm uses a time step dt, particle mass m, and caps the maximum particle displacement to d max at each Euler integration step; convergence is judged based on the energy tolerance etol, as already explained in Sect. 2.3. The DETO process is considered as converged when new i − old i < 0.004 for all particles. The stresses per particle ab,i (with ab = xx, xy, or yy) can be computed from the interaction forces using the virial approach [29] (see Part 1 for more details [8]): r a and F b indicate respectively the a and b component (x, y, or z) of the particle position and interaction force vectors. Here we use V i = V tot N , where V tot is the total volume of the design domain and N is the number of particles in the system. From these stresses we can then compute Von Mises deviatoric stresses per particle: having dropped the subscript i to simplify the notation.

Results
This section presents three optimization problems where material nonlinearity may significant impact the resulting topologies. Results are first obtained for the Lin and Weak interactions in Table 1, which are analogous to the linear elastic and elastoplastic materials considered in Ref. [20]. Two of the systems are 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 stiffening in compression (Weak-T), and asymmetric material hardening in compression and stiffening in tension (Weak-C) also tested using the asymmetric potentials Weak-C and Weak-T in Table 1. 3.1 Three-support system with imposed displacement from the top The system is shown in Fig. 4; it is analogous to one originally analysed in Ref. [20] using FEM-based TO. Figure 5a shows the topology resulting from DETO when the material is linear elastic.
Most of the structure gets concentrated into a central pillar, which provides the shortest and stiffest path to transfer the load from the point A down to the central support. The benefit of increasing the cross section of the pillar is limited by the size of the support, to the extent that for the target solid fraction used here, f = 0.3 , additional stiffness is gained by creating diagonal branches that reach for the lateral support, despite such branches are longer than the central pillar and thus contribute less efficiently to the overall stiffness. A similar result was obtained in Ref. [20] using FEM-based TO; in that work, however, the optimum structure did not feature the diagonal branches. This difference may be due to the difference between an FE-based description and our DE-based one. Another possible explanation lies in different optimization procedures, e.g. the  Table 1 for some of the parameters used in this manuscript: k 0 = 100 kN, D = 1 mm, a = 400 mm − 1 ; b Corresponding force-elongation curves from Table 1, which are proportional to the stress-strain behaviors of the materials  Fig. 4 Optimization problem for a beam on three-supports with imposed displacement at point A at midspan. The value of imp has been fine-tuned to obtain an appreciable impact of material nonlinearity different updating schemes for or parameters such as the maximum change of i allowed between subsequent optimization steps. We found that the lateral branches appear also when imposing much smaller displacements, which excludes that they result from geometric nonlinearity and the fact that DE analyses compute forces in the deformed configuration. In any case, additional simulations not presented here have shown that the overall stiffness changes only very slightly when the mass is all concentrated into the central pillar, rather than being partly distributed to the thin diagonal branches in Fig. 5a. Figure 5b shows the optimization result for the symmetrically nonlinear material. The limiting factor for U tot in this case is that some pairs of particles may reach the maximum asymptotic value of their interaction force (see Fig. 3), thus entering into the analogous of a plastic flow regime. This happens near the supports and under the plate applying the imposed displacement, as shown by the sharp diagonal fronts of large deviatoric stress in Fig. 5b. The magnitude of such stresses is consistent with the quantitative estimation in Sect. 2. Because of these mechanisms, the response of the system is controlled by the thinnest cross section across which the load is transferred. Therefore, if the thick central pillar in Fig. 5a was retained, all its mass in excess to its smallest cross section, viz the size of the support below it, would not contribute to the maximum U tot . Therefore, when the material nonlinearity is considered in the optimization process, the excess mass is removed from the central pillar and used to thicken the diagonal branches, exploiting as much additional area from the lateral supports as possible. The result is analogous to that obtained in Ref. [20]. Figure 5c shows the evolution of the objective function, U tot , during the optimization process. As expected, the weaker nonlinear material ends up with significantly lower U tot . The snapshots within the figure show how the systems in Fig. 5a and b appear after 8 optimization steps only. Both systems then feature thick diagonal branches, but with the key Fig. 5 Optimization results for the three-support system in Fig. 4 with target solid fraction f = 0.3 , assuming a linear elastic and b symmetrically strain-hardening material, as per Lin and Weak expressions in Table 1. The colors represent the intensity of local deviatoric von Mises stress; c evolution of the objective function U tot , viz the total strain energy of the systems during the optimization. The base case with inputs in Table 2 is compared with cases with no filtering and with smaller ̇u imp ; d evolution of force-displacement curves during the optimization difference that the Weak system is already clearly utilizing the branches (light color meaning intense von Mises stresses in them), whereas the Lin system is not utilizing them significantly (dark color meaning little stress). As a result, at this step during the optimization mass tends to move away from the branches in the Lin case, whereas it tends to move towards the branches in the Weak case. Figure 5d shows the force-displacement curves for the Lin and Weak materials, evolving during the optimization process. Clearly the final solutions are much stronger than the initial ones, where all particles had 3 . At small displacements the two systems feature similar stiffness, whereas the nonlinearity caused by the material in the Weak system becomes evident at larger u imp .
For both types of material, Fig. 5c and d compare results for three different cases: the base case with input data in Table 2, the base case but without filtering, and the base case but with a smaller loading rate ̇u imp = 10 −5 (instead of 10 −4 mm s −1 in the base case).
For the linear material all cases give identical result. For the nonlinear material, instead, the case without filtering reaches a less optimal solution with lower U tot , whereas the other two cases returns the same evolution of U tot . A close scrutiny of the force-displacement curves for the Weak system indicates that the curves for the base case are the highest, suggesting a more optimum outcome. However, when reaching the target u imp , the base systems continues to minimize its strain energy which causes a drop of force while u imp remains fixed at 0.4 mm. The case with lower loading rate features a lower curve but with no further relaxation at u imp = 0.4 mm. As a result, both the base case and the slower one attain the same final value of force, and thus of U tot , at u imp = 0.4 mm; this explains why their U tot are identical in Fig. 5c. By contrast, the force-displacement curves for the case without filtering are intermediate between the base and slower cases but, when u imp = 0.4 mm is reached, a significant further relaxation sees the force dropping below those of the other cases (see the thin vertical lines at u imp = 0.4 in Fig. 5d). This explains why U tot in Fig. 5c is smaller in this case than for the others.
In terms of geometry evolution, what limits the unfiltered case is that the system rapidly converges to a configuration with all i ≈ 0 or 1, getting effectively stuck into a local energy minimum. This hinders a full transfer of mass towards the lateral branches, hence a full exploitation of the supports. A similar effect of unfiltered systems getting stuck into sub-optimal local minima was already observed and discussed in Part 1 for systems with linear material [8]. Figure 6 shows results that are particularly relevant for structural design.The Lin from Weak series explores how the structure in Fig. 5b, optimized for a nonlinear material (viz for best performance approaching failure), behaves in the linear elastic range. The results show that the stiffness of the structure is lower than that in the Lin structure, which was originally optimized assuming a linear elastic material. The loss in stiffness is 11%, from a gradient of 53.5 kN/mm in the Lin case to 47.75 kN/mm in the Lin from Weak case. In the same figure, the Weak from Lin series explores how the structure in Fig. 5a, optimized for a linear material (viz for maximum stiffness in service conditions) behaves when approaching failure. The results show that the maximum force and the strain energy at u imp = 0.4 mm are both substantially smaller than in the Weak structure, which was originally optimized assuming nonlinear material. The maximum force and strain energy go from 7.83 and 2.34 kN mm for the Weak case, to 6.12 kN and 2.11 kN mm for the Weak from Lin case, decreasing by 22% and 10% respectively. A 15% loss in maximum force was obtained in Ref. [20] for a system with same geometry, but using FEM-based TO and elastoplastic material. An 11% loss of stiffness in service conditions is likely to be less problematic than a 22% loss of strength approaching failure. Therefore, the designer should use TO with linear elastic materials carefully and favour optimization using realistic material behaviors when addressing strength and structural failure.
3.2 Three-support system with mid support settlement The system in Fig. 7 has very similar geometry as the previous one in Fig. 4. The differences are that the displacement is imposed at the mid support instead of above the beam, and that the lateral supports are only half as wide as before. This problem was also originally addressed in Ref. [20], there using FEM-based TO with elastoplastic material. Figure 8 shows the optimum geometries obtained from DETO.
In all cases, the resisting mechanism is akin to that in the seminal work of Michell [3], where the central ties connect the settling plate to the compressed arch above, which transfers the load to the stable lateral supports. The linear elastic Lin case produces a structure that is very similar to that in Ref. [20], despite the already mentioned methodological differences. A material that is strain-hardening both in tension and in compression produces the Weak structure in Fig. 8b, with a flattening of the arch at is its top and with fewer thicker ties linking the settling mid support with the compressed arch. Another important detail is that the Weak structure concentrates more mass near the later supports, which are instead not fully utilized in the Lin case. An analogous tendency to fully exploit the supports has been already discussed in the previous section, and was also observed in Ref. [20] for this case study.
The structures in Fig. 5, in the previous section, were fully under compression when loaded, hence considering asymmetric materials in tension and compression was not useful then. Here instead, Fig. 8 shows how asymmetric material behaviors lead to different optimum structures. In particular, Fig. 8c shows that a material that is weak, viz strain-hardening, in compression and strong, viz strain-stiffening, in tension produces a structure with thin central ties under high stress, and a thicker compressed arch that fully utilizes the lateral supports. By contrast, in Fig. 8d, a material that is weak in tension and strong in compression creates thick central ties and a shallower and thinner compressed arch which utilizes the lateral supports only in part. In this latter case, the limiting factor is the size of the settling central support, which controls the maximum cross section in tension and thus the maximum force that the structure can carry. Figure 9 shows the force-displacement curves for the four systems in Fig. 8.
As expected, all curves start with the same gradient in the initial linear regime. The Weak-T system displays an initially increasing gradient, due to strain-stiffening in the compressed arch. At displacements over 0.5 mm, however, the strain-hardening behavior of the central ties takes over and plastic flow caps the maximum force. The Weak system features the smallest strength, but the Weak-C is only marginally better, as opposed to the significantly stronger Weak-T system. This happens because the strengthcontrolling element in the Weak and Weak-C system is the compressed arch. The Weak-C system can transfer a bit more mass from the central ties into the arch, but eventually the minimum cross sectional area of the arch is limited by the size of the lateral supports, Fig. 6 Force-displacement curves for the configuration in Fig. 5b assuming linear elastic material (Lin from Weak), and for the configuration in Fig. 5a assuming nonlinear material (Weak from Lin). The curves are compared with the base cases for linear and nonlinear materials already shown in Fig. 5.d (solid curves). All curves here were obtained using loading rate ̇u imp = 10 −5 mm s −1 Fig. 7 Optimization problem for a beam on three supports with imposed settlement of the mid support which both the Weak and Weak-C systems utilize in full or almost. In the Weak-T system, instead, strength is controlled by the central ties and therefore the system has more freedom to move mass away from the compressed arch and alter the overall geometry to maximize its strain energy.
The different optimum solutions in Fig. 8 raise the question of how much an incorrect assumption of material behavior in the TO process may affect the structural performance. As an example, consider a structure where the elements under compression are confined using FRP to induce strain-stiffening in a material that would otherwise be symmetrically strain-hardening. In our model, this means turning a Weak system into a Weak-T one. If optimized assuming Weak-T behavior, the geometry in Fig. 8d would be obtained. However, if the FRP system failed in the actual structure, the material behavior would go back to Weak , for which the optimum geometry would be that in Fig. 8b instead. This raises two questions: how much strength loss may be caused by an incorrect assumption of material behavior? Which of the four material behaviors considered here would produce the most robust structure, in case the material ends up behaving differently? The results in Fig. 10 address these questions.  Fig. 7, with solid fraction f = 0.25 and for different material behaviors as per Table 1: a linear elastic, b symmetrically strain-hardening, c strain-hardening in compression and strain-stiffening in ten-sion, and d strain-hardening in tension and stiffening in compression. The colors represent the intensity of the deviatoric von Mises stress Each subfigure in Fig. 10 shows how one of the optimized structures in Fig. 8 would behave for any of the four material types in this manuscript. A first take is that all four structures, irrespective of the material assumption underlying them, feature a similar force-displacement curve when the material behaves linearly (compare the black solid Lin curves across the four subfigures in Fig. 10). This means that, for the structural system considered here, stiffness is not sensitive to the geometric details and the risk of losing service performance due to an incorrect material assumption is low.
To address the question on strength loss from unexpected material behavior, consider Fig. 10d. Going back to our example with the FRP, a structure optimized assuming Weak-T material should resist a force of ca. 5 kN, if the material behaves as predicted. However, if the FRP system fails and the material ends up behaving as Weak , the maximum force drops to 1 kN, with an 80% strength loss that would likely entail collapse. An analogous loss of strength would occur for structures optimized assuming Lin or Weak-C materials, in Fig. 10a and c, albeit less pronounced in the latter case due to the aforementioned, similar resisting mechanisms in the Weak and Weak-C cases. The only case not involving strength loss is that of a structure optimized assuming Weak material, in Fig. 10b. At first sight, this may be simply reduced to a "design for the worst-case scenario" message. However, designing for the worst case is a way to define suitably large cross sections for the various structural elements. Here the problem is different, as optimization with fixed f implies that any increase in cross section at one place requires a reduction of cross section elsewhere. Under this constraint, it is a nontrivial finding that the geometry optimized assuming Weak material gives the most robust structure with respect to other possible material behaviors.

Doubly fixed beam
In the previous section, the load was transferred to the lateral supports via a serial arrangement of ties working in tension, followed by the arch working in compression. In this section, a problem is devised to obtain elements in tension working in parallel with elements in compression. The system in Fig. 11 is proposed to this end; it is similar, but not identical, to the system that was used in Part 1 to highlight the impact of geometric nonlinearity [8]. Figure 12 shows the optimum geometries obtained with different assumptions on material behavior. Despite the symmetry of the system in Fig. 11, the structure optimized assuming Lin material is asymmetric with respect to the horizontal axis: see Fig. 12a. The asymmetry stems from geometric nonlinearity, which generates additional tensile stresses and thus favors concentration of material in the lower half of the structure (see Part 1 for more discussion on this point [8]). Figure 12b shows the optimum structure for a symmetrically strain-hardening material, Weak . The material nonlinarity enhances the asymmetry caused by the geometric nonlinearity, while mass is more concentrated in the main compressed arch and lower deck in tension, removing some of the diagonal struts that were present in Fig. 12a. Figure 12c shows the optimum structure for a material that is weak in compression only, Weak-C . This case features further concentrates mass in the lower deck, which is now fully exploited in tension, whereas the weaker compressed arch is significantly reduced in size. An almost specular geometry, except for a slight asymmetry due to geometric nonlinariy, is obtained for the Weak-T material, as shown in Fig. 12d. Figure 13 shows the force-displacement curves for the four structures in Fig. 12. The curve for the Weak case shows that the imposed displacement of 0.6 mm is triggering significant nonlinearity. Indeed, an upper bound for the strain in the structure can be estimated in , considering the diagonal struts and ties (assumed at 45 • ) in the optimum structures immediately below and above the center of the beam, and assuming that the very top and bottom rows of particles do not move vertically at all. This leads to an upper bound strain of 2.1%, which is indeed well in the nonlinear regime as per Fig. 3, while still far from strain levels that would require consideration of diverging energy and force upon strong compression. Fig. 11 Optimization problem for a beam partly fixed on both ends and with imposed settlement at the mid point Fig. 12 Optimum geometries for the problem in Fig. 11, with solid fraction f = 0.4 and for different material behaviors as per Table 1: a linear elastic, b symmetrically strain-hardening, c strain-hardening in compression and strain-stiffening in com-pression, and d strain-hardening in tension and stiffening in compression. The colors represent the intensity of the deviatoric von Mises stress The results in Fig. 13 agree conceptually with those in the previous section, with all materials providing similar stiffness at small deformations, and with significant differences emerging at larger u imp . As expected, the Weak material results in the lowest strength. The Weak-C and Weak-T materials lead to very similar force-displacement curves, which well reflect their almost specular geometries, combined with their specular material behaviors (see Fig. 3 in Sect. 2). Both structures with Weak-C and Weak-T materials overshoot the Lin curve at u imp < 0.9 mm; this is due to the strain-stiffening behavior of the Weak-C and Weak-T materials respectively in tension and in compression, which is eventually overtaken by strain-hardening in compression and tension. Figure 14 explores how robust the structure in Fig. 12 are with respect to wrong assumptions of material behavior. The results in Fig. 14 corroborate those in Fig. 10 in the previous section. Namely, all structures feature a similar stiffness, meaning comparable performance in service conditions. By contrast, strength is sensitive to material behavior and geometry. Out of the structures considered here, only the structure assuming Weak material preserves a similar strength if the material ends up behaving differently: see Fig. 14.b. Instead, structures optimized assuming Weak-C or Weak-T materials, in Fig. 14c and d, would end up with as little as half their design strength if the material turns out to feature a different type of nonlinearity. This means that, also for the parallel tension-compression system considered here, assuming the weakest material behavior for the optimization leads to the structure that is most robust against other unexpected material behaviors approaching failure.

Conclusion
The manuscript has presented a first application of DE-based TO to materials with nonlinear stress-strain behavior. In the DEM, the material behavior is captured by the interaction potentials between particles. The optimization problem under imposed displacement has been formulated for complex interaction potentials, and then particularized for potentials that model linear, strain-hardening, and strain-stiffening materials under monotonically increasing strain.
Four different material behaviors have been analysed, exploring their impact on three structural problems. General findings are that: • Structural optimization with linear or strainstiffening materials may lead to a partial utilization of the supports, whereas full utilization is to be sought when materials are expected to enter a strain-hardening regime near the supports themselves. This confirms previous literature results from FEM-based TO with elastoplastic materials [20]. • Structural optimization with the frequently used assumption of linear material behavior, is likely to generate significantly sub-optimum structures toward failure, when material nonlinearity takes over. On the other hand, structures optimized Fig. 13 Preliminary force displacment curves for the doubly fixed beam case assuming nonlinear material have been shown to produce similar stiffness in service conditions as structures that have been optimized assuming linear material. Therefore the safest approach in design is to optimize structures for nonlinear material behaviours and large deformations. • If the material ends up behaving differently from what is assumed during TO, the structure may end up being significantly weaker than predicted in design, causing significant risk of failure. Out of the material behaviors considered here, for both structures with serial and parallel tension-compression load paths, assuming the weakest material for TO has produced the most robust structural geometry against unexpected material behaviors.
Overall, this manuscript has shown how the newly proposed DETO method can incorporate material nonlinearity, and has provided some design-relevant insights from several case studies. So far, material nonlinarity has only been considered via reversible interactions between particles. Our results are therefore representative of structures that approach failure in a quasi-static regime of monotonically increasing strain. Despite this limitation, the inclusion of material nonlinearity is an important contribution towards future incorporation of irreversible behaviors, including for example inter-particle collisions and bond failure. This outlines a path where DETO will eventually enable the optimization of systems undergoing discontinuous and post-failure processes, which are challenging the current optimization methods based on the FEM or on other continuum descriptions.