Evolution of adaptive force chains in reconfigurable granular metamaterials

Abstract

Under an externally applied load, granular packings form force chains that depend on the contact network and moduli of the grains. In this work, we investigate packings of variable modulus (VM) particles, where we can direct force chains by changing the Young's modulus of individual particles within the packing on demand. Each VM particle is made of a silicone shell that encapsulates a core made of a low-melting-point metallic alloy (Field's metal). By sending an electric current through a co-located copper heater, the Field's metal internal to each particle can be melted via Joule heating, which softens the particle. As the particle cools to room temperature, the alloy solidifies and the particle recovers its original modulus. To optimize the mechanical response of granular packings containing both soft and stiff particles, we employ an evolutionary algorithm coupled with discrete element method simulations to predict the patterns of particle moduli that will yield specific force outputs on the assembly boundaries. The predicted patterns of particle moduli from the simulations were realized in experiments using quasi-2D assemblies of VM particles and the force outputs on the assembly boundaries were measured using photoelastic techniques. These studies represent a step towards making robotic granular metamaterials that can dynamically adapt their mechanical properties in response to different environmental conditions or perform specific tasks on demand.

---

1 Introduction

Mechanical metamaterials process mechanical inputs, such as forces, pressures, or waves, to achieve programmable stiffness, shape transformations, and force propagation.^1–3 Many current mechanical metamaterial approaches employ continuum solids or linkages/mechanisms with a fixed structure and therefore demonstrate only fixed responses. We are interested in mechanical metamaterials that can exhibit increased dynamic plasticity, enabling adaptation to different environmental inputs or task demands by reconfiguring their physical structure. Granular metamaterials—consisting of discrete, macroscale particles—offer an advantageous platform for such dynamic programmability, as individual particle properties can be tuned to achieve different bulk responses.⁴ For example, one can imagine a future granular metamaterial that uses changes in individual particle properties to route forces around a damaged region of the material thereby maintaining its function.

Existing granular metamaterials most often incorporate ordered particle assemblies, as it is difficult to predict the material properties of disordered granular assemblies. Ordered granular metamaterials have been used for vibration mitigation,^5–9 acoustic switches,^10,11 energy absorption,¹² elastic modulus and density tuning,^13,14 and non-reciprocal behaviors.¹⁵ Most of these prior studies have confined their approach to granular systems with inert, rigid grains and fixed elastic moduli. In this work, we expand the parameter space by including variable modulus particles, thereby enabling access to a wider set of possible material responses.

In previous studies, we developed discrete element method (DEM) simulations to study the vibrational response of jammed packings composed of binary mixtures of spherical particles with different masses. For example, in one recent study,¹¹ we showed that 2D and 3D mixed-mass granular assemblies can be used to either transmit or block vibrations with particular frequencies. Expanding on that result, we coupled evolutionary algorithms (EA) with DEM simulations to create more complex logic gates.^16–18 Using similar granular systems, in this work we seek to control the transmission of force chains, in both simulations and experiments.

In general, granular packings are collections of macroscopic particles that are densely packed, forming contact networks through which interparticle forces are transmitted. It is well-known that a large fraction of the external load applied to granular packings is often carried by a small number of particles,^19–24 resulting in the emergence of stable force chains.²⁵ Due to the large number of degrees of freedom, nonlinear response, and out-of-equilibrium behavior of granular systems, force chain dynamics are challenging to predict analytically. However, there have been many computational studies of force chain networks, such as contact dynamics simulations,²⁶ DEM simulations,^27,28 random matrix theory, and statistical mechanics-based models,²⁹ as well as experimental studies of force distributions.^30–32 Previous studies have emphasized that individual particle properties, such as the elastic moduli and density, influence how forces propagate through granular assemblies.³³ We will exploit this relationship by switching the elastic moduli of individual particles to achieve specific force outputs from selected particles.

We realize granular metamaterials with adaptable force chains using an assembly of variable modulus (VM) particles. While there are many potential approaches to particle modulus modulation, we adapt the approach previously used by Pashine, et al.³⁴ to make variable stiffness bonds in an allosteric metamaterial. Our VM particles are fabricated by incorporating a low-melting-point alloy (Field's metal) in soft silicone shells. A VM particle possesses a larger modulus when the Field's metal core is solid (at low temperature) and a smaller modulus when the core is liquid (at high temperature).

Using discrete element method (DEM) simulations combined with evolutionary algorithms, we identify particle configurations that optimize specific force outputs on the boundaries of a granular assembly. These optimized configurations are then fabricated and tested in experiments, where we demonstrate that the force chain networks within the packing can be dynamically reconfigured by modulating the modulus of individual particles. To explore practical trade-offs, we implement a multi-objective optimization pipeline that simultaneously maximizes boundary force outputs and minimizes power consumption, based on the assumption that power is expended when switching particle stiffness between soft and stiff states. Across all cases studied, we find strong agreement between the simulation predictions and experimental measurements of output forces. The novelty of this work lies in the integration of three components: experiments with variable modulus particles, DEM simulations that quantitatively capture stiffness-dependent interactions, and evolutionary algorithms that discover particle arrangements producing optimal force networks consistent with experimental results. Taken together, these results represent a step toward dynamic granular metamaterials capable of adapting their mechanical responses to evolving environmental or task-specific demands.

2 Methods

To controllably modify the force chain network of a packing, we create a pipeline that takes a set of objectives (e.g., maximize the force between a specific particle and the boundary) and hardware constraints as inputs, and outputs the configuration of large- and small-modulus particles that best achieves the objectives. Our pipeline has several steps. First, we characterize the contact mechanics and interparticle force law of the VM particles. Second, we develop DEM simulations to represent the physical properties of the granular packings. Finally, we combine DEM simulations and evolutionary algorithms to identify the grain configurations that achieve the objectives (Fig. 1).

Fig. 1 Description of the pipeline to generate granular packings with specified force outputs. (A) We start with DEM simulations that can generate force distributions that mimic those of the VM particle packings in experiments. We apply an input force F_i to the top boundary of the packing and measure output forces exerted on the bottom boundary from particles (blue outlines) in the bottom row. We then define optimization objectives to target specific force outputs to satisfy possible hardware constraints. (B) Using multi-objective optimization, the evolutionary algorithm identifies a set of particle configurations that best satisfies the objectives. These solutions represent the Pareto front, where each solution represents tradeoffs between the objectives. The algorithm starts with a set of randomly generated configurations in the first generation (Gen 0). During the optimization process, particle configurations that better satisfy the objectives are replaced in the solution set until the stopping criteria are reached (Gen n). (C) An example optimal particle configuration from the Pareto front. Stiff and soft grains are shaded dark and light green, respectively. Black lines indicate inter-particle contacts with the line thickness proportional to the force magnitude. Output forces are maximized for the odd-numbered particles with blue outlines. (D) We can then realize the selected granular packing in experiments and compare the force outputs to those predicted by the DEM simulations. The color scale indicates the particle temperature T, which in turn determines the particle moduli (particles with T > 62 °C are soft and particles with T < 62 °C are stiff).

The inverse design problem—designing particle configurations to match a desired force chain output—is an arduous task. Without a systematic methodology, we would have to evaluate an exponentially large number of configurations to decide which configuration gives the desired force output. While it is theoretically feasible to find packings for complex objectives analytically, this task is nontrivial given our nonlinear force model and inability to control the interparticle contacts that occur in granular packings. Therefore, we use DEM simulations coupled with an evolutionary algorithm to identify the optimal grain configurations. In particular, we implement a multi-objective optimization algorithm to search for solutions that satisfy the objectives, subject to experimental constraints. We can then make the grain configurations obtained from the optimization pipeline in experiments and evaluate their force networks, ensuring that the objectives have been met.

We focus on granular packings composed of monodisperse cylindrical particles (with diameter D) confined between two flat walls, such that their cylindrical axes are perpendicular to the walls, and arranged on a 2D triangular lattice. The experimental system is modeled in 2D as collection of monodipserse disks as shown in Fig. 1A. The packing consists of N_r = 5 rows, each with either 4 or 5 particles per row, yielding a total number of particles, N = 23. We index each particle p between 1 and N, increasing from left to right and from bottom to top in the packing, as shown in Fig. 1A. For a specific configuration denoted as C, we assign each grain with an elastic modulus chosen from one of two values: soft k_soft or stiff k_stiff with k_soft < k_stiff. We define the input force F_i as the total force applied downwards on the top boundary of the system (Fig. 1A) After applying F_i, we measure the output forces exerted on the boundary from the pth particle in the bottom row, F_o,p.

In the remainder of this section, we describe in detail the characterization of the Young's modulus of the particles and interparticle contact force law, the experimental setup, the DEM simulations, and the optimization approaches.

2.1 Experimental design

2.1.1 Variable modulus particles. We make the cores of the VM particles using Field's metal, an alloy of bismuth, indium, and tin with a melting temperature of T₀ = 62 °C. We encapsulate a solid cylinder of Field's metal inside a soft silicone shell made from Smooth-on Dragonskin™ 10 (Fig. 2A). We place a small copper heater inside the shell that allows us to melt the Field's metal via Joule heating. The high-resistance copper heater is separated from the Field's metal by a thin layer of silicone to avoid shorting the heater. By running a current (≈1 A) through the resistive heater we induce Joule heating in the particles. This heating melts the Field's metal inside at around 3.5 W, which changes its Young's modulus. While solid, Field's metal has a Young's modulus of 9.25 GPa,³⁵ which reduces to near zero at the melting temperature.

Fig. 2 Variable modulus (VM) particle characterization. (A) Illustration of the components of a VM particle (left) and an image of a fully assembled VM particle (right). Conductive wires (not shown) are connected to the leads of the copper heater. (B) Force F on VM particle from the two walls in the Instron plotted versus the symmetric displacement d of the walls (relative to the location of the walls at contact) in the soft (blue solid lines) and stiff (magenta solid lines) states. The dashed lines indicate fits to eqn (1) for the Hertzian contact law. Data for a particle with the same silicone geometry, but without Field's metal inside, is labeled as an “empty particle.” Shaded areas represent ±1 standard deviation from six trials. (C) Particle modulus k plotted as a function of temperature T. The melting temperature of the Field's metal is highlighted by the vertical dashed line. k for the “empty particle” is also indicated by the horizontal dashed line.

Next, we describe measurements of the Young's modulus of the VM particles. We compress single VM particles between two steel walls using a material testing system (Instron™ 3365). We find that the relationship between the symmetric displacement d of the walls (relative to their location at contact with the particle) and the compressive force F is best fit by assuming a Hertzian contact law:³⁶

whereE_s and E are the Young's moduli of steel and the VM particle, ν_s and ν give the Poisson ratios of steel and the VM particle,and R_s and R are the radii of the steel walls and the VM particle. We assume that E_s/E → ∞ and R_s/R → ∞. Note that eqn (1) gives the contact force law between two elastic spheres, not between two elastic cylinders with parallel cylindrical axes. The force law between the cylinders in the experiments can mimic the force law between contacting spheres when the cylinders have shape imperfections and the cylindrical axes are not parallel.^37,38 The force F versus displacement d is shown in Fig. 2B for the silicone shell alone, the VM particle when heated (soft state), and the VM particle at room temperature (stiff state). Fitting the force versus the displacement curve yields the VM particle modulus:

In Fig. 2C, we plot the dependence of k on the temperature T of the Field's metal. The particle modulus decreases rapidly when T > T₀. We find k_soft ∼ 0.34 MPa for the soft state (T > T₀) and k_stiff ∼ 1 MPa for the stiff state (T < T₀). Thus, the VM particle design possesses a three-fold change in modulus. k_soft is roughly 1.5 times the modulus of the empty silicone shell. For the empty silicone shell, we report k ≈ 0.2 MPa, which is below that for a solid silicone material (∼0.53 MPa). Also note that since the Field's metal becomes a liquid when heated above the melting temperature, an isolated Field's metal particle can in principle yield an infinite-fold modulus change between the liquid and solid states. The three-fold change for the VM particle modulus highlights the substantial contribution of the outer silicone shell to the VM particle modulus. Since the silicone shell encapsulates the metal core, compressive strains predominantly affect the silicone rather than the metal, reducing the modulus ratio between the stiff and soft states compared to the case where only the metal is compressed in the absence of the silicone shell.

2.1.2 Particle assembly. We assemble the VM particles into a triangular lattice where each particle can be actuated (heated) individually using a microcontroller (Arduino™). However, heat from actuated particles eventually dissipates into the neighboring non-actuated particles, causing thermal crosstalk. To mitigate thermal crosstalk, we established limits on the actuation power of each VS particle (see ESI,† S1).

We apply a uniform load F_i to the top layer of particles using the same materials testing system employed to characterize the moduli of the VM particles. The force propagates through the packing to create a particular force chain network that is a function of the positions of the particles and their individual moduli. Based on the locations of the stiff and soft particles, the load is transmitted differently to the bottom row of particles, exerting different amounts of force at various locations on the bottom boundary.

The assembly boundaries are made of a photoelastic material (ClearFlex™ 95) placed between two circular-crossed polarizers. We measure the output forces F_o,p on the bottom boundary by analyzing the stress-induced birefringence patterns.^23,39,40 When the boundary is not under stress, no light passes through the polarizers. Any stress on the boundary rotates the direction of the incoming electric field, allowing light to pass through the polarizer and appear as a fringe pattern. Because different wavelengths of light will form different patterns on the photoelastic boundary, we use a filter to only observe one wavelength of light (λ = 530 nm).

To measure F_o,p, we consider the intensity I_out that occurs through the application of a point force F_o,p on an elastic half-plane through a polariscope.^41,42 The light intensity I_out at any given point (x, y) within the half-plane is

where t is the thickness of the material, I₀ is the maximum intensity, K is the material-dependent stress-optic coefficient, and σ₁ − σ₂ is the principal stress difference, which in our system corresponds to the normal stress measured radially outward from the particle–wall contact. Therefore,where r is the distance from the point force to the point (x, y), θ is the angle of the vector pointing from the point force to (x, y) relative to the x-axis (parallel to the bottom boundary), and ϕ is the angle that determines the components of the force in the x–y plane (Fig. 3A).

Fig. 3 Assessment of the optimization of the VM particle assembly. (A) Schematic of the particle assembly. The optimization objective [scr O, script letter O]₀ is to maximize the force output F_o,3 of the middle particle in the bottom row given an input force applied to the top boundary F_i. (B) [scr O, script letter O]₀ for the best solution and averaged over the whole population, plotted as a function of the number of generations. The results are averaged over three independent trials with random initialization. (C) The probability distribution of [scr O, script letter O]₀ for 5000 randomly generated configurations with vertical lines that indicate [scr O, script letter O]₀ for the average (blue dashed line) and best solutions (red solid line) at G = 250. We also show [scr O, script letter O]₀ for two random configurations (random 1: cyan dotted line and random 2: magenta dashed-dotted line) pictured in panel D. (D) The first column shows three configurations with the predicted inter-particle contact networks from the DEM simulations. The soft and stiff particles are shown in light and dark shades of green, respectively. The best solution found by EA is shown in the first row, along with two random configurations (random 1 and random 2) on the second and third rows. The second column shows thermal images of the experimental packings with temperature increasing from dark green to white. The third column illustrates the birefringence patterns on the photoelastic boundaries. The black cords emanating from the bottom right of the packings provide electrical current to each particle. (E) F_o,3 plotted as a function of F_i obtained from experiments for the same configurations in (D). Shaded areas represent ±1 standard deviation from three independent trials. In the inset to (E), we plot the objective function from the experiments [scr O, script letter O]^exp₀ versus that in the the DEM simulations [scr O, script letter O]^sim₀, where the dashed line indicates [scr O, script letter O]^exp₀ = [scr O, script letter O]^sim₀.

Using eqn (5) and (6), we can construct an image for a given force on a half-plane. By comparing this image to the experimentally obtained image, we can fit for 2tKF_o,p/λ using minimum chi-square estimation.⁴¹ This fit yields F_o,p to within a proportionality factor, 2tK/λ. However, we can eliminate the proportionality constant by summing all of the forces on the bottom boundary: . See ESI,† S2 for more details.

2.2 DEM simulations

We employ discrete element method (DEM) simulations in two dimensions that mimic the experimental setup. We assemble a collection of monodisperse, frictionless disks into a triangular lattice (see Fig. 1A). We apply a constant downward force F_i/(k_softD²) = 0.01 to the top boundary. Each particle interacts with its neighbors and the walls, assuming purely repulsive Hertzian interactions, as in eqn (1). We integrate Newton's equations of motion for each disk with damping forces proportional to the particle velocities using a modified velocity-Verlet integration scheme until the net force on all disks . Determining which disks are in contact allows us to construct the network of interparticle contacts, as shown in Fig. 1C. In ESI,† S3, we show that F_o,p is proportional to F_i over a broad range of F_i. Hence, the specific choice of F_i in the DEM simulations does not impact the optimal solutions predicted by the evolutionary algorithm, as long as F_i is sufficiently small.

We note that there are some differences between the DEM simulations and the experiments. For example, the DEM simulations do not include particle–substrate and inter-particle friction. We mitigate the effects of friction by applying a thin layer of cornstarch over the surfaces of each VS particle in the experiments, thus decreasing the effective friction in the experiments. In addition, the particle sizes are slightly polydisperse (<2%) in the experiments, while the simulations assume monodisperse disks. Despite these differences, we find strong qualitative agreement between the results of the DEM simulations and the experiments.

2.3 Optimization setup

To solve the inverse design problem of finding particle configurations that achieve a desired output force pattern, we employ evolutionary algorithms (EAs) in our design pipeline (Fig. 1). EAs are a class of population-based gradient-free global optimization methods inspired by concepts from biology, genetics, and natural evolution.⁴³ In the general case, the algorithm starts with a set of randomly generated solutions (called a population). In each subsequent optimization step (i.e., generation), possible solutions are evaluated using a fitness function. The solutions with the highest fitness are chosen to survive to the next generation and produce alternate solutions through a process of mutation and selection that can potentially create solutions with higher fitness. EAs have been useful in myriad problems in materials science,⁴⁴ robotics,⁴⁵ and data retrieval applications.⁴⁶

In this study, to find an optimal configuration of soft and stiff particles to achieve a desired force output, we utilize a particular EA called age-fitness Pareto optimization (AFPO).^47,48 AFPO is a multi-objective evolutionary algorithm that balances the diversity of solutions with their fitness values, allowing for an extensive search in configuration space without converging to local optima in a rugged fitness landscape. When the objectives are not mutually exclusive and no explicit coupling is imposed during optimization, AFPO converges to a set of optimal solutions—commonly referred to as the Pareto front—rather than a single best solution. Note that AFPO can also be applied to single-objective optimization, in which case it does not produce a Pareto front, but instead converges to the single best solution. On the Pareto front, no solution is better than the others (and each solution optimizes the multiple objectives differently) unless additional trade-offs are specified.

To select a single best solution from the Pareto front, one can rank Pareto-front solutions with methods used for multiple criteria decision-making, such as the technique for order of preference by similarity to ideal solution (TOPSIS) method.⁴⁹ TOPSIS implements trade-offs between objectives that oppose each other by assigning a weight vector W that indicates how much each given objective should be weighted while traversing the Pareto front. TOPSIS then chooses the solution with the shortest distance from the ideal solution and the greatest distance from the negative ideal (worst) solution. Note that TOPSIS is applied to the Pareto front only after it has been identified by the AFPO algorithm.

The AFPO algorithm has several components including a genome representation, a variation operator, and a fitness function. In this study, we used a direct encoding scheme for the genome, where the genotype is a binary string (0 for k_soft and 1 for k_stiff) that represents the modulus values for all of the particles. The initial genotype is a string of N random bits, [scr U, script letter U](0,1). The variation operator is a bit-flip mutation with a given mutation probability. The bit-flip mutation changes the modulus of a particle from k_soft to k_stiff or vice versa. We consider both single- and multi-objective fitness functions. Additional details concerning the implementation of the AFPO algorithm are provided in Table 1.

Table 1 Details of the AFPO algorithm for the single- and multi-objective optimizations carried out in this study

Parameter	Value
Genome representation (C)	Binary string C = [x1, x2, …, xp, …, x23], xp ∈ {0,1}
Phenotype ([scr K, script letter K])	[scr K, script letter K] = [k1, k2, …, kp, …, k23],
	kp = xpksoft + (1 − xp)kstiff,
	kstiff/ksoft = 3.0
Variation operator	Bit-flip mutation
Mutation probability	0.05
Objective functions	Feasibility check	[scr O, script letter O]0 = Fo,3/Fi
	Case I
	Case II
	Case III
	Case IV
Population size (P)	Feasibility check	P = 50
	Cases I–IV	P = 100
Number of generations (G)	Feasibility check	G = 250
	Cases I–IV	G = 500
Initialization	{C^n}, n = 1…P
	For each n, [xp] = [scr U, script letter U](0,1)

---

We explore five optimization objectives as outlined in Table 1. In the feasibility check, we consider a single-objective optimization problem to maximize the force output from a single particle on the bottom row (Section 3.1). We then consider two other single-objective optimization problems to further test the efficacy of our optimization pipeline (cases I and II). In these two cases, the objective function is defined so that the force output from the odd- or even-numbered particles on the bottom boundary is maximized and has units of force before normalization (Section 3.2). In Secttion 3.3, we consider a multi-objective optimization problem (cases III and IV) where we can change from the optimal configuration in case I to that in case II (or vice versa) with an additional hardware constraint that limits the number of changes in the particle modulus as the second objective. The formulation of the objective functions and optimization results for each case are presented in the next section.

3 Results

3.1 Feasibility check: maximizing the output force from a single particle on the bottom row

To verify that the DEM simulations are sufficiently accurate and that the EA generates optimal solutions, we first study a simple case of maximizing the output force exerted by the middle particle in the bottom row on the boundary (F_o,3 in Fig. 3A). For this case, we define the objective function,and the goal is to find a configuration C* that maximizes [scr O, script letter O]₀. As shown in Fig. 3B, using EA, the optimal configuration emerges after ∼100 generations and has higher [scr O, script letter O]₀ than that generated using a random search with over 5000 configurations (Fig. 3C).

In Fig. 3D, we show that C* contains two lines of stiff particles stretching from the top corners to the central particle in the bottom row. This solution is intuitive, as larger forces will be carried by the stiffer particles. Thus, arranging stiff particles such that they form a connected path from the top to the bottom boundary without splitting will result in the largest force carried by the central bottom particle.

To validate the results obtained by the EA, we first construct the three configurations (best solution, random 1, and random 2) in the experiments (Fig. 3D) and measure F_o,3 as a function of F_i for each configuration. In Fig. 3E, we show that the output force F_o,3 is approximately linear with F_i for F_i > 4 N, which is consistent with the results from the DEM simulations in ESI,† S3. For F_i < 4 N, F_o,3 is nearly constant, which is caused by the low signal-to-noise values for the photoelastic material at small forces. We also note that the top plate provides comparable forces to each particle in the top row (e.g., ∼10 N/5 yields ∼2 N per particle) and these force magnitudes are similar to those used to measure the stiffness of the VM particles (∼1.5 N, Fig. 2B).

To de-emphasize the uncertainties of the photoelastic stress measurements at small F_o,3, we determine the objective function in the experiments, [scr O, script letter O]^exp₃, as the slope of F_o,3 versus F_i (namely ∂F_o,3/∂F_i) for F_i > 4 N, instead of evaluating F_o,3/F_i at a given F_i. Note that ∂F_o,3/∂F_i can be different from F_o,3/F_i, and F_o,3/F_i can vary with F_i in the linear regime of the experimental curve. This is because F_o,3(F_i) has a non-zero intercept at F_i = 0, presumably due to errors in the photoelastic measurements.

As shown in Fig. 3E, the best solution ([scr O, script letter O]^exp₀ ≈ 0.5) substantially outperforms the two random configurations ([scr O, script letter O]^exp₀ ≈ 0.32 and 0.13 for random 1 and random 2, respectively).

We observe quantitative agreement for [scr O, script letter O]₀ between the experiments and DEM simulations for the best solution and random 1 configurations: [scr O, script letter O]^sim₀ ≈ [scr O, script letter O]^exp₀ ≈ 0.46 for the best solution and [scr O, script letter O]^sim₀ ≈ [scr O, script letter O]^exp₀ ≈ 0.28 for random 1. However, [scr O, script letter O]^sim₀ ≈ 0.27 for random 2 is approximately twice that measured in experiments. (Note that this discrepancy will not affect our results since the [scr O, script letter O]^sim₀ value for random 2 is much less than that for the best solution.) These results demonstrate our ability to experimentally construct configurations with maximal [scr O, script letter O]₀ that have been predicted by EA.

3.2 Maximizing force outputs from multiple particles

We next consider objective functions that yield configurations with maximal output forces from mutliple particles in the bottom row, instead of just one particle as described in Section 3.1. Specifically, we will consider the following two cases:

• In case I, we seek a configuration that maximizes the force outputs from the odd-numbered particles on the bottom row (Fig. 4A).

Fig. 4 Assessment of the optimization for cases I and II. (A) For case I, we seek the optimal solution that maximizes the output forces from the odd-numbered particles on the bottom row using the objective function [scr O, script letter O]₁ in eqn (8). (B) [scr O, script letter O]₁ for the best solution and averaged over the whole population plotted versus the number of generations. The results are averaged over three independent trials with random initialization. (C) The probability distribution of [scr O, script letter O]₁ for 5000 randomly generated configurations with vertical lines that indicate [scr O, script letter O]₁ for the average (blue dashed line) and best solution (red solid line) at G = 500. (D) We show one of the two optimal asymmetric solutions for [scr O, script letter O]₁, . The other optimal solution is a mirror image of . We plot F_o,p versus the location p of the particles on the bottom row from the DEM simulations and experiments. The blue shading indicates the odd-numbered particles. (E) For case II, we seek the optimal solution that maximizes the output forces from the even-numbered particles on the bottom row using the objective function [scr O, script letter O]₂ in eqn (9). (F) [scr O, script letter O]₂ for the best solution and averaged over the whole population, plotted versus the number of generations. The results are averaged over three independent trials with random initialization. (G) The probability distribution of [scr O, script letter O]₂ for 5000 randomly generated configurations with vertical lines that indicate [scr O, script letter O]₂ for the average (blue dashed line) and best solution (red solid line) at G = 500. (H) is the optimal solution for [scr O, script letter O]₂. We plot F_o,p versus p from the DEM simulations and experiments. The blue shading indicates the even-numbered particles. F_o,p/F_i is determined in experiments by measuring the slope of F_o,p versus F_i curve (namely, ∂F_o,p/∂F_i) in panels (D) and (H).

• In case II, we seek a configuration that maximizes the force outputs from the even-numbered particles on the bottom row (Fig. 4E).

For case I, we define the objective function,

which is the geometric mean of the force outputs of the odd-numbered particles on the bottom row. Similarly, for case II, we definewhich is the geometric mean of the force outputs of the even-numbered particles on the bottom row. Both [scr O, script letter O]₁ and [scr O, script letter O]₂ are normalized by the input force F_i. We defined [scr O, script letter O]₁ and [scr O, script letter O]₂ so that they have the same units—force—before normalizing by the input force F_i. For both cases I and II, we use the optimization pipeline described in Table 1. When calculating [scr O, script letter O]₁ and [scr O, script letter O]₂ for the experiments, we assume that [scr O, script letter O]^exp₁ = ((∂F_o,1/∂F_i)(∂F_o,3/∂F_i)(∂F_o,5/∂F_i))^1/3 and [scr O, script letter O]^exp₂ = ((∂F_o,2/∂F_i)(∂F_o,4/∂F_i))^1/2.

Similar to the case of the feasibility check, the optimal configuration that maximizes [scr O, script letter O]₁ for case I emerges after ∼80 generations (Fig. 4B) and has a much larger [scr O, script letter O]₁ than that generated using a random search with over 5000 configurations (Fig. 4C). possesses lines of stiff particles that extend from the force input on the top row of particles to the odd-numbered particles in the bottom row (Fig. 4D), similar to the stiffness pattern in C*. Note that the stiffness pattern in is not symmetric about the y-axis, hence there are two solutions: shown in Fig. 4D, and a mirror image of about the y-axis with different genome representations (Table 1), but the same [scr O, script letter O]₁. In addition, in Fig. 4D we show that the F_o,p/F_i values in the experiments, which are determined by ∂F_o,p/∂F_i, are close to those predicted by the DEM simulations for all particles p. In addition, all odd-numbered particles have higher output forces than those for the even-numbered particles on the bottom row, which is consistent with maximizing [scr O, script letter O]₁.

We observe similar results for cases I and II. We obtain the optimal configuration that maximizes [scr O, script letter O]₂ after ∼80 generations (Fig. 4F). It has a much larger [scr O, script letter O]₂ than that generated using a random search (Fig. 4G). is symmetric about the y-axis and the lines of stiff particles are reflected at the side walls, which directs forces from the top plate to the even-numbered particles and avoids the propagation of large forces to the odd-numbered particles (Fig. 4H). We again find that F_o,p/F_i in the experiments are very close to those predicted by the DEM simulations for all p, as shown in Fig. 4H. In , all even-numbered particles have higher output forces than those for the odd-numbered particles on the bottom row, which is consistent with maximizing [scr O, script letter O]₂.

3.3 Optimal solutions with two objectives

Given that there are different force outputs from the particles on the bottom row in , and , one can envision an application where the system needs to switch from (large force outputs from odd-numbered particles) to (large force outputs from even-numbered particles) or vice versa. The VM particles provide the capability to switch the packing from one configuration to another, i.e., by local heating, one can change k for each particle from soft to stiff and vice versa. In contrast, to switch from to with inert particles, one would need to create , disassemble it, and then build . In the movie in the ESI,† we show that using VM particles we can change the stiffness network from to without disassembling the packing by changing the temperatures of the appropriate particles.

To change from to (and vice versa) requires N_s = 19 particle modulus changes. However, what if there is a hardware constraint for the VM particles where significant power is consumed to change the modulus, but not to maintain the particle modulus as soft or stiff? In this case, minimization of the power consumption corresponds to minimizing the number of switches in modulus between the two configurations and a multi-objective EA is required to identify the optimal solution.

For case III, we consider the following multi-objective optimization problem. Starting from , we seek to maximize [scr O, script letter O]₂ and minimize

where is the number of particle modulus changes from to C₂, δ(x,y) = 1 when x = y and δ(x,y) = 0 when x ≠ y, and k^C,p is the modulus of the pth particle in configuration C.

We ran the multi-objective optimization using AFPO with parameters outlined in Table 1. The algorithm converged to a Pareto front with N_s + 1 = 20 optimal solutions at G = 500, as shown in Fig. 5A. There is no “best” solution among the 20 solutions on the Pareto front if we do not consider trade-offs between maximizing [scr O, script letter O]₂ and minimizing [scr O, script letter O]₃. To implement trade-offs, we include a TOPSIS weight vector: W = [w,1 − w] with 0 ≤ w ≤ 1. The first and second components of W correspond to the weight for [scr O, script letter O]₃ and [scr O, script letter O]₂, respectively. We find that the best solution given by TOPSIS traverses the Pareto front from top-right to bottom-left in the [scr O, script letter O]₃–[scr O, script letter O]₂ plane. Namely, [scr O, script letter O]₃ and [scr O, script letter O]₂ decrease monotonically with increasing w, as shown in Fig. 5B. We also observe that the dependence of both [scr O, script letter O]₃ and [scr O, script letter O]₂ on w on the Pareto front is non-linear.

Fig. 5 Multi-objective optimization to find a configuration that maximizes the output forces on the even-numbered particles while minimizing the number of particle modulus changes starting from . (A) Solutions in the [scr O, script letter O]₃–[scr O, script letter O]₂ plane at G = 500 (yellow circles). We also highlight the Pareto front (black squares). The orange diamond, magenta star, and purple triangle indicate the best solutions based on the TOPSIS method for weight vectors W = [w,1 − w], where the first and second elements indicate the weights for [scr O, script letter O]₃ and [scr O, script letter O]₂, respectively. (B) [scr O, script letter O]₂ (left axis, solid blue line) and [scr O, script letter O]₃ (right axis, dashed red line) for the best solution selected from the Pareto front using TOPSIS plotted versus w. Black squares connected by dotted lines indicate solutions from the Pareto front in panel (A). (C)–(E) Starting from (left), we show the best solution (right) from the Pareto front with a given TOPSIS weight vector W: (C) W = [1,0], C^[1,0]₂, (D) W = [0.5,0.5], C^[0.5,0.5]₂, and (E) W = [0,1], C^[0,1]₂.

Next, we present several examples of the best solutions selected by TOPSIS for different values of w. For instance, if we set w = 1, the best solution is with no particle modulus changes (Fig. 5C). Similarly, if we set w = 0, only objective [scr O, script letter O]₂ is considered, and the best solution is (Fig. 5E). These two TOPSIS weight vectors have obvious solutions, and confirm that and are on the Pareto front. However, if we set 0 < w < 1 (e.g., w = 0.5), the best solution is not or , and is instead a compromise between maximizing [scr O, script letter O]₃ and minimizing [scr O, script letter O]₂, as shown in Fig. 5D.

For case IV, we consider the following multi-objective optimization problem. Starting from , we seek to maximize [scr O, script letter O]₁ and minimize

We again observe a Pareto front with N_s + 1 = 20 optimal solutions at G = 500, as shown in Fig. 6A. We implement a similar TOPSIS weight vector: W = [w,1 − w], with the first and second components corresponding to the weight for [scr O, script letter O]₄ and [scr O, script letter O]₁, respectively. Similar to case III, we find that [scr O, script letter O]₄ and [scr O, script letter O]₁ decrease monotonically with increasing w for the best solution given by TOPSIS, as shown in Fig. 6B. We again observe that the dependence of both [scr O, script letter O]₄ and [scr O, script letter O]₁ on w on the Pareto front is non-linear. We also present several examples of the best solutions selected by TOPSIS for values of w. For instance, if we set w = 1, the best solution is with no particle modulus changes (Fig. 6C). Similarly, if we set w = 0, the best solution is , where only objective [scr O, script letter O]₁ is considered (Fig. 6E). If we set 0 < w < 1 (e.g., w = 0.5), the best solution is not or , and is instead a compromise between maximizing [scr O, script letter O]₃ and minimizing [scr O, script letter O]₂, as shown in Fig. 6D.

Fig. 6 Multi-objective optimization to find a configuration that maximizes the output forces on the odd-numbered particles while minimizing the number of particle modulus changes starting from . (A) Solutions in the [scr O, script letter O]₄–[scr O, script letter O]₁ plane at G = 500 (yellow circles). We also highlight the Pareto front (black squares). The orange diamond, magenta star, and purple triangle indicate the best solutions based on the TOPSIS method for weight vectors W = [w,1 − w], where the first and second elements indicate the weights for [scr O, script letter O]₄ and [scr O, script letter O]₁, respectively. (B) [scr O, script letter O]₁ (left axis, solid blue line) and [scr O, script letter O]₄ (right axis, dashed red line) for the best solution selected from the Pareto front using TOPSIS plotted versus w. Black squares connected by dotted lines indicate solutions from the Pareto front in panel (A). (C) and (D) Starting from (left), we show the best solution (right) from the Pareto front with a given TOPSIS weight vector W: (C) W = [1,0], C^[1,0]₁, (D) W = [0.5,0.5], C^[0.5,0.5]₁, and (E) W = [0,1], C^[0,1]₁.

4 Conclusions

In this article, we developed a pipeline that combines experiments using variable modulus (VM) particles, discrete element method (DEM) simulations, and evolutionary algorithms (EAs) to design granular packings with adaptive mechanical responses. Specifically, we use the pipeline to design granular packings composed of N = 23 same-sized VM particles arranged in a triangular lattice with different patterns of particle moduli that maximize the output forces from particles on the bottom row of the system. We realized VM particles in experiments by embedding a Field's metal core within a silicone shell. By running current through a resistive heater co-located with the Field's metal, we achieved a modulus ratio of k_stiff/k_soft ∼ 3 when the Field's metal is heated above the melting temperature. The system boundaries were made of photoelastic material to enable measurements of the output forces. We characterized the inter-particle contact forces as a function of compressive strain in experiments to calibrate the inter-particle forces in the DEM simulations. We carried out evolutionary algorithms (EA) to identify the optimal configurations of soft and stiff particles for specific force outputs on the bottom boundary without enumerating all possible particle modulus combinations.

We considered five optimization cases in this study. For the feasibility check, we found the configuration C* of particle moduli that maximizes the output force from the middle particle in the bottom row. In case I (case II), we identified the configuration that maximizes the output forces from the odd-numbered (even-numbered) particles in the bottom row. In case III (case IV), we found the configuration C₂ (C₁) that maximizes the output forces from the even-numbered (odd-numbered) particles in the bottom row, while minimizing the number of particle modulus changes starting from . The first three cases involve single-objective optimization, and the remaining two involve multi-objective optimization. We have shown that EA converges to the optimal solution for the cases involving single-objective optimization and a set of solutions on the Pareto front for the cases involving multi-objective optimization. We verify that the force outputs for the best solutions match in experiments and DEM simulations for the cases of single-objective optimization. We further show that we can identify the best solution in the cases of multi-objective optimization when trade-offs between the objectives are implemented using a TOPSIS weight vector. This study has demonstrated the ability to search efficiently through a large ensemble of granular packings for those that satisfy specific objectives concerning the mechanical response. Although there have been several numerical studies aimed at identifying optimal properties in granular packings using EA,¹¹ few have validated the EA methods using experiments.⁴⁴

We envision several interesting directions for future studies. First, we can examine all of the solutions on the Pareto front in cases III and IV, implement the optimal solutions in experiments, and determine the power consumption when generating them. Second, our results can serve as a proof-of-concept for mechanical logic gates in granular materials, where continuous force values are discretized into binary inputs (“0” or “1”) and routed through different logic gate architectures. Previous work has shown logic capabilities within granular materials,⁵⁰ but these studies relied on the embedded capabilities of the granular material rather than optimizing for a desired logic gate. Although no traditional form of logic has been investigated in the current study, we believe that the proposed pipeline can serve as a tool to optimize a material to perform logical operations. Third, the present study can be readily extended to bulk granular materials in three dimensions (3D). We anticipate that the greater number of contacts per particle in 3D packings will increase the complexity of the configurational search space for identifying optimal force networks.

The heterogeneity of force chain networks in granular packings is often controlled by positional disorder. Previous studies of force chain networks in granular materials have analyzed force propagation in disordered networks using methods from statistical physics and mathematics.^29,51–54 In this work, we tuned the force chain networks through particle modulus variations (k_stiff/k_soft ∼ 3) on a fixed triangular lattice, rather than through positional disorder. Although this modulus ratio is substantial and has allowed us to significantly vary the morphology of force chain networks, an optimized particle design can induce a more dramatic particle modulus ratio and thus more dramatic changes in the force chain networks in VM particle packings. To efficiently sample and quantify force chain networks, we implemented DEM simulations with a simplified contact model for VM particles and the system boundaries. However, the experimental system and the DEM simulations can have differing force outputs for some of the sub-optimal configurations. Increasing the complexity of the inter-particle contact model in the DEM simulations, e.g., using deformable particle simulations in 3D,⁵⁵ can improve the agreement between the predicted and experimentally measured force outputs.

Author contributions

S. W.: methodology, validation, formal analysis, investigation, data curation, writing, review and editing, visualization, and project administration. A. P.: methodology, formal analysis, validation, writing, review and editing, and visualization. D. W.: methodology, formal analysis, writing, review and editing, and visualization. N. P.: methodology, investigation, supervision, review, and editing. J. Z.: methodology, formal analysis. A. K. M.: methodology. M. D. S.: investigation, resources. R. K.-B.: conceptualization, supervision, project administration, review and editing, resources, and funding acquisition. C. S. O.: conceptualization, supervision, project administration, review and editing, resources, and funding acquisition. J. B.: conceptualization, supervision, project administration, review and editing, resources, and funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

We have deposited all relevant data for this paper in the repository https://github.com/AtoosaParsa/AdaptiveForceChains.

Acknowledgements

This material is based upon work supported by the National Science Foundation under the Designing Materials to Revolutionize and Engineer our Future (DMREF) program (award no. 2118988).

Notes and references

1. C. Jiang, F. Rist, H. Wang, J. Wallner and H. Pottmann, Comput.-Aided Des., 2022, 143, 103146.
2. W. Yang, Z.-M. Li, W. Shi, B.-H. Xie and M.-B. Yang, J. Mater. Sci., 2004, 39, 3269–3279.
3. Y. Wang, L. Li, D. Hofmann, J. Andrade and C. Daraio, Nature, 2021, 596, 238–243.
4. B. Jenett, C. Cameron, F. Tourlomousis, A. P. Rubio, M. Ochalek and N. Gershenfeld, Sci. Adv., 2020, 6, eabc9943.
5. G. Gantzounis, M. Serra-Garcia, K. Homma, J. Mendoza and C. Daraio, J. Appl. Phys., 2013, 114, 093514.
6. A. Leonard, L. Ponson and C. Daraio, J. Mech. Phys. Solids, 2014, 73, 103–117.
7. L. Bonanomi, G. Theocharis and C. Daraio, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2015, 91, 033208.
8. N. Nejadsadeghi, L. Placidi, M. Romeo and A. Misra, Mech. Res. Commun., 2019, 95, 96–103.
9. E. Kim and J. Yang, Funct. Compos. Struct., 2019, 1, 012002.
10. F. Li, P. Anzel, J. Yang, P. G. Kevrekidis and C. Daraio, Nat. Commun., 2014, 5, 1–6.
11. Q. Wu, C. Cui, T. Bertrand, M. D. Shattuck and C. S. O'Hern, Phys. Rev. E, 2019, 99, 062901.
12. K. Fu, Z. Zhao and L. Jin, Adv. Funct. Mater., 2019, 29, 1901258.
13. M. A. Porter, P. G. Kevrekidis and C. Daraio, Phys. Today, 2015, 68, 44–50.
14. N. Nejadsadeghi and A. Misra, Int. J. Mech. Sci., 2019, 161, 105042.
15. A. Merkel and J. Christensen, Commun. Phys., 2019, 2, 1–6.
16. A. Parsa, D. Wang, C. S. O’Hern, M. D. Shattuck, R. Kramer-Bottiglio and J. Bongard, International Conference on the Applications of Evolutionary Computation (Part of EvoStar), 2022, pp. 93–109.
17. A. Parsa, D. Wang, C. S. O’Hern, M. D. Shattuck, R. Kramer-Bottiglio and J. Bongard, Proceedings of the Genetic and Evolutionary Computation Conference, 2022, pp. 122–129.
18. A. Parsa, S. Witthaus, N. Pashine, C. S. O’Hern, R. Kramer-Bottiglio and J. Bongard, Proceedings of the Genetic and Evolutionary Computation Conference, 2023, pp. 193–201.
19. H. Herrmann, J. Hovi and S. Luding, Physics of Dry Granular Media, Springer, Netherlands, 2014.
20. H. M. Jaeger, S. R. Nagel and R. P. Behringer, Rev. Mod. Phys., 1996, 68, 1259–1273.
21. R. Behringer, D. Howell, L. Kondic, S. Tennakoon and C. Veje, Phys. D, 1999, 133, 1–17.
22. D. W. Howell, R. P. Behringer and C. T. Veje, Chaos, 1999, 9, 559–572.
23. B. Utter and R. Behringer, Eur. Phys. J. E:Soft Matter Biol. Phys., 2004, 14, 373–380.
24. T. S. Majmudar and R. P. Behringer, Nature, 2005, 435, 1079–1082.
25. R. Behringer, H. Jaeger and S. Nagel, Chaos, 1999, 9, 509–510.
26. F. Radjai, M. Jean, J.-J. Moreau and S. Roux, Phys. Rev. Lett., 1996, 77, 274–277.
27. S. Luding, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 55, 4720–4729.
28. J. F. Peters, M. Muthuswamy, J. Wibowo and A. Tordesillas, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2005, 72, 041307.
29. S. F. Edwards and D. V. Grinev, Granular Matter, 2003, 4, 147–153.
30. A. H. Clark, A. J. Petersen, L. Kondic and R. P. Behringer, Phys. Rev. Lett., 2015, 114, 144502.
31. L. Papadopoulos, J. G. Puckett, K. E. Daniels and D. S. Bassett, Phys. Rev. E, 2016, 94, 032908.
32. D. Wang, J. Ren, J. A. Dijksman, H. Zheng and R. P. Behringer, Phys. Rev. Lett., 2018, 120, 208004.
33. C. S. Campbell, Powder Technol., 2006, 162, 208–229.
34. N. Pashine, A. M. Nasab and R. Kramer-Bottiglio, Soft Matter, 2023, 19, 1617–1623.
35. W. Shan, T. Lu and C. Majidi, Smart Mater. Struct., 2013, 22, 085005.
36. M. J. Puttock and E. G. Thwaite, Elastic compression of spheres and cylinders at point and line contact, by M. J. Puttock and E. G. Thwaite, Commonwealth Scientific and Industrial Research Organization Melbourne, 1969, p. 64.
37. S. V. Franklin and M. D. Shattuck, Handbook of granular materials, CRC Press, 2016.
38. Y. Zhao, Y. Zhao, D. Wang, H. Zheng, B. Chakraborty and J. E. S. Socolar, Phys. Rev. X, 2022, 12, 031021.
39. G. D. J. de Jong, in Étude photo-Élastique d'un empilement de disques, ed. R. J. Schotting, H. C. J. van Duijn and A. Verruijt, Springer, Netherlands, Dordrecht, 2006, pp. 148–161.
40. T. Wakabayashi, J. Phys. Soc. Jpn., 1950, 5, 383–385.
41. A. Resnick, Contemp. Phys., 2016, 58, 1.
42. A. Flamant, C. R. Acad. Sci., 1892, 114, 1465–1468.
43. D. E. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley Longman Publishing Co., Inc., USA, 1st edn, 1989.
44. M. Z. Miskin and H. M. Jaeger, Soft Matter, 2014, 10, 3708–3715.
45. D. Shah, B. Yang, S. Kriegman, M. Levin, J. Bongard and R. Kramer-Bottiglio, Adv. Mater., 2021, 33, 2002882.
46. R. Chiong, T. Weise and Z. Michalewicz, Variants of evolutionary algorithms for real-world applications, Springer, 2012, vol. 2.
47. M. D. Schmidt and H. Lipson, Proceedings of the 12th annual conference on Genetic and evolutionary computation, 2010, pp. 543–544.
48. M. Schmidt and H. Lipson, in Age-Fitness Pareto Optimization, ed. R. Riolo, T. McConaghy and E. Vladislavleva, Springer New York, New York, NY, 2011, pp. 129–146.
49. C.-L. Hwang, K. Yoon, C.-L. Hwang and K. Yoon, Multiple attribute decision making: methods and applications a state-of-the-art survey, 1981, pp. 58–191.
50. F. Li, P. Anzel, J. Yang, P. G. Kevrekidis and C. Daraio, Nat. Commun., 2014, 5, 5311.
51. M. Kramar, A. Goullet, L. Kondic and K. Mischaikow, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2013, 87, 042207.
52. E. Stanifer, P. K. Morse, A. A. Middleton and M. L. Manning, Phys. Rev. E, 2018, 98, 042908.
53. M. Kramar, A. Goullet, L. Kondic and K. Mischaikow, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2013, 87, 042207.
54. L. Kondic, A. Goullet, C. S. O'Hern, M. Kramar, K. Mischaikow and R. P. Behringer, Europhys. Lett., 2012, 97, 54001.
55. D. Wang, J. D. Treado, A. Boromand, B. Norwick, M. P. Murrell, M. D. Shattuck and C. S. O’Hern, Soft Matter, 2021, 17, 9901–9915.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sm00965g

‡ These authors contributed equally to this work.

---