Tuning the size and stiffness of inflatable particles

Abstract

We describe size-varying cylindrical particles made from silicone elastomers that can serve as building blocks for granular materials with tunable structural and mechanical properties. The particle size variation, which is achieved by inflation, gives rise to changes in stiffness under compression. We design and fabricate inflatable particles that can become stiffer or softer during inflation, depending on key parameters of the particle geometry, such as the ratio of the fillet radius to the wall thickness, r/t. We also conduct numerical simulations of the inflatable particles and show that they only soften during inflation when localization of large strains occurs in the regime r/t → 0. This work introduces novel particle systems with tunable size and stiffness that can be implemented in particle packings for soft robotic applications.

---

1. Introduction

Granular materials composed of discrete frictional particles possess a host of complex yet tunable mechanical properties that have been the subject of study for decades.^1–3 Granular metamaterials can be designed to exhibit targeted mechanical responses, such as phononic band gaps and the ability to propagate solitons and edge modes.^4–6 Granular metamaterials are typically composed of uniform particles with fixed mechanical properties.^7–10 However, recent work has expanded this paradigm by introducing mixtures of dissimilar particles^11,12 and by designing particles with tunable properties that enable dynamic reconfiguration of the assembly post-fabrication.¹³ These advances point toward a new class of “robotic” granular materials: systems in which each particle can modulate its own size, shape, stiffness, or other intrinsic properties to influence the global behavior of the assembly. Robotic particles could enable packings that adapt their bulk mechanical properties on demand, recover from damage by rerouting force chains, or implement reprogrammable mechanical logic.

In this work, we present the design and fabrication of inflatable particles with independently tunable size and stiffness, advancing the toolkit for robotic granular materials. Prior studies have demonstrated stiffness modulation and particle-scale actuation;^13,14 a unique contribution of the current work is the ability to control the stiffness-inflation pressure relationship of an inflatable particle. By tuning specific geometric parameters, we achieve either stiffening or softening behavior under inflation, providing a new level of control over the mechanical response of individual particles and their assemblies.

This study focuses on two types of inflatable cylindrical particles: type I: cylindrical shells and type II: toroidal shells. Type I particles are empty inside and expand radially, as well as along the cylindrical axis of the particle, when inflated. In contrast, type II particles mostly expand radially, with minimal deformation along the cylindrical axis. Images of the two particle types are shown in Fig. 1. We limit the expansion of type II particles along the cylindrical axis by connecting the top and bottom surfaces with a central pillar.

Fig. 1 Schematic of the two types of particles (left) with images in the uninflated (middle) and inflated (right) states. Side views of (A) type I and (B) type II particles. (C) Top view of a type II particle. The type I and II particles have been inflated to ∼3.4 kPa and ∼8.6 kPa, respectively. All scale bars correspond to 10 mm.

The particle geometries are specified by several dimensionless ratios. In this work, we vary the shapes of both types of particles mainly by changing the ratio of the cylinder's corner radius (r in Fig. 1) to its thickness (t in Fig. 1). We measure the force response as a function of pressure when these particles are compressed radially. Studying these two particle designs allows us to determine how the particle deformation during inflation affects its compressive stiffness.

For both types of particles, we find that the compressive stiffness can increase, decrease, or stay the same during inflation. With finite element method simulations, we show that the slope of the stiffness versus inflation pressure curve is negative only when substantial strain localization occurs near the cylinder's “corners” (where the top wall meets the side wall) during inflation, resulting in strain softening (see Section 2.3). We also demonstrate that the mode of particle deformation—whether cylindrical or radial—during inflation has little impact on stiffness variation, as both type I and type II particles exhibit similar stiffness behavior across different inflation pressures. These findings are used to design particles with specific variation of the compressive stiffness with inflation pressure, i.e., increasing, decreasing, or constant stiffness. We also generate example packings of inflatable particles to illustrate that they can possess different compressive stiffness even at the same packing fraction.

The remainder of the article is organized as follows. In Section 2, we describe the manufacturing processes for type I and II particles. We also describe the finite element method (FEM) simulations conducted to understand the compressive stiffness for type I particles during inflation. In Section 3, we discuss experimental measurements of the particle stiffness versus inflation pressure for type II particles, as well as measurements of the collective stiffness of packings of type II particles. In Section 4, we summarize the conclusions and propose promising directions for future research. In Appendix A, we show that the FEM simulation results for the particle compressive stiffness do not depend strongly on the friction coefficient between the particles and compression plates. In Appendix B, we include the experimental results for the relation between the change in radius of type II particles and inflation pressure.

2. Methods

2.1. Particle fabrication

The inflatable particles are fabricated from EcoFlex™ 50, a two-part elastomer that reaches a 00–50 shore hardness upon curing. Shore hardness is a scale that measures the resistance to localized deformation of a material, and Shore 00–50 corresponds to a relatively soft material with a Young's modulus of ∼0.1 MPa.^15,16 A schematic of the fabrication process for both particle types is shown in Fig. 2. The particles are molded as two halves and then attached to each other. The molds for each of the parts and for the final assembly are 3D-printed in PLA on a Prusa™ MK3S+ printer, with a precision δ = 0.1 mm. The fabrication process involves first molding the top and bottom halves (mold A in Fig. 2) and then assembling the two parts using mold B in Fig. 2 to make the particle. During the molding process, a silicone tube is attached to the top part of the particle using Silpoxy™, which serves as an air inlet. After the two parts of the particle are cured, they are glued using a thin layer of uncured EcoFlex™ 50 and pressed in place using the second part of the mold. This particle assembly process ensures sufficient contact between the two particle layers so that leaks are avoided.

Fig. 2 Summary of the particle fabrication process. Each type of particle requires two sets of molds: one for molding the two halves (mold A) and one for assembling the two halves together (mold B). The components of the molds are 3D-printed out of Polyvinyl Lactic Acid (PLA). The blue (green) molds are used for type I (type II) particles. The fabrication process for both particle types is similar, with differences in the mold design leading to differences in the internal structure of the particles. The seven-step fabrication process for each particle type is shown. The particles are molded from EcoFlex™ 50.

The geometry of the uninflated type I particles is characterized by four particle shape parameters: the height H, particle radius R, wall thickness t, and fillet radius r, as shown in Fig. 1. In addition to these parameters, uninflated type II particles are also parametrized by the inner annular radius R_a. The structural and mechanical properties of type I particles are governed by three dimensionless particle shape parameters – H/t, R/t, and r/t – which are derived from the system's four geometric parameters. Type II particles with five defining particle shape parameters are characterized by four dimensionless particle shape parameters: the same three as type I particles (H/t, R/t, and r/t) and R_a/t. The experiments select H = 10 mm, R = 10 mm, t = 2 mm, and δ ≤ r ≤ 4 mm for type I particles and H = 10 mm, R = 20 mm, 6 mm ≤ R_a ≤ 12 mm, 1 mm < t < 4 mm, and δ < r < 3 mm for type II particles. For type I particles, we fix H/t = 5, R/t = 5, and vary r/t, with the smallest r/t denoted as (r/t)_min = δ/t = 0.05. For type II particles, we fix zero, one, two, or three of the dimensionless particle shape parameters at a time, while varying the rest, and set (r/t)_min = 0.1.

The air inlet in type I particles allows airflow into the cavity of the particle, as shown in Fig. 1A. The central pillar in type II particles has channels connecting the inlet tube to the outer cavity to allow airflow into the particles, as shown in Fig. 1B. The air inlets in the type I and II particles enable changes in the particle's size and stiffness.

2.2. Experimental measurements of particle stiffness

In this work, we consider the compressive stiffness, which is defined as the slope of the normal force versus displacement curve of a particle under uniaxial compression. The stiffness of the particles is measured under uniaxial compression in the direction perpendicular to the cylindrical axis using a materials tester (Instron™ 3365) over a range of inflation pressures p < 20 kPa. The particles are inflated by connecting them to a high-pressure air line whose output pressure can be regulated to within an error of Δp = ±0.345 kPa. The inflation at a given p will change the particle radius, which we denote as R(p), as well as all of the other particle shape parameters (i.e., r, t, H, and R_a). The compressive stiffness k of the particles during inflation is determined from the slope of the force versus displacement curve obtained from the materials tester. Specifically, it is calculated based on the compressive force versus the displacement of the compression plate, measured from the point of initial contact with the particle at each fixed inflation pressure. The maximum displacement is 2 mm, which corresponds to a compressive strain of ∼5%. We also image the particles during the compression tests at each inflation pressure to determine their size and shape, enabling comparisons to the FEM simulations of particle deformation.

2.3. FEM simulations for type I particles

To understand the particle compressive stiffness as a function of inflation pressure, we employ two FEM models to simulate the deformation of type I particles: a hyperelastic model implemented in ABAQUS, and a more general spring network model that allows us to tune the pressure-dependent Young's modulus. We focus on modeling type I particles due to their simple geometry compared to type II particles.

For the first FEM approach, we use an incompressible neo-Hookean model¹⁷ in ABAQUS to simulate inflatable particles made from EcoFlex™ 50. The strain energy density for a neo-Hookean material is

U = C₁(I₁ − 3) (1)

where C₁ = 30 kPa is proportional to the Young's modulus, I₁ = λ₁² + λ₂² + λ₃² is the first invariant of the right Cauchy-Green deformation tensor [scr F, script letter F], and λ₁, λ₂, and λ₃ are the eigenvalues of [scr F, script letter F] in descending order. For an incompressible material, C₁ is half of the shear modulus in the limit of zero deformation, which yields a Young's modulus of 6C₁ = 180 kPa at zero strain for EcoFlex™ 50.^18,19 In addition, we set the rigid parallel plates that apply the compressive strain to be rough and exert frictional forces once in contact with the particle (Fig. 3A). Previous studies have shown that the friction coefficient μ between Ecoflex and steel varies over a wide range , depending on the contact pressure and duration of the contact.^20,21 In this study, we set μ = 0.5, which is near the lower end of the reported experimental range.^20,21 We show in Appendix A that the results for the compressive stiffness do not depend strongly on μ.

Fig. 3 (A) (left panel) Sketch of the spring network model for a type I particle (cyan) with r/t = 4 compressed between two parallel rough, rigid plates (black) in a direction perpendicular to the cylindrical axis. (center panel) Cross-sectional view of the particle in the left panel, which shows the connections between nodes of the spring network (blue lines) and rough surface of the walls modeled by spheres on a square lattice. (right panel) The inflation pressure is applied through forces F⃑_p,f split evenly on each node of each triangle f that tessellate the interior surface of the cylindrical shell of the particle. (B) Spring network meshes (blue lines) for type I particles with different ratios of the fillet radius to the wall thickness r/t (from left to right and from top to bottom): r/t = 0, 0.5, 1, and 2. The Cartesian coordinate systems in both panels illustrate the particle orientations.

Unlike linear elastic materials, hyperelastic materials possess a nonlinear stress versus strain relation with a slope that can decrease at large strains. The decrease in slope, known as strain softening, can be described by several phenomenological nonlinear material models.²² In our studies, we seek to develop particle geometries that give rise to strain softening during particle inflation at increased pressure. Under uniaxial extension, , where λ is the extensional deformation, and eqn (1) becomes

We can then calculate the second derivative of U with respect to λ, which gives the Young's modulus:

For this model, Y decreases with increasing λ. Note that Y = 6C₁ in the limit of zero deformation with λ = 1.

The degree of strain softening, i.e., how Y varies with λ, is determined by the ratio of the constant term and the coefficient of the 1/λ³ term in eqn (3). The ratio is fixed at 0.5 for the neo-Hookean model. To vary the degree of strain softening, we developed another FEM-based model for type I particles. We meshed the particles using N_t ∼ 2 × 10⁴ Delaunay tetrahedra that include a total of N ∼ 6600 nodes connected by N_l ∼ 37 000 springs as shown in Fig. 3A and B. The potential energy for spring i is

where u_l is the characteristic energy scale of the spring, λ_i = l_i/l_i,0 is the stretch of the spring with instantaneous length l_i and rest length l_i,0, and 0 ≤ a ≤1. The cylindrical shells are meshed such that the average equilibrium rest length is 〈l_i,0〉 = 0.75 mm, which yields a number density of n_s ∼ 25 mm⁻³ for the spring networks of type I particles. This mesh density is sufficiently large such that the results do not change when we increase the mesh density above this value. We set u_l = 0.0072 J so that the effective Young's modulus of the spring network for the type I particle is Y ∼ n_su_l ∼ 180 kPa, which is the same as that for the ABAQUS model. The second derivative of U_i with respect to λ_i iswhich resembles eqn (3), except that we can tune the local Young's modulus by changing a. Here, a = 1 corresponds to a linear spring, and a < 1 corresponds to a spring that softens upon stretching.

To simulate purely repulsive, frictional contacts between the inflatable particles and compression plates, we model each rough, rigid parallel plate using a collection of M monodisperse spheres fixed to maintain a square lattice with spacing d = d_s, where d_s is the sphere diameter, as shown in Fig. 3A. Each node of the spring network within the inflatable particle also includes a sphere with diameter d_s that only interacts with the spheres within the compression plate through the interaction potential:

where u_p the repulsive energy scale, r^s_nm is the distance between the nth sphere in the inflatable particle and the mth sphere in the sth wall, and Θ(·) is the Heaviside step function, which ensures that the node spheres and spheres that make up the walls do not interact when they are not in contact. We set d_s = 0.75 〈l_i,0〉 = 0.75 mm and u_p = 5u_l. To inflate the particles to pressure p, we first identify all N_f triangles in the mesh that cover the interior of the particle. For the fth triangle, we assume that the total force resulting from the inflation pressure is

F⃑_p,f = pA⃑_f (7)

where A⃑_f is a vector normal to the triangle pointing outwards with a magnitude equal to the triangle area. We then evenly distribute F⃑_p,f on the three nodes associated with the fth triangle.

For both FEM-based modeling approaches, we first apply a given inflation pressure, followed by compression tests with the inflation pressure held fixed. In the ABAQUS model, at each p, we perform a static analysis with a compression step size of 0.35 mm and record the reaction force on one of the compression plates. In the spring network model, we apply the compression quasi-statically. That is, we apply successive compression steps with size 0.35 mm to reach a desired total compression. After each compression step, we use the FIRE algorithm²³ to minimize the total potential energy at fixed p until the system reaches a force-balanced state, such that the net force magnitude on each node is less than 10⁻¹²u_l/〈l_i,0〉. In both ABAQUS and the spring network model, we only measure the perpendicular component of the force F⃑ on one of the compression plates in response to each compression step. In the spring network model, , where r⃑^s_m is the position of the mth sphere on the compression plate with s = 1. We determine the particle stiffness k by calculating the slope of a linear fit of F versus the compressive strain up to strains of 10%.

3. Results

In this section, we present experimental measurements of the compressive stiffness k for type I and II particles, along with numerical simulation results for k in type I particles, all analyzed as a function of inflation pressure. First, we show that the slope of k versus inflation pressure p can be tuned from positive to negative as we change the internal structure of the type I particles. In particular, the slope of k versus p decreases with decreasing r/t, becoming negative for r/t ⪅ 0.5 for type I particles. We implement a neo-Hookean model for type I particles in ABAQUS and show that it recapitulates k versus p over a wide range of r/t. For particles with negative k versus p slopes, we observe localization of large strains near the corners of the cylindrical particles during compression. In contrast, when the slope of k versus p is positive, localization of large strains does not occur. Using a more general spring network model, we show that strain softening, arising from the localization of large strains, gives rise to the negative k versus p slopes. Finally, we find that the dependence of the slope of k versus p on r/t is similar for both type I and type II particles.

3.1. Compressive stiffness for type I particles

While it was expected that the compressive stiffness k of inflatable particles depends on the inflation pressure p, a key finding of this work is that modifying the particle geometry allows us to tune this pressure dependence. In Fig. 4A, we show k(p) for type I particles for several values of r/t. We find that k increases with pressure for large values of r/t, k is nearly independent of p for r/t∼ 0.5, and k decreases with p over a large range of p for r/t = (r/t)_min. This behavior is in contrast to linearly elastic particles, where k is independent of pressure.

Fig. 4 (A) Compressive stiffness k of type I particles with different ratios of the fillet radius to the wall thickness r/t (indicated by color) normalized by k₀ at zero inflation pressure plotted versus the inflation pressure p (in kPa). Both the experimental (solid symbols with dashed lines) and ABAQUS simulation results using a neo-Hookean model (open symbols with solid lines) are indicated. Each data point in experiments is an average over three samples measured three times each, with the error bars on k/k₀ indicating the standard deviation. Error bars on p are the precision Δp in the pressure gauge. (B) k/k₀ versus p for type I particles simulated in ABAQUS with the same r/t = 2, R/t = 5, and H/t = 5, but different t:t = 1 mm (upward triangles), 1.5 mm (squares), 2 mm (circles), 2.5 mm (diamonds), 3 mm (leftward triangles), and 4 mm (stars). (C) Spatial distribution of the maximum principal strain, λ₁ − 1, at p = 6.89 kPa inflation pressure from the ABAQUS simulations. The particles from left to right and top to bottom have r/t = (r/t)_min, 0.5, 1, and 2. The cross section is through the particle center and perpendicular to both compression plates and the particle's long axis. The Cartesian coordinate system illustrates the orientation of the particles. The color bar indicates values of λ₁ − 1 for all panels. (D) Mean of the largest 1% of 〈λ₁〉 − 1 values in the particle as a function of the inflation pressure for type I particles with different r/t from ABAQUS simulations.

The ABAQUS simulations of the neo-Hookean model recapitulate the experimental results to within an average normalized root-mean-square error of , where k^E and k^A are compressive stiffness values from the experiments and ABAQUS simulations and n > 10 is the number of measurements. We observe that the largest deviation between experimental measurements and ABAQUS simulations are from the particles with r/t = (r/t)_min at p > 4 kPa, where k/k₀ from experimental measurements reaches a plateau while k/k₀ from the ABAQUS simulations reaches a minimum at p ∼ 4 kPa and starts to increase with p. We postulate that the increase in k/k₀ for p > 4 kPa is caused by the overestimation of stress in the neo-Hookean model at large tensile strains (between 50% and 200%).²⁴ Note that ABAQUS simulations using the Mooney-Rivlin and Yeoh hyperelastic models can match the experimental results for k(p) at r/t = (r/t)_min when p > 4 kPa since they describe the stress–strain relationship better at large strains.²⁴ We further demonstrate that scaling all of the geometrical parameters by the same factor does not change k(p). Specifically, in Fig. 4B, k(p) does not vary significantly when we change r, t, H, and R, while keeping the ratios r/t = 2, H/t = 5, and R/t = 5 fixed.

3.2. Numerical simulations of the compressive stiffness of type I particles

To explain the compressive stiffness k versus pressure p behavior for type I particles for different values of r/t in Fig. 4A, we calculated the spatial distribution of the principal strains at each inflation pressure. In Fig. 4C, we show how the spatial distribution of the maximum principal strain λ₁ − 1 depends on r/t. For large r/t, the maximal strain distribution is spatially uniform. As r/t decreases, λ₁ − 1 becomes non-uniform and the largest maximal strains occur near the fillet. In Fig. 4D, we also show that the mean value of the largest 1% of maximal strains in the particle (〈λ₁〉 − 1) increases with decreasing r/t. We calculate 〈λ₁〉 − 1 by finding the largest 1% of λ₁ from all of the elements (tetrahedra) that compose the particle in the ABAQUS simulations at each p. In the previous section, we showed that the most significant strain softening occurs in the r/t→0 limit. The results in Fig. 4D suggest that the localization of the largest strains near the fillet controls the strain softening (i.e., decreases in k with increasing p) in the r/t → 0 limit.

To further understand the pressure-dependent compressive stiffness, we also carried out numerical studies of the spring network model for type I particles, which allowed us to tune the degree of strain softening during deformation. Since the smallest element that describes the spring-network-modelled particles is a single spring or bond, we investigate the deformation of elements using the spring stretch λ_i, as opposed to λ₁ from the ABAQUS simulations. In Fig. 5A and B, we plot the mean of the largest 1% of strains, 〈λ_i〉 − 1, as a function of p for a = 1 and 0, where the parameter a that tunes the degree of strain softening is defined in eqn (4). The largest 1% is identified at each p. We find that 〈λ_i〉 − 1 decreases with increasing r/t, for all values of a. Therefore, the most substantial strain softening occurs in the r/t→0 limit. In addition, we show in Fig. 5C that 〈λ_i〉 − 1 increases with decreasing a, which confirms that the most significant strain softening occurs as a → 0, as well as r/t → 0.

Fig. 5 The mean of the largest 1% of strains 〈λ_i〉 − 1 plotted as a function of the inflation pressure p for type I particles with various r/t modeled using spring networks with (A) a = 1 and (B) 0. (C) Mean of the largest 1% of strains 〈λ_i〉 − 1 plotted versus a at p = 6.89 kPa for type I particles with various r/t modeled using spring networks. The vertical dashed lines in (A) and (B) indicate the pressure used for the data in (C).

In Fig. 6, we plot k(p) from the spring network model for three values of a and four values of r/t in panels (A)–(D) with comparisons to the experimental and ABAQUS simulation data. We find that the spring network model with a = 0 is most similar to the experimental data for k(p) when considering all values of r/t for type I particles. Note that for linear spring networks with a = 1, the type I particles exhibit only increasing stiffness with increasing inflation pressure for all r/t. These results emphasize that strain softening is necessary to achieve decreasing compressive stiffness k with increasing p. a = 0.33 for the spring network model gives the same form for the strain-dependent Young's modulus Y(λ) for both the spring network model and neo-Hookean model implemented in ABAQUS, yet the results for k(p) for these two models are slightly different. The deviations are caused by the fact that ABAQUS enforces local incompressibility, whereas the spring network model does not.

Fig. 6 Compressive stiffness k for type I particles plotted versus the inflation pressure p, normalized by the stiffness k₀ at p = 0, at various fillet radius-wall thickness ratios r/t: (A) r/t = (r/t)_min, (B) 0.5, (C) 1, and (D) 2. Each panel shows results for the spring network model with a = 1 (blue triangles), 0.33 (magenta diamonds), and 0 (green left triangles), experimental measurements (black solid circles with dashed lines), and ABAQUS simulations with a neo-Hookean model (black open circles with dash-dotted lines). The error bars for k/k₀ are calculated using the standard deviation from N_r independent runs, where N_r = 9 for the experiments and 10 for the spring network simulations. The error bars for p in experiments are the precision Δp in the pressure gauge.

3.3. Compressive stiffness of type II particles

Type I particles expand both radially and along the cylinder axis during inflation, whereas type II particles mainly expand radially. Purely radial expansion is desirable for studies of quasi-2D packings of inflatable particles that are easily visualized during applied deformation. In this section, we determine the compressive stiffness k(p) for type II particles. In this case, four ratios, r/t, H/t, R/t, and R_a/t, specify the geometry of the particle. We experimentally measure the compressive stiffness versus the inflation pressure for different combinations of the initial values of r/t, H/t, R/t, and R_a/t.

We first demonstrate that scaling all of the geometrical parameters that describe the type II particle shape by the same factor does not change k(p)/k₀ significantly. Specifically, in Fig. 7A, there is <4% deviation in k(p)/k₀ over the full range in pressure when we vary r, t, H, R, and R_a, while keeping the ratios r/t = 1.5, H/t = 5, R/t = 10, and R_a/t = 6 unchanged.

Fig. 7 Relative stiffness k/k₀ as a function of the inflation pressure p (in kPa) for type II particles with different combinations of r/t, H/t, R/t, and R_a/t. In (A), we show k(p) when (circles) all parameters are scaled proportionally at fixed r/t = 1.5, as well as for (pluses) r/t = 1.0, (crosses) 0.5, and (diamonds) (r/t)_min all at fixed R/t = 10 and R_a/t = 6. In (B), we show (cyan) k(p) for R_a/t = 6 (circles), 5 (diamonds), 4 (leftward triangles), and 3 (rightward triangles) at fixed r/t = 1.5 and R/t = 10. We also show k(p) when two ratios are varied: R_a/t = 6 and R/t = 10 (blue pluses) and R_a/t = 4 and R/t = 6.67 (red squares) both at fixed r/t = 1.0. In (C), we show k(p) when all four ratios are varied. In all plots, each data point is an average over three samples measured three times each. The error bars indicate the standard deviation and the dotted lines are added as guides for the data.

Next, we vary one parameter, r/t, while keeping the other three parameters unchanged. Similar to the results for type I particles in Fig. 4A, the sign of the slope of k(p) is determined by r/t for type II particles as shown in Fig. 7A. At low pressures (⪅6 kPa), k(p) has a strong dependence on r/t. For large r/t (≥1), the slope of k(p) is positive; and for small r/t (≤0.5), the slope is negative. At larger pressures (≳7 kPa), k(p) reaches a plateau for most type II particles. We suspect that restricting type II particles to inflate only in the radial direction results in greater radial strains at the top and bottom regions compared to type I particles, leading to increased strain localization.

In Fig. 7B, we keep r/t fixed, while varying the rest of the dimensionless particle shape parameters. We find similar behavior for k(p) as that shown in Fig. 7A, i.e., k first increases with p and then reaches a plateau. If we vary R_a/t between 3 and 6, while keeping r/t, R/t, and H/t unchanged, we find no significant difference in k(p). We further observe that while keeping r/t = 1 fixed, varying the other three ratios alters the behavior of k(p). Although the slope remains positive in both cases, it decreases for smaller values of R_a/t, R/t, and H/t.

Finally, we examine the case where all four ratios, r/t, H/t, R/t, and R_a/t, are varied. As shown in Fig. 7C, k(p) increases with p at small p and the slope of k(p) at small p is mainly determined by r/t. However, the four dimensionless ratios that prescribe the initial shape of the particle are important for determining the p-dependent compressive stiffness, especially at large p. Similar results are found for the evolution of the particle radius R with p, as shown in Appendix B.

3.4. Packings of inflatable particles

As a proof of principle that our particle designs can enable adaptable granular packings, we fabricated an assembly composed of two types of type II inflatable particles. A schematic of the assembly is shown in Fig. 8A, and the stiffness k/k₀ versus pressure p for each particle type is shown in Fig. 8B. Under inflation, particle A (r/t = 3.0, R_a/t = 12, R/t = 20, H/t = 10) increases in stiffness, whereas particle B (r/t = (r/t)_min, R_a/t = 6, R/t = 10, H/t = 5) decreases in stiffness. Note that both particle types have the same absolute height H, ensuring the packing remains quasi-2D. We inflated the two particle types to different pressures to create two packings with distinct mechanical properties, while keeping the packing fraction fixed (see Fig. 8C). Fig. 8D shows that the two packings exhibit different responses under uniaxial compression. The force versus displacement curves are somewhat noisy at small displacements but become approximately linear at larger displacements. This behavior arises because the compression plate does not contact both particle types simultaneously at low displacement. We determine the compressive stiffness by fitting a line to the curves between 1 and 3 mm displacement. Packing I yields a slope of 144 ± 3 N m⁻¹, and Packing II yields a slope of 95 ± 2 N m⁻¹. This change corresponds to an increase in stiffness of over 50% when switching from the softer to the stiffer packing at a fixed packing fraction. These results demonstrate that adaptive granular assemblies composed of inflatable particles can achieve a wide range of bulk mechanical properties, e.g., uniaxial compressive stiffness, at constant packing fraction, a capability inaccessible with conventional particles of fixed stiffness.

Fig. 8 A tunable particle assembly with two different force responses at the same packing fraction (ϕ = 0.77 ± 0.01). (A) Schematic of a four-particle packing with two variations of type II particles: A (white) and B (cyan). (B) Normalized stiffness k/k₀ plotted versus pressure p for particles A (dashed grey line) and B (solid cyan line). Red triangles (green circles) label the state of particles A and B in packing I (packing II). (C) Images of packings I and II labelling the inflation pressure and inflated radius of each particle. (D) Force versus displacement for packings I and II. The three colors in the scatter plots for each packing correspond to three independent trials, and the dashed line is a linear fit to the data for displacements >1 mm.

4. Conclusions and future directions

In this article, we combined experimental and numerical studies to characterize the properties of two types of inflatable particles with initial cylindrical shapes. We investigated how the initial geometry of the particle affects the pressure dependence of the compressive stiffness k(p). The initial ratio of the fillet radius r to the shell thickness t plays a dominant role in determining k(p). We show that k can decrease, remain the same, or increase with the inflation pressure p, for both types of particles, for different values of r/t. In particular, k(p) decreases with pressure for small r/t(⪅0.8) and increases with pressure for r/t ≳ 0.8. We find that at small r/t, large strain localization near the “corners” of the cylindrical particles leads to strain softening, resulting in a decrease in compressive stiffness with increasing pressure.

This work introduces a novel method for fabricating particles with controllable size and stiffness. Our results suggest several promising future directions. First, further investigations of the dependence of k(p) on all dimensionless ratios that characterize the particle shape are necessary to describe the pressure-dependent stiffness, especially at large p. Second, it will be interesting to further study the structural and mechanical properties of packings of the inflatable particles presented in this work. For example, size and stiffness changes in the particle packings can change the topology of the interparticle contact networks, as well as the force distributions. Tuning the particle size and stiffness of even a single particle can have a significant impact on the contact networks, force networks, and mechanical properties of the packing. One potential application of this study is to use packings of inflatable particles to build quasi-2D robotic granular materials, whose size, shape, and stiffness can change independently and can be optimized to yield packings with specified and adaptable properties. Finally, frictional behavior in inflatable particle packings remains an open question. Prior studies^20,21 show that the friction coefficient between soft materials can vary with contact pressure. Similar effects likely arise in Ecoflex–Ecoflex contacts, potentially influencing tangential forces and history-dependent contact behavior. For instance, interparticle forces may differ depending on whether contact is formed before or after inflation. These complexities suggest that inflatable particle contacts may exhibit highly nonlinear and path-dependent dynamics.

This work also raises intriguing questions about the general behavior of stretchable and inflatable systems. It may be useful to conceptualize stress and strain distributions in stretchable sheets through a few key design parameters, which could, in turn, enable inverse design of complex inflatable structures. Balloon-like inflation systems are particularly interesting due to their non-monotonic pressure response to inflation, a phenomenon that has been extensively studied experimentally and theoretically.^25–29 However, the dependence of compressive stiffness on geometric parameters, as demonstrated in this work, represents a distinct effect occurring at smaller strains, separate from previously observed non-monotonic behaviors. Our findings suggest that further research could explore the role of geometry in governing bi-stability and non-monotonic responses in inflatable systems.

Inflatable structures play a central role in soft robotics, serving as actuators. However, current design strategies primarily focus on the final inflated shape rather than the internal stress distribution within the material. Stiffness modulation in soft robotic systems is typically achieved through decoupled mechanisms such as phase change or particle/layer jamming.^30–34 Our results suggest a coupled approach, where both shape transformation and stiffness modulation are integrated into a single system through thoughtful geometric design of inflatable components, opening new possibilities for adaptive and multifunctional soft robots.

Appendices

Appendix A

In this Appendix, we show how the friction coefficient μ between Ecoflex and the rigid, rough compression plates affects the compressive stiffness of type I particles in the ABAQUS simulations. In Fig. 9, we find that k/k₀ is insenstive to μ over the full range of pressure p. In particular, the largest change in the stiffness over this range of pressure satisfies k/k₀ < 2%.

Fig. 9 The normalized compressive stiffness k/k₀ from ABAQUS simulations of type I particles with a ratio of the fillet radius to the wall thickness r/t = 2 plotted versus the inflation pressure p (in kPa). The friction coefficient between the Ecoflex particles and steel walls is varied over the range 0.15 ≤ μ ≤ 1.

Appendix B

In this Appendix, we describe the relationship between the particle radius and inflation pressure p for all of the type II particles studied in the experiments in Section 7. Since type II particles mainly expand in the radial direction when inflated, we use the ratio R(p)/R to quantify the change in particle size. As shown in Fig. 10A, R(p)/R increases linearly with p for only one type II particle (r/t = 3, R_a/t = 12, R/t = 20, and H/t = 10). For the rest of the type II particles, R(p)/R is concave upward, consistent with strain softening. Fig. 10B shows this data plotted on a log–log scale. Most of the curves display two regimes: an approximately linear regime at lower inflation pressures and a second non-linear regime at higher inflation pressures. R(p)/R is is most strongly influenced by R/t and H/t, however, all four dimensionless particle shape parameters affect R(p), especially at large p.

Fig. 10 (A) Radius of the inflatable particles relative to the uninflated values (R(p)/R) plotted as a function of the inflation pressure p for the type II particles studied in Fig. 7 for different combinations of r/t, H/t, R/t, and R_a/t. Each data point is an average over three samples measured three times each. The error bars indicate the standard deviation and the dotted lines are added as guides to the data. (B) Data in (A) replotted on a log–log scale. The dashed and dotted lines have slopes 1 and 2, respectively.

Author contributions

N. P.: methodology, experimental fabrication and measuremnts, data analysis, writing, review and editing, visualization, and project administration. D. W.: methodology, numerical and finite element simulations, data analysis, writing, review and editing, visualization, and project administration. R. B.: finite element simulations, review and editing. M. G.: experimental fabrication, review and editing. M. D. S.: conceptualization, supervision, and resources. C. S. O.: conceptualization, supervision, project administration, review and editing, resources, and funding acquisition. R. K-B.: conceptualization, supervision, project administration, review and editing, resources, and funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data presented in this article, including data from experiments and FEA simulations are available on GitHub at https://github.com/nidhipashine/Inflatable-particles.

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). This work was also supported by the High Performance Computing facilities operated by Yale's Center for Research Computing.

References

1. R. P. Behringer, C. R. Phys., 2015, 16, 10–25.
2. A. J. Liu and S. R. Nagel, Annu. Rev. Condens. Matter Phys., 2010, 1, 347–369.
3. R. P. Behringer and B. Chakraborty, Rep. Prog. Phys., 2018, 82, 012601.
4. J. Zhang, D. Wang, W. Jin, A. Xia, N. Pashine, R. Kramer-Bottiglio, M. D. Shattuck and C. S. O'Hern, Phys. Rev. E, 2023, 108, 034901.
5. N. Gaspar, Mech. Mater., 2010, 42, 673–677.
6. M. Z. Miskin and H. M. Jaeger, Nat. Mater., 2013, 12, 326–331.
7. K. Fu, Z. Zhao and L. Jin, Adv. Funct. Mater., 2019, 29, 1901258.
8. N. Pashine, D. Wang, J. Zhang, S. K. Patiballa, S. Witthaus, M. D. Shattuck, C. S. O’Hern and R. Kramer-Bottiglio, Extreme Mech. Lett., 2023, 63, 102055.
9. H. Kocharyan and N. Karanjgaokar, Extreme Mech. Lett., 2023, 58, 101943.
10. D. Haver, D. Acuña, S. Janbaz, E. Lerner, G. Düring and C. Coulais, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2317915121.
11. Q. Wu, C. Cui, T. Bertrand, M. D. Shattuck and C. S. O'Hern, Phys. Rev. E, 2019, 99, 062901.
12. A. Parsa, D. Wang, C. S. O’Hern, M. D. Shattuck, R. Kramer-Bottiglio and J. Bongard, Applications of Evolutionary Computation, Cham, 2022, pp. 93–109.
13. S. Witthaus, A. Parsa, D. Wang, N. Pashine, J. Zhang, A. MacKeith, M. D. Shattuck, J. Bongard, C. S. O’Hern and R. Kramer-Bottiglio, Soft Matter, 2025, 21, 6088–6099.
14. M. S. Li, B. H. Do, C. L. Le, C. S. O'Hern and R. Kramer-Bottiglio, 2025 8th IEEE-RAS International Conference on Soft Robotics (RoboSoft), 2025.
15. A. N. Gent, Rubber Chem. Technol., 1958, 31, 896–906.
16. K. Larson, Can you estimate modulus from durometer hardness for silicones?, Dow Corning Corporation, 2016, 1–6.
17. L. R. G. Treloar, Trans. Faraday Soc., 1943, 39, 241–246.
18. Z. Cui, W. Wang, L. Guo, Z. Liu, P. Cai, Y. Cui, T. Wang, C. Wang, M. Zhu, Y. Zhou, W. Liu, Y. Zheng, G. Deng, C. Xu and X. Chen, Adv. Mater., 2022, 34, 2104078.
19. F. Pineda, F. Bottausci, B. Icard, L. Malaquin and Y. Fouillet, Microelectron. Eng., 2015, 144, 27–31.
20. Y. Yu, H. Yuk, G. A. Parada, Y. Wu, X. Liu, C. S. Nabzdyk, K. Youcef-Toumi, J. Zang and X. Zhao, Adv. Mater., 2019, 31, 1807101.
21. R. Berthold, J. Burgner-Kahrs, M. Wangenheim and S. Kahms, Meccanica, 2023, 58, 2165–2176.
22. L. A. Mihai and A. Goriely, Proc. R. Soc. A, 2017, 473, 20170607.
23. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
24. L. Marechal, P. Balland, L. Lindenroth, F. Petrou, C. Kontovounisios and F. Bello, Soft Robot., 2021, 8, 284–297.
25. R. S. Stein, J. Chem. Educ., 1958, 35, 203.
26. A. Needleman, Int. J. Solids Struct., 1977, 13, 409–421.
27. I. Müller and H. Struchtrup, Math. Mech. Solids, 2002, 7, 569–577.
28. G. Muhaxheri and C. D. Santangelo, Phys. Rev. E, 2024, 110, 024209.
29. E. Verron and G. Marckmann, Thin Wall Struct., 2003, 41, 731–746.
30. T. L. Buckner, M. C. Yuen, S. Y. Kim and R. Kramer-Bottiglio, Adv. Funct. Mater., 2019, 29, 1903368.
31. Y. Hao, J. Gao, Y. Lv and J. Liu, Adv. Funct. Mater., 2022, 32, 2201942.
32. E. Brown, N. Rodenberg, J. Amend, A. Mozeika, E. Steltz, M. R. Zakin, H. Lipson and H. M. Jaeger, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 18809–18814.
33. Y. S. Narang, J. J. Vlassak and R. D. Howe, Adv. Funct. Mater., 2018, 28, 1707136.
34. B. Yang, R. Baines, D. Shah, S. Patiballa, E. Thomas, M. Venkadesan and R. Kramer-Bottiglio, Sci. Adv., 2021, 7, eabh2073.

---

Footnote

† These authors contributed equally to this work.

---