Measuring and upscaling micromechanical interactions in a cohesive granular material

Abstract

The mechanical properties of a disordered heterogeneous medium depend, in general, on a complex interplay between multiple length scales. Connecting local interactions to macroscopic observables, such as stiffness or fracture, is thus challenging in this type of material. Here, we study the properties of a cohesive granular material composed of glass beads held together by soft polymer bridges. We characterise the mechanical response of single bridges under traction and shear, using a setup based on the deflection of flexible micropipettes. These measurements, along with information from X-ray microtomograms of the granular packings, then inform large-scale discrete element model (DEM) simulations. Although simple, these simulations are constrained in every way by empirical measurement and accurately predict mechanical responses of the aggregates, including details on their compressive failure, and how the material's stiffness depends on the stiffness and geometry of its parts. By demonstrating how to accurately relate microscopic information to macroscopic properties, these results provide new perspectives for predicting the behaviour of complex disordered materials, such as porous rock, snow, or foam.

---

1 Introduction

Disordered porous materials are ubiquitous in nature and industry. The mechanical behaviour of these materials depends on the local interactions between their components,^1–4 which can be of various shapes, compositions, and length scales, from nanometres in colloidal glasses to millimetres in sedimentary rocks. Understanding and predicting the relationship between their microscopic structure and macroscopic mechanical response is of importance for applications in engineering, geoscience and materials science, and is the topic of active research in both experimental and numerical domains.^2–11 As an original example, it has been shown how birds such as swallows build strong nests using their saliva to form cohesive bonds between mud granules, in an analogy with the mechanical properties observed in artificial cohesive granular materials.⁴

Linking the mechanical properties of a disordered solid to its local structure presents difficulties. For example, the failure of a disordered material can suddenly change from ductile to brittle, depending on the density, strength, and plasticity of its contacts.^12–15 Experimentally, it is hard to study such effects in realistic systems: these are complex by nature and have intrinsic characteristics that cannot be systematically varied. While different bottom-up approaches using model systems have emerged within the last ten years,^16–20 they are still rare. Thus, most progress in this field has been made in numerical modelling,^5–9,21 with the inherent difficulty that a too simplistic model might not capture the complexity observed experimentally, while an over-detailed one may not give any useful general law.

A tunable experimental model of a cohesive granular material has been recently introduced,^19,20 which we will use here to help overcome these challenges. It consists of glass beads held together by elastic bridges of polydimethylsiloxane (PDMS), a solidified elastomer, whose Young's modulus can be easily varied by changing its composition (see Fig. 1). Despite its apparent simplicity this material shows a rich range of behaviours, allowing investigation of the mechanical properties of cohesive granulates in particular and of disordered solids in general. Its elastic¹⁹ and fracture²⁰ properties can be finely controlled by varying the stiffness, size, and density of the PDMS bonds linking the beads. The wide range of tunable parameters makes this clean and well-defined system an appropriate candidate to investigate scaling laws in the mechanics of disordered structures including, for examples, avalanches in snow, onsets of earthquakes, fracture of rocks and building materials, or poroelasticity and fluid–rock interactions in rocks under large hydrostatic pressure.

Fig. 1 Structure of a cohesive granular material. (a) A micrograph of a cohesive granular sample shows the capillary bridges of solid PDMS formed between glass beads. (b) We model each pair of connected beads as spheres of diameter D linked by a truncated cylinder of diameter d and height l and separated by a gap of size δ⁰. (c) X-ray tomogram of a sample compressed in situ between two pistons. The inset shows a cross-section of the 3D data. The sample is a cylinder of 4.25 mm diameter × 4.05 mm height, made of beads of diameter D = 200.9 μm.

This cohesive granular material also shows intriguing mechanical properties. For instance, it fails under uniaxial compression by sliding along a shear plane (as expected for this type of material) but only after a remarkably large yield strain.¹⁹ Such a feature is the result of an interplay between microscopic properties at the scale of individual beads and the large-scale (re-)organisation of the packing under mechanical load. A full understanding of the macroscopic properties of the material requires a good characterisation of the local contacts and a proper upscaling of microscopic laws towards larger assemblies of beads.

To this end, we report here direct measurements of the constitutive relationship of the bonds between individual particles in this type of cohesive granular material, and then demonstrate how to use these measurements to accurately predict the elastic properties of such a material. To link the microscopic and macroscopic length scales we use a discrete element model (DEM), developed by some of us to investigate the failure mechanisms of cohesive granulates, and which incorporates a Griffith-like energy-based failure criterion for the bonds.¹⁵ Both the experiments and model involve hard particles, connected by softer bridges, and in this limit we show that the bridge geometry is a key parameter in upscaling mechanical properties. The experimental observations will inform the parameters used in the model and we will then show that this approach can reproduce, both qualitatively and quantitatively, key features observed in the experiments, such as the results of uniaxial compression tests.

2 Micromechanics of cohesive granulates

We first characterised the micromechanical interactions between particles in a cohesive granular medium. If two particles or grains are mechanically bound to each other, they may interact via normal and tangential displacements, and by rolling and twisting. These mechanisms can be modelled by elastic springs with associated spring constants for the normal and tangential deformations (k^bond_n and k^bond_t, respectively), and bending stiffnesses for the rolling and twisting ones. Aiming for a micromechanical model which reproduces the essential properties of this type of material with as few parameters as necessary, we developed a protocol for measuring the normal and tangential spring constants over a single pair of beads, while assuming contributions from bending and twisting to be less relevant. Here, we describe these methods and give direct observations of the elasticity of the bonds between particles taken from a cohesive granular material.

2.1 Experimental methods

Micromechanical testing was conducted using flexible glass micropipettes (see Fig. 2), in an approach adapted from one used to measure the forces exerted by cells and other micro-organisms.²²

Fig. 2 Schematics of the micromechanical tests in the normal (a) and tangential (b) configurations. The dotted frames correspond to the images shown as inserts in the example force–displacement measurements given in (c) and (d). Linear fits (red lines) to the force–displacement data give the spring constants k^bond_n or k^bond_t. Panels (e) and (f) summarise the spring constant measurements made of different bridge diameters d and at various speeds. The lines are the results of simulations of bridges in similar conditions, assuming an initial bridge length, δ⁰, of either 1 μm (dotted line) or 2 μm (solid line).

Samples of cohesive granular material were prepared, following published methods.¹⁹ The elastic bridges consist of polydimethylsiloxane (PDMS, Sylgard 184, Dow Corning), a curable elastomer, using a mass ratio of base to cross-linker of 40 : 1. Uncured PDMS was mixed with glass beads of 365 μm diameter (Sigmund Lindner GmbH, 5% polydispersity) and at a polymer volume fraction W = 2.3%, poured into a cylindrical mould and cured at 90 °C overnight. The Young's modulus of the PDMS itself, E_p = 50 ± 5 kPa, was measured by unconfined uniaxial compression of a cylinder of cured polymer. The high base to cross-linker ratio, and low E_p, allows the bridges to be soft enough to give measurable deformations during the micromechanical tests, while still behaving rigidly and elastically. Finally, like many other elastomers, PDMS is known to be nearly incompressible, with a Poisson ratio ν ≥ 0.49 measured over a wide range of curing conditions.^23,24

For each test a single pair of beads, connected by a PDMS bridge, was detached from a sample of cohesive granular material using a scalpel under an inverted microscope (IX-73, Olympus), taking care not to stretch or bend the bridge. The beads were moved onto a glass slide, allowing for half of one bead to protrude over the edge of the slide. A droplet of epoxy was deposited onto the end of the dangling bead, which was then quickly put into contact with a flat silicon wafer. After the epoxy was cured the glass slide was removed, leaving the pair of beads attached normally to the wafer. Finally, a glass micropipette was glued onto the other bead, resulting in configurations as shown in Fig. 2.

The micropipettes are narrow glass filaments (borosilicate glass capillaries, World Precision Instruments, TW100-6), pulled by an horizontal pipette puller (P-97 Micropipette Puller, Sutter Instruments) and bent into a characteristic U-shape in a microforge (MF-900, Narishige). This shape gives them an elastic response in extension, as the glass would break before any plastic deformation would occur. It also allows pure normal (Fⁿ) or shear (F^t) forces to be independently directed across the bridge, minimising any off-axis contributions. A pipette’s spring constant is fixed by its geometry and measured using a pre-calibrated AFM cantilever (Veeco).²² For example, one typical pipette had a diameter of 0.3 mm and side lengths, as defined in Fig. 2, of a₁ = 2.5 mm, a₂ = 13 mm, a₃ = 3.5 mm, a₄ = 8.5 mm and a₅ = 117 mm, resulting in a spring constant of 40 ± 2 N m⁻¹. During an experiment, one end (top of a₅) of the pipette was kept fixed, while the wafer was mounted on and moved by a motorised translation stage (Newport, Conex LTA-HS).

Images were recorded of experiments at 10 frames per second with a digital camera (FLIR Systems, Grasshopper 3, GS3-U3-41C6M-C) and a 4× (1.37 μm per pixel) or 10× objective (0.55 μm per pixel). Bead positions were measured using cross-correlation image analysis with a sub-pixel resolution.²² Glass spikes made on the pipette with the microforge allowed for independent tracking of the beads and pipette. The normal/extensional (δ_n) or tangential (δ_t) deformation of the bridge was deduced from the relative displacement of the two beads, while the force exerted on the pipette was calculated from its deflection and spring constant.²² Finally, the bead diameters were found by fitting circles to their shapes in the microscopy images, while the bridge diameters, d, were measured across their middle. For both diameters, the measurement uncertainty is estimated to be ± 2 μm. The initial bridge spanning distance, δ⁰, i.e. the initial surface-to-surface separation between the beads, was always found to be lower than 2 μm, but could not be measured more accurately.

2.2 Results of micromechanical tests

Tests were performed on bridges with diameters d between 50 and 270 μm (see Movies S1 and S2 in the ESI† for examples of normal and shear tests, respectively). In each case a spring constant was calculated from a linear fit to the force–displacement curve in the linear regime of deformation, as demonstrated in Fig. 2(c and d). The restoring forces were linear over a wide range of displacements, especially compared to the initial particle separation, δ⁰. This confirms that we can treat the bonds as Hookean springs. We further checked that the response was in the quasi-static limit by repeatedly testing some bridges at increasing deformation speeds of 1 to 10 μm s⁻¹, with minimal differences in the spring constants.

As shown in Fig. 2(e and f), values for the normal spring constants, k^bond_n, lay in the range of 200–800 N m⁻¹, while tangentially, k^bond_t was between 5–30 N m⁻¹. Despite the dispersion of results for bridges prepared in a similar manner, it is clear that k^bond_n ≫ k^bond_t, meaning that bridges are strong in tension, but easy to shear.

2.3 Scaling spring constants: FEM simulation

The micromechanical measurements were made for beads of one particular size, D, and a fixed polymer content W and stiffness E_p. In order to explore how a bridge's spring constants depend on these variables, we performed simulations replicating our micromechanical tests, but for various bond geometries.

To this end, we used COMSOL Multiphysics to build a finite-element model (FEM) simulation of an elastic bridge. The bridge was modelled as a cylinder of diameter d, truncated by two spherical caps of diameter D and surface separation δ⁰, as in Fig. 1(b). It was treated as a linear elastic material, with Young's modulus E_p and Poisson ratio ν = 0.49, consistent with literature measurements of PDMS as a practically incompressible material.^23,24 A normal or tangential displacement was imposed on one of the sphere–bridge interfaces, while the other was kept fixed, using no-slip boundary conditions for both. The total reaction force exerted on either interface was measured as a function of the displacement, from which follow the spring constants k^bond_n and k^bond_t of the simulated bridge.

The key result is shown in Fig. 2(e and f). Namely, if we use the experimental values of D, d and E_p, and take δ⁰ = 2 μm, then the simulations accurately predict the magnitudes of k^bond_n and k^bond_t along with their relatively weak dependence on d. The simulation's sensitivity to δ⁰ can also potentially explain why there is such a large variation in the experimentally measured spring constants. For example, by adjusting δ⁰ between 1 and 2 μm (which is below our experimental resolution), the simulated spring constants span the range of most experimental values.

Otherwise, the scaling of the spring constants with the bridge geometry is not trivial, although it is clear that there should be a linear dependence on E_p. In the following we will therefore use the COMSOL simulations to predict representative values for k^bond_n and k^bond_t, given specific bead and bridge sizes. In particular, we are interested in mimicking experimental uniaxial compression tests with a discrete element model of large-scale bead assemblies, and for these cases have simulated bridges with microscopic parameters matching experimental observations (see Table 1).

Table 1 Key parameters used in the micromechanical modelling. Values are given for the materials used in the micromechanical tests (Section 2), as well as the two samples, A and B, for which the initial bead positions are known from X-ray microtomography (Section 3.1). In all cases E_p, d and D are measured experimentally. The bond stiffnesses are directly measured for the micromechanical test samples, and these results used to validate a FEM simulation of an elastic bridge, which is then used to predict k^bond_n and k^bond_t for the other samples

Property	Symbol	Micromechanical tests	Sample A	Sample B
Bond Young's modulus	E p	50 kPa	1.5 MPa	0.64 MPa
Bond diameter	d	50–270 μm	73.9 μm	72.1 μm
Bead diameter	D	365 μm	200.9 μm	200.9 μm
Bond stiffness, normal	k ^bondn	200–800 N m^−1	5760 N m^−1	2345 N m^−1
Bond stiffness, tangential	k ^bondt	5–30 N m^−1	325 N m^−1	131 N m^−1

---

3 Upscaling micromechanical interactions

Cohesive granular materials can behave as linear elastic solids, with well-defined failure conditions.^{2,8,11,16,19,20,25,26} However, predicting the homogeneous elastic response of members of this class of material from knowledge of their micro-structure, i.e. upscaling a microscopic model, is made difficult by their discrete and disordered nature. Even more challenging is the prediction of how, and when, such a system will fail, because of the importance of nonlinear effects such as force chains and strain localisation when approaching failure.^8,15,26–30

Here, we pursue this goal of upscaling the micromechanical properties of a cohesive granular material. Experimentally, we imaged such materials using X-ray microtomography, and extracted the positions of the constituent beads along with statistical information about the bridge networks. These measurements, combined with the micromechanical characterisation of single bridges (see Section 2), are then used as inputs for a minimal model of cohesive granular materials.¹⁵ This model is thus constrained in all its parameters by observations at the microscopic scale. In Section 4, we will then simulate unconfined uniaxial testing of the modelled material, and find that the results compare favourably to in situ experimental compression tests.

3.1 X-ray microtomography: bead and bridge geometry

X-ray microtomography data obtained in a previous study¹⁹ were further analysed in the present work (sample A, see Fig. 1(c)), and complemented by data from a similar sample which featured softer bonds (sample B). Briefly, unconfined uniaxial compression tests of cohesive granular samples were performed in situ within an X-ray micro-computed tomography system (GE Nanotom). The samples were composed of monodisperse beads of diameter D = 200.9 ± 1.9 μm (Whitehouse Scientific), mixed with PDMS of Young's modulus E_p = 1.5 ± 0.15 MPa (sample A) or E_p = 0.64 ± 0.05 MPa (sample B). These materials were shaped into cylinders of 4.25 mm diameter × 4.05 mm height for sample A, and 4.25 mm diameter × 4.82 mm height for sample B. Tomography scans were performed before compression, and consisted of sets of 2000 projections with a resolution of 1132 × 1132 pixels and a voxel size of 5 μm. The samples were then compressed at a speed of 5 μm min⁻¹, and their stress–strain curves will be compared with simulations in Section 4, and Fig. 3.

Fig. 3 Comparison of experimental (a and b) and simulated (c and d) stress–strain curves. In each case the Young's modulus, E, is obtained by a least-squares fit to the linear region (black lines). Each simulated curve corresponds to a different initial realisation of the particle diameters, but using the same particle positions. Panels (a) and (b) share the same y-axis, as do (c) and (d). The grey area in (c) and (d) gives the number of bonds that break, N, in each increment of strain.

Beads were detected within the reconstructed sample volumes using Matlab, and particle positions determined by finding the centroid of each individual connected volume after segmentation, with a precision of about one voxel (see Tables S1 and S2 in the ESI† for particle positions and Movies S3 and S4 in the ESI† for illustrating their detection). However, the X-ray contrast between the PDMS bridges and the glass beads was insufficient for the automatic detection of which particles were linked by bridges. Instead, we manually counted bridges on beads within the interior of the sample. For example, in sample B the coordination number Z, here defined as the average number of bridges per bead, was found to be 7.4 ± 0.1 (standard error), which is consistent with previous results from a similar sample.²⁰ The distributions of bridge diameters were also extracted by manual inspection of the microtomograms. For example, for sample A measurements from 100 bridges were well-approximated by a Gaussian distribution with a mean of 73.9 μm and standard deviation of 12.7 μm.

The characteristic properties of the two experimental samples are given in Table 1, along with the bond stiffnesses estimated through the micromechanical testing and COMSOL bridge model.

3.2 Discrete element model of cohesive granulates

To link the microscopic measurements to the bulk properties of a cohesive granular material, we used a discrete element model (DEM) that simulates the interactions between particles and their bonds. As the development of this model has been detailed elsewhere,¹⁵ we focus only on its main features here, including the modifications made for this study.

3.2.1 Model setup and compression test.
Bead positions. The initial bead locations are taken from the relative positions of particles as measured in either of the two microtomograms. The natural dispersion in bead sizes and voxel resolution of the tomography data lead to a small uncertainty in these positions, which could result in slightly overlapping beads. We avoided this issue by randomly varying the diameter of each particle with the Box–Muller method³¹ using a Gaussian distribution with a mean diameter of D = 200.9 μm and a standard deviation of 5%, while constraining the sizes of the particles to prevent overlap. Ensembles of at least 5 different realisations of particle sizes were used, and the error bars in Fig. 4 and 5 give the standard deviations of measurements on such ensembles.

Fig. 4 Comparison of how the Young's modulus of a cohesive granular sample, E, depends on the Young's modulus of the bridges composing it, E_p. The simulations (open symbols) reproduce the experimental relationship (closed squares) seen between E and E_p, and show the importance of the microscopic bridge length δ⁰ in determining the macroscopic stiffness of the sample. The dashed line shows E = E_p, for comparison.

Fig. 5 The Young's modulus of the material, E, increases with increasing content of PDMS, W, in a similar manner in both the experiments (closed squares) and the simulations (open circles), up to W = 2.7% (dashed line). Experimentally, this value corresponds to the appearance of trimer structures in the system,¹⁹ while only pendular bridges are observed for lower W. Here, E_p = 250 kPa in both experiments and simulations.

Bridge network. To match the experimental value of the average coordination number (i.e. bonds per particle), nearby particles were considered to share a bond at the start of the simulation only if their initial surface-to-surface distance was smaller than 0.064D; this condition leads to Z = 7.4 for particles in the interior of the simulated packing. Aiming for a model with as few parameters as possible, we assign all bridges in each sample the same values of k^bond_n and k^bond_t, estimated from the average diameters D and d and using δ⁰ = 2 μm, unless otherwise mentioned. In other words, for the purposes of setting up the bond stiffnesses we treat all bonds as identical, rather than introducing a function that would depend on their specific initial geometries. Spring constants corresponding to samples A and B are given in Table 1. These were evaluated by a look-up-table of discrete values from the FEM calculations, giving a precision of about 2%.
Compression test. Once a model sample has been configured it is confined between two flat, rigid walls. Uniaxial compression tests are then simulated by keeping the top wall fixed and moving the lower one upwards at a constant speed. For example, sample A was compressed at 46.4 mm s⁻¹, corresponding to a speed of 10⁻⁴ in the non-dimensionalised terms of the simulation,¹⁵ and the same dimensionless speed was used for all simulations. This velocity is a compromise between requiring the dynamics to be slow enough to reproduce the quasi-static experiments, while keeping a reasonable computation time.

3.2.2 Model for the interactions. The DEM simulates the dynamics of a collection of spherical particles that can interact with each other through contact forces (glass–glass) and cohesive bonds (PDMS), modelled as linear, Hookean springs. Each bead, i, is represented by a sphere with mass m_i, centre position r_i, diameter D_i, angular velocity ω_i and moment of inertia I_i. The particle can interact with its neighbours through normal forces and torques. The resulting translational and rotational equations of motion arewhere n_ij = (r_i − r_j)/|r_i − r_j| is the normal unit vector between particles i and j and where their interaction force is decomposed into a normal component, Fⁿ_ij, and any tangential forces, F^t_ij.

The normal force between two particles, Fⁿ_ij, has three possible contributions: a contact force due to particle overlap, F^c_ij, an elastic force due to a bond, F^bond_ij, and a dissipation force, F^diss_ij. The first two of these are conservative, and expressed as

where δⁿ_ij = |r_i − r_j| − (D_i + D_j)/2 is the surface separation of the particles. The terms involving k^bond_n are only applied if there is a bond between the two particles. For bead overlap we set the spring constant k^glass_n = 10k^bond_n, which is high enough to represent a system of stiff beads connected by softer springs, but low enough to avoid numerical instabilities; we have confirmed¹⁵ that simulation results are robust to changes in k^glass_n. Finally, we include dissipation due to particle deformation as F^diss_ij = −ζ(v_ij·n_ij), where v_ij is the relative velocity between beads i and j, and ζ is a dissipation rate.¹⁵

Tangential forces can also arise between particles which share a bond, or which overlap. For this, a linear restoring force is generated by any relative tangential displacement of two beads around their first point of contact (i.e. initial positions, for a bond). A cohesive bond is treated as a spring with a tangential spring constant k^bond_t. For particle overlap, as with normal forces, we use k^glass_t = 10k^bond_t. In either case the tangential force is limited by a Coulomb friction law with a friction coefficient μ = 0.5.³²

The interaction of a particle with either the top or bottom wall is treated the same way as with another particle, although no bonds need to be considered, only overlap. The total compressive stress in the system, σ, is then the sum of all the forces acting on a wall, divided by the initial cross-sectional area of the sample.

Finally, under significant enough tension or shear, a PDMS bond will fail; fracture testing has shown that the dominant mode of failure is by bonds detaching, or peeling away from the glass beads.²⁰ For the corresponding mechanism in our model we adopt a Griffith-like failure criterion, by comparing the strain and surface energies of a bond. As with the DEM simulations, the bond is here treated as elastic and incompressible. This leads to the failure condition¹⁵

For this, the bond is treated as a cylinder of height l_ij. Heterogeneity in the failure condition is provided through the distribution of values of l_ij. The normal displacement of the bond is Δr_n = δⁿ_ij − δ_ij⁰, and G = 7 J m⁻² is the measured interfacial energy between PDMS and glass.^20,33 Whenever a bond meets this criterion, it is removed from the simulation.

4 Elasticity of cohesive granulates

Having measured the constitutive relationship of the individual bonds between the particles of a cohesive granular material, and outlined a model of this class of material that is designed around these microscopic measurements, we now test the predictions of this model. Specifically, we will compare the stress–strain curves of simulations of uniaxial compression against experiments, and then explore how the elastic response changes with the stiffness and size of the bridges. A related manuscript¹⁵ makes use of the same model, but focuses on the macroscopic failure mechanisms, and how these depend on the packing density of the material.

4.1 Stress–strain curves

First, we compare the results of uniaxial compression in matched experimental and simulated geometries. The experiments are described in Section 3.1, with a typical sample shown in Fig. 1(c), and the compression tests follow protocols described in more detail elsewhere.¹⁹ We show in Fig. 3(a and b) the normal stress, σ, measured as a function of the compressive strain, ε, for samples A and B. The corresponding results from simulations in the same geometries (i.e. using values from Table 1 and particle positions taken from X-ray tomograms) are given in Fig. 3(c and d). Videos of the simulations are given as Movies S5 and S6 in the ESI.†

For both samples the model accurately captures the shape of the stress–strain curves, including the yielding behaviour at large ε, although the modelled system behaves about twice as stiffly as the experiments. Different realisations of the simulations show little variation in the predicted response. All curves start with a gradual build-up of stress, with a concave shape over a strain range corresponding to a displacement of ∼100 μm, or one bead radius. This regime corresponds to the progressive contact between the piston and the sample, whose surface is uneven on this scale, rather than any inherent non-linear elasticity.¹⁹ A linear response follows, from which the Young's modulus is measured by a least-squares fit. Here, we have used an extrapolation of the linear regime of deformation to define our zero-point of strain (i.e. ε = 0). With this convention, the linear regime extends up to about 4% in the experiments as well as in the simulations, after which point failure sets in.

In the simulations, we also tracked the number of bonds broken in each time step, binned into strain windows of size Δε = 4.2 × 10⁻⁴. These data, shown in Fig. 3(c and d) for representative cases, give a measure of the microcrack activity, which remains low until the end of the linear regime. This confirms that the response is essentially elastic until the peak stress or failure point is reached. Our companion paper¹⁵ explores the failure processes of cohesive granular materials in more detail; here, we simply note that failure proceeds by a shear band forming across the sample (see Movies S5 and S6 in the ESI†), as in the experiments.¹⁹

4.2 Variation with polymer stiffness

The elastic responses of diverse types of cohesive granular materials are known to be controlled by the cohesive component (i.e. the matrix), despite the fact that it represents only a small minority phase.^4,11,16,19 This is particularly true when the matrix material is softer than the particles or grains it binds together. For our samples the Young's modulus of the material, E, increases with, and is about an order of magnitude higher than that of, the stiffness of the PDMS composing its cohesive bridges, E_p.^19,20

In order to compare our model predictions to the experimental results, we systematically varied E_p in simulations of sample A, keeping the other parameters fixed. The results are shown in Fig. 4 alongside those from similar experiments¹⁹ (note that these used slightly larger beads with higher polydispersity: D = 210 μm ± 5%). The simulations correctly capture the scaling of E with E_p, including the deviation from a linear dependency of E ∝ E_p for stiff aggregates, although they again predict stiffer samples than seen experimentally. Investigating this discrepancy in more detail, we repeated the simulations using spring constants calculated assuming initial particle separations of δ⁰ = 1 μm (rather than 2 μm, see Section 2.3). The two different sets of spring constants, which span the range of behaviours observed by the micromechanical testing, had a significant impact on the stiffness of the aggregate, with E changing by about 60%. This increase is comparable to the increase in k^bond_n, and demonstrates the sensitivity of the mechanical properties of this class of materials to the exact geometry and nature of the bonds between particles.

4.3 Variation with polymer content

Finally, we focus on the stiffness of the aggregate when varying the volume fraction of polymer in the material, W. Experimentally, E increases approximately linearly with W up to a critical value of W* = 2.7%. After this point a much slower increase of the stiffness with polymer content is observed (see Fig. 5).¹⁹ This change coincides with the point where the polymer bridges coalesce into more complex structures like trimers (i.e. the pendular–funicular transition^19,34). In other words, above W* the addition of more matrix material does not strengthen the bonds between particles, but rather fills in the larger pore spaces.

To test the model's response in the pendular or capillary regime, we systematically varied the bridge diameter d in simulations of sample A, while fixing all other parameters. The spring constants k^bond_n and k^bond_t were calculated as in Section 2.3, for each value of d. Empirically,¹⁹ the volume fraction W ∝ d² (i.e. the bridge cross-section) and we used this to estimate W from d.

We compare the moduli E from the DEM simulations with experiments on similar samples (again, using D = 210 μm beads) in Fig. 5. The two sets of results agree up to the end of the pendular regime (dashed line in Fig. 5), although the simulations are again consistently stiffer than the experiments. As expected, the simulations diverge from the experiments in the funicular regime, W > W*, where the polymer in the experiments starts to form larger suprastructures, rather than simple independent bridges between particles. In this regime, different assumptions about the elasticity of the particle bonds would be required, as only a fraction of the matrix material will be supporting the stress.

5 Discussion & conclusion

This work has presented direct measurements of the micromechanical properties of cohesive granular media, where hard particles are connected by softer bonds. These measurements complement previous characterisation of the bulk elasticity and failure criteria of this type of material.^19,20 A simple discrete element model,¹⁵ was then used to mediate between the two limits, by upscaling the microscopic interactions in order to make predictions about larger assemblies of cohesive particles.

Despite their minimal set of ingredients, the DEM simulations correctly predict the relationships between the microscopic nature of the bridges and the macroscopic properties of a sample, and provide insights on the underlying mechanisms. Thus, they confirm how crucial are the details of the bond or matrix geometry to the elastic and fracture properties of cohesive granular materials. A small difference in the bridge length can significantly change, for example, the stiffness of the whole aggregate.

In particular, the simulations reproduce two major experimental results: the scaling of the material stiffness, E, with the matrix stiffness, E_p, over several orders of magnitude; and the scaling of E with the matrix content, W, within the pendular regime. The fact that the simulations diverge from the experimental results precisely at W = W*, which corresponds to the transition to the funicular regime, shows the role of suprastructures, such as trimers, in the softening of the samples at large W. Wang et al.¹¹ observed a similar transition in stiffness with increasing matrix content using millimetric beads and a harder cement, suggesting that a cross-over of stiffness with W is a fairly general result.

Nevertheless, the material appeared to be consistently stiffer in the simulations than experiments. This is arguably a result of our choice of a single value for the bridge length δ⁰. Changing δ⁰ from 1 to 2 μm while calculating the bond spring constants significantly decreased the stiffness of the material, as shown in Fig. 4, and it is clear that a realistic granular material would show a distribution of values for d and δ⁰. Indeed, a close inspection of tomography scans indicates the presence of bridges with a spanning distance significantly higher than 2 μm,²⁰ which would presumably contribute to a lower stiffness of the bridges and of the ensemble, although this would need to be confirmed in practice.

Interestingly, the model also correctly predicts the strain at which a sample ceases to have a linear response to compression and starts to fail (Fig. 3), which demonstrates the accuracy of the criterion used for breaking a bond in the simulations. Regarding eqn (4), we see that bonds are more likely to break in tension than in shear, and that the breakage depends on physical parameters, such as the Young's modulus and the interfacial energy. This is a key result for future modelling of failure in such a cohesive granular medium, or for the design of this type of material.

Finally, the numerical simulations also give access to quantities difficult to measure in experimental systems, such as the spatio-temporal dynamics of microscopic failures within the sample, and a companion article¹⁵ uses these signals to further explore failure in this model. In cohesive granular materials, brittle failure is expected to be abrupt with minimal precursory damage, while a transition into a ductile mode involves significant activity prior to catastrophic breakdown, which can be analysed for precursor signals. Studies of these materials will benefit from the present relatively simple minimal model, informed in all aspects by explicit microscale quantification, which accurately predicts elastic bulk response.

Conflicts of interest

A. H. and L. G. are inventors on a patent regarding the cohesive granular material described in this paper (WO/2018/014935).

Acknowledgements

We thank Matthias Schröter and Ran Holtzman for helpful discussions in the early development of this work and Soumyajyoti Biswas for generous input into the model design. Y. Y. was supported by the Leading Graduate Course for Frontiers of Mathematical Sciences and Physics, The University of Tokyo, MEXT.

Notes and references

1. J. Dvorkin, G. Mavko and A. Nur, Mech. Mater., 1991, 12, 207–217.
2. E. D. Cubuk, R. J. S. Ivancic, S. S. Schoenholz, D. J. Strickland, A. Basu, Z. S. Davidson, J. Fontaine, J. L. Hor, Y.-R. Huang, Y. Jiang, N. C. Keim, K. D. Koshigan, J. A. Lefever, T. Liu, X.-G. Ma, D. J. Magagnosc, E. Morrow, C. P. Ortiz, J. M. Rieser, A. Shavit, T. Still, Y. Xu, Y. Zhang, K. N. Nordstrom, P. E. Arratia, R. W. Carpick, D. J. Durian, Z. Fakhraai, D. J. Jerolmack, D. Lee, J. Li, R. Riggleman, K. T. Turner, A. G. Yodh, D. S. Gianola and A. J. Liu, Science, 2017, 358, 1033–1037.
3. M. Harrington and D. J. Durian, Phys. Rev. E, 2018, 97, 012904.
4. Y. Jung, S. Jung, S. Lee, W. Kim and H.-Y. Kim, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2018509118.
5. F. Kun, I. Varga, S. Lennartz-Sassinek and I. G. Main, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 062207.
6. F. Kun, I. Varga, S. Lennartz-Sassinek and I. G. Main, Phys. Rev. Lett., 2014, 112, 065501.
7. A. Neveu, R. Artoni, P. Richard and Y. Descantes, J. Mech. Phys. Solids, 2016, 95, 308–319.
8. J. Gaume, H. Löwe, S. Tan and L. Tsang, Phys. Rev. E, 2017, 96, 032914.
9. Y. Yamaguchi, S. Takada and T. Hatano, J. Phys. Soc. Jpn., 2018, 87, 094802.
10. M. Ozawa, L. Berthier, G. Biroli, A. Rosso and G. Tarjus, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 6656–6661.
11. W. Wang, J. Pan and F. Jin, J. Mater. Civ. Eng., 2019, 31, 04019047.
12. A. Shekhawat, S. Zapperi and J. P. Sethna, Phys. Rev. Lett., 2013, 110, 185505.
13. M. M. Driscoll, B. G.-g. Chen, T. H. Beuman, S. Ulrich, S. R. Nagel and V. Vitelli, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 10813–10817.
14. E. Berthier, J. E. Kollmer, S. E. Henkes, K. Liu, J. M. Schwarz and K. E. Daniels, Phys. Rev. Mater., 2019, 3, 075602.
15. Y. Yamaguchi, S. Biswas, T. Hatano and L. Goehring, Phys. Rev. E, 2020, 102, 052903.
16. J.-Y. Delenne, V. Topin and F. Radjai, Acta Mech., 2009, 205, 9–21.
17. S. Arif, J.-C. Tsai and S. Hilgenfeldt, J. Rheol., 2012, 56, 485–499.
18. W. Li, J. M. Rieser, A. J. Liu, D. J. Durian and J. Li, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 062212.
19. A. Hemmerle, M. Schröter and L. Goehring, Sci. Rep., 2016, 6, 35650.
20. A. Schmeink, L. Goehring and A. Hemmerle, Soft Matter, 2017, 13, 1040–1047.
21. L. Jing, Int. J. Rock Mech. Min., 2003, 40, 283–353.
22. M. Backholm and O. Bäumchen, Nat. Protoc., 2019, 14, 594–615.
23. A. Müller, M. C. Wapler and U. Wallrabe, Soft Matter, 2019, 15, 779–784.
24. R. H. Pritchard, P. Lava, D. Debruyne and E. M. Terentjev, Soft Matter, 2013, 9, 6037–6045.
25. M. Jiang, W. Zhang, Y. Sun and S. Utili, Granular Matter, 2013, 15, 65–84.
26. L. Brendel, J. Török, R. Kirsch and U. Bröckel, Granular Matter, 2011, 13, 777–786.
27. F. A. Gilabert, J.-N. Roux and A. Castellanos, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 031305.
28. P. Baud, T.-F. Wong and W. Zhu, Int. J. Rock Mech. Min., 2014, 67, 202–211.
29. S. Roy and M. S. Tirumkudulu, J. Rheol., 2016, 60, 559.
30. J. McBeck, K. Mair and F. Renard, J. Geophys. Res.: Solid Earth, 2019, 124, 9920–9939.
31. G. E. P. Box and M. E. Muller, Ann. Math. Stat., 1958, 29, 610–611.
32. I. Penskiy, A. P. Gerratt and S. Bergbreiter, J. Micromech. Microeng., 2011, 21, 105013.
33. J. Chopin, A. Prevost, A. Boudaoud and M. Adda-Bedia, Phys. Rev. Lett., 2011, 107, 144301.
34. M. Scheel, R. Seemann, M. Brinkmann, M. Di Michiel, A. Sheppard, B. Breidenbach and S. Herminghaus, Nat. Mater., 2008, 7, 189–193.

---

Footnotes

† Electronic supplementary information (ESI) available: 6 movie files, 2 data tables and a pdf document with captions for ESI files. See DOI: 10.1039/d1sm00458a

‡ These authors contributed equally to this study.

---