# Stability of cohesive granular heaps: Contact Dynamics simulations

Lydie Staron, Fanshuo Ma, Pierre-Yves Lagrée

► To cite this version:

Lydie Staron, Fanshuo Ma, Pierre-Yves Lagrée. Stability of cohesive granular heaps: Contact Dynamics simulations. Powders & Grains 2025, Dec 2025, Candolim, Goa, India. pp.02023, 10.1051/epj-conf/202534002023 . hal-05442297

HAL Id: hal-05442297

https://hal.science/hal-05442297v1

Submitted on 5 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Stability of cohesive granular heaps: Contact Dynamics simulations

Lydie Staron¹,*, Fanshuo Ma¹, and Pierre-Yves Lagrée¹

¹Sorbonne Université, CNRS - UMR 7190, Institut Jean Le Rond d'Alembert, F-75005 Paris, France

**Abstract.** Using Contact Dynamic simulations, we study how a rupture length in the adhesive interaction at contact, inducing hysteresis, reflects on the stability of cohesive systems. We examine the combined effect of the adhesive range and the algorithm performances. We observe how long-range adhesive interactions increase the stability of a pile under gravity, compared to short-range adhesive interactions, and compare their increased sensitivity to algorithm precision with the short-range case. The results suggest that long-range adhesive interactions are more prone to propagate local computational imprecisions associated to erroneous displacements, the way some material propagate dislocations, thereby recalling a kind of numerically-induced ductility.

## 1 Introduction

In nature or in the lab, material cohesion can have different origins. Beside the intensity of adhesion (van der Waals forces being for instance generally stronger than capillary bridges), they mostly differ in the range of the force at contact, and the existence of an hysteresis. While van der Waals forces are short-range and reversible, capillary bridges are long-range (namely acting up to a distance of typically $10^{-1}d$) and characterised by an hysteresis. This means that, in the latter case, an adhesive contact forms when two grains touch, but continue existing as long as the separation between the grains does not exceed a finite distance, called debonding or rupture length.

The size of debonding lengths increases significantly the shear strength of the material [1]. In this contribution, a numerical illustration of this effect is given for gravity-induced stability: "wet-like" long-range adhesion leads to more stable pile than "dry-like" short-range adhesion does. Moreover, considering displacement, volume variations, the sensitivity of the gravity-driven slumping to the Contact Dynamics algorithm performance is assessed. In particular, it is shown how long-range adhesion is more prone to propagate numerical error than short-range adhesion is.

## 2 Simulations of gravity-driven slumps

### 2.1 Numerical method of Contact Dynamics in 2D

The granular systems are simulated applying a Contact Dynamics algorithm in two dimensions, following the protocol described in [2]. The details of the algorithm will be found in Radjai & Richefeu (2009) [3].

Systems are made cohesive by the introduction of an adhesive contact force threshold $F_{adh}$, allowing contacts to

outlive tensile forces up to this limit (Fig. 1). The threshold $F_{adh}$ is expressed as dependent on the grains weight through the granular Bond number $B_{og}$ [4]:

$$F_{adh} = B_{og}m_{ij}g, \quad (1)$$

with $m_{ij} = 2 \times (1/m_i + 1/m_j)^{-1}$, where i and j are the two grains at contact.

In addition, an hysteresis is introduced in the contact interaction in the shape of a rupture or debonding length: while an adhesive contact forms when two grains touch, it does not vanish before the two grains are separated of a finite distance $l_r$ (Fig. 1). This mimics the existence of a capillary bound stretching until it breaks. The rupture length is simply given as a fraction of the contacting grains mean diameter $l_r = k_r(d_i + d_j)/2$.

In the following, the intensity of contact adhesion is varied by tuning $B_{og}$; typically, $B_{og} = F_{adh}/mg$ varies in the interval [30:360]. Meanwhile, the range of the adhesive force will be varied by tuning $k_r$, taking successively the value $k_r = 0$ and $k_r = 10^{-1}$.

The contacting grains also interact through Coulomb friction, involving the coefficient of friction $\mu$, set to $\mu = 0.2$ and not varied in this paper. The coefficient of restitution is set to zero.

### 2.2 Algorithm performance

The Contact Dynamics algorithm is implicit, relying on an iterative process to compute contact forces, and requiring a convergence criteria. For each time step, the convergence criteria need be met for the iterative process to end, thus signalling the success of the computation. Meanwhile, a maximum number of iterations is imposed for the force computation process to stop in case of the convergence criteria not being met, to avoid unduly large computation. Let's call `itmx` this number. Obviously,

<footnote>*e-mail: lydie.staron@upmc.fr</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** left: Signorini's graph shifted toward negative values to account for adhesion and right: schematic representation of the rupture length $l_r$ implemented through the local definition of an effective radius $d(1 + k_r)$, with $l_r = k_r d$.</figcaption></figure>

the greater this number, the higher the precision of the algorithm. In the following, we consider `itmx=300` and `itmx=1000` and study its effect on systems stability.

## 2.3 Experiments principle

The stability of squat cohesive granular columns is tested under gravity by allowing them to stand or collapse unchecked onto a rough horizontal surface (as in Fig. 2). We first generate rectangular granular steps by letting 7955 grains rain and relax in a container of width $W = 130d$ and height $H = 60d$. Grains have a diameter uniformly distributed in the range [0.004:0.006] to disrupt crystal-like ordering; $d$ denotes the mean diameter ($d = 0.005m$). The contact adhesion is set to zero during raining and ensuing relaxing to achieve a dense state (volume fraction $\phi \approx 0.82$). When rest is achieved, the adhesion is set to $B_{og}=F_c/mg=100$ during several computation steps to stick grains together.

The right-hand-side wall closing the container is removed only at the start of the collapse simulation, when the Bond number is set to the chosen value, between 30 and 360. The first time step of computation immediately equilibrates adhesion forces, smoothing the initial discontinuity in adhesion intensity. The influence of the initial state on the ensuing failure/stability, thus minimised, remains nevertheless a question of interest, yet beyond our scope. Once the right-hand-side wall removed, the system evolves as a result of gravity during 10000 time-steps of $dt=2 \cdot 10^{-4}s$, corresponding to a duration $T \approx 11.5 \sqrt{H/g}$, which is sufficient for a failure to occur or stability to manifest itself.

In addition to the Bond number $B_{og}$ taking values between 30 and 360, the rupture length $l_r$ is also varied and set to 0 (i.e. no hysteresis) and 0.1 (corresponding to typical capillary bridges).

## 3 Adhesion range and packing stability

In this section, we consider a precision-wise demanding algorithm `itmx=1000`. We increase the contact adhesion embedded in the granular Bond number from $B_{og} = 30$ to $B_{og} = 360$, considering both short-range $l_r/d = 0$ and long-range $l_r/d = 0.1$ contacts interactions. For each value of $B_{og}$, the system is left free to adapt the gravity field, either yielding of standing, during the duration $T$. Fig. 2 shows

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Snapshots of the final state of the system for seven values of the Bond number $B_{og} = 30, 100, 120, 140, 200, 260$ and 360 with zero rupture length $k_r = l_r/d = 0$ (left), and a finite rupture length $k_r = l_r/d = 0.1$ (right). Grains displaced a distance larger than $1d$ are colored in black. `itmx = 1000`.</figcaption></figure>

seven snapshots of final states for $B_{og} = 30$ to 360, corresponding to unstable or stable states, for both short-range and long-range contact adhesion.

Different possible signatures of the systems stability are now discussed for both adhesion ranges, including the effect of algorithm performance.

### 3.1 The mass centre shift

A blatant evidence that a system has yielded is the displacement of the material. In very unstable cases (e.g. $B_{og} = 30$ in Fig. 2), detached grains are flowing sideway. For larger cohesion, the mean deformation may be invisible to the naked eye. We also observe motion in the bulk (black areas) without clear failure localisation.

We hence quantify in the simplest manner the stability of the system by computing the total displacement of the centre of mass $G$ over the simulation duration $T$:

$$ \frac{R_G}{d} = \frac{|r_G(T) - r_G(0)|}{d}. \qquad (2) $$

The outcome is displayed in Fig. 3(a). In the case of short-range adhesion, we observe a sharp transition between unstable piles ($R_g/d \ge 4$) and the seemingly more stable group ($R_g/d \le 2$). The latter however, in the light of Fig. 2, includes some badly damaged cases. In the case of long-range adhesion, seemingly stable piles ($R_g/d \le 1$) coincid-
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Total displacement of the systems centre of mass $R_G/d$ and (b) total volume variation $\epsilon = (V(T) - V_0)/V_0$, as a function of the granular Bond number $B_{og}$ for short-range adhesion $\ell_r/d=0$ and long-range adhesion $\ell_r/d=0.1$ (itmx=1000).</figcaption></figure>

ing with a damaged state are also observed (for instance for $B_{og} = 120$). The mass centre shift thus seems an insufficient proxy of systems stability. We can nevertheless conclude that long-range interaction have a stabilising effect on the systems. A second conclusion is that additional measures of pile's stability are needed.

## 3.2 The volume variation

Another measure of piles stability can be searched for in volume variations over the simulation duration $T$:

$$ \epsilon = \frac{V(T) - V(0)}{V(0)}, \qquad (3) $$

where the volume $V$ is estimated from the profile of the deposit given by the grains forming the top line.

The outcome is displayed in Fig. 3(b). For both short- and long-range adhesion, transitions between dilating ($\epsilon \ge 0$) and shrinking ($\epsilon < 0$) systems are observed: unstable piles tend to expand, while more stable systems tend to condense. The transition between the two trends does not coincide with transitions in the behaviour of $R_g$: $\epsilon$ designates a larger group of dilating systems than simply including those recognised as unstable by the shift of $R_g$. When compared with Fig. 2, we conclude that $\epsilon$ rather renders the level of damage of the piles.

At this stage, it is however apparent that stability of a collection of grains under gravity is an elusive mechanism [2].

## 3.3 Algorithm performance and grains overlap

Because the systems are numerical, strictly numerical measures can give valuable indications on their evolution. We chose to evaluate the mean grains overlap at contact $\langle \delta \rangle$, which reveals how efficient the program is in guaranteeing the rigidity of the grains: it provides a picture of the algorithm's difficulties reflecting the system physical evolution (large stress, rapid changes etc).

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Mean grains overlap at contact computed over the total durations $T$ of the simulations, for $\ell_r/d = 0$ and $\ell_r/d = 0.1$ (itmx=1000).</figcaption></figure>

The mean overlap is computed over the whole duration $T$ of the simulations. The outcome is displayed Fig. 4. The fluctuations of $\langle \delta \rangle$ with the adhesion strength $B_{og}$ make it difficult to draw definite conclusions on the pile stability. However, it clearly appears that larger adhesion induces larger overlap. In other words, the existence of adhesion is challenging the algorithm efficiency.

A spontaneous question thus arises: how does the algorithm efficiency challenge the stability of the cohesive piles? To address this question, we perform the same simulation series as discussed in Sec. 3, only changing the exit condition of the resolution loop: the number of iteration is lowered from itmx=1000 to itmx=300. In other words, the demand on algorithmic convergence is lowered.

## 4 Algorithm performance and systems stability

Fig. 5 and Fig. 6 show the changes operated when lowering the efficiency of the solver. Changes are more significant for $\ell_r/d = 0.1$. In this case, the systems mass centre is less displaced for small adhesion, as reflected by $R_g$ (Fig. 5), while no damage is seen in the packing (Fig. 6). However, for large adhesion, $R_g$ displays a non-negligible displacement coinciding with a marked shrinking of the bulk (Fig. 6(a) and (b)). Fig. 5 shows the corresponding diffuse motion propagating in the bulk while adhesion increases. In the extreme case $B_{og} = 360$, the bulging of the pile's right side becomes apparent.

The corresponding algorithm's performance in terms of grains overlap at contact is displayed in Fig. 7. We observe how increasing contact adhesion affects the overlap significantly for $\ell_r/d = 0.1$, although its value remains comparable with the case itmx=1000. Interestingly, for all systems, the negative volume variations (shrinking) cannot be accounted for by the volume absorbed in overlapping contacts. The latter, evaluated from the mean overlap $\langle \delta \rangle$ and the number of contact (between 13000 and 15000), falls one order of magnitude under the figures given by the evolution of $\epsilon$ (not shown here). This means that the differences of the two behaviours, i.e. between
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Snapshots of the final state of the system for a less efficient algorithm (itmx = 300 instead of itmx = 1000), for the same values of the Bond number $B_{og}$ as displayed in Fig. 2, with rupture length $k_r = \ell_r/d = 0$ (left) and $k_r = \ell_r/d = 0.1$ (right). Grains displaced a distance larger than $1d$ are colored in black.</figcaption></figure>

itmx=300 and itmx=1000, is the result of an intricate reorganisation of grains position, obeying a different force network.

In this respect, the case of $\ell_r/d = 0.1$ for large adhesion $B_{og} = 120$ to $B_{og} = 360$, is particularly intriguing. While the pile does not fail, in the case of itmx=300, the pile is like "collapsing inwardly" (Fig 5). Important rearrangements are seen in the bulk, and propagate all the more that the contact adhesion is strong. This occurs only for a lax algorithm with itmx=300, for which poorly resolved forces occur more often than in the case itmx=1000. It is noteworthy that bulging appears at very large adhesion, when the grains displacement map shows that most of the bulk undergoes some motion. Interestingly, rather than showing only in the degradation of grains rigidity (reflected by $\langle\delta\rangle$), it shows in the lessen rigidity of the whole packing who slumps and seems to melt. This alerts us on the fact that poor computation performances do not only show where one would expect, but may lead to a whole different collective behaviour.

## 5 Conclusion

We have compared the stability of cohesive piles varying the intensity of contact adhesion, and its range. We observe how longer-range adhesion induces larger stability. Meanwhile, the outcome of two algorithms differing only

<figure><img src="image_7.png" /><figcaption>**Figure 6.** (a) Total displacement of the systems centre of mass $R_G/d$ and (b) total volume variation $\epsilon = (V(T) - V_0)/V_0$, as a function of the granular Bond number $B_{og}$ for short-range adhesion $\ell_r/d = 0$ and long-range adhesion $\ell_r/d = 0.1$, for itmx=300. Corresponding evolutions for itmx=1000 are shown in gray.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Mean grains overlap at contact computed over the total durations $T$ of the simulations, for $\ell_r/d = 0$ and $\ell_r/d = 0.1$, for itmx=300. Corresponding evolutions for itmx=1000 are shown in gray.</figcaption></figure>

by their performance was discussed. We observe that the stability of the system deteriorates with lower algorithm performances, and that long-range adhesion is more sensitive to the latter. In particular, we see how displacements propagate in the bulk without failure opening, inducing the whole system to shrink and slump. This suggests that long-range adhesive contact networks are more prone to propagate local computational imprecisions, the way some material propagate dislocations, thereby recalling a kind of numerically-induced ductility.

## References

[1] Khamseh, S., J.-N. Roux, F. Chevoir, Phys. Rev. E **92**, 022201 (2015).

[2] F. Ma, P.-Y. Lagrée, L. Staron, Phys. Rev. E **111**, 055401 (2025).

[3] F. Radjai, V. Richefeu, Mech. Mater. **41**, 715-728 (2009).

[4] P.G. Rognon, J.-N. Roux, M. Naaïm, F. Chevoir, J. Fluid Mech **596**, 21-47 (2008).