# How vorticity and agglomeration control shear strength in soft cohesive granular flows

O. Bouillanne, Guilhem Mollon, Aurélien Saulot, Sylvie Descartes, Nathalie Serres, Guillaume Chassaing, Karim Demmou

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

## To cite this version:

O. Bouillanne, Guilhem Mollon, Aurélien Saulot, Sylvie Descartes, Nathalie Serres, et al., How vorticity and agglomeration control shear strength in soft cohesive granular flows. GRANULAR MATTER, 2022, 24 (2), 10.1007/s10035-022-01216-8. [hal-03659784](https://hal.science/hal-03659784)

HAL Id: hal-03659784

https://hal.science/hal-03659784

Submitted on 12 May 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# How vorticity and agglomeration control shear strength in soft cohesive granular flows

Olivier Bouillanne<fnref n="1" />,<fnref n="2" />, Guilhem Mollon<fnref n="1" />, Aurélien Saulot<fnref n="1" />, Sylvie Descartes<fnref n="1" />, Nathalie Serres<fnref n="2" />, Guillaume Chassaing<fnref n="2" />, Karim Demmou<fnref n="2" />

## Abstract

Deformable granular flows present complex kinematics. These materials can have various flow regimes: plastic, agglomerated, rigid-like granular flow, etc. In this paper, a multibody meshfree model is used to investigate the consequences of cohesion, stiffness, and viscosity of the particles on their collective sheared flows in tribological contacts. An approach derived from fluid mechanics postprocessing tools, based on vortex detection, is employed to understand the links between these parameters and the emerging friction coefficient of the sheared interface. These results pave the way to complete kinematic studies of third body simulations in tribological contacts.

**Keywords** Soft grains · Sheared flow · Agglomeration · Vorticity · Tribological contact · Third body

# 1 Introduction

## 1.1 Tribology and fretting

Tribology is the science that describes the interaction of the surfaces of two solid bodies in contact and relative motion. This contact generates friction and can lead to various mechanical phenomena, such as wear and fatigue. Fretting is at the center of many practical topics, from moving cut stones by the ancient Egyptians [1] to brushing teeth [2]. Mechanical systems can be lubricated or dry. While the understanding of lubricated contacts phenomenology is now fairly accurate, that of dry contacts remains fragmentary.

An approach popularized by Godet [3] and then Berthier [4, 5] is to consider that, at the interface between two rubbing bodies, there is a third body. This third body, composed of a combination of wear debris from the first bodies and particles from the outside, is located between the surfaces and allows loads to be transmitted while accommodating the relative velocity of the two rubbing bodies. This view

is supported by countless experimental characterizations of worn surfaces. Hereinafter, these two rubbing bodies will be called “first bodies”, following the conventional tribological terminology.

Fretting is a tribological phenomenon in which two surfaces in contact undergo oscillatory movements of amplitude smaller than the size of the contact. This fretting can generate fatigue or wear, or a mixture of both. Many mechanical parts undergo fretting: splines, bolted or riveted assemblies [6], etc. This affects many fields of engineering: aeronautics, nuclear, biomedical, civil engineering, etc. The concept of third body has led to a better understanding of fretting, especially following advances made by the fretting maps [7] and the study of the different fretting regimes (stick condition, stick-slip condition, and gross-slip condition) [8].

The third body may protect or damage surfaces and can show various aspects, as can be seen in Fig. 1. These SEM views are obtained as a result of two titanium alloy (Ti-6Al-4 V) surfaces rubbing against each other in a fretting-test tribometer, i.e. following a small amplitude reciprocating sliding motion during a large number of cycles. As in any such test, a third body is visible. On the left-hand micrograph, various third body behaviors are detectable: granular on the left, and rather plastic in the center. On the right-hand micrograph, the third body is granular and large agglomerated wear debris are visible.

<footnote>☑ Guilhem Mollon
Guilhem.mollon@insa-lyon.fr</footnote>

<footnote n="1">¹ LaMCoS, INSA-Lyon, CNRS UMR5259, Université de Lyon, 69621 Villeurbanne, France</footnote>

<footnote n="2">² Safran Aircraft Engines, Moissy-Cramayel, France</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1 SEM views of different areas of mechanical contact after friction tests, with various third body accommodation regimes</figcaption></figure>

Experimental approaches can reveal this third body after the test, but traditional mechanical models tend to neglect the impact that the third body can have on wear. The main reason for this is the fact that state-of-the-art measurement techniques are not able to resolve in time and space the variability of the local third body behaviors and their influence on *local* friction and wear. There is therefore a lack of knowledge of such local effects, and the best way to circumvent this experimental bottleneck is to employ numerical simulations.

Conventional numerical methods of modeling wear are descriptive and hardly predictive. This is due to the multitude of factors influencing interface behavior: the amplitude of the movement between the two surfaces, the pressure distribution on these surfaces, the number of accumulated cycles and their frequency, the properties of the materials and the condition of their surfaces, the contact temperature, the atmospheric environment... Collins [9] thus counts about fifty variables that play a role in the fretting process. The relationship between the third body and wear is established but not yet totally predictable [10, 11].

Surface degradation leads to wear after a while, especially if the particles that have become detached from the surface come out of contact permanently. If these particles remain within the contact, they form the third body and can help to limit this wear. The stresses applied by the flow of the third body in the interface are probably one of the causes of surface degradation. It is therefore necessary to characterize in detail the stresses applied to the first bodies by the third body. Such forces are obviously related to the flow regime of the third body in the interface, especially on its shear strength and on its spatial and temporal variations. Kinematic indicators will allow us to characterize these regimes

and to understand the link between the third body and the stresses, then with the degradation, and finally with the wear.

The kinematics of third body flows must therefore be better understood, to grasp the interaction between the third body and stresses. Numerical [12] and experimental [13–15] approaches to calculate the vorticity of a granular flow have already been proposed, but they were limited to rigid granular materials. The terms *agglomerate* and *vortex* are complementary. The term *agglomerate* refers to a set of particles moving in a coordinated manner, while a *vortex* refers to a set of particles rotating around a certain geometrical point. Hereafter, in this study, the term agglomerates or structures will be used, because the vortices studied here have much longer duration than those described in the existing granular physics literature.

## 1.2 Numerical tribology

Different numerical approaches are used to describe the tribological contacts. The first one is based on finite elements. To model the wear track, one way to do this is to calculate the stresses undergone by the surfaces and to deduce a displacement of the nodes composing the mesh of these solids. [16]. Another way, more radical, consists of recalculating a complete mesh as the wear progresses [17]. The "wearbox" numerical tool, on the other hand, relies on a measurement of the amount of energy dissipated in contact to apply a deformation to the mesh [18, 19]. The shortcoming of these methods is that they do not consider the debris bed present in the interface, i.e. the third body. Ding [20] tried to take this into account in a new finite element model (without representing this third body), but the wear tracks that can be simulated are shallow. Semi-analytical approaches are
also used [21, 22], which allow reducing the computation
time and simulate the full system. The reader can refer to the
article [23] for more information on the numerous methods
used in finite elements.

Although finite element models allow modeling the stress
and displacement fields of bodies in contact, they do not
allow considering the discontinuities of the material. The
discrete element approaches used in tribology have their ori-
gins in the numerical modeling of discontinua of Cundall
and Hart [24]. These simulations, while quite varied, are
based on three key principles:

* The equations of motion of mechanics;

* The detection of contacts between particles;

* The calculation of the forces between the particles.

Discrete element models (DEM) are used in tribology
to model different aspects of contact. For example, DEM
can be used to model the contact and the degradation of
materials between the wheel of a train and the rail [25].
DEM simulations are also used to simulate the flow of the
third body—as well as its degradation—in different con-
texts: thermo-mechanical response [26], two-phase materials
[27], or electrical contacts [28, 29].

Some authors tried to combine discrete elements with
finite elements to combine the advantages of both methods
[30–32]. In particular, the particles of a discrete model are a
priori rigid, which is not the case when combining these two
methods. The field of deformable granular materials is of
growing interest, both experimentally [33, 34] and numeri-
cally [35–38], and extends well beyond the purely tribologi-
cal applications.

Finally, the last type of approach concerns nanometric
models. Several techniques exist, based on the displace-
ment of atoms. The Monte Carlo method and the molecular
dynamics method are the most widespread. These simu-
lations are based on atomic interaction laws derived from
quantum mechanics [39]. Some authors focus on modeling
adhesive wear and the formation of agglomerates in con-
tact [40, 41]. These simulations are informative but time-
consuming, and only allow representing a tiny fraction of
the contact. For more information on these simulations, the
reader may refer to [42].

**2 Methods**

**2.1 Numerical model**

This work is largely based on a previous work [43], where an
extensive simulation campaign was performed, with various
mechanical parameters. In what follows, we revisit these
stored numerical results with dedicated techniques to extract

relevant new kinematic data. It is therefore worth provid-
ing some general information about the model. The reader
interested in more details may refer to [43]. In that paper,
the author describes a simulation of a flow of deformable
particles in a dry sliding interface. These simulations are
performed on the open-source code MELODY, a software
developed in LaMCoS, INSA Lyon, France [44, 45]. This
software generalizes the conventional Discrete Element
Method (DEM) and allows simulating complex granular
materials, with either rigid or highly deformable grains with
arbitrary shapes, within the framework of large deformation
hyperelasticity in 2D. The versatility of this software allows
it to be used in several fields, for example in mechanical
engineering [46, 47] or geomechanics [48].

2.1.1 Principles

Figure 2 shows a sketch of the numerical model, which
is purely dimensionless. Two rough rigid walls (length
$L = 100$) are placed on either side of a collection of 2000
deformable particles with a unit average diameter, to repre-
sent two bodies rubbing against each other with a third body.
The two first bodies have a self-affine roughness, with a $R_α$
equal to 1 space unit. They are rigid and non-degradable.
Their cohesion with particles is very high, in order to pre-
vent wall-slip effects, and to focus on velocity accommoda-
tion within the third body only.

The deformable particles follow a hyperelastic law. This
makes it possible to apply elasticity to bodies that undergo
large deformations and rotations. The neo-Hookean law con-
siders two main parameters, which are Young's modulus *E*
and Poisson's ratio *v*. The Poisson's ratio is fixed at 0.495
to impose a quasi-incompressibility. Finally, the density is
equal to 1 (dimensionless model). To make the simulation
stable, damping is added via the $\alpha$ Rayleigh parameter: the
damping matrix is proportional to the stiffness matrix via
this parameter. The result is a viscoelastic-like behavior of
the Kelvin-Voigt type. The contact law governing the inter-
actions between the particles and with the walls is governed
by a cohesive frictionless model with only one parameter
$c$, which is the strength per length unit of contact between
two surfaces. By using deformable particles that can move
freely in space, the model combines the laws of dynamics
with those of continuum mechanics.

A gap (noted *D*) between the two first bodies is imposed. Finally, a relative tangential velocity $V = 1$ is imposed on the two bodies: the upper body goes rightwards at $V/2$, and the lower one, leftwards at $-V/2$. This generates a shear within the collection of particles that will accommodate velocities while transmitting loads throughout the simulation, like the third body within a tribological contact. Periodic boundary conditions are applied to the extremities of the simulation. For each
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Sketch of the numerical model with two rigid first-bodies and a third-body represented by deformable grains. Colors are arbitrary and only serve to distinguish grains on the top and the left figure. Velocity field magnitude is plotted in the lower-right corner</figcaption></figure>

simulation, the simulated time corresponds to three sheared spatial periods, i.e. a time of 300-time units (called Δt and equal to 1 in this dimensionless context). To ensure that the flow regime of the third body is well established, only the last 200Δt will be considered for quantitative interpretation.

The normal and tangential forces experienced by the upper and lower bodies are measured over time and noted $F_N(t)$ and $F_T(t)$ respectively. Their time average is noted $\overline{F_N}$ and $\overline{F_T}$. The stresses, which are the ratio of the forces to the length L (2D simulation), are noted $\overline{\sigma}_N$ and $\overline{\sigma}_T$. The ratio of these two quantities provides an average coefficient of friction at the scale of the whole interface, as an output quantity of the simulation. This approach is a good conceptualization of the experimental reality where the friction coefficient of a given interface is a measured quantity that emerges from a large number of local phenomena and events. The elastic behavior of each grain can certainly be seen as a limitation when applied to plastic and ductile materials (e.g. metals), but constitutes an improvement with respect to perfectly rigid grains. Besides, irreversible deformations occur within the flow because of the cohesive character of the contact laws. It means that plasticity is present in the simulations, albeit at the scale of the flow instead of that of the grains.

### 2.1.2 Parametric space

Three dimensionless parameters are studied and lead to various flows of the third body.

$$ \tilde{E} = \log_{10} \left( \frac{E}{\overline{\sigma}_n} \right) \qquad (1) $$

$$ \tilde{c} = \log_{10} \left( \frac{c}{\overline{\sigma}_n} \right) \qquad (2) $$

$$ \tilde{\alpha} = \log_{10} \left( \alpha \cdot V/D \right) \qquad (3) $$

These three parameters are a normalized stiffness $\tilde{E}$ (between 0.15 and 1.2), a normalized cohesion $\tilde{c}$ (between -1 and 1.4) and normalized viscosity $\tilde{\alpha}$ (between -3.3 and -1.5). The first two parameters are normalized by the average normal stress experienced by the first bodies, $\overline{\sigma}_n$. The third parameter is normalized by the mean strain rate in the interface, $V/D$. The normalized stiffness characterizes the resistance to deformation of particles under stress. The normalized cohesion quantifies the cohesive force between the particles in the system. This cohesion therefore varies roughly between one order of magnitude below the normal stress (negligible cohesion) and one order
of magnitude above it (dominating cohesion). Finally, the normalized viscosity characterizes the response of the grains to their rate of deformation. As demonstrated in [43], the chosen range for $\tilde{\alpha}$ is consistent with the solid viscosity of steel (for example) under typical tribological conditions.

### 2.1.3 Consequences of stiffness, cohesion, and viscosity of third body

The coefficient of friction, noted $\mu$, represents the ratio of the measured tangential force to the normal force. This is an output value, not an input value. It characterizes the shear strength of this confined soft flowing granular material.

$$ \mu = \frac{\overline{F_T}}{\overline{F_N}} \qquad (4) $$

The evolution of the friction coefficient as a function of the values of $\tilde{E}$, $\tilde{c}$ and $\tilde{\alpha}$ can be seen in Fig. 3. The coefficient of friction has a wide range, from 0.1 to over 1.0, but remains in the typical range of observed values for dry sliding. The increase in the $\tilde{\alpha}$ parameter leads to a general increase in the coefficient of friction. The two maps of Fig. 3 show the coefficient of friction at imposed $\tilde{\alpha}$. In particular, when $\tilde{\alpha} \approx -1.8$, a maximum coefficient of friction is observed, at $\tilde{E} = 0.7$ and $\tilde{c} = 0.7$ (that is, for a Young modulus of the grains and an interparticle cohesion roughly one order of magnitude larger than the applied normal stress).

These three parameters greatly influence the flow and the microstructures of the third body. Figure 4 shows different kinematic patterns of the third body as a function of the parameters $\tilde{E}$ and $\tilde{c}$, at imposed $\tilde{\alpha} = -1.8$. For example, when $\tilde{E}$ and $\tilde{c}$ are low, the observed behavior is close to

picture A. The particles deform a lot and fill the pore space. They adopt a structured pattern with a preferred horizontal direction which promotes a perfectly laminar flow: because of the low cohesion, each layer of particles slides over the others. On the other hand, in the case of B, the particles deform little. The observed regime can be described as granular. Cases D and F are also particularly interesting. In these cases, the particles have sufficient cohesion to form agglomerates, similar to those seen in Fig. 1. These third body agglomerates flow, shear and roll into the interface. The other cases (C, E, H, and I) are intermediate to those described above. These simulations are quantitatively consistent with the experimental variety of third body visual aspects [49].

## 2.2 Coherence

In some simulations, agglomerates are observed. Grains form various sizes of cohesive clusters, in which the particles move in a motion similar to that of a rigid body. An example can be seen in Fig. 5. The agglomerates are difficult to observe in a succession of still images but are clearly visible in animations.

The characterization of the size of these agglomerates would make it possible to analyze the consequences of the rheology of the third body particles on the first bodies. Tools developed to characterize this size of agglomerates, and derived from fluid mechanics are documented and adapted hereafter.

The agglomerate identification function is derived from vortex identification functions. These mathematical descriptors were proposed in the context of vortex studies

<figure><img src="image_5.png" /><figcaption>Fig. 3 Coefficient of friction in function of $\tilde{E}$ and $\tilde{c}$, for $\tilde{\alpha} = -1.8$ and $\tilde{\alpha} = -3.2$</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4 Microstructure of third body as a function of $\tilde{E}$ and $\tilde{c}$ (from [43])</figcaption></figure>

in fluid mechanics and are called $\Gamma_1$ and $\Gamma_2$ [50, 51]. These functions compare the observed motion of points near a central point with that of particles forming an ideal vortex. The function proposed here, called “coherence”, is based on the same concept but compared with the motion of an ideal rigid body. To do this at a given time step, a velocity field is first interpolated on a regular Eulerian grid from the motion of the particles, as can be seen in Fig. 6. The interpolation used is cubic, which allows $C^2$ continuity. The velocity is interpolated with respect to the particle

center-of-mass velocity, because the deformation-related velocity field in particles is small compared to the rigid body velocity of these particles, as can be seen in Fig. 2. The time step used to extract grains positions and to compute velocities is equal to one time unit, and corresponds to an average strain level of 0.05 between two extractions. This is a good compromise between too small (producing a large amount of redundant data) and too large (overlooking complex trajectories of the grains) time steps.
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Snapshots of a simulation at three successive times, for $\bar{E} = 0.24$, $\bar{c} = 0.94$, $\bar{\alpha} = -1.65$. Several rigid body movements can be observed</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6** Velocity vector of each particle in the simulation at time $T = t + 2\Delta t$ (black arrows), and interpolated velocity field from the particle velocities (colored field)</figcaption></figure>
The velocity field is discretized into a square mesh, with a spatial step of $\Delta x = \Delta y = 1$, and each point of this mesh has a velocity $\vec{V}$. We then define a local quantity for each point, called “coherence”, which quantifies to what extent the velocity field in the neighborhood of this point corresponds to a rigid-body motion. Coherence at a given point $P$ is therefore calculated by comparing, at each point $P$, the velocity of each point $M$ located in a disk of radius $R$ and center $P$, with the velocity that an ideal rigid body of center $P$ and radius $R$ would have. The complete formula is as follows:

$$C(x, y, R) = \frac{\int_{M \in S} (V - V_{ave}) \cdot V_{circ} dS}{\sqrt{\int_{M \in S} (V - V_{ave})^2 dS} \cdot \sqrt{\int_{M \in S} (V_{circ})^2 dS}} \quad (5)$$

where $\overline{PM} \le R$, $V_{circ}$ is an ideal circumferential relative velocity of point $M$ with respect to point $P$ in a rigid-body rotation (the angular velocity is arbitrarily chosen to be equal to 1), and $V_{ave}$ is the average velocity of all points $M$ within the radius $R$ (cf. Fig. 7). By construction, $-1 \le C < 1$. A value of $C = 1$ means the particles rotate in a coordinated rigid-body motion in the positive direction, and vice versa for $C = -1$. The radius $R$ of the interrogation circle is a user-defined parameter. It seems that a radius of two to three times the typical size of a particle allows the detection of structures properly while keeping the computation time reasonable.

<figure><img src="image_9.png" /><figcaption>Fig. 7 Coherence calculation scheme</figcaption></figure>

By repeating the operation on the whole simulation, it is thus possible to determine a coherence field, as can be seen in Fig. 8A, computed with the velocity field extracted from Fig. 5. Red zones correspond to a coherence close to 1, and blue ones to -1. Expectedly, the coherence map is essentially negative because the natural direction of rotation due to the shearing of the interface is negative in the present sign convention. By applying a filter to select only the areas where $|C| > 0.8$, agglomerates are highlighted (Fig. 8B). It is finally possible to stack these detected areas to observe their evolution over time (Fig. 8C). Independent structures are easily identified and labeled based on connectivity relationships in the space-time domain.

Each detected agglomerate is characterized by its coherence. An additional relevant property to be calculated is the average angular velocity, noted $\omega$. Each agglomerate detected at an instant $T$ can be associated with a certain number of grains belonging to it. The average angular velocity is calculated by computing the average of the circumferential component of the velocity of each grain to the barycenter of the agglomerate.

$$\omega(x, y, R) = \sum \frac{V_{circ}}{R_i} / n \quad (6)$$

With $R_i$ the distance between each point in the structure to the barycenter and $n$ the number of points in. To consider the difference in distance between the first two bodies according to the simulations, this quantity is normalized by the large-scale shearing of the third body layer, $V/D$, and noted $\tilde{\omega}$.

$$\tilde{\omega} = \frac{\omega}{V/D} \quad (7)$$

In this case, the $V/D$ shear rate is negative. Structures rotating in a positive direction will therefore have a positive value of $\tilde{\omega}$. Normalization makes it easy to compare the angular velocity of the structures with the shear rate. If $\tilde{\omega} = 1$, the structure rotates as a disk in a rolling movement without sliding in the contact. If $\tilde{\omega} > 1$, it rotates faster, and if $0 < \tilde{\omega} < 1$, the structure rotates slower. Finally, if $\tilde{\omega} < 0$, the structure rotates in the opposite direction to the shear.

## 3 Results

75 simulations were performed with different values of $\tilde{E}$, $\tilde{c}$ and $\tilde{\alpha}$ in [43], and we revisit their stored numerical results. Coherence is applied to this series of simulations used to evaluate the coefficient of friction as a function of the three parameters identified above. Seven notable simulations are retained and are the same as those used in the previous study. They are noted from A to G and located in the parametric space in Fig. 4.
<figure><img src="image_10.png" /><figcaption>**Fig. 8** A Coherence field extracted from the velocity field (Fig. 6); B Filtered coherence field at the same time-step at $|C| > 0.8$; C Stacked filtered coherence over time</figcaption></figure>

**3.1 Stacked coherence**

In Fig. 9, the coherence is plotted over time and stacked, following the methodology shown in Fig. 8C.

In case A, the particles are soft and hardly cohesive. There are very few agglomerates detected, and they all rotate in the negative direction. Because of the shear generated by the two upper and lower bodies, this is naturally the preferred direction of rotation, while positive rotations correspond to agglomerates which rotates against the main shear direction. The lifespan of the structures is very short.

Cases B, C, and D follow a path of decreasing stiffness. C is located at the maximum of the friction coefficient. The other parameters do not vary very much. First of all, for case B, structures are not very spatially extensive but persist for a fairly long time. No observed structure rotates in the positive direction. For case C, on the other hand, the structures, although numerous, are not very extensive spatially and temporally. Some structures rotating in the positive direction are to be noted, unlike case A. Finally, in case D, the negatively

rotating structures are massive, both spatially and temporally. Structures that rotate in the positive direction are also observed, and their spatial extent is notable. We thus observe for these three cases, with decreasing $\tilde{E}$, first many structures, then a drastic reduction in their size, and finally the appearance of large, numerous, and persistent agglomerates.

Simulations E, C, and F follow an increasing cohesion path. Case E is similar to case A, except that there are more structures, and they have a slightly larger spatial extent. Case F is very interesting and shows the impact of an important cohesion. The structures are very large, whatever their direction of rotation. In these three cases, the increase in cohesion results in an increase in structures size. Cases E and C are similar, while regime F has very large structures, but is seemingly shorter-lived than in case D.

Case G, visible in Fig. 10, is similar to case C, except that the viscosity is much lower. The alpha viscosity has an important role in whether or not structures are formed since the detected structures have a longer lifetime and a larger spatial extent.
<figure><img src="image_11.png" /></figure>

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

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

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

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

<figure><img src="image_16.png" /><figcaption>Fig. 9 A–F Stacked filtered coherence over time for six characteristic simulations. The structures in red rotate in the positive direction, and vice versa for the blue ones</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 10** Stacked filtered coherence over time for G case. The structures in red rotate in the positive direction, and vice versa for the blue ones</figcaption></figure>

## 3.2 Size of structures

Figure 11 shows the proportion of grains belonging to a structure compared to all the grains in the simulation. This number of grains is very dependent on the cohesion $\tilde{c}$, and to a lesser extent, on the viscosity $\tilde{\alpha}$. On the other hand, the stiffness $\tilde{E}$ seems to have very little influence. The maximum grain content in a structure, 11%, is rather low. This is due to the choice of a rather high detection threshold ($C > 0.8$). The data should be analyzed relatively rather than absolutely. The proportion of grains in the structures is an important parameter since it gives information on the extent of the agglomeration phenomenon.

## 3.3 Angular velocity of structures

The probability density function (PDF) is used to extract statistical data from all the simulations, including the angular velocity of the agglomerates. These PDF are in terms of equivalent grains, i.e. they represent the density of probability for a given grain to belong to a structure with a certain diameter or a certain angular velocity. They are visible in Fig. 12, left column. The seven cases analyzed above have very different profiles.

Path B (purple) → C (medium blue) → D (green), which follows a decreasing stiffness, is interesting because we observe that case C (maximum peak of the coefficient of friction) corresponds to a minimum angular velocity and a minimum of grains in the structures. For path E (dark blue) → C (medium blue) → F (light blue), which follows an increasing cohesion, the case C is intermediate between E and F. Case G does not have very different characteristics from case C, which has the same stiffness and cohesion, but not the same viscosity. Case A, in orange, which has very few structures, shows only structures rotating in the shear direction, with very high dispersion. Thus, many structures rotate nearly twice as fast as the shear rate.

To each of the 75 simulations can be attributed an angular velocity mode $\tilde{\omega}$ corresponding to the peak of the angular velocity of the PDF function. Some correlations are interesting to observe, for example, between $\tilde{E}$ and $\tilde{\omega}$ in Fig. 13A. No correlation (Pearson's coefficient of 0.22, Pearson's coefficient is relevant to quantify the quality of linear regression) is noted between $\tilde{\omega}$ and $\tilde{E}$ for the whole simulations, but a correlation is visible when the data are binned according to the value of the normalized viscosity. The linear relationship is much stronger at a low viscosity (Pearson's coefficient of 0.51 for $\tilde{\alpha} = -2.2$ and 0.55 for $\tilde{\alpha} = -3.2$); at $\tilde{\alpha} = -1.8$, the correlation is weaker.

<figure><img src="image_18.png" /><figcaption>**Fig. 11** Proportion of grains in structures for A $\tilde{\alpha} = -1.8$, B $\tilde{\alpha} = -2.2$, C $\tilde{\alpha} = -3.2$</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 12 Probability density function of structures for various cases of normalized angular velocity of structures $\tilde{\omega}$ and diameter of structures $D$</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. 13** A $\tilde{\omega}$ as a function of normalized rigidity $\tilde{E}$; B $\tilde{\omega}$ as a function of solid fraction $S_f$. The size of the dots is relative to the proportion of grains in agglomerates in the simulation (in %)</figcaption></figure>

Another remarkable correlation, visible in Fig. 13B, is that between the average angular velocity of the structures $\tilde{\omega}$ and the solid fraction $S_f$. The solid fraction is the ratio of the surface occupied by the particles in the simulation to the surface of the space that contains them. A solid fraction close to 1 indicates that the grains are very tightly packed and fill all the pore space, while a value close to 0 indicates the opposite.

## 3.4 Diameter of structures

One way to quantify the size of the structures is to use a cylindrical representation. For this purpose, each structure is associated with a lifetime $L$ and a diameter $D$. This diameter $D$ is based on the average area of the structure in each time step along the lifetime. It is thus possible to approximately represent each structure as a cylinder whose height is the lifetime and whose diameter characterizes the average spatial extent over time. An example of cylindrical representation can be found in Fig. 14.

In the same way as for the normalized angular velocity, statistical data of the structure diameters can be extracted from the studied simulations. The probability density function and the cumulative distribution function of the diameter of structures can be seen in Fig. 12, right column. In path B→C→D, case C corresponds to a minimum size of agglomerates. In the same way, as for the angular velocities,

<figure><img src="image_21.png" /><figcaption>**Fig. 14** Cylindrical representation of detected structures for case F (cf. Fig. 9). The structures rotating in the positive direction are in red, and in blue for those rotating in the negative direction. The intensity of the color depends on the angular velocity</figcaption></figure>

case C corresponds to an intermediate step in the E→C→F paths.

An interesting relationship can be observed between the cohesion and the characteristic (i.e. the mode of the PDF)
<figure><img src="image_22.png" /><figcaption>Fig. 15 Average diameter of the structures as a function of the cohesion parameter</figcaption></figure>

diameter of the structures (cf. Fig. 15). The Spearman's correlation coefficient (relevant for non-linear correlations) is quite high, over 0.67. This indicates that the larger the cohesion, the larger the diameter of the structures, which agrees with the observations that can be made in Fig. 9. The relation between these quantities is non-linear, with a sharp increase of the characteristic diameter of the agglomerates when the normalized cohesion goes beyond 0.5.

## 3.5 Coefficient of friction

A correlation can be found between the normalized angular velocity structures $\tilde{\omega}$, and the coefficient of friction and can be visualized in Fig. 16. The Pearson correlation coefficient between these two data is -0.70. When a large friction coefficient is observed, the average angular velocity of the field is very low. In contrast, the friction coefficient is found to be very low when the normalized angular velocity gets closer to 1, i.e. when the agglomerates rotate following the natural shear rate of the flow.

As described in [43], based on energetic arguments, coefficient of friction $\mu$ can be expressed as a sum of two partial friction coefficients: a “surface-related” coefficient of friction $\mu_S$, and a “bulk-related” coefficient of friction $\mu_B$. The former describes the part of the energy dissipated during the creation and destruction of surfaces in the granular sample during shear, and the latter the energy dissipated by inelastic deformations within the deformable particles.

<figure><img src="image_23.png" /><figcaption>Fig. 16 Coefficient of friction as a function of mean normalized angular velocity</figcaption></figure>

$$ \mu = \mu_S + \mu_B \qquad (8) $$

These components are particularly complex and difficult to correlate with other parameters. However, two significant correlations can be observed. In Fig. 17A, a link can be made between the normalized angular velocity of structures $\tilde{\omega}$ (cf. Sect. 3.2) with the bulk-related coefficient of friction, with a Pearson coefficient of -0.58. As the structures rotate more, there is less energy loss through dissipation in the bodies.

In Fig. 17B, a clear link can be made between the size of the structures (cf. Sect. 3.4) and the bulk-related coefficient of friction (Spearman coefficient of -0.83). This relationship is monotonic but non-linear. The higher the coefficient of friction associated with the interparticle energy dissipations, the smaller structures there are. In contrast, the appearance of large agglomerates (i.e. with $D > 4$) correlates well with a large drop of the surface-related energy dissipation in the system.

## 4 Discussion and conclusion

Based on the quantitative results presented in the previous section, we can draw a general picture of the causes and consequences of the agglomeration phenomena in soft cohesive granular flows. It appears that the existence of persistent rigid-like agglomerates in this class of flows is promoted by the interparticle cohesion and to a lesser extent
<figure><img src="image_24.png" /><figcaption>**Fig. 17** A Bulk-related coefficient of friction as a function of mean normalized angular velocity; B Surface-related coefficient of friction in function of mean diameter</figcaption></figure>

by the viscous response to shear deformation (cf. Fig. 11). Agglomeration is indeed negligible when the quantities $\tilde{c}$ and $\tilde{\alpha}$ are low. The size of the agglomerates is mostly related to cohesion (cf. Fig. 15), which is simply explained by the fact that more cohesive particles tend to pull more surrounding particles in the same motion. The normalized cohesion value $\tilde{c} = 0.5$ seems to be a threshold above which very large structures can develop. Viscosity also plays a role in the creation of large structures (cf. Fig. 11).

The stiffness and the viscosity parameters play an important role in the angular velocity of the structures (cf. Fig. 13A). The softer and less viscous the grains are, the more they tend to form structures with an angular velocity corresponding to the natural shear rate of the interface (although a few outlying cases are noticed). When the stiffness is important, the structures are more reluctant to follow the natural shear-related rotation. The solid fraction also plays a role in the angular velocity of the structures (cf. Fig. 13B). When the solid fraction is high, the grains are closer to each other, their contacting areas are larger, and their motions are more coordinated, which seems to increase the angular velocity of agglomerates. Since a high stiffness and a low viscosity (to a lesser extent) decrease the solid fraction, we can speculate that this decrease in solid fraction provides less incentive for structures to rotate. Hence, a simplistic view of this class of granular flows could be that cohesion and viscosity favor the formation and the growth of agglomerates, and that stiffness and viscosity restrict

their natural rotation within the sheared flow. The proposed results hence provide a contribution to the growing topic of the rheology of cohesive granular flows [52, 53], with the added complexity of a cohesive strength increasing with the contact area between the grains (and hence with their deformability and the applied stress).

A link can be made with the surface-related coefficient of friction (cf. Fig. 17B). At a diameter greater than 4, the value of the surface-related coefficient of friction drops and becomes very low. This could be due to a geometrical effect linked to a specific surface limit beyond which the behavior changes radically. This value is approximate and could be dependent on the parameters chosen in the simulations. It is at the same value that a behavior change is observed between the cohesion and the average diameter of structures. This decrease of the surface-related dissipation can be explained by the fact that the agglomerates become larger. The specific surface area thus decreases in the sheared sample and the contact area mobilized between the agglomerates diminishes. The energy dissipation therefore no longer occurs at the interfaces between the agglomerates, but within them, via the inelastic component of their deformation. Likewise, as structures rotate faster (i.e. closer and closer to the natural angular velocity associated with rolling in the sheared flow) they tend to endure less deformation, resulting in a decrease in the bulk-related coefficient of friction (cf. Fig. 17A).

Finally, the friction coefficient strongly depends on the angular velocity of the agglomerates (cf. Fig. 16). It can
be assumed that the agglomerates roll in the contact, and
decrease the coefficient of friction, like balls in bearings.
This link between high cohesion, the generation of roll-
ing agglomerates, and a limitation of the resulting fric-
tion coefficient at the scale of the interface is consistent
with the hypothesis reported in a previous study [46], but
is now supported by more quantitative data. The accom-
modation of the relative velocity of the two walls is more
likely to occur via rolling at low stiffness and high cohe-
sion, while it is more likely to occur via particle deforma-
tion when the cohesion is low.

The methodological framework proposed in this paper for the detection and quantification of agglomerated structures is still imperfect. A blind spot exists for structures with zero angular velocities, although it seems minimal. The presence of large pore spaces is also ignored by the interpolation algorithm, and this could be corrected in future work. The main limitation, however, seems to be the rather short duration of the simulations, which are imposed by the large computational cost of the multibody meshfree approach but limit the statistical accuracy of the detection and characterization of the agglomerates. Much longer simulations could be made possible in future studies by the use of the Soft Discrete Element Method [54]. The thickness of the sheared layer, which provides an upper bound for the agglomerates diameters, might also play an important role and should be investigated in future studies.

From a more applicative perspective, the tool developed
here will be applied to fretting simulations in a tribological
context. By applying this approach to realistic mechanical
contacts, we expect to establish a link between the third body
rheological properties, the formation and characteristics of
agglomerates, the stress fluctuation patterns undergone by
the surfaces bounding the third body flow, and finally their
damage and wear.

**Acknowledgements** The authors gratefully acknowledge Safran Aircraft Engines and the French National Research and Technology Agency (ANRT) for financially supporting this research project. The authors acknowledge that this study contains original material. Its publication has been approved tacitly by the responsible authorities at the institutes where the work has been carried out.

Declarations

**Conflict of interest** The authors declare that they have no conflict of interest.

References

1. Dowson, D.: History of Tribology, 2nd edn. Wiley, Hoboken (1998)

2. Descartes, S., Courtieux, L., Berthier, Y., Peditto, F.: Tribological study of oral care silica. Tribol. Int. **82**, 551–560 (2015). https://doi.org/10.1016/j.triboint.2014.02.023

3. Godet, M.: The third-body approach: a mechanical view of wear. Wear **100**, 437–452 (1984). https://doi.org/10.1016/0043-1648(84)90025-5

4. Berthier, Y., Vincent, L., Godet, M.: Velocity accommodation in fretting. Wear **125**, 25–38 (1988). https://doi.org/10.1016/0043-1648(88)90191-3

5. Berthier, Y., Vincent, L., Godet, M.: Fretting fatigue and fretting wear. Tribol. Int. **22**, 235–242 (1989). https://doi.org/10.1016/0301-679X(89)90081-9

6. Kounoudji, K.A., Renouf, M., Mollon, G., Berthier, Y.: Role of third body on bolted joints’ self-loosening. Tribol. Lett. **61**, 25 (2016). https://doi.org/10.1007/s11249-016-0640-8

7. Vingsbo, O., Söderberg, S.: On fretting maps. Wear **126**, 131–147 (1988). https://doi.org/10.1016/0043-1648(88)90134-2

8. Vincent, L., Berthier, Y., Dubourg, M.-C., Godet, M.: Mechanics and materials in fretting. Wear **153**, 135–148 (1992). https://doi.org/10.1016/0043-1648(92)90266-B

9. Collins, J.A.: Failure of Materials in Mechanical Design, 2nd edn. Wiley, Honoken (1993)

10. Arnaud, P., Baydoun, S., Fouvry, S.: Modeling adhesive and abrasive wear phenomena in fretting interfaces: a multiphysics approach coupling friction energy, third body and contact oxygenation concepts. Tribol. Int. **161**, 107077 (2021). https://doi.org/10.1016/j.triboint.2021.107077

11. Baydoun, S., Fouvry, S., Descartes, S.: Modeling contact size effect on fretting wear: a combined contact oxygenation-third body approach. Wear (2021). https://doi.org/10.1016/j.wear.2021.204168

12. Wang, D., Dijksman, J.A., Barés, J., Zheng, H.: Strain dependent vorticity in sheared granular media. EPJ Web Conf. **249**, 02010 (2021). https://doi.org/10.1051/epjconf/202124902010

13. Forterre, Y., Pouliquen, O.: Stability analysis of rapid granular chute flows: formation of longitudinal vortices. J. Fluid Mech. **467**, 361–387 (2002). https://doi.org/10.1017/S0022112002001581

14. Abedi, S., Rechenmacher, A.L., Orlando, A.D.: Vortex formation and dissolution in sheared sands. Granul. Matter. **14**, 695–705 (2012). https://doi.org/10.1007/s10035-012-0369-5

15. Rognon, P., Miller, T., Einav, I.: A circulation-based method for detecting vortices in granular materials. Granul. Matter. **17**, 177–188 (2015). https://doi.org/10.1007/s10035-015-0549-1

16. Ding, J., Leen, S.B., McColl, I.R.: The effect of slip regime on fretting wear-induced stress evolution. Int. J. Fatigue **26**, 521–531 (2004). https://doi.org/10.1016/j.ijfatigue.2003.09.001

17. Basseville, S., Cailletaud, G.: An evaluation of the competition between wear and crack initiation in fretting conditions for Ti–6Al–4V alloy. Wear **328–329**, 443–455 (2015). https://doi.org/10.1016/j.wear.2015.03.010

18. Mary, C., Fouvry, S.: Numerical prediction of fretting contact durability using energy wear approach: optimisation of finite-element model. Wear **263**, 444–450 (2007). https://doi.org/10.1016/j.wear.2007.01.116

19. Paulin, C., Fouvry, S., Meunier, C.: Finite element modelling of fretting wear surface evolution: application to a Ti–6Al–4V contact. Wear **264**, 26–36 (2008). https://doi.org/10.1016/j.wear.2007.01.037

20. Ding, J., McColl, I.R., Leen, S.B., Shipway, P.H.: A finite element based approach to simulating the effects of debris on fretting wear. Wear **263**, 481–491 (2007). https://doi.org/10.1016/j.wear.2006.12.056

21. Gallego, L., Nélias, D., Deyber, S.: A fast and efficient contact algorithm for fretting problems applied to fretting modes I, II
and III. Wear **268**, 208–222 (2010). https://doi.org/10.1016/j.wear.2009.07.019

22. Gallego, L., Fulleringer, B., Deyber, S., Nélias, D.: Multiscale computation of fretting wear at the blade/disk interface. Tribol. Int. **43**, 708–718 (2010). https://doi.org/10.1016/j.triboint.2009.10.011

23. Renouf, M., Massi, F., Fillot, N., Saulot, A.: Numerical tribology of a dry contact. Tribol. Int. **44**, 834–844 (2011). https://doi.org/10.1016/j.triboint.2011.02.008

24. Cundall, P.A., Hart, R.D.: Numerical modelling of dicontinua. Eng. Comput. **9**, 101–113 (1992). https://doi.org/10.1108/eb023851

25. Chapteuil, É., Renouf, M., Zeng, C., Berthier, Y.: Influence of copper/graphite properties on the tribological and electrical behavior of copper-graphite third body layer. Lubricants **6**, 109 (2018). https://doi.org/10.3390/lubricants6040109

26. Rivière, J., Renouf, M., Berthier, Y.: Thermo-mechanical investigations of a tribological interface. Tribol. Lett. **58**, 48 (2015). https://doi.org/10.1007/s11249-015-0523-4

27. Champagne, M., Renouf, M., Berthier, Y.: Modeling wear for heterogeneous bi-phasic materials using discrete elements approach. J. Tribol. **136**, 1–11 (2014). https://doi.org/10.1115/1.4026053

28. Renouf, M., Fillot, N.: Coupling electrical and mechanical effects in discrete element simulations. Int. J. Numer. Methods Eng. **74**, 238–254 (2008). https://doi.org/10.1002/nme.2157

29. Descartes, S., Renouf, M., Fillot, N., Gautier, B., Descamps, A., Berthier, Y., Demanche, P.: A new mechanical–electrical approach to the wheel-rail contact. Wear **265**, 1408–1416 (2008). https://doi.org/10.1016/j.wear.2008.02.040

30. Renouf, M., Nhu, V.-H., Saulot, A., Massi, F.: First-body versus third-nody: dialogue between an experiment and a combined discrete and finite element approach. J. Tribol. **136**, 1–9 (2014). https://doi.org/10.1115/1.4026062

31. Li, W., Huang, Y., Fu, B., Cui, Y., Dong, S.: Fretting damage modeling of liner-bearing interaction by combined finite element–discrete element method. Tribol. Int. **61**, 19–31 (2013). https://doi.org/10.1016/j.triboint.2012.11.019

32. Leonard, B.D., Ghosh, A., Sadeghi, F., Shinde, S., Mittelbach, M.: Third body modeling in fretting using the combined finite–discrete element method. Int. J. Solids Struct. **51**, 1375–1389 (2014). https://doi.org/10.1016/j.ijsolstr.2013.12.036

33. Dijksman, J.A., Brodu, N., Behringer, R.P.: Refractive index matched scanning and detection of soft particles. Rev. Sci. Instrum. **88**, 051807 (2017). https://doi.org/10.1063/1.4983047

34. Vu, T.-L., Barés, J.: Soft-grain compression: beyond the jamming point. Phys. Rev. E. **100**, 042907 (2019). https://doi.org/10.1103/PhysRevE.100.042907

35. Cantor, D., Cárdenas-Barrantes, M., Preechawuttipong, I., Renouf, M., Azéma, E.: Compaction model for highly deformable particle assemblies. Phys. Rev. Lett. **124**, 208003 (2020). https://doi.org/10.1103/PhysRevLett.124.208003

36. Nguyen, T.H., Nezamabadi, S., Delenne, J.-Y., Radjai, F.: Compaction of granular materials composed of deformable particles. EPJ Web Conf. **140**, 05013 (2017). https://doi.org/10.1051/epjconf/201714005013

37. Harthong, B., Jérier, J.-F., Richefeu, V., Chareyre, B., Dorémus, P., Imbault, D., Donzé, F.-V.: Contact impingement in packings of elastic–plastic spheres, application to powder compaction. Int. J. Mech. Sci. **61**, 32–43 (2012). https://doi.org/10.1016/j.ijmecsci.2012.04.013

38. Favier de Coulomb, A., Bouzid, M., Claudin, P., Clément, E., Andreotti, B.: Rheology of granular flows across the transition

from soft to rigid particles. Phys. Rev. Fluids. **2**, 102301 (2017). https://doi.org/10.1103/PhysRevFluids.2.102301

39. Doucet, J.-P., Weber, J.: Computer-Aided Molecular Design: Theory and Applications. Academic Press, Cambridge (1996)

40. Aghababaei, R., Warner, D.H., Molinari, J.-F.: Critical length scale controls adhesive wear mechanisms. Nat. Commun. **7**, 11816 (2016). https://doi.org/10.1038/ncomms11816

41. Molinari, J.-F., Aghababaei, R., Brink, T., Frérot, L., Milanese, E.: Adhesive wear mechanisms uncovered by atomistic simulations. Friction. **6**, 245–259 (2018). https://doi.org/10.1007/s40544-018-0234-6

42. Robbins, M.O., Müser, M.H.: Computer simulations of friction, lubrication, and wear. In: Modern Tribology Handbook, Two Volume Set, pp. 717–765 (2000). https://doi.org/10.1201/9780849377877.ch20

43. Mollon, G.: Solid flow regimes within dry sliding contacts. Tribol. Lett. **67**, 120 (2019). https://doi.org/10.1007/s11249-019-1233-0

44. Mollon, G.: A multibody meshfree strategy for the simulation of highly deformable granular materials. Int. J. Numer. Methods Eng. **108**, 1477–1497 (2016). https://doi.org/10.1002/nme.5258

45. Mollon, G.: A unified numerical framework for rigid and compliant granular materials. Comput. Part. Mech. **5**, 517–527 (2018). https://doi.org/10.1007/s40571-018-0187-6

46. Zhang, Y., Mollon, G., Descartes, S.: Significance of third body rheology in friction at a dry sliding interface observed by a multi-body meshfree model: influence of cohesion between particles. Tribol. Int. **145**, 106188 (2020). https://doi.org/10.1016/j.triboint.2020.106188

47. Quacquarelli, A., Mollon, G., Commeau, T., Fillot, N.: A dual numerical-experimental approach for modeling wear of Diamond Impregnated Tools. Wear **478–479**, 203763 (2021). https://doi.org/10.1016/j.wear.2021.203763

48. Casas, N., Mollon, G., Daouadji, A.: DEM analyses of cemented granular fault gouges at the onset of seismic sliding: peak strength, development of shear zones and kinematics. Earth Sp. Sci. Open Arch. (2021). https://doi.org/10.1002/essoar.105071281.1

49. Jaza, R., Mollon, G., Descartes, S., Paquet, A., Berthier, Y.: Lessons learned using machine learning to link third body particles morphology to interface rheology. Tribol. Int. **153**, 106630 (2021). https://doi.org/10.1016/j.triboint.2020.106630

50. Graftieaux, L., Michard, M., Grosjean, N.: Combining PIV, POD and vortex identification algorithms for the study of unsteady turbulent swirling flows. Meas. Sci. Technol. **12**, 1422–1429 (2001). https://doi.org/10.1088/0957-0233/12/9/307

51. Berson, A., Michard, M., Blanc-Benon, P.: Vortex identification and tracking in unsteady flows. C. R. Mec. **337**, 61–67 (2009). https://doi.org/10.1016/j.crme.2009.03.006

52. Macaulay, M., Rognon, P.: Viscosity of cohesive granular flows. Soft Matter **17**, 165–173 (2021). https://doi.org/10.1039/D0SM01456G

53. Mandal, S., Nicolas, M., Pouliquen, O.: Rheology of cohesive granular media: shear banding, hysteresis, and nonlocal effects. Phys. Rev. X **11**, 021017 (2021). https://doi.org/10.1103/PhysRevX.11.021017

54. Mollon, G.: The soft discrete element method. Granul. Matter. **24**, 11 (2022). https://doi.org/10.1007/s10035-021-01172-9