The role of deformability in determining the structural and mechanical properties of bubbles and emulsions

TL;DR

This study compares the structural and mechanical properties of jammed particle packings using soft particle and deformable particle models, revealing how deformability influences jamming behavior and void space in emulsions and foams.

Contribution

It introduces a deformable particle model for compressed bubbles and emulsions and compares its predictions with the soft particle model, highlighting differences in jamming and void structure.

Findings

Both models show $ ho$-scaling of contact number and shear modulus near jamming.

Deformability affects void space and Plateau border formation.

Results align with experimental data on foams and emulsions.

Abstract

We perform computational studies of jammed particle packings in two dimensions undergoing isotropic compression using the well-characterized soft particle (SP) model and the deformable particle (DP) model that we developed for compressed bubbles and emulsions. In the SP model, circular particles are allowed to overlap, generating purely repulsive forces. In the DP model, particles minimize their perimeter, while deforming at fixed area to avoid overlap during compression. We directly compare the structural and mechanical properties of jammed particle packings generated using the SP and DP models as a function of the true packing fraction $\rho$, instead of the reduced number density $\phi$. We show that near jamming onset the excess contact number $\Delta z=z-z_J$ and shear modulus ${\cal G}$ scale as $\Delta \rho^{0.5}$ in the large system limit for both the SP and DP models, where $\Delta \rho = \rho-\rho_J$ and $z_J \approx 4$ and $\rho_J \approx 0.842$ are the values at jamming onset. $\Delta z$ and ${\cal G}$ for the SP and DP models begin to differ for $\rho \gtrsim 0.88$. In this regime, $\Delta z \sim {\cal G}$ can be described by a sum of two power-laws in $\Delta \rho$, i.e. $\Delta z \sim {\cal G} \sim C_0\Delta \rho^{0.5} +C_1\Delta \rho^{1.0}$ to lowest order. We show that the ratio $C_1/C_0$ is much larger for the DP model compared to to that for the SP model. We also characterize the void space in jammed packings as a function of $\rho$. We find that, unlike the SP model, the DP model is able to describe the formation of Plateau borders as the system approaches $\rho = 1$. We further show that the results for $z$ and the shape factor ${\cal A}$ versus $\rho$ for the DP model agree with recent experimental studies of compressed foams and emulsions.

Figures

Tables

Table 1. Table 1: Parameters for the scaling forms for the excess contact number Δ ​ z ​ ( Δ ​ ϕ ) Δ 𝑧 Δ italic-ϕ \Delta z(\Delta\phi) and Δ ​ z ​ ( Δ ​ ρ ) Δ 𝑧 Δ 𝜌 \Delta z(\Delta\rho) (Eqs.  1 and  13 ) for the SP and DP models.

Model	$z_J\pm 0.005$	${\varphi }_J,{\rho }_J\pm 0.005$	$z_0^{\varphi ,\rho }$	${\alpha }_0^{\varphi ,\rho }$	$z_1^{\varphi ,\rho }$	${\alpha }_1^{\varphi ,\rho }$
SP($\varphi$)	3.97	0.84	3.7	0.5	0	1.0
SP($\rho$)	3.97	0.84	2.6	0.5	4.0	1.0
DP($\rho$)	3.97	0.835	3.3	0.5	7.1	1.0

Table 2. Table 2: Parameters for the scaling forms for the shear modulus 𝒢 ​ ( Δ ​ ϕ ) 𝒢 Δ italic-ϕ {\cal G}(\Delta\phi) and 𝒢 ​ ( Δ ​ ρ ) 𝒢 Δ 𝜌 {\cal G}(\Delta\rho) (Eqs.  2 and  14 ) for SP and DP models.

Model	${\varphi }_J,{\rho }_J\pm 0.005$	${𝒢}_0^{\varphi ,\rho }$	${\beta }_0^{\varphi ,\rho }$	${𝒢}_1^{\varphi ,\rho }$	${\beta }_1^{\varphi ,\rho }$
SP($\varphi$)	0.84	0.42	0.5	0	1.0
SP($\rho$)	0.84	0.38	0.5	0.35	1.0
DP($\rho$)	0.835	0.07	0.5	1.7	1.0

Equations

$$z-z_J\approx z_0^{\varphi }(\varphi -{\varphi }_J{)}^{{\alpha }_0^{\varphi }},$$

$$\mathcal{G}\approx {\mathcal{G}}_0^{\varphi }(\varphi -{\varphi }_J{)}^{{\beta }_0^{\varphi }},$$

$$z-z_J\sim (\rho -{\rho }_J{)}^{{\alpha }^{\prime }},$$

$$\mathcal{G}\sim (\rho -{\rho }_J{)}^{{\beta }^{\prime }},$$

$$U_{SP}(r_{ij})$$

$${\Sigma }_{\mu \nu }=A^{-1}\sum _{i>j}^N{f}_{ij\mu }{r}_{ij\nu },$$

$$U_{DP}$$

$$U_{\mathrm{int}}$$

$$\rho =\sum _{i=1}^N\frac{\pi {\sigma }_i^2}{4A}-\frac{1}{A}\sum _{i>j}^N{a}_{ij}^{\mathrm{ov}},$$

$$a_{ij}^{\mathrm{ov}}=\frac{\sqrt{(-r_{ij}+{\sigma }_{ij})(r_{ij}-{\overline{\sigma }}_{ij})(r_{ij}+{\overline{\sigma }}_{ij})(r_{ij}+{\sigma }_{ij})}}{2}$$

$$\Delta \rho \approx C(\Delta \varphi -\Delta {\varphi }^{1.5}),$$

$$\rho =\sum _{i=1}^N\frac{\pi {\sigma }_i^2}{4A}-\frac{1}{A}\sum _{i>j}^N{a}_{ij}^{\mathrm{ov}}+\frac{2}{A}\sum _{k>i,j}^N{a}_{ijk}^{\mathrm{ov}},$$

$$z-z_J\approx z_0^{\rho }(\rho -{\rho }_J{)}^{{\alpha }_0^{\rho }}+z_1^{\rho }(\rho -{\rho }_J{)}^{{\alpha }_1^{\rho }}$$

$$\mathcal{G}\approx {\mathcal{G}}_0^{\rho }(\rho -{\rho }_J{)}^{{\beta }_0^{\rho }}+{\mathcal{G}}_1^{\rho }(\rho -{\rho }_J{)}^{{\beta }_1^{\rho }}$$

Full text

\makeFNbottom

[TABLE]

††footnotetext: a* Department of Mechanical Engineering and Materials Science, Yale University, New Haven, Connecticut, 06520, USA; E-mail: arman.boromand@yale.edu. ††footnotetext: b Beijing National Laboratory for Condensed Matter Physics and CAS Key Laboratory of Soft Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing, China.††footnotetext: c Program in Computational Biology and Bioinformatics, Yale University, New Haven, Connecticut, 06520, USA*††footnotetext: d* Department of Physics, Emory University, Atlanta, Georgia 30322, USA.††footnotetext: e School of Physical Sciences, University of Chinese Academy of Sciences, Beijing, China.††footnotetext: f Benjamin Levich Institute and Physics Department, The City College of New York, New York, New York 10031, USA.††footnotetext: g Department of Physics, Yale University, New Haven, Connecticut, 06520, USA. ††footnotetext: h Department of Applied Physics, Yale University, New Haven, Connecticut, 06520, USA; E-mail: corey.ohern@yale.edu. *

1 Introduction

Soft materials, such as grafted core-shell particles, dendrimers, star polymers, emulsions, foams, and hydrogels, are a class of materials for which their microstructure can be altered by external fields, applied deformation, and thermal fluctuations at room temperature.. 1, 2 The ability to vary particle microstructure enables the design of soft materials with novel functional properties and processing capabilities. Molecular architecture, surface interactions, and deformability of soft particles can be harnessed to develop novel soft composites with optimized energy absorption, self-healing behavior, high mechanical strength, and other desirable properties.3, 4, 5, 6 In addition, many biological systems such as biofilms 7, cell aggregates 8, and tissues 9 can be considered as collections of soft and deformable particles.

The interactions between soft particles, e.g. the softness, range, and strength of the attraction and repulsion between soft particles is controlled by their composition and microstructure. In turn, the interactions between soft particles determine the collective mechanical and rheological properties of packings of soft particles. Significant challenges remain in understanding the influence of particle microstructure and interactions on the macroscopic properties of soft matter systems. In this article, we study the role of particle deformability in determining the structural and mechanical properties of packings of quasi-2D emulsions, modeled as collections of purely repulsive, deformable particles at and above the jamming transition.

Systems composed of soft, frictionless particles, such as foams and emulsions, can jam, or develop a non-zero static shear modulus ${\cal G}$, when they are isotropically compressed to packing fractions that approach random close packing $\phi_{J}$. 10, 11, 12 For $\phi<\phi_{J}$, packings of purely repulsive, frictionless spherical particles have an insufficient number of interparticle contacts for them to be mechanically stable. As a result, the packings exist at zero pressure ($p=0$) and are fluid-like, and particle rearrangements cost zero energy. 13 When compressed to $\phi_{J}$, the packings develop a connected interparticle contact network with an isostatic number of contacts per particle $z_{J}=N_{c}/N$, where $N_{c}=2N^{\prime}-1$, $N^{\prime}=N-N_{r}$, $N$ is the total number of particles, and $N_{r}$ is the number of rattler particles with less than $2$ contacts, for frictionless circular particles in two spatial dimensions (2D) with periodic boundary conditions. 14, 15 For $\phi>\phi_{J}$, the packings become solid-like with $z>z_{J}$ and a nonzero shear modulus ${\cal G}>0$. 16, 15, 17

A number of computational studies have been performed to investigate the structural and mechanical properties of compressed foams and emulsions in the solid-like regime $\phi>\phi_{J}$. 18, 19, 20, 21 One of the most frequently used models for characterizing their structural and mechanical properties is the soft particle (SP) model 18, 16, for which there is a potential energy cost proportional to the square of the overlap between pairs of spherical particles, and no energy cost when the particles do not overlap. These studies find that the excess contact number above the isostatic number,

[TABLE]

and the shear modulus,

[TABLE]

obey power-law scaling relations with $\phi-\phi_{J}$, where the scaling exponents $\alpha_{0}^{\phi}=\beta_{0}^{\phi}=0.5$ in the large-system limit. 16, 15, 22 Further, these studies have found that the exponents do not vary with the shape of the purely repulsive interaction potential and are the same in 2D and 3D. 15, 22, 19, 17.

In Eqs. 1 and 2, we define the packing fraction (or reduced number density) in 2D for packings of circular disks as $\phi=\sum_{i=1}^{N}\pi\sigma_{i}^{2}/4A$, where $\sigma_{i}$ is the diameter of disk $i$ and $A=L_{x}L_{y}$ is the area of the simulation box with edge lengths $L_{x}=L_{y}$ in the $x$- and $y$-directions. Note that when using this definition of $\phi$, the area of particle overlaps is multiply-counted. A positive feature of the SP model is its simplicity, however, a negative aspect is that the particles do not conserve area when the packings are compressed above jamming onset.

There have also been a number of experimental 23, 24, 25, 26 studies (as well as computational studies 20, 21) of the structural and mechanical properties of compressed foams and emulsions in 2D and 3D. These studies also find power-law scaling of the excess contact number,

[TABLE]

where $\alpha^{\prime}$ occurs in the range $0.4$ to $1$ depending on the particular study. The power-law scaling of $z-z_{J}$ is measured versus $\rho-\rho_{J}$, not $\phi-\phi_{J}$, where $\rho$ is the true packing fraction of the system. For systems at jamming onset, $\rho_{J}=\phi_{J}$, but $\phi>\rho$ when particles overlap in the SP model. Also, $\phi\geq 1.0$ is allowed, whereas $\rho\leq 1.0$ is a hard constraint. In addition, experimental studies of compressed emulsions in 3D have shown that the shear modulus obeys power-law scaling in $\rho-\rho_{J}$,

[TABLE]

but the scaling exponent $\beta^{\prime}>\beta^{\phi}_{0}$. 12, 27, 28

In light of the discrepancies between the power-law scaling exponents found in the experimental studies of compressed foams and emulsions and those obtained from the computational studies of the SP model, we employ the recently developed deformable particle (DP) model 29for foams and emulsions to understand how particle deformability affects the packing fraction dependence of the structural and mechanical properties of jammed particle packings.

Other computational methods have been employed to model particle deformability in soft matter systems, such as foams and emulsions. There are two main classes of methods for modeling particle deformability: Lattice Boltzmann 30, 31 and particle-based methods 20, 32, 33, 34, 35. The lattice-based methods have typically focused on two- or multi-phase modeling, whereas the DP model focuses only on the shape degrees of freedom of the particles (i.e. bubbles or droplets). Our work on the DP model differs from the previous studies using particle-based methods. First, the work by Rognon, et al. 35 has been limited to small systems composed of $2-5$ particles. Second, the froth model by Kern, et al. 33 is limited to the dry regime, where the true packing fraction approaches unity. In contrast, the DP model can be used to study a wide range of packing fractions, from values where the particles are out of contact to confluent systems. The model proposed by Kahara et al. 32 is most similar to the DP model. However, by modeling the pressure of the carrier fluid, their study is limited to the wet regime. In addition, their studies have focused on the rheological properties of bubbles during shear. In contrast, this article will focus on the structural and mechanical properties of jammed deformable particles generated during isotropic compression.

A central assumption of the SP model is that the particles remain spherical as particle overlap increases when the system is compressed above jamming onset. For this reason, studies that employ the SP model typically quantify the system properties as a function of $\phi$ instead of the true packing fraction $\rho$ 23, 24.

Even though the SP model does not conserve particle area as the system is compressed, one can also measure the structural and mechanical properties as a function of the true packing fraction $\rho$ when using the SP model by attributing half of an overlap between particles $i$ and $j$ to particle $i$ and the other half to particle $j$. (See Fig. 1.)

A key feature in jammed packings of foams and emulsions is that the particles maintain their area (volume) during compression over the full range of packing fraction. Bubbles in foams and emulsion droplets can deform, become non-spherical, and form additional contacts that do not occur in the SP model at a comparable value of packing fraction as that shown in Fig. 2 (a) and (b). In this article, we show that the soft particle and deformable particle models show similar results for the scaling of the excess contact number and shear modulus versus $\rho-\rho_{J}$ for packing fractions close to jamming onset. However, for larger $\rho$, we find that $z-z_{J}$ and ${\cal G}$ for the SP and DP models begin to differ significantly. In this regime, $z(\rho)$ for the deformable particle model is similar to that found for experimental studies of compressed emulsions and foams in 2D. We also study the geometric properties of the void space of jammed packings as a function of $\rho$. We show that unlike the SP model, the DP model is able to recapitulate the formation of Plateau borders 36, 37, 38, 39, where bubble edges have a relative orientation of $120^{\circ}$ and form a void with shape factor ${\cal A}\approx 4.87$, near confluence. (The shape factor ${\cal A}=p^{2}/4\pi a$, where $p$ is the perimeter and $a$ is the area of the void space 29, 26.)

The remainder of the article is organized as follows. In Sec. 2, we describe the soft particle and deformable particle models for 2D compresed foams and emulations, and the isotropic compression protocol that we employ to numerically generate jammed packings. In Sec. 3, we compare the results for the structural and mechanical properties of jammed packings using the SP and DP models. We show the variation of $\rho$ with $\phi$ above jamming onset and quantify the decrease in the average area of the particles as a function of increasing packing fraction for the SP model. We measure $\rho$ versus the shape factor ${\cal A}$ to determine at what shape factor the SP and DP models reach confluence. We then show the power-law scaling results for $z-z_{J}$ and ${\cal G}$ versus $\rho-\rho_{J}$ for the SP and DP models. We also characterize the connected void regions, by measuring the number, size, and shape of the voids as a function of packing fraction for the SP and DP models. In Sec. 4, we summarize the results for the SP and DP models, compare the results for the DP model to those from recent experiments on compressed foams and emulsions, and discuss future research directions.

2 Simulation Methods

In this article, we study the structural and mechanical properties of isotropically compressed jammed packings of $N$ purely repulsive, frictionless bidisperse particles in 2D using the soft particle and deformable particle models. For both models, the simulation cell is square with periodic boundaries in both directions.

2.1 Soft Particle Model

For the soft particle model, pairs of circular disks $i$ and $j$ interact via the purely repulsive pairwise potential:

[TABLE]

where $\epsilon$ is the characteristic energy scale of the interaction, $\sigma_{ij}$=$(\sigma_{i}+\sigma_{j})/2$ is the average diameter and $r_{ij}$ is the center-to-center separation between disks $i$ and $j$, and $\Theta(.)$ is the Heaviside step function that sets the interaction potential to zero when disks $i$ and $j$ do not overlap. We focused on systems composed of $N/2$ large and $N/2$ small disks with equal mass $m$ and diameter ratio, $\sigma_{L}/\sigma_{S}=1.4$ to avoid crystallization. The total potential energy for the SP model is $U_{SP}=\sum_{i>j}U_{SP}(r_{ij})$ and the stress tensor is given by

[TABLE]

where $\mu$, $\nu=x$ and $y$ and ${\vec{f}}_{ij}=-{\vec{\nabla}}_{r_{ij}}U_{SP}(r_{ij})$. To measure the shear modulus ${\cal G}$, we apply an infinitesimal affine shear strain $\gamma$ to the $x$-positions of the particle centers, $x^{\prime}_{i}=x_{i}+\gamma y_{i}$, and measure the resulting shear stress $\Sigma_{xy}$. We then calculate the shear modulus ${\cal G}=-d\Sigma_{xy}/d\gamma$ (at fixed area). For the SP model, we measure energy, length, and stress in units of $\epsilon$, $\sigma_{LS}$, and $\epsilon/\sigma^{2}_{LS}$.

2.2 Deformable Particle Model

To model bubbles and droplets, we consider the deformable particle model 29 for foams and emulsions with a potential energy that includes the following three terms:

[TABLE]

Each deformable “particle” (indexed by $m=1,\ldots,N$) is modeled as a polygon with $N_{v}$ circulo-line edges to represent $N_{v}-1$ shape degrees of freedom. The circulo-lines have width $\delta$ 40 and are indexed by $i=1\ldots,N_{v}$. (See Fig. 2c).

We consider $N/2$ large particles with $N_{v}=17$ and $N/2$ small particles with $N_{v}=12$, and $a_{L0}/a_{S0}=(17/12)^{2}\approx 2.0$, which is similar to the area ratio of the large and small disks in the SP model. We have also studied DP packings with larger numbers of vertices (while maintaining $a_{L0}/a_{S0}\approx 2.0$), and the structural and mechanical properties are similar to those for $N_{v}=17$ and $12$ for the large and small particles, respectively. The location of the $i$th circulo-line in particle $m$ is ${\vec{v}}_{m,i}$, and the bond vector ${\vec{l}}_{m,i}={\vec{v}}_{m,i+1}-{\vec{v}}_{m,i}=l_{m,i}{\hat{l}}_{m,i}$ connects circulo-lines $i+1$ and $i$.

The first term in $U_{DP}$ is proportional to the total length of the interface, i.e. the perimeter, $p_{m}=\sum_{i=1}^{N_{v}}{l_{m,i}}$ of the $m$th particle with a proportionality constant equal to the line tension $\gamma$. The second term is quadratic in $a_{m}$ with a minumum at $a_{L,S0}$, which penalizes deviations in area from the reference value $a_{L,S0}$. Here, we study $k_{a}a_{L,S0}^{2}>10^{3}$, which implies that the fluctuations in the particle areas are $\lesssim 10^{-3}$. We characterize the shape of the deformable particles by calculating the particle shape parameter ${\cal A}_{m}=p_{m}^{2}/4\pi a_{m}$, which equals $1$ for circular disks and is greater than $1$ for all non-spherical shapes 41, 29.

Note that for the DP model for foams and emulsions, we remove the constraint on the elastic interface 29, i.e. the preferred bond length, $l_{m,i}=0$. As a result, the spacing between the vertices on a given deformable particle can change as they interact with vertices on neighboring particles. This allows us to correctly model the formation of elongated edges when deformable particles make contacts (such that $l_{m,i}>0$), as well as model the formation of Plateau borders 37, 36 with $l_{m,i}\rightarrow 0$.

The third term, $U_{\rm int}$, penalizes overlaps between deformable particles by including purely repulsive interactions between pairs of contacting circulo-lines on neighboring deformable particles:

[TABLE]

where $\epsilon_{r}$ gives the strength of the repulsive interactions, $d_{\rm min}$ is the minimum distance between circulo-lines $j$ and $k$ on contacting deformable particles $m$ and $n$, and $\Theta(.)$ ensures that there is no interaction when the circulo-lines on different particles are out of contact. The stress tensor for packings of deformable particles is obtained using $\Sigma_{\mu\nu}=A^{-1}\sum_{i=1}^{N}f_{i\mu}^{\rm ext}r^{c}_{i\nu}$, where ${\vec{f}}^{\rm ext}_{i}=-{\vec{\nabla}}_{r_{i}}U_{\rm int}$ is the force on particle $i$ arising from $U_{\rm int}$ and ${\vec{r}}_{i}$ is the position of the centroid of particle $i$.

To measure the shear modulus ${\cal G}$, we apply an infinitesimal affine shear strain $\gamma$ to the $x$-positions of the $i=1,\ldots,N_{v}$ circulo-lines on each particle $m$, $v_{xm,i}^{\prime}=v_{xm,i}+\gamma v_{ym,i}$, and measure the resulting shear stress $\Sigma_{xy}$. We then calculate the shear modulus ${\cal G}=-d\Sigma_{xy}/d\gamma$ (at fixed area). For the DP model, we measure energy, length, and stress in units of $\epsilon_{r}$, $l$, and $\epsilon_{r}/l^{2}$, where $l=(\sqrt{a_{S0}}+\sqrt{a_{L0}})/2$. The structural and mechanical properties of DP packings at jamming onset do not depend on the parameters $\gamma$, $k_{a}$, and $\epsilon_{r}$. However, above jamming onset, the properties can depend on these parameters. We focused on the parameter regime, $\epsilon_{r}>k_{a}(a_{L,S0})^{2}>\gamma\langle p_{m}\rangle$, which is typical for foams and emulsions.

2.3 Isotropic Compression Packing Protocol

The protocol to generate jammed packings is similar for the SP and DP models. The protocol proceeds in two stages. For each initial condition, we first identify the packing fraction $\rho_{J}=\phi_{J}$ at jamming onset. For the SP model, the system is initialized using random locations for the disks at $\rho=0.20$. For the DP model, we place the particle centers randomly at $\rho=0.20$ and then position the $N_{v}$ circulo-lines equally spaced around each particle center. We successively isotropically compress the system (by decreasing the size of the simulation cell) using small packing fraction increments ($\delta\phi$=$10^{-4}$ for the SP model and $\delta\rho=10^{-4}$ for the DP model) and minimize the total potential energy per particle, $U_{SP}/N$, for the SP model (or $U_{DP}/(NN_{v})$ for the DP model) after each compression step using over-damped molecular dynamics simulations until the kinetic energy per particle satisfies $K/N<10^{-20}$ for the SP model and the total kinetic energy per circulo-line $K/(NN_{v})<10^{-20}$ for the DP model. If $U_{SP}/N$ or $U_{DP}/(NN_{v})$ is “zero” (i.e. $U_{SP}/N<10^{-15}$ or $U_{DP}/(NN_{v})<10^{-15}$) after minimization, the system is subsequently compressed. If $U_{SP}/N$ or $U_{DP}/(NN_{v})$ is nonzero, i.e. there are finite particle overlaps and $U_{SP}/N>10^{-13}$ or $U_{DP}/(NN_{v})>10^{-13}$, after minimization, the system is subsequently decompressed. The increment by which the packing fraction is changed at each compression or decompression step is gradually decreased.

We terminate the process of finding the onset of jamming at $\rho_{J}$ or $\phi_{J}$ when the system satisfies $10^{-13}<U_{SP}/N<10^{-16}$ and $K/N<10^{-20}$ for the SP model or $10^{-13}<U_{DP}/(NN_{v})<10^{-16}$ and $K/(NN_{v})<10^{-20}$ for the DP model. This process yields mechanically stable packings at jamming onset.

The second stage of the protocol involves sampling the system at set of packing fractions $\rho-\rho_{J}>0$ above jamming onset with adjacent values separated by $\delta\rho=10^{-4}$ for the DP model or a set of $\phi-\phi_{J}>0$ with adjacent values separated by $\delta\phi=10^{-4}$ for the SP model. Ensemble averages are obtained by averaging over systems at fixed $\phi-\phi_{J}$ or $\rho-\rho_{J}$ for the SP and DP models, respectively, where $\phi_{J}$ and $\rho_{J}$ are determined separately for each initial condition. The distribution of $\phi_{J}$ and $\rho_{J}$ for the SP and DP models are shown for several system sizes in the Appendix.

3 Results

Below, we compare the results for the structural and mechanical properties of particle packings as a function of packing fraction above jamming onset obtained using the soft particle and deformable particle models.

These studies allow us to investigate the effect of particle deformability on the structural and mechanical properties of jammed solids. As discussed in Sec. 1, there are several key differences between the SP and DP models. For example, the SP model allows overlap between particles and concomitant decreases in the particle area as the system is compressed above jamming onset. In contrast, the particles in the DP model deform to prevent interparticle overalps, and thus they maintain their areas, and do not remain circular in shape.

A method to minimize the effects of the loss of particle area for the SP model is to quantify the structural and mechanical properties of jammed packings generated using the SP model as a function of the true packing fraction $\rho$, not $\phi$. To measure $\rho$ at each $\phi>\phi_{J}$, we need to subtract from $\phi A$ the muliply-counted areas of overlapping disks. For $\phi<1.2$, we only need to consider overlaps between pairs of disks, i.e. subtract off the area of each lens between pairs of overlapping disks. (See Fig. 1.) In this case, the true packing fraction is

[TABLE]

where

[TABLE]

is the area of the lens between overlapping disks $i$ and $j$ and ${\overline{\sigma}}_{ij}=(\sigma_{i}-\sigma_{j})/2$. For $\phi>1.2$, the lens between overlapping disks $i$ and $j$ can overlap with the lens of other overlapping pairs of disks, which modifies Eq. 9.

In Fig. 3 (a), we plot the deviation in the true packing fraction from that at jamming onset, $\Delta\rho=\rho-\rho_{J}$, versus the deviation in the reduced number density $\phi$ from the value at jamming onset, $\Delta\phi=\phi-\phi_{J}$, for jammed packings generated using the SP model. On a linear scale, $\rho\approx\phi$ for $\phi\lesssim 0.88$. More generally, we find

[TABLE]

where $C$ is weakly dependent on $\phi$. (See Figs. 10 (a) and (b) in the Appendix.) Note that Eq. 11 is exact when higher-order overlaps (i.e. the lens from overlapping disks $i$ and $j$ overlaps with the lens from other overlapping disks) do not occur. To study the structural properties of the SP model near confluence, $\rho\rightarrow 1$, we also considered cases where three disks mutually overlap. The true packing fraction becomes

[TABLE]

where $a_{ijk}^{\rm ov}$ is the area of the Reuleaux triangles that form when three disks mutually overlap. Using this approximation, we find that $\rho\rightarrow 1$ near $\phi\approx 1.24$.

As discussed in the introduction, for packings generated using the SP model, the area of the particles decreases with increasing packing fraction above jamming onset. We calculate the average area of the particles (normalized by the average at jamming onset) $\langle a\rangle/a_{J}$ versus $\Delta\rho$ for packings generated using the SP and DP models in Fig. 3 (b). On a linear scale, $\langle a\rangle$ begins deviating significantly from $a_{J}$ for $\Delta\rho\gtrsim 0.04$ for the SP model. In the Appendix, we show that $a/a_{J}-1\propto\Delta\rho^{\zeta}$ with an exponent $\zeta\approx 1.5$ at small $\Delta\rho$ and $\zeta\approx 2.5$ at large $\Delta\rho$ for the SP model. In contrast, $\langle a\rangle\approx a_{J}$ over the full range of $\Delta\rho$ for the DP model.

In Fig. 3 (c), we quantify how the particles deform during isotropic compression above jamming onset. In general, the packing fraction increases with the shape parameter ${\cal A}-1$. In the Appendix, we show that $\rho-\rho_{J}$ grows as a power-law in the deviation of the shape parameter from that at jamming onset, $({\cal A}-1)^{\omega}$. At small ${\cal A}-1$, $\omega\approx 0.66$ and $\approx 1.0$ for the SP and DP models, respectively. At large ${\cal A}-1$, $\omega\approx 0.3$ for both the SP and DP models. $\phi-\phi_{J}$ for the SP model also grows as a power-law with ${\cal A}-1$, but much faster than $\rho-\rho_{J}$.

At what shape parameter do 2D foams and emulsions reach confluence? Fig. 3 (c) shows that the DP and SP models reach confluence ($\Delta\rho\approx 0.16$) at different values of the shape parameter, ${\cal A}\approx 1.07$ for the DP model and ${\cal A}>1.10$ for the SP model. Thus, we find similarites and differences between the shapes of the particles for packings generated using the SP and DP models as they are compressed above jamming onset. An interesting similarity is that the packing fraction for both the SP and DP models scales as $\rho-\rho_{J}\sim({\cal A}-1)^{\omega}$, where $\omega\approx 0.3$ at large $\rho$.

Next, we compare the contact number $z$ and shear modulus ${\cal G}$ versus the packing fraction for the SP and DP models. On a linear scale, which emphasizes the values at large packing fraction, we find weak system-size dependence for $z$ and ${\cal G}$ for packings generated via the SP and DP models. Note that for the deformable particle model, multiple circulo-lines on one deformable particle can be in contact with multiple circulo-lines on another deformable particle. These multiple circulo-line contacts are treated as a single contact between two deformable particles.

In Fig. 4, for both $z$ and ${\cal G}$, measured in packings of $N=256$ particles, the results for the SP and DP models are similar near jamming onset $\phi_{J}\approx\rho_{J}$. For $z$ and ${\cal G}$, the results for the SP and DP models begin to deviate near $\rho\approx 0.88$. We find that more interparticle contacts form as the packings are compressed above jamming onset for the DP model, compared to that for the SP model. As a result, the shear modulus grows more rapidly with $\rho$ for packings generated using the DP model. We also show the best fits of $z$ and ${\cal G}$ for the SP model to the power-law scaling form with $\Delta\phi$ in Eqs. 1 and 2. As found previously, the scaling exponents $\alpha_{0}^{\phi}=\beta_{0}^{\phi}\approx 0.5$ for the SP model.15, 13 By plugging Eq. 11 into Eqs. 1 and 2, we can convert $z(\Delta\phi)$ and ${\cal G}(\Delta\phi)$ to $z(\Delta\rho)$ and ${\cal G}(\Delta\rho)$. To lowest order in $\Delta\rho$, we find

[TABLE]

and

[TABLE]

for the SP model, where $\alpha_{0}^{\rho}\approx\beta_{0}^{\rho}\approx 0.5$ and $\alpha_{1}^{\rho}\approx\beta_{1}^{\rho}\approx 1.0$. We show fits of $z(\Delta\rho)$ and ${\cal G}(\Delta\rho)$ for the SP model to Eqs. 13 and 14 as dashed-dotted lines in Fig. 4. The combination of the two power-laws in $\Delta\rho$ with exponents $0.5$ and $1.0$ accurately describes the data for $\Delta z$ and ${\cal G}$. However, $\Delta z$ begins deviating from Eq. 13 for $\rho\gtrsim 0.98$ near confluence. Moreover, we find that Eqs. 13 and  14 can be used to fit the data for the DP model as well. The parameters for the fitting functions are shown in Tables 1 and  2.

In Fig. 5 (a) and (d), we present the excess contact number $z-z_{J}$ and shear modulus ${\cal G}$ versus $\Delta\phi$ on logarithmic axes for jammed packings using the SP model for system sizes ranging from $N=32$ to $256$. By plotting the data on logarithmic scales, we can identify several different regimes in $\Delta\phi$: 1) $\Delta\phi<10^{-3}$, 2) $10^{-3}<\Delta\phi<0.2$, and 3) $\Delta\phi>0.2$. Regime $1$, where $\phi\approx\phi_{J}$, is difficult to see on the linear scales shown in Fig. 4 (a). In this regime, there is strong system-size dependence and $\Delta z\sim\Delta\phi$ and ${\cal G}\sim\Delta\phi^{0}$. 22, 42 As found previously, at intermediate $\Delta\phi$, $\Delta z\sim\Delta\phi^{\alpha_{1}^{\phi}}$ and ${\cal G}\sim\Delta\phi^{\beta_{1}^{\phi}}$, where $\alpha_{1}^{\phi}\approx\beta_{1}^{\phi}\approx 0.5$. The characteristic $\Delta\phi^{*}$ at which the $\alpha$ and $\beta$ exponents cross-over from $1$ to $0.5$ and [math] and $0.5$, respectively, scales as $\Delta\phi^{*}\sim N^{-2}$. (See Tables 1 and 2.) Thus, regime $2$ extends to smaller $\Delta\phi$ as the system size increases. In regime $3$, at large $\Delta\phi$, we find that the power-law scaling behavior of $\Delta z$ and ${\cal G}$ with $\Delta\phi$ breaks down.

Using Eq. 11, we can convert the data from functions of $\Delta\phi$ to functions of $\Delta\rho$ for the SP model. We show $\Delta z(\Delta\rho)$ and ${\cal G}(\Delta\rho)$ in Fig. 5 (b) and (d). We find two different regimes in $\Delta\rho$. At low $\Delta\rho$, $z-z_{J}\sim\Delta\rho$ and ${\cal G}\sim\Delta\rho^{0}$. At intermediate and large $\Delta\rho$, the forms for $z-z_{J}$ and ${\cal G}$ include a sum of two power-laws as shown in Eqs. 13 and 14. The characteristic packing fraction that separates the small and large $\Delta\rho$ regimes also scales as $\Delta\rho^{*}\sim N^{-2}$.

In Fig. 5 (c) and (f), we show $\Delta z$ and ${\cal G}$ versus $\Delta\rho$ for the DP model. The data also show two regimes in $\Delta\rho$. At small $\Delta\rho$, $\Delta z\sim\Delta\rho$ and ${\cal G}\sim\Delta\rho^{0}$. At intermediate and large $\Delta\rho$, we fit $\Delta z$ and ${\cal G}$ to a sum of two power laws with exponents $0.5$ and $1$. (See Tables 1 and 2.) Note that the scaling of $\Delta z$ and ${\cal G}$ for the DP model is similar to that for the SP model (Eqs.13 and 14), but the ratio of the coefficient for the linear term in $\Delta\rho$ to that for the $\Delta\rho^{0.5}$ term is much larger than the corresponding ratio for the SP model. If we use only a single power law exponent to fit the data for $\Delta z$ and ${\cal G}$ at large $\Delta\rho$, we find a best fit exponent of $0.8$ for both $\Delta z$ and ${\cal G}$. Winklemann, et al. 20 recently studied packings of bubbles and showed that $z-z_{J}\sim(\rho-\rho_{J})^{\alpha^{\prime}}$ with an exponent $\alpha^{\prime}=1.0$, which is larger than the value we find if we use a single power-law to fit the data for $\Delta z$ for $\Delta\rho>10^{-2}$.

In Fig. 5 (a) and (b), the data for $\Delta z$ begin to plateau at large $\Delta\phi$ or $\Delta\rho$, indicating that the packings are reaching confluence with $z_{\rm conf}=6$ at $\rho=\rho^{*}$, where $\rho^{*}\rightarrow 1$ in the large system limit. To understand the scaling of $z$ near confluence, we plot $z_{\rm conf}-z$ versus $\rho^{*}-\rho$. For the SP model, we use Eq. 12 to determine $\rho$. We find that the data for $z_{\rm conf}-z$ versus $\rho^{*}-\rho$ for different system sizes can be collapsed onto a master curve when we scale $z_{\rm conf}-z$ and $\rho^{*}-\rho$ by $N$ as shown in the inset to Fig. 6 (a). However, as $\rho\rightarrow\rho^{*}$, $(z_{\rm conf}-z)N$ does not approach zero for the SP model. Instead, $(z_{\rm conf}-z)N\rightarrow 2$. $(z_{\rm conf}-z)N$ approaches a finite value because Eq. 12 does not account for all multiply-counted overlaps between mutually overlapping pairs of disks (i.e. beyond three mutually overlapping disks). For the DP model, $z_{\rm conf}-z$ verus $\rho^{*}-\rho$ can also be collapsed onto a master curve when both $z_{\rm conf}-z$ and $\rho^{*}-\rho$ are scaled by $N$ (Fig 6 (b)). Unlike the SP model, for the DP model we find that $z_{\rm conf}-z$ scales as $\rho^{*}-\rho$ in the large system limit over several orders of magnitude in $\rho^{*}-\rho$.

In addition to the contact number and shear modulus, we also studied the variation of the structure of the void space of jammed packings as a function of packing fraction $\rho$. To do this, for each jammed packing, we identify all of the void space that is not occupied by particles. We then determine the connected void regions (i.e. one can reach any part of a connected void region from any point in the region). The topology of each connected void can be characterized by the number of edges, or the smallest number of particles that form a loop on the perimeter of the void region (using the depth first search algorithm). 43 In Fig. 7 (a), we compare the probability $P_{l}(\rho)$ to have a void with $l$ sides as a function of $\rho$ for the SP and DP models. Near $\rho_{J}$, the SP and DP models are identical and we find that the distributions $P_{l}(\rho_{J})$ are the same for the two models. In this regime, the probability of $3$- and $4$-sided voids are similar ($\sim 0.4$), while the probabilities for $5$- and $6$-sided voids are much smaller ($\sim 0.18$ and $\sim 0.02$). For $\rho\gtrsim 0.88$, $P_{l}(\rho)$ for the DP model begins to deviate from that for the SP model. Similar behavior was found for the deviation in the contact number $\Delta z(\rho)$ for $\rho\gtrsim 0.88$. For the DP model, we find an increase in the probability of $3$-sided voids over that for the SP model and a comparable decrease in the probability of $4$-sided voids relative to the SP model. For the DP model, we find that $P_{3}=1$ for $\rho\gtrsim 0.97$. In contrast, $P_{3}=1$ only in the limit $\rho\rightarrow 1$ for the SP model. Further, we find that the shape parameter ${\cal A}={\cal A}_{p}=\pi/(4\sqrt{3}-2\pi)\approx 4.87$ of the $3$-sided voids for the DP model is independent of $\rho$, indicating that the DP model correctly captures the structure of the Plateau borders that form as $\rho\rightarrow 1$. However, the shape parameter of the $3$-sided voids varies from ${\cal A}\approx{\cal A}_{p}$ to less than $2$ as $\rho$ increases from $\rho_{J}$ to $1$ for the SP model, which indicates that the void structure for the SP model differs significantly from that for jammed packings of foams and emulsions near confluence.

In Fig. 8, we show preliminary studies comparing the results for jammed packings generated using the DP model to the results from optical microscopy experiments of quasi-2D compressed emulsion droplets 23. We find that the contact number $z$ versus $\rho$ for the DP model closely matches that from the experiments (Fig. 8 (c)). In addition, we show that the evolution of the shape factor of the particles with packing fraction above jamming onset is similar for the DP model and experiments (Fig. 8 (d)). However, we encourage additional experiments on compressed foams and emulsions to be performed with packing fraction $\rho\gtrsim 0.95$ to determine whether the DP model can recapitulate the structural and mechanical properties of compressed emulsions near confluence. For example, new experimental studies can test the prediction for the DP model that $z_{\rm conf}-z$ scales linearly with $1-\rho$ near $\rho=1$.

4 Discussion and Conclusions

In this article, we investigated the structural and mechanical properties of jammed packings undergoing isotropic compression in 2D using the soft particle (SP) model 18, 44 and the new deformable particle model 29 that we developed for bubbles and emulsions. The SP model has been widely used to characterize the structural, mechanical, and rheological properties of jammed particulate systems including granular materials 15, dense colloidal suspensions 45, foams 18, and emulsions 46. The key difference between the two models is that in the SP model, the particles decrease in area as the system is compressed above jamming onset, while the DP model conserves particle area during compression. Studies that have employed the SP model typically characterize the properties of jammed packings as a function of the reduced number density $\phi$ 16, 13, 17, rather than using the true packing fraction $\rho$. In this work, we provide direct comparisons of the structural and mechanical properties of packings generated using the SP and DP models as a function of the packing fraction $\rho$.

First, we showed explicitly that the SP and DP models give the same results near jamming onset, where the disks are undeformed. In particular, we showed that the probability distribution $P(\rho_{J})$ of jamming onsets $\rho_{J}\approx\phi_{J}$ are nearly identical for the SP and DP models with $\rho_{J}\approx 0.842$ in the large system limit. (For the detailed discussion of this point, see the Appendix.) In addition, we find similar scaling behavior for the excess contact number $z-z_{J}$ and shear modulus ${\cal G}$ versus $\rho-\rho_{J}$. Near jamming onset, for both the SP and DP models, $z-z_{J}\sim(\rho-\rho_{J})^{1.0}$ and ${\cal G}\sim(\rho-\rho_{J})^{0}$ for small systems. 22, 42 This scaling behavior occurs for $\rho-\rho_{J}<\Delta\rho^{*}$, where $\Delta\rho^{*}\sim N^{-2}$. In the large system limit, $z-z_{J}\sim{\cal G}\sim(\rho-\rho_{J})^{0.5}$ near jamming onset for both the SP and DP models.

For packings that are compressed above jamming onset, we determined the relation between the packing fraction $\rho$ and reduced number density $\phi$ for the SP model. Using this relation, we showed that for the SP model $\Delta z(\Delta\rho)\sim{\cal G}(\Delta\rho)$ can be represented as a sum of power-laws in $\Delta\rho$ (not a single power-law), with $\Delta\rho^{0.5}$ and $\Delta\rho^{1.0}$ as the lowest order terms. The scaling of $\Delta z$ and ${\cal G}$ is similar for the DP model, but the ratio of the coefficient for the linear term in $\Delta\rho$ to that for the $\Delta\rho^{0.5}$ term is much larger than the corresponding ratio for the SP model. As a result, we find that $z(\rho)$ and ${\cal G}(\rho)$ are larger for the DP model compared to the SP model for $\rho\gtrsim 0.88$.

In addition, we characterized the void space in jammed packings 26, 47 as a function of packing fraction for both models. We showed that the DP model can recapitulate the formation of Plateau borders as the deformable particles become tightly packed. The probability to obtain $3$-sided voids becomes unity for $\rho\gtrsim 0.97$ and the shape factor of the $3$-sided voids (${\cal A}\approx{\cal A}_{p}$) is independent of $\rho$. In contrast, for the SP model, the probability of $3$-sided voids only becomes unity in the limit $\rho\rightarrow 1$ and ${\cal A}<{\cal A}_{p}$ over the full range of packing fraction. Thus, the SP model does not capture the topological features of the void space of packings of compressible bubbles near confluence. We believe that the results from this article will inspire additional experimental and theoretical studies of the collective behavior of droplets and bubbles in emulsions and foams. For example, the DP model can be extended to include attractive interactions to investigate the mechanical response of attractive emulsions in both 2D and 3D 48. In addition, active forces can be added to each circulo-line in the DP model to simulate collective motion, such as swarming and migration, in cell aggregates as well as living tissues.

5 Appendix

In this Appendix, we provide additional calculations to support the conclusions described in the main text. In Fig. 9, we show the distribution of the reduced number densities $P(\phi_{J})$ for the SP model (panel (a)) and the distribution of packing fractions $P(\rho_{J})$ for the DP model (panel (b)) at jamming onset using the isotropic compression protocol discussed in Sec. 2.3. (We also include the fraction $f_{J}$ of packings with $\phi_{J}$ (or $\rho_{J}$) at or below a given value in the insets to Fig. 9 (a) and (b).) The distributions for the SP and DP models were obtained using the same initial random particle positions. At jamming onset $\phi_{J}=\rho_{J}$, and thus the distributions $P(\phi_{J})$ and $P(\rho_{J})$ for the SP and DP models are nearly identical. Similar to previous studies, we also find that the root-mean-square deviation in the packing fraction at jamming onset scales as $\Delta\phi_{J}\sim N^{-\theta}$, with $\theta\approx 0.55$, and $\phi_{J}\rightarrow 0.842$ in the large system limit.

In Fig. 10, we show a series of scaling laws for several physical quantities and compare the behavior for packings generated using the SP and DP models. In Fig. 10 (a), we show that $\rho-\phi\sim(U_{SP}/N)^{\kappa}$, where $\kappa\approx 0.75$, for the SP model. The scaling exponent $\kappa$ can be obtained by assuming $\phi-\rho\propto\Delta\phi^{1.5}$ from Eq. 11 and using $\Delta\phi\sim(U_{SP}/N)^{0.5}$ from Eq. 5.

In Fig. 10 (b), we show $\Delta\rho=\rho-\rho_{J}$ versus $\Delta\phi=\phi-\phi_{J}$ for packings generated using the SP model. The dashed line obeys Eq. 11 with $C\approx 1.2$. As discussed in Sec. 3, the relation between $\Delta\rho$ and $\Delta\phi$ only holds when multiply-counted overlaps (i.e. beyond pairwise overlaps) are absent. Combining the results from Fig. 10 (b) and 3 (a), we see that Eq. 11 holds over nearly $10$ orders of magnitude in $\Delta\phi$. In Fig. 10 (c), we show the data for the relative deviation in the area of the particles $1-\langle a\rangle/a_{J}$ versus $\rho-\rho_{J}$ from Fig. 3 (b) on logarithmic scales for the SP model. In contrast to the DP model, the SP model does not conserve particle area when ovecompressed above jamming onset. We find two power-law scaling regimes. $1-\langle a\rangle/a_{J}\sim\Delta\rho^{\zeta}$ with $\zeta\approx 1.5$ and $2.5$ at small and large $\Delta\rho$, respectively.

In Figs. 10 (c) and (d), we characterize the change in particle shape as the packings are compressed above jamming onset for the SP and DP models. We show that $\Delta\rho\sim({\cal A-}1)^{\omega}$, but the scaling exponent takes on different values in the for small and large values of ${\cal A}-1$. At small ${\cal A}-1$, $\omega\approx 0.66$ and $\approx 1.0$ for the SP and DP models, respectively. At large ${\cal A}-1$, $\omega\approx 0.3$ for both the SP and DP models.

6 Acknowledgments

We acknowledge support from NSF Grants No. PHY-1522467 (A.B., A.S., and C.S. O’Hern), No. CBET-1336401 (J.L., C.S. Orellana, and E.R.W), No. CBET-1804186 (E.R.W), No. CMMI-1462439 (C.S. O’Hern) and No. CMMI-1463455 (M.D.S.), the President’s International Fellowship Initiative (PIFI) and Hundred Talents Program of the Chinese Academy of Sciences (A.B. and F.Y.), and the Raymond and Beverly Sackler Institute for Biological, Physical, and Engineering Sciences (A.B.). This work was also supported by the High Performance Computing facilities operated by Yale’s Center for Research Computing.

Conflicts of interest

There are no conflicts to declare.