Efficient Simulation of Deformable Particles

Deformable particles refer to macroscopic elastic particles that can undergo deformations such as stretching, bending, and twisting in response to external forces or interactions with other particles. Examples of deformable particles include foams [1], bubbles, emulsions [2], and biological cells and tissues [3, 4, 5]. Efficient simulation of deformable particles is crucial in various fields, including computer graphics, engineering, biological systems, and materials science. Unlike systems made up of hard, rigid particles where the arrangement of particles determines the mechanical properties, in soft granular systems, the deformability of the particles significantly influences the mechanical properties. Simulating deformable particles efficiently is complex and computationally intensive, especially with large numbers of particles. Many numerical studies on athermal solids simplify deformable particles as spheres interacting via purely repulsive pairwise potentials, such as Hertzian, Harmonic, or the repulsive part of the Lennard-Jones model [6, 7, 8, 9, 10, 11, 12]. These models often neglect the change in particle shape during interaction, considering the overlap ($\delta$) between particles as effective deformation. This simplification limits the applicability of soft potential models to real-life situations. Soft granular systems are often modeled by Hertzian potential, where the normal contact force between two elastic spherical beads is given by Hertz contact force, $F\propto\delta^{3/2}$ where $\delta$ is the overlap [13]. However, this approximation is only valid for small deformations within the elastic limit. Experiments on elastic homogeneous soft cylinders have shown different non-linear behavior at higher values of deformation [14].

There have also been several studies that try to incorporate the shape of the deformed particle into the model. For instance, studies on foams consider the shape of the foam during deformation [1, 15, 16], where the bubbles are generated using PLAT software [17]. However, the bubbles generated by PLAT show different structural properties compared to well-studied soft particle models. For example, in bubbles, the variation of the excess coordination number ($\delta z$) with excess packing fraction ($\delta\phi$) above jamming show linear behavior whereas soft particles models show a square-root variation [16]. There are also confluent models, such as the Vertex or Voronoi models, extensively used in biological systems such as epithelial tissues and the eye facets of Drosophila [18, 19, 20]. The non-confluent vertex model is also used to study the dynamics of embryonic tissues [21].

Only a few numerical approaches for modeling deformable particles have been developed over the past two decades, with most of them being based on either the discrete element method (DEM) or the finite element method (FEM)  [22, 23, 24, 25]. One such DEM model is the deformable polygon (DP) model, which is effective in both lower packing fractions as well as confluent systems [26]. The DP model accounts for particle shape during deformation and successfully reproduces the structural properties observed near the jamming transition [27]. It is widely used in the study of biological systems and soft granular packings [28, 29]. Even though it is an elegant model, it is computationally expensive as each deformable particle is defined by several degrees of freedom. The aim of this article is to reduce the computational cost of creating deformable particles with properties similar to those of the DP model using a geometric approach. Our model significantly reduces the number of degrees of freedom, making it more accessible for larger system sizes and potentially faster for both static packings through energy minimization and for studying dynamical properties at finite temperatures.

In Section 2, we present a brief overview of the DP model. In Section 3, we introduce a geometric method for determining the shape of a deformable particle by following certain assumptions on the deformed particle shape. In Section 4, we verify the accuracy of the analytical shape of a single deformable particle compressed by polygonal confinements by comparing it with an equivalent deformable polygon. In Section 5, we extend our methods to jammed packing multiple deformable particles in a rectangular confinement. We provide a comparison between the numerically generated packing of deformable polygons and the corresponding packing generated from the center of mass and radius of the particles using the geometric shape functions defined in this study. In Section 6, we provide a possible protocol for obtaining the force on each particle and particle position update. Finally, in Section 7, we conclude by suggesting possible limitations and future directions.

In this study, we derive a theoretical form for the shape function of a defomable particle and use the deformable polygon model as a template to verify our theoretical description. In the deformable polygon model [26, 27], the deformable particle is represented by a circular ring of soft beads connected by springs where any deviation from the circular shape has the following energetic cost:

where $k_{L}$ is the spring constant between two consecutive beads in the ring and $L^{0}_{m}$ is the rest length of the springs in the $m$-th ring whereas the $L_{m,i}$ is the $i$-th spring length in the $m$-th ring after the deformation. Similarly, $k_{A}$ is the spring constant associated with the area fluctuation with respect to the undeformed ring. $N_{v}$ is the number of beads (or vertices) on each ring and $N$ is the number of DP particles in the system. Since an undeformed ring is a regular polygon, we can obtain the rest length between two vertices and the area of the DP particle as $L_{m}^{0}=2R_{0}\sin{(\pi/N_{v})}$ and $A_{m}^{0}=N_{v}R_{0}^{2}\sin{(\pi/N_{v})}$ respectively. In a deformed ring, we can obtain $L_{m,i}$ by measuring the distance between the vertex $\vec{v}_{m,i}$ and $\vec{v}_{m,i+1}$ and $A_{m}$ using the Surveyor’s Area Formula (Shoelace Formula) [30]. Here $U_{int}$ is the interaction energy between two neighboring DP particles which is calculated by measuring the overlap between beads from the neighboring particles. Here each bead situated on the vertices of the ring is considered as a harmonically repulsive circular particle of radius $d$.

Let us consider a two dimensional system of circular deformable particles of radius $R_{m}$ where the center of mass of the undeformed particle can be represented by $\vec{r}_{m}^{\,cm}$. The locus of the outer perimeter of the deformable particle with respect to this center can be parameterised as

where $x_{m}^{\,cm},y_{m}^{\,cm}$ correspond to the Cartesian coordinates of the particle center in the undeformed state. The aim of this article is to derive the shape function $g(\theta)$ given the particle positions and their radii. When the particle is undeformed, it assumes a circular shape to minimize its area, resulting in $g_{m}(\theta)=R_{m}$. To determine the shape function in the deformed state, we make the following assumptions:

• 1.

  The contact region between two neighboring particles is a straight line segment ($l$) whose length depends on the distance between the neighbors as well as the magnitudes of the parameters $k_{L}$ (perimeter stiffness) and $k_{A}$ (area stiffness).

• 2.

  In the area between two contact regions, the particle forms a circular arc that is tangent to both contact line segments.

• 3.

  We consider fully deformable particles with no overlap between neighboring particles. Consequently, the total potential energy of the system is determined solely by changes in the particle shapes, excluding any interaction potential $U_{int}$.

• 4.

  Instead of a polygon, we treat each particle as a continuous deformable ring. Thus, the total potential energy can be expressed as a sum of an area contribution and a perimeter contribution:

  		$\displaystyle U=\sum_{m=1}^{N}(U^{m}_{P}+U^{m}_{A}),$		(3)
  		$\displaystyle U^{m}_{P}=\frac{k_{P}}{2}(P_{m}-P^{0}_{m})^{2},\text{ and, }U^{m}_{A}=\frac{k_{A}}{2}(A_{m}-A^{0}_{m})^{2},$

  where $A^{0}_{m}=\pi R_{m}^{2}$ and $P_{m}^{0}=2\pi R_{m}$ are the initial area and perimeter of the undeformed particle, respectively, and $k_{P}=k_{L}/N_{v}$. The perimeter and area of the deformed particle can be expressed as:

  	$\displaystyle P_{m}=$	$\displaystyle\frac{1}{2}\int_{0}^{2\pi}\sqrt{g_{m}^{2}+\left(\frac{\partial g_{m}}{\partial\theta}\right)^{2}}\mathrm{d}\theta,$		(4)
  	$\displaystyle A_{m}=$	$\displaystyle\frac{1}{2}\int_{0}^{2\pi}g_{m}^{2}\text{d}\theta.$

• 5.

  We assume that the deformable particles do not slip on top of each other to mimic the behaviour of the frictional deformable polygons.

Under these constraints, we consider a single deformable particle of undeformed radius $R_{0}$ ($m=0$) situated at $\vec{r}_{0}^{\,cm}$ which is surrounded by $z$ neighbors with radii $R_{i}$ and positioned at $\vec{r}_{i}^{\,cm}$ where $i$ ranges from $1$ to $z$. The condition for a particle $i$ to be a neighbor can be represented as: $\left|\vec{r}_{i}^{\,cm}-\vec{r}_{0}^{\,cm}\right|\leq(R_{i}+R_{0})$. The distance between the particle centre $\vec{r}_{0}^{\,cm}$ and the contact surface with the $i$-th neighbor can be written as $R_{0}-a_{i}$ where,

The angle between two consecutive contacts ($i$ and $i+1$) (see Fig. 2) can be written as,

Let $l_{i}$ be the part of the perimeter of the particle that is in contact with the $i$-th neighbor. Then the region between two such consecutive contacts can be represented by two line segments of length $l_{i}$, $l_{i+1}$ and a circular arc with a radius $q_{i}$ and an arclength of $\phi_{i}q_{i}$. To derive a relation between these quantities we compare the energy associated with area fluctuation to the energy associated with the perimeter fluctuation between two such contact regions. For a fixed value of $k_{P}$ and $k_{A}$, we can write the energy change due to perimeter fluctuation and area fluctuation between two contact regions as

In our study, we consider $k_{P}\gg k_{A}$, which implies a very high energy cost for any change in the particle perimeter compared to energy cost for change in the particle area. Thus, the major contribution to deformation comes from the area fluctuation. Considering the non-slip contact constraint, we can compare the change in the energy due to area ($U_{P}^{i,i+1}$) and perimeter ($U_{A}^{i,i+1}$) change between two contacts for $k_{A}/k_{P}\to 0$, which gives

This relation indicates that the perimeter of the deformable particle between two contacts remains constant before and after the deformation. We can also obtain an additional relationship between $l_{i},l_{i+1}$ and $q_{i}$ that is independent of Eq. (8), by using the geometric constraint that the circular arc is tangent to both the contact line segments (see Fig. 2), which can be represented as:

A detailed derivation of the above geometric relations are given in the A. Next, using the perimeter constraint in Eq. (8) and the geometric constrain in Eq. (9), we can obtain,

The position of the center of the circular arc with respect to the particle center $\vec{r}_{0}^{\,cm}$ can be represented as

Next using the formulae given in Eqs. (5), (6),  (10) and  (11), we can obtain the functional form of the contour of the deformed particle situated at $\vec{r}_{0}^{\,cm}$ in polar coordinate system and the angular variation of the $g_{j}(\theta)$ can be written as:

where $H$ corresponds to the Heaviside theta function. $\theta_{i}^{s}$ and $\theta_{i}^{e}$ correspond to the start and end polar angle of the circular arc region with respect to the particle center $\vec{r}_{0}^{\,cm}$ and can be represented as

In an undeformed particle, $a_{i}=l_{i}=p_{i}=0$, which makes $g_{0}(\theta)=q_{i}=R_{0}$. Here, the particle shape is determined using the constraint that the perimeter of the particle is always constant at any level of deformation. However, different types of constraints can also be used, such as fixed area deformable particles or mixed situations where both the perimeter and the area are allowed to fluctuate. For instance, considering deformable particles where $K_{A}\gg K_{P}$, we can obtain a different constraint on the particle shape i.e.

The above constraint corresponds to a situation where the particle can change its shape to fit within the confines of the surrounding geometry while maintaing a constant area during deformation. Therefore, by using the constraint for constant area of the particle as given in Eq. (14) and the geometrical constraints for the contact region and the circular arc described in Eq. (9), we can obtain the parameters associated with the shape of the deformed particle. For example, we can obtain the lengths of the contact regions within an edge between two consecutive contacts $i$ and $i+1$ as,

where $a_{i},\phi_{i}$ can be obtained using Eq. (5) and Eq. (6). For any values of $k_{P}$ and $k_{A}$, the given functional forms of $a_{i},\phi_{i},p_{i},\theta_{i}^{p}$ are all same. The only differece lies in the contact length and the radius of the circular arc which varies with the different stiffness constants.

To verify the validity of the shape function defined in Eq. (12), we perform simulations with a single deformable particle modeled using a spring-mass ring system, as detailed in Section 2. We place this particle within a closed confinement shaped as a $n$-sided polygon (with $n$ chosen to be 3, 4, 5, or 6). Initially, the particle is positioned at the center of the confinement, with the distance from the center to each side of the polygon greater than the radius of the particle. We then isotropically compress the polygonal box. By the end of the compression, the distance from the box sides to the center is less than the radius of the particle, causing the its shape to deform gradually. Here we assume non-slip boundary to mimic interaction between two deformable polygons i.e. the DP is not allowed to slip on the surface of the walls. The final shape of the deformable particle is achieved by minimizing the total shape energy using the FIRE algorithm [31]. Throughout the simulation, we maintain the spring constant $k_{A}$ at $1$. We vary the spring constant $k_{L}$ and the degree of overcompression, denoted by $a$. For each value of $k_{L}$, we incrementally increase the compression by steps of $0.005$ (i.e., $a/R_{0}=0.005,0.01,...,0.05$), and minimize the energy at each step to obtain the final shape of the particle. The energy-minimized shape in the deformable polygon (DP) model is represented by the bead positions $\vec{v}_{0,i}$. In Fig. 3, we have presented the shape of a DP particle compressed by polygonal boxes with different numbers of sides. We have also plotted the theoretical shape of these particles defined in the following paragraph on top of the DP particles (denoted by black lines in Fig. 3).

Below we give the expression for the shape of the deformed particle in a regular polygonal confinement. Given the center of mass of the particle, we express the shape function as outlined in Eq. (12). This shape function is designed to describe the deformed geometry of the particle as it adapts to the constraints of the surrounding confinement. Let $R^{\prime}$ denote the perpendicular distance from any side of the regular polygonal confinement to the center of the polygon, and $R_{0}$ denote the radius of the particle when it is in an undeformed, relaxed state. As we compress the confinement, the distance $R^{\prime}$ becomes smaller than $R_{0}$, forcing the particle to deform to fit within the new boundaries. To understand how the particle shape changes in response to this compression, we leverage the discrete rotational symmetry of the confinement. Regular polygons, by definition, have sides and angles that are uniformly distributed. For instance, in a regular $n$-sided polygon, the angle between any two adjacent sides is $2\pi/n$. This symmetry allows us to focus on a single section of the particle, knowing that the behavior in this section will be replicated around the entire perimeter. Using this symmetry we define the contact regions between the particle and the confinement as:

Similarly, other parameters such as the position and radius of the deformed circular arc can also be obtained as

Now, using the parameters defined in Eq. (16) and Eq. (17), we can define the shape of the deformed particle using Eq. (12). Similarly, the expression for the shape of a deformed particle in an irregular confinement can also be obtained. In Fig. 4 (a), we have shown the deformed particle obtained from simulation (represented by red points) and compared it with the perimeter-constrained geometric shape (represented by the black dashed line) for $k_{L}=10^{6}$, $N_{v}=314$ and overcompression, $a/R_{0}=0.02$. For a given bead position $\vec{v}_{0,i}$, we calculate the corresponding angle $\theta_{i}$ using the equation:

This angle $\theta_{i}$ allows us to compare the simulated shape with the theoretical shape. Next, we measure the difference between the deformable polygon (DP) model and the theoretical shape as $\left|\vec{v}_{0,i}\right|-g_{0}(\theta_{i})$. In Fig. 4 (b), we have shown the difference between two models with increasing magnitude of $k_{L}$. We observe that as $k_{L}$ increases, the difference between the two models decreases, indicating that the theoretical model becomes a better approximation of the actual deformed shape. To numerically quantify this, we measure the net difference between the two shapes using the summation:

We measure $I$ for different values $k_{L}$ and observe that at higher values of compression, the analytical shape deviates from the DP shape more. However, $I$ decreases monotonically with $k_{L}$ for all the values of overcompression $a$. This suggest that for $k_{A}/k_{L}\to 0$, $I\to 0$ for all $a/R_{0}$. This monotonic decrease in $I$ implies that the theoretical shape function $g_{0}(\theta)$, which is derived by considering the perimeter of the particle to be fixed, accurately predicts the shape of the deformable polygon when the stiffness ratio $k_{A}/k_{L}$ is very small. Therefore, for sufficiently large $k_{L}$, the analytical model $g_{0}(\theta)$ becomes an excellent approximation of the shape of a deformable particle confined in a polygonal box. We also measure the variation in $I$ with increasing number of beads on the DP for a fixed value of perimeter spring constant $k_{P}$ which is presented in Fig. 3 (e). We observe that the $I$ decreases with increasing $N_{v}$ and eventually saturates to a fixed value for any levels of overcompression $a/R_{0}$. This saturation point is dependent on $k_{P}$, with higher values of $k_{P}$ leading to a lower saturation value of $I$. This outcome further validates the accuracy of the geometric method used to predict the shape of the DP, especially in the limit where $N_{v}\to\infty$ and $k_{P}\to\infty$ both approach infinity.

To check the validity of this geometrical method for jammed packing of deformable polygons in a rectangular box, we consider a bidisperse system of soft particles where the ratio of the radius between the larger and smaller particles is $R^{l}/R^{s}=1.4$. We begin by distributing these soft particles randomly within a large square box. Subsequently, we compress the box isotropically while minimizing the total energy of the system using the FIRE algorithm. The compression is halted when the system becomes jammed, which is characterized by a finite positive stress. At this point, the particles are in close contact, and no further compression can occur without significantly increasing the energy of the system. Next, we replace each soft particle with a $N_{v}$-sided deformable polygon. The number of vertices for each particle is chosen based on its size i.e. the larger particle is assigned a greater number of vertices compared to the smaller particle, such that:

where $N_{v}^{l}$ denotes the number of vertices on the larger particle, $N_{v}^{s}$ denotes the number of vertices on the smaller particle, $R^{l}$ is the radius of the larger particle, $R^{s}$ is the radius of the smaller particle, and $p$ is the radius of each bead on the vertices of the deformable polygons. We then minimize the energy of this packing where the energy is defined by Eq. (1).

Next, we isotropically compress the system with an incremental compressive strain of magnitude $\delta\epsilon$. At each strain increment, the system is relaxed to its energy minimum, ensuring the protocol can be considered quasistatic. At the jamming point, where the particles are in their undeformed states, each polygon can be approximated as a circular particle. However, as we further compress the system, the relaxation of the net energy occurs through two primary mechanisms: (1) particle rearrangements and (2) deformation of the particles. For this study, we have taken a system of $N=128$ deformable particles where the radii of smaller and larger particles are $R^{l}=0.7$ and $R^{s}=0.5$ respectively. We choose the radius of each bead on the vertex of the deformable polygons as $p=0.05$ resulting in a ratio of vertices $N_{v}^{l}/N_{v}^{s}=1.44$. Specifically, we select $N_{v}^{l}=72$ vertices for the larger particles and $N_{v}^{s}=50$ vertices for the smaller particles.

Next, we extend the geometrical method described in Section 3 to multiparticle jammed packings within a rectangular box. Similar to a single particle in confinement, the contact region between two particles acts as a boundary. As we compress the box, the average distance between particles decreases, effectively bringing the boundary closer to the particle centers. We implement techniques detailed in Section 3 to obtain the shape of the individual deformable particle in the system. However, the deformed particle shape can also overlap with other neighboring deformed particles which we refer as secondary overlap. These secondary overlaps are treated similarly to primary overlaps between two circles with the radii of the circles being the radii of the circular arc. The same procedure is followed to resolve secondary deformations, and this process can be iterated for higher-order deformations. For most configurations the final deformed shape can be obtained just performing three iterations. In Fig. 5, we have shown a jammed packing of deformable polygons (red circles) for an istropic strain above the jamming with an amplitude, $\epsilon=0.01$. From this packing, we obtain the center of mass and radius of each particle. Using the extracted parameters and the constrained geometric shape method described by Eq. (12), we generate the shapes of the particles in the packing. These shapes are represented by black solid lines in Fig. 5.

Using the functional form $g_{0}(\theta)$ in Eq. (12), we can obtain the area of the deformed particle as:

Using the above formula for area we can determine the total energy ($U$) of the packing using Eq. (3). We can also determine the force on each deformable particle by taking the positional derivative ($-\nabla U$) of the total potential energy i.e.

Since the potential energy only contains the area fluctuation term, the final positions of the particles in a jammed packing are highly dependent on their deformed area. However, other terms related to perimeter fluctuation and interaction between the neighbors can be included to force for higher values of $k_{P}$ and $k_{r}$.

The above formulation is invalid when the particle has only one neighbor as an edge can only be defined when coordination number $z\geq 2$. In such a scenario, we consider both the particles repel each other harmonically without deformation just like in soft particle systems via potential: $U_{ij}=\frac{k_{A}}{2}\left(1-\frac{r_{ij}^{\,cm}}{R_{i}+R_{j}}\right)^{2}$. We displace the particle in the direction normal to the contact in using the force $-\partial_{r_{i}}U_{ij}$ until it forms a new contact or the present contact breaks to form a rattler. Given the forces, we can update the particle positions at each time step using a verlet algorithm for dynamical systems and can achieve energy minimized packing through damped dynamics/ FIRE minimization.

In this study, we have proposed an efficient method to obtain jammed packings of deformable particles at higher packing fractions. Similar to soft particle models, our geometric formulation only requires the particle positions and their radii to determine the shape of individual particles after deformation, which in turn can be used to determine the total force acting on each particle. In the deformable polygon (DP) model, the total number of degrees of freedom in the system is $N_{v}N$, whereas our proposed model has only $N$ degrees of freedom for the same packing. This reduction in complexity makes our model more efficient and scalable, enabling us to create larger system packings more effectively. By employing this geometric method to calculate the force between two deformable particles, we can develop a faster algorithm for generating static packings through energy minimization as well as for studying dynamical properties at finite temperatures.

One significant limitation of the proposed method arises when dealing with systems composed of multiple deformable particles with varying degrees of deformability. Specifically, if the particles in the system have different spring constants $k_{P}$ and $k_{A}$. The primary issue stems from the assumption that the contact regions between particles are perfectly flat. In reality, when deformable particles with different spring constants interact, the contact regions can exhibit complex shapes, including concave or convex curvatures, depending on the relative deformability of the particles. This geometric method, in its current form, does not account for these complex contact geometries, potentially leading to inaccuracies when predicting the behavior of systems with heterogeneous deformability. Despite this limitation, the proposed geometric approach can be extended to incorporate the shape of the contact regions, taking into account the specific deformabilities of the interacting particles which would enhance the accuracy of the method.

Using this geometric method, various structural and mechanical properties of the jammed packings can be studied and verified. For example, the average coordination number can be analyzed to understand the connectivity and stability of the packing. Additionally, the elastic moduli, which describe the stiffness of the material, can be measured above the jamming transition to explore how the mechanical response of the system changes with increased packing density. Our model also allows for direct studies on the effect of particle deformability on order-disorder transitions and amorphization transitions. Understanding these transitions is key to exploring how systems evolve from ordered to disordered states and how they become amorphous under varying conditions.

I would like to thank Kabir Ramola, Stephy Jose, Surajit Chakraborty, Rashmiranjan Bhutia, and Shilpa Prakash for useful discussions. This project was funded by intramural funds at TIFR Hyderabad from the Department of Atomic Energy (DAE). We acknowledge support of the Department of Atomic Energy, Government of India, under Project Identification No. RTI4007.