Maximized response by structural optimization of soft elastic composite systems

Abstract Soft actuators triggered in a wire—and contactless way advance soft robotics, for instance, concerning microsurgical perspectives. For optimal performance in this and other contexts, maximized stimuli-responsiveness is frequently desirable. We demonstrate on the example of soft magnetoelastic systems how analytical theoretical measures in combination with computer simulations provide tools to develop optimized components. To enhance the overall macroscopic response, we adjust microstructural properties. Our strategy guides us towards ideally structured soft materials that can be fabricated using modern technologies.

To optimize the deformational or mechanical response of our magnetoelastic actuator systems, we utilize simulated annealing (SA) with adaptive cooling rates.For this purpose, the algorithm in Ref. 1 is adapted to our situation, see also the main text.
At the beginning of each process of optimization, we initialize the positions of the inclusions at random, unless indicated otherwise.One inclusion is randomly selected at a time.To generate a new trial configuration, its position is modified by a displacement randomly chosen in the range of [−δ/2, δ/2) for each direction.This marks one Monte Carlo (MC) step [2], where we measure lengths in units of the radius R of the elastic sphere when optimizing the magnetostrictive effect or in units of side length a of the elastic cube when addressing the magnetorheological (MR) effect.The positional constraints mentioned in the main text are obeyed, that is, minimal distances between the inclusions and from the boundary of the elastic sphere/cube are ensured.Every n jump MC steps, we instead generate a completely new trial position for the selected inclusion, still maintaining the positional constraints.Individual positional changes favorable to the optimization are always accepted.Those implying an opposite trend are only accepted with a certain probability, namely e −∆u/T , if ∆u > 0. Here, T is a temperature-like variable and ∆u marks the difference between the new trial configuration and the previous configuration, where u = ±u ⊥ 20 /R or u = ±u ⊥ 00 /R in our case, depending on the deformational mode to be minimized or maximized in the magnetostrictive case.For the optimization of the MR effect, we instead insert u = ±∆µ rel .Since u is dimensionless by definition, T is also introduced in a dimensionless way.Afterwards, we move on to the next MC step, testing for a new trial configuration.
As our optimization is based on SA, we decrease the temperature T during the course of optimization [3].Our "cooling schedule" starts from an initial temperature T = T i .We equilibrate the system for a period of N eq MC steps at the current temperature T .Afterwards, we use N prod steps to sample a specific heat-like variable where u 2 and ⟨u⟩ mark the mean squared and mean values of u, see above.Still, the temperature remains fixed during these N prod MC steps.During the next N cool MC steps of our optimization, we reduce the temperature each time before generating a new trial configuration.We decrease the temperature geometrically, that is we set T = (1 − k)T old , where T old denotes temperature in the previous MC step and k determines the cooling rate.
In our algorithm, the cooling rate for this period is set depending on the specific heat that was evaluated in the previous N prod steps: Here, we introduced the notations k s and k f for slow and fast cooling rates (k s ≤ k f ), respectively, and C * V as a characteristic threshold value of the specific heat.Consequently, our idea is to cool the system more slowly when more changes are occurring during the process of measuring the specific heat, thus allowing the system more opportunities to keep changing.Otherwise, the system is cooled faster to speed up the optimization.Together, these N eq + N prod + N cool MC steps denote one period of optimization.The same steps are repeated during the next period.Finally, our routine is stopped when no positional changes have been accepted for a number of patience periods (during sampling), or when the temperature has decreased to a predetermined final value of T f .
To summarize, our scheme makes use of the parameters δ, n jump , T i , N eq , N prod , N cool , k s , k f , C * V , patience, and T f .All parameter settings that were utilized for the optimization, leading to the structures depicted in Figs.2A-C -7 of the main text, we usually keep most of the parameter settings as in Tab.I.It is, however, necessary to rescale the temperature-like variable in the exponents e −∆u/T , because the typical magnitudes of deformation entering ∆u change with ν and N .For simplicity, we rescale initial and final temperatures equally.Therefore, we define T i (ν, N ) = T i α(ν, N ) and , where α(ν, N ) takes care of the required rescaling.We set α(0.5, 1000) = 1 when addressing ±u ⊥ 20 , α(0.3, 1000) = 1 when addressing ±u ⊥ 00 , and α(0.5, 200) = 1 when addressing the MR effect.The employed values of α and all other adjusted parameter settings, if different from the values shown in Tab.I, when varying ν and N are listed in Tabs.II-VIII.For elevated numbers of inclusions N , it is necessary to adjust the initial input configuration to fit all inclusions into the sphere when optimizing the magnetostrictive effect.In that case, we use a face-centered cubic (fcc) initial configuration of hexagonal layers oriented perpendicular to the magnetization direction in an ABC stacking configuration.One inclusion is positioned at the center of the elastic sphere.From there, we construct all layers.All nearest-neighbor distances in the fcc configuration approximate the minimal allowed distance.Then, we randomly delete inclusions from the full fcc configuration until the imposed number of inclusions N is reached.For optimizations of the MR effect, we also start with a simple stacking of hexagonal layers of inclusions oriented perpendicular to the magnetization direction.The bottom layer is placed at a distance of 0.061a from the lower boundary of the cube and one inclusion is positioned at the center of the layer, with inclusions in different layers placed above each other (denoted by hexagonal in Tab.VIII).Besides, we also utilize fcc arrangements as starting configurations (denoted by fcc in Tab.VIII), similarly to the simple stacking of hexagonal layers, however, with fcc featuring an ABC stacking.(As a remark we add that we never observed a decrease in performance after applying our scheme of optimization when compared to the performance of the initial configuration.)Furthermore, all optimizations are performed starting from 8 different initializations of the employed random number generator Mersenne Twister [4], leading to comparable results.Deviations in the results are at most of the order of the symbol sizes in the presented figures.Among these, we always selected those configurations of most competitive performance when plotting Figs.2-7 of the main text.

II. DETAILS ON THE REGULAR LATTICES TO WHICH WE COMPARE OUR OPTIMIZED STRUCTURES
In the main text, we compared the performance of our optimized configurations to those resulting from more regular arrangements.The latter were constructed artificially, guided by the optimized textures.In the following, we include details on the employed regular structures, first for the studies on magnetostrictive effects, then for the studies on magnetorheological (MR) effects.

A. Regular lattice structures for comparison with optimized magnetostrictive effects
We begin by minimizing u ⊥ 20 .The performance of our optimized structures in Figs.2D,E of the main text is compared to results obtained for more regular textures.Specifically, these were regular structures consisting of hexagonally arranged chain-like aggregates.To construct them, we start from layers of hexagonally positioned inclusions.Their normal vectors coincide with the magnetization direction (z-direction).To maximize the magnetic repulsion within each hexagonal layer, we approach the minimal allowed inclusion distance by setting the lattice constant to 0.121R.The distance between neighboring layers is set to an identical value of 0.121R.Inclusions within nearest-neighboring layers are located above and below each other, so that they form chains.For simplicity, we anchor the lattices by placing one first inclusion at the center of the sphere and by building up the lattice from there.For a regular arrangement of N = 1001 inclusions (to be compared to the optimized structure for N = 1000), we confine ourselves to the 11 innermost layers.Moreover, we only consider inclusions that are located within a distance of ρ = 0.62R from that center axis of the sphere oriented parallel to the magnetization direction.The resulting deformations characterized by u ⊥ 20 for this configuration and for different values of the Poisson ratio ν are displayed in Fig. 2D of the main text.
To accommodate increasing numbers of inclusions N x ŷ ẑ m FIG.S1.To best approximate the deformational response of the optimized structure for maximal relative elongation u ⊥ 20 , we identified a simple cubic lattice.Here, we illustrate its orientation.The normal vector of the plane indicated in red coincides with the magnetization direction m (by construction).
we increase ρ.Moreover, we slightly shift the regular arrangements of inclusions homogeneously along the magnetization direction by a distance ∆z.We observe that the optimized structures do not feature a layer exactly in the middle of the sphere.Indeed, through numerical trials with different values of ∆z the performance was increased further.Appropriate numbers n + and n − of layers above and below the center of the sphere, respectively, are determined by numerical trials as well to further optimize the deformational response.Next, when maximizing u ⊥ 20 (see Fig. 3 of the main text), we identified a simple cubic (sc) lattice to best approximate the optimized structures.As in the previous case, a lattice constant of 0.121R is employed.The orientation of this regular structure provides additional degrees of freedom to improve its performance.We start from an sc lattice for which the lattice vectors coincide with the coordinate axes of our Cartesian coordinate system, the z-direction again indicating the magnetization direction.Then, a global rotation is applied to this lattice.We found a rotation by a rotation matrix to approximately lead to best performance.It sets a rotation by a rotation angle of about 56.6 • around the

axis (
The resulting orientation is illustrated in Fig. S1. We again place the first inclusion at the center of the sphere.Furthermore, to set the number of inclusions, we only include lattice positions within a distance r from the center.We use r/R = 0.747, 0.81, 0.82, 0.83, 0.85, 0.86, 0.88, 0.9, 0.94 to obtain N = 1021, 1237, 1309, 1357, 1419, 1503, 1575, 1743, 1935, respectively.The resulting performance for u ⊥ 20 is depicted in Figs.3D,E of the main text.In the latter figure, we also display results for a facecentered cubic (fcc) configuration of r/R = 0.84.There, the edges of the cubic unit cells are oriented along the axes of our Cartesian coordinate system, leading to N = 1985.In that case, we likewise use a nearest-neighbor distance of 0.121R.
Third, we optimize u ⊥ 00 , starting with its minimization, see Fig. 4C of the main text.In that case, we use regular chain-like structures with a vertical shift of neighboring chains by half the inter-inclusion distance along each chain, as was similarly observed for the optimized structures, see Fig. S2.The lateral arrangement of these chains in the optimized structure is mostly irregular.For our regular structures, we arrange the chains in a square lattice.The corresponding lattice constant is set to approximately 0.108R.We construct the configuration starting from a chain that runs through the center of the sphere.The center is located symmetrically between the two innermost inclusions on that chain.As in the previous case, we use a radial cut-off distance r from the center of the sphere to adjust the number of inclusions N .We set r/R = 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.84, 0.878, 0.89, 0.907, 0. Next, we compare the results for textures of maximized u ⊥ 00 , see Fig. 5C of the main text, to the performance of regular configurations of inclusions.To this end, we arrange the inclusions in 11 layers of orientation normal to the magnetization direction.Within each layer, the inclusions are hexagonally organized, with a distance between nearest neighbors of 0.121R.The textures are constructed starting from one inclusion at the center of the sphere.In this case, however, the inter-layer spacing l is always chosen larger than 0.121R.Thus, nearest neighbors are found within the layers, not along the magnetization direction.By adjusting l and only keeping inclusions that are located within our constraints inside the sphere, we change the number of inclusions N .Setting l/R = 0.27, 0.22, 0.18, 0.155, 0.126 leads to N = 1021, 1239, 1517, 1745, 1973, respectively.

B. Regular lattice structures for comparison with optimized magnetorheological effects
Similarly, we construct by hand regular lattice structures, the performance of which we compare to our optimized ones that maximize the MR effect as described in the main text for an elastic cubical system.In general, the optimized structures guide us to the particular parameters used for our constructed regular lattice structures.In this way, the latter are adapted to the former.
For the positive MR effect, maximizing ∆µ rel for uniaxial elongation, see Figs. 6A,B of the main text, we utilize an fcc arrangement (ABC stacking of hexagonal layers) with a lattice constant of 0.121a, where a is the side length of the elastic cube.The bottommost hexagonal layer when viewed along the magnetization direction is always inserted at 0.1a from the bottom surface of the cube.The lattice is anchored by placing an inclusion at the center of the cube.A cut-off value ρ for the lateral distance of an inclusion from the center axis of the cube oriented parallel to the magnetization direction is set to tune the number of inclusions.Here, the value of ρ is chosen per layer, adapted from the optimized structures, generally decreasing the numbers of inclusions in the topmost and bottommost layers when compared to the other layers.For example, for N = 201 inclusions (for comparison to the optimized structures of N = 200), we use ρ = 0.2a for the topmost and bottommost layers and ρ = 0.35a for another 6 layers in between, resulting in a total of 8 layers.We also generate fcc configurations of 255, 300, 354, 399, 449, 497 inclusions, here always using 9 layers of varying ρ.The layer distance is chosen such that an fcc lattice of lattice constant 0.121a results.
When instead maximizing the positive MR effect for simple shear deformations, see Figs. 6C,D of the main text, textures consisting of hexagonally structured layers with their normals oriented along the direction of the shear displacements are introduced.The same lattice constant as above is used.This time, for each layer an inclusion is placed on that center axis of the cube that is oriented along the direction of shear displacements.The number of layers and locations of layers are adapted from the optimized configurations.We alternately select opposite sides from one layer to the next with respect to the central shear plane.Approximately, our procedure consists of cutting inclusions from these sides, starting from maximal distance from the central shear plane.We observe such asymmetric features in the optimized structures, see Fig. 6D of the main text.This procedure leads to N = 187, 291, 347, 402, 448, 499 inclusions organized in 4, 5, 6, 6, 7, 8 layers, respectively.In the last two cases, we use an alternating AB stacking of the layers to fit all inclusions into the cube.
Turning to the less frequently considered negative MR effect, we first consider the case of uniaxial elongation of the cube in the main text, see Figs. 7A,B.Neighboring chains forming when optimizing the structures are vertically shifted with respect to each other, similarly to what is displayed in Fig. S2.We generally use a square lattice arrangement of these chains for simplicity, as its performance shows good agreement in the resulting negative MR effect.For comparison, for N = 180 we also arrange 24 chains on a circle of radius 0.42a, showing a magnitude of the negative MR effect in line with the trend of corresponding square lattice arrangements of the chains.We shift the inclusions in neighboring chains vertically by half of the nearest-neighbor distance in a single chain (0.121a).Their lateral distance is chosen such that neighboring inclusions from different chains have the same total distance from each other of 0.121a.Therefore, each inner inclusion on a chain has two nearest neighbors on the same chain and two more per neighboring chain, for a total of up to 10 nearest neighbors.To tune the inclusion numbers of the square lattice arrangements of the chains, we selectively delete chains, leading to arrangements of N = 150, 203, 248, 301, 346, 398, 450, 503 with 20, 27, 33, 40, 46, 53, 60, 67 chains, respectively.The chains for deletion are selected by comparison to the optimized structures, which are displayed as the insets in Fig. 7B of the main text, approximately such that the square lattice of chains has a gap in the middle.Results for these square lattice arrangements as well as the one circular arrangement of the chains mentioned above for N = 180 are presented in Fig. 7B of the main text (chain-like config. ).The data point for N = 180 of circular arrangement of the chains is highlighted by the black circular mark.
Lastly, we turn to the negative MR effect under simple shear deformations, see Figs. 7C,D of the main text.The resulting optimized structures in that case are similar to the ones for the positive MR effect under simple shear.However, the layers here are oriented parallel to the plane of simple shear.Furthermore, as can be inferred from Fig. 7C of the main text, the layers do not consist of regular hexagons anymore.Instead, they resemble a centered rectangular lattice, with the longer edge along the magnetization direction.We choose 0.146a for the edge length along the direction of shear displacements and 0.194a for the edge length along the magnetization direction.Again, we set the numbers and locations of the layers as inferred from the optimized structures.To restrict the number of inclusions in each layer, we use a cut-off value ρ for the lateral distances of the inclusions from that center axis that is oriented perpendicular to the shear plane, similar to the maximization of ∆µ rel for uniaxial elongation (in contrast to the rather diagonal reduction in Fig. 7C of the main text).Here, ρ is cho-sen for each layer again by comparison to the optimized structures.For example, we observe that the optimized structures for N = 200 inclusions feature a center layer (see the leftmost inset of Fig. 7D of the main text) with fewer inclusions compared to the other layers.Therefore, we choose a smaller value for ρ for that layer.In this way, we obtain N = 203, 248, 297, 347, 399, 450, 497 with 5, 6, 7, 7, 7, 8, 9 layers, respectively.For N = 497, we use an alternating AB stacking of the layers.
, 3A-C of the main text for N = 1000 and ν = 0.5, Figs.4A,B, 5A,B of the main text for N = 1000 and ν = 0.3, Figs.6A,C and 7A,C of the main text for N = 200, are listed in Tab.I.When turning to other values of the Poisson ratio ν or the number of inclusions N , see subfigures B, C, D, or E of Figs. 2

FIG. S2 .
FIG. S2.Configurations that maximize the decrease in volume, minimizing u ⊥ 00 , see Fig.4Aof the main text, contain chain-like elements along the magnetization direction.Here, we illustrate the typical nearest-neighbor structure of these chains, where inclusions in neighboring chains are vertically shifted relative to each other by about one half (vertical) nearest-neighbor distance.

TABLE I .
Parameter settings in our SA algorithm to minimize or maximize u ⊥ 20 , u ⊥ 00 , or ∆µ rel , see the first column.To minimize u ⊥ 00 , the parameter C * V is arbitrary because we set ks = k f , see Eq. (S2).To maximize u ⊥ 00 , the parameters in the third row from the bottom are used.Yet, when varying N (see Tab. VII) we use the alternative parameter settings (alt.) of the second row from the bottom.

TABLE II .
Parameter settings in our SA algorithm when minimizing u ⊥ 20 for different values of the Poisson ratio ν, leading to the results in Fig. 2D of the main text.

TABLE IV .
Parameter settings in our SA algorithm when maximizing u ⊥ 20 for different values of the Poisson ratio ν, leading to the results in Fig.3Dof the main text.

TABLE VII .
Parameter settings in our SA algorithm when maximizing u ⊥ 00 for different inclusion numbers N , leading to the results in Fig.5Cof the main text.

TABLE VIII .
Parameter settings in our SA algorithm when maximizing the magnitude of the MR effect ∆µ rel for different inclusion numbers N , leading to the results in Figs.6B,D and Figs.7B,D of the main text.