# Effect of interlocking on the compressive strength of agglomerates composed of cohesive nonconvex particles

Trieu-Duy Tran, Saeid Nezamabadi, Jean-Philippe Bayle, Lhassan Amarsid, Farhang Radjai

► To cite this version:

Trieu-Duy Tran, Saeid Nezamabadi, Jean-Philippe Bayle, Lhassan Amarsid, Farhang Radjai. Effect of interlocking on the compressive strength of agglomerates composed of cohesive nonconvex particles. Advanced Powder Technology, 2025, 36 (2), pp.104780. 10.1016/j.apt.2025.104780. hal-04916992

HAL Id: hal-04916992

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

Submitted on 31 Jan 2025

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 Attribution 4.0 International License
# Effect of interlocking on the compressive strength of agglomerates composed of cohesive nonconvex particles

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

Trieu-Duy Tran<sup>a,b,e</sup>, Saeid Nezamabadi<sup>a,d</sup>, Jean-Philippe Bayle<sup>b</sup>, Lhassan Amarsid<sup>c</sup>, Farhang Radjai<sup>a,*</sup>

## ARTICLE INFO

**Article history:**

Received 6 September 2024

Received in revised form 11 December 2024

Accepted 5 January 2025

Available online xxxx

**Keywords:**

Granular materials

Cohesive agglomerates

Hexapod

Interlocking

Compressive strength

Geometrical cohesion

## ABSTRACT

Despite the important role of fine particle agglomerates in nature and powder technology, the physical mechanisms underlying their strength are still not well understood. In particular, the effects of particle shape in association with cohesive-frictional interactions between particles remain to be elucidated on a quantitative basis. We use particle dynamics simulations to create agglomerates composed of hexapod-shaped particles and analyze their mechanical behavior under diametral compression to show-case the effect of non-convex particle shape and interlocking on their compressive strength. Two different regimes are identified as a function of hexapod aspect ratio. In the first regime, where the hexapods are featured by their rough surface, the strength of the agglomerate is primarily controlled by cohesion at the contact level. In the second regime, where the hexapods are characterized by their long arms that can interlock, the tensile strength at the local level is scaled up by orders of magnitude to yield a high compressive strength at the scale of the agglomerate. We demonstrate that this amplifying effect of interlocking is enhanced by friction coefficient between hexapods, which hinders their disentanglement under the action of the external load.

© 2024 Published by Elsevier B.V. on behalf of The Society of Powder Technology Japan. All rights are reserved, including those for text and data mining, AI training, and similar technologies.

## 1. Introduction

Many natural materials and industrial products are agglomerates of fine particles bonded together by the action of surface forces or a binding phase [1,2]. Well-known examples are metal powder compacts [3], Li-ion battery electrodes [4], pharmaceuticals [5], nuclear fuel pellets [6,7], fertilizers [8], foodstuff [9], and iron ore granules [10]. A key property of agglomerates for most applications is their compressive strength, which depends on the interactions of the primary particles and microstructural features such as porosity, disorder, and distribution of the binding phase. Irrespective of the chemical and physical origins of adhesion in agglomerates, dimensional analysis implies that compressive strength is proportional to the adhesion force between constituents [11,12]. Simple theoretical expressions account also for the connectivity of the bonding network and effects of particle size

polydispersity [12]. There has been, however, amazingly little work on the role of particle shape with respect to the microstructure and strength of agglomerates. The combined effects of non-spherical particle shape and adhesion force on the compressive strength of agglomerates represent therefore an issue of both fundamental nature and practical importance for cohesive granular materials.

In cohesionless granular materials, aspherical particle shape can enhance shear strength through geometrical constraints such as rotation hinderance by face-face contacts between polyhedra and interlocking between nonconvex particles. There are several examples revealing nonlinear and sometimes counterintuitive effects of particle shape. For example, the non-monotonic variation of packing fraction with the level of shape asphericity (any deviation from spherical shape) leads to lower packing fraction but higher frictional strength at large level of asphericity [13–15]. Simulations have also shown that the influence of interparticle friction on the shear strength is enhanced by aspherical shape [16–18]. In the case of cohesive agglomerates, one therefore expects similar effects of even higher complexity due to the lar-

<footnote><sup>a</sup> LMGC, CNRS, University of Montpellier, Montpellier F-34090, France</footnote>

<footnote><sup>b</sup> CEA, ISEC, DMRC, University of Montpellier, Marcoule F-30207, France</footnote>

<footnote><sup>c</sup> CEA, DES, IRESNE, DEC, Cadarache F-13108, France</footnote>

<footnote><sup>d</sup> IATE, CIRAD, INRAE, Montpellier SupAgro, University of Montpellier, Montpellier F-34060, France</footnote>

<footnote><sup>e</sup> Center for Environmental Intelligence, VinUniversity, Hanoi V-100000, Viet Nam</footnote>

<footnote>* Corresponding author.
E-mail address: franck.radjai@umontpellier.fr (F. Radjai).</footnote>
ger variability of their microstructure and the fact that they may fail in tension or shear. These remarks raise at least three general questions: (1) To what extent is the effect of adhesion between primary particles amplified by particle shape asphericity? (2) What is the role of interparticle friction? (3) How does the agglomerate fail?

In this paper, we use particle dynamics simulations to investigate the compressive strength of agglomerates composed of hexapods. Hexapod is a nonconvex particle shape composed of three mutually orthogonal rounded-cap cylinders. It can represent an ideal example of dendritic crystallites or similar nonconvex shapes that can interlock geometrically [19–25]. Interlocking or entanglement between particles arises as a result of contacts between the recesses and protrusions of the particles. To quantify the effect of interlocking, we consider hexapods of increasing aspect ratio, defined as the thickness-to-length ratio of the six arms of the hexapods. With increasing aspect ratio, particle shape varies from sphere to overlapped multi-sphere shape, and then towards thin hexapods resembling a 3D cross. Hence, the level of interlocking increases with aspect ratio.

Agglomerates are built by radial compaction with zero friction between hexapods to obtain dense nearly spherical agglomerates. Then, adhesion and friction forces are activated at all contact points inside the agglomerate and its strength is measured by simulating diametral (Brazilian) compression of the agglomerates between two platens [26,27]. Diametral compression is commonly used in powder technology and rock mechanics to characterize the strength of agglomerates and aggregates [28–31,27,32,33]. The average stress obtained by dividing the force at failure by the largest section of the agglomerate represents its overall compressive strength. Its value reflects the nucleation and propagation of cracks at locations where the shear stress or tensile stress reach a yield threshold. For this reason, the strength obtained by diametral compression reflects not only the adhesion force but also the friction coefficient between primary particles and the microstructure of the agglomerate.

We perform a detailed parametric analysis to quantify the effects of both hexapod aspect ratio and friction coefficient on the compressive stress and bond network. As we shall see, particle shape has a dramatic effect on the agglomerate strength when the aspect ratio of hexapods is large enough to allow for effective interlocking between hexapods. We also show that the friction coefficient between hexapods plays a crucial role for this amplifying effect of interlocking. A detailed analysis of the microstructure at the end of the radial compaction was recently published by the same authors for the same range of aspect ratios as in this paper [25]. The microstructures of the aggregates are very similar at the end of radial compaction and after the removal of the confining pressure and the addition of cohesive forces between particles. Therefore, this paper will be mainly focused on the mechanical behavior of cohesive agglomerates under compression while useful microstructural aspects for the physical interpretation of the observed behavior will also be extensively covered.

In the following, we first introduce in Section 2 the numerical procedures including the numerical method, definition of particle shape, the method used to create aggregates, and the parameters that we applied to diametral compression test. In Section 3, we characterize salient features of the microstructure. In Section 4, we investigate the effects of interlocking on the compressive strength of the aggregates. Then, in Section 5, we focus on the effect of friction coefficient on the scale-up from adhesion between particles to the overall strength of the agglomerate and its micromechanical origins. We conclude in Section 6 with a brief summary of the most important results and several perspectives of this work.

## 2. Model description and numerical procedures

### 2.1. Simulation method

The simulations were carried out by means of an in-house code named Rockable ([34,35]), which is based on Discrete Element Method (DEM) ([36–38]). Particles of arbitrary shape can be built and simulated by means of this code, which efficiently handles frictional contacts, adhesion laws, various time-stepping schemes, and boundary conditions. The basic elements are points, lines, and planes, which are transformed by Minkowsky sum with a sphere of radius $h$ ([39]). As a result, the spheres, cylinders and plates are treated as sphero-points, sphero-segments, and sphero-planes, respectively. All other particle shapes are built by assembling these basic sub-elements.

Contact detection is based on the search of neighbors and overlaps between sub-elements. The contact points are determined with their overlaps and contact normal directions. The total interaction force $\vec{f}$ at a contact point between two particles is the sum of normal and tangential components $f_n$ and $\vec{f}_t$, respectively:

$$ \vec{f} = f_n \vec{n} + \vec{f}_t, \qquad (1) $$

where $\vec{n}$ is the contact normal. The interaction force is a function of the overlap $\delta_n$ and cumulative tangential displacement $\vec{\delta}_t$.

The total normal force $f_n$ is the sum of a repulsive elastic force $f_n^e$ and the adhesion force $-f_c$:

$$ f_n = f_n^e - f_c \qquad (2) $$

The adhesion force $-f_c$ is equivalent to the pull-off force required to break a cohesive bond between two particles. In our simulations we used a constant value of $f_c$ irrespective of local curvature. Indeed, depending on the nature of cohesion, $f_c$ may depend on the surface tension, particle size, local curvature at the contact points or the amount of the binding material [40]. We neglect this variability from one contact to another since the focus of this paper is on the effects of particle shape and microstructure. Furthermore, and in order to focus on the overall effects of the adhesion force, the value of $f_c$ is defined from a reference cohesive stress $\sigma_c$ as discussed in Section 2.1.

We used a linear elastic contact law, in which $f_n^e = -k_n \delta_n$, where $k_n$ is the normal stiffness. The normal force law can be formulated as follows:

$$ f_n = \begin{cases} -k_n \delta_n - f_c & \text{for } \delta_n < 0, \\ 0 & \text{for } \delta_n \ge 0, \end{cases} \qquad (3) $$

The contact is in equilibrium when the adhesion force is balanced by the elastic repulsive force so that $f_n = 0$ and thus $\delta_c = -f_c/k_n$. When $\delta_n$ decreases in absolute value, a tensile force is mobilized ($f_n < 0$) and increases in absolute value up to $-f_n = f_c$ at $\delta = 0$, where the cohesive bond fails. In this work, it is assumed that cohesive bonds are fully irreversible. Hence, the failed bonds transform into cohesionless contacts. A graphical representation of the normal force law is shown in Fig. 1(a).

The tangential force is governed by a linear elastic law together with a Coulomb friction criterion:

$$ \vec{f}_t = \begin{cases} -k_t \vec{\delta}_t & \text{for } |\vec{f}_t| \le \mu f_n^e, \\ -\mu f_n^e \frac{\vec{\delta}_t}{|\vec{\delta}_t|} & \text{otherwise,} \end{cases} \qquad (4) $$

where $\mu$ is interparticle friction coefficient and $k_t$ is tangential stiffness. Note that, for a contact in equilibrium ($f_n = 0$), the friction force threshold is given by $-k_t \delta_c = k_t f_c / k_n$. Compared to the Coulomb criterion $|\vec{f}_t| \le \mu f_n$ for cohesionless contacts, the Coulomb
<figure><img src="image_3.png" /><figcaption>Fig. 1. Graphical representation of the normal force law (a) and tangential elasto-frictional force law (b). Dashed line in (a) represents the force law for a cohesionless contact.</figcaption></figure>

cone for cohesive contacts is shifted to account for the adhesion force added to the normal force, as shown in Fig. 1(b) ([37,38]).

In addition to frictional sliding and bond breaking, which are natural mechanisms of dissipation, we assume that the collisions are inelastic with a restitution coefficient $e_n < 1$. In cohesionless contacts, the value of $e_n$ is controlled by adding a viscous normal force $f_n^v$ to the contact force:

$$f_n^v = -2\alpha_n \sqrt{k_n m} \delta_n, \quad (5)$$

where $m$ is reduced particle mass. With this expression of the damping force, the coefficient $\alpha_n$ is given by ([41–43]):

$$\alpha_n = - \frac{\ln(e_n)}{\sqrt{\ln^2(e_n) + \pi^2}}. \quad (6)$$

This expression holds in the absence of adhesion force. The effective restitution coefficient for cohesive contacts is lower due to the presence of cohesive bond. It can be calculated as a function of $\alpha_n$ and $f_c$ ([44,45]). In all simulations, we set $e_n = 10^{-5}$, corresponding to $\alpha_n \approx 0.96$. Energy dissipation in dense granular systems is a multi-contact process involving elastic wave propagation across the contact network. While high-frequency modes are efficiently damped by viscous terms, low-frequency vibration modes are damped slowly. This effect can, however, be mitigated by setting the restitution coefficient to small values [46].

## 2.2. Particle shape

We use symmetric hexapods, which are composed of three mutually orthogonal cylinders with rounded caps. Hexapod can also be described as the Minkowsky sum of a sphere of diameter $h$ with three mutually orthogonal segments of the same center and length $L$ ([39]). Each arm has a diameter $h$ and total length $d = L + h$ ([23]). The aspect ratio $\alpha$ of the hexapod is the same as the aspect ratio of each arm:

$$\alpha = \frac{L + h}{h}. \quad (7)$$

Fig. 2 shows the shapes of hexapods with increasing $\alpha$. We kept $L$ at a constant value and reduced $h$ to obtain different values of $\alpha$. The case $\alpha = 1$ corresponds to a sphere defined either by setting $L = 0$ for arbitrary value of $h$ or by setting $h$ to infinity for arbitrary value of $L$. It is noteworthy that for $\alpha \le 2$ the cylindrical shapes of the three arms are ‘hidden’ behind their semi-spherical caps. Only when $\alpha > 2$ the cylindrical arms emerge. At large values of $\alpha$, the hexapod approaches the limit of three mutually orthogonal segments.

Hexapod shape may also be described through a parameter $\eta$ which was defined in ([15]) as $\eta = \Delta r/r$, where $\Delta r$ is the difference

between the radius $r$ of the smallest circumscribing sphere and the radius $r - \Delta r$ of the largest inscribed sphere. This definition for hexapods implies

$$\eta = \frac{L}{L + h} = 1 - \frac{1}{\alpha}. \quad (8)$$

This parameter varies from 0 for spheres to 1 for infinitely thin hexapods. Hence, it represents the degree of deviation from spherical shape. In the case of hexapods, this parameter corresponds to the degree of non-convexity. For $\alpha \le 2$ ($\eta \le 0.5$), $\eta$ plays the role of a roughness parameter whereas for $\alpha > 2$ ($\eta > 0.5$), this parameter represents more precisely a shape parameter.

## 2.3. Simulation parameters

We use the hexapod length $L$ as the reference length and we express all lengths as a function of $L$ and the aspect ratio $\alpha$. In particular, we have

$$h = \frac{L}{\alpha - 1}, \quad (9)$$

and

$$d = L + h = L \frac{\alpha}{\alpha - 1} = \frac{L}{\eta}. \quad (10)$$

In the same way, we define a cohesive stress $\sigma_c$, which can be set to an arbitrary value without changing the simulation results or loosening generality provided all stresses are normalized by $\sigma_c$. Accordingly, we set the adhesion force between hexapods to

$$f_c = \sigma_c L^2. \quad (11)$$

In the same way, the normal stiffness $k_n$ in the simulations was set to

$$k_n = 4 \times 10^4 L \sigma_c. \quad (12)$$

The prefactor $4 \times 10^4$ was chosen to make the particles stiff enough to keep the maximum overlap $\delta_{max}$ between particles to a low value (below $10^{-5}L$).

In our simulations, we set $\sigma_c = 0.5$ MPa and $L = 10^{-3}$ m. However, as stated previously, the results expressed in terms of the normalized stresses, forces and lengths do not depend on this choice. We also set $k_t = 0.8k_n$. Five different values of $\mu$ from 0 to 0.8 were used in the simulations. The values and ranges of all numerical and physical parameters used in this paper are given in Table 1.

## 2.4. Radial compaction

To create isotropic agglomerates of hexapods, we implemented the radial compaction process. The particles are placed at random locations inside a spherical ball of radius $R$. The latter is free to change according to the equation of motion:

$$\sum_i f_i + p(4\pi R^2) = m \frac{d^2 R}{dt^2}, \quad (13)$$

where $p$ is the confining pressure, $f_i$ is the normal force exerted by particle $i$ on the ball, and $m$ is the fictitious mass of the ball. Before the application of radial compaction, a small randomly-oriented velocity is added to each particle. The particles move around and a nearly uniform distribution of particle centers is reached inside the ball. The ball is then allowed to contract by applying the pressure $p$ on the ball. The hexapods inside the ball are swept by the walls and pushed inward. The contraction slows down due to energy dissipation by inelastic collisions and friction. All particle
<figure><img src="image_4.png" /><figcaption>Fig. 2. Hexapods of increasing aspect ratio $\alpha$.</figcaption></figure>

<table><caption>**Table 1**

Numerical and physical parameters.</caption><thead><tr><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>Number of particles, N<sub>p</sub></td><td>8385</td></tr><tr><td>Particle aspect ratio, α</td><td>[1–9]</td></tr><tr><td>Reference cohesive stress, σ<sub>c</sub></td><td>5 MPa</td></tr><tr><td>Reference hexapod length, L</td><td>10<sup>-3</sup> m</td></tr><tr><td>Adhesion force, f<sub>c</sub></td><td>σ<sub>c</sub>L<sup>2</sup></td></tr><tr><td>Particle diameter, d</td><td>L<sup>2/3</sup>σ<sub>c</sub></td></tr><tr><td>Particle density, ρ</td><td>2500 kg m<sup>-3</sup></td></tr><tr><td>Normal stiffness, k<sub>n</sub></td><td>4 × 10<sup>4</sup>Lσ<sub>c</sub></td></tr><tr><td>Tangential stiffness, k<sub>t</sub></td><td>0.8 k<sub>n</sub></td></tr><tr><td>Normal damping, x<sub>n</sub></td><td>0.96</td></tr><tr><td>Normal resitution coefficient, e<sub>n</sub></td><td>10<sup>-5</sup></td></tr><tr><td>Friction coefficient, μ</td><td>[0, 0.2, 0.4, 0.6, 0.8]</td></tr></tbody></table>

velocities and contraction rate tend asymptotically to zero and the hexapods become fully jammed in a static configuration with $\sum \vec{f}_i = -p(4\pi R^2)$ in finite time.

Since we are interested in the effects of particle shape on the overall strength of the agglomerates, it is important to avoid variability arising from the preparation protocol. For this reason, during radial compaction the adhesion force and friction coefficient were set to zero. This leads to well-defined isostatic jammed configurations. At low values of confining pressure $p$, the theoretical contact coordination number (average number of contacts per particle) is $Z_c \approx 6$ for spherical particles and $Z_c \approx 12$ for hexapods [25,47–49].

Once the jammed configuration is stable, the spherical ball used for radial compaction is deleted and the adhesion force is activated between particles. To avoid elastic rebound, the adhesion force $f_c$ must be larger than the repulsive force $pd^2$ induced by the compaction pressure. For all values of $\alpha$, this force is large enough to hinder the explosion of the agglomerate despite the removal of the confining ball. The agglomerate is allowed to relax to a stable state under the action of adhesion forces between particles before the application of diametral compression.

We generated hexapods with $\alpha$ varying from 1 to 9 corresponding to the shapes displayed in Fig. 2. The parameter $\eta$ ranges from 0 to 0.9. Radial compaction was used according to the aforementioned procedure to build stable cohesive agglomerates composed of $N_p = 8385$ mono-sized hexapods for each value of $\alpha$. The unconfinement agglomerate is self-equilibrated under the action of the elastic repulsive forces and attractive adhesion forces between particles. Fig. 3(a) shows the initial state of the agglomerate of hexapods with $\alpha = 9$ in which the particles are represented in two different colors depending on the compressive or tensile nature of the force to which they are subjected. Each of these two subsets of particles percolates across the agglomerate and the two subsets define

together a bi-percolating network. A similar bi-percolating structure was previously observed in simulated samples of cohesive spherical particles ([50]). Fig. 3(b) shows the connectivity network defined by the segments joining the centers of touching hexapods. The color depends on whether the total force along the segment between two hexapods is tensile or compressive. Since there can be several contacts between two hexapods, we define the compressive or tensile nature of a link from the sign of the projection of the vector sum of all force vectors between two hexapods on their center-to-center line. Here also we observe a bi-percolating bond network.

## 2.5. Diametral compression

The agglomerates created by radial compaction for given values of aspect ratio $\alpha$, adhesion force $f_c$, and friction coefficient $\mu$, are subjected to diametral compression between two rigid plates. The top plate moves slowly downward while the bottom plate is kept at a fixed position. The top plate speed and time step are set to small values such that the contact deflections $\delta = f_n/k_n$ are well resolved throughout the system. The contact deflections can be compared with the length $L$ of the hexapods. The ratio $\sigma_e = k_n/L$ defines the order of magnitude of the elastic modulus of the agglomerate. Given the values of system parameters in Table 1, we have $\sigma_c/\sigma_e = f_c/k_nL = 2.5 \times 10^{-5}$. This is also the order of magnitude of the relative elastic deflection $\delta/L$ in our simulations.

Fig. 4 displays a sideview snapshot of diametral compression for the agglomerate composed of hexapods of aspect ratio $\alpha = 9$. Initially, no hexapods are in contact with either of the two plates. As the top plate moves downward, the number of contacts between the agglomerate and the two plates grows and the contact forces between hexapods increase along the compression axis. Fig. 5 displays the top and side views of the contact force network shortly after failure. Only the centers of the hexapods are shown and they are interconnected by segments whose thickness is proportional to the total force acting between them and whose color represents either tensile or compressive contacts. We see that strongest forces occur along the central axis and they are compressive whereas tensile contacts are mostly perpendicular to the axis. At the onset of breakage, an increasing number of bonds break along the ortho-radial direction, leading to the opening of a transverse crack, as observed in Fig. 5(b). The broken bonds either disappear from the contact list or turn into cohesionless contacts.

The vertical deformation $\epsilon$ of the agglomerate is given by:

$$ \epsilon = \ln \left( 1 + \frac{s}{D_v} \right) \approx \frac{s}{D}, \qquad (14) $$
<figure><img src="image_5.png" /><figcaption>Fig. 3. (a) Side view of the agglomerate composed of hexapods of aspect ratio 9. The two colors represent hexapods subjected to compressive and tensile pressures. (b) Connectivity network in the same agglomerate defined by segments joining the centers of touching hexapods. The tensile and compressive links are colored in red and white, respectively.</figcaption></figure>

where s is the total downward displacement of the top plate, $D_v$ is the distance between the top and bottom plates, and D is the initial diameter of the agglomerate as obtained by radial compaction and subsequent free relaxation under the action of the internal adhesion forces between the particles. We also measure the vertical force component F between the agglomerate and the top plate by summing the normal forces between primary particles and the plate. Let us note that, as the agglomerate is in static equilibrium, the forces at each horizontal layer of the agglomerate are balanced so that the total vertical force acting between all horizontal layers is equal to F. Let $D_h$ be the diameter of the largest section perpendicular to the compression axis and passing through the center of the agglomerate. By dividing F by the largest sectional area $\pi D_h^2/4$ we obtain the average vertical stress in this section:

<figure><img src="image_6.png" /><figcaption>Fig. 4. Sideview snapshot of an agglomerate composed of hexapods of aspect ratio 9 during diametral compression. The colors encode the horizontal velocity components of the particle centers of mass.</figcaption></figure>

$$ \sigma = \frac{4F}{\pi D_h^2}. \qquad (15) $$

The value of $D_h$ during diametral compression is estimated from the positions of the particles located at the boundary of the central section of the agglomerate. This estimate of $\sigma$ is objective as long as the agglomerate is not damaged or broken. When the cracks appear and propagate across the agglomerate, the sectional area becomes ill-defined. The evolution of $\sigma$ as measured using this procedure is therefore reliable mainly up to the breakage of the agglomerate although the force F is correctly evaluated.

It is worth noting that the stress tensor can also be calculated from the values of contact forces and branch vectors (vectors joining particle centers) ([51,52,50,53]). The vertical component is given by

$$ \sigma_{zz} = \frac{1}{V} \sum_{k=1}^{N_c} f_z^k \ell_z^k = n_c \langle f_z \ell_z \rangle, \qquad (16) $$

where V is the volume of the agglomerate, $N_c$ is the number of contacts, $n_c = N_c/V$ is the number density of contacts, $f_z^k$ and $\ell_z^k$ are the z-components of the contact force vector and branch vector, respectively, at the contact k, including all internal contacts and the contacts with the plates. In contrast to $\sigma$ given by Eq. (15), which corresponds to the average stress component at the center of the sphere, the expression (16) corresponds to the average vertical stress component in the whole volume of the agglomerate and its value is lower by a factor 1/2. In the following, we will use $\sigma$ according to (15) for the sake of consistency with experiments in which the vertical stress component is calculated from Eq. (15). Eq. (16) can also be applied to each individual particle to calculate its average pressure from its contacts with neighboring particles. We will use this equation to evaluate the probability density function of particle pressures.
<figure><img src="image_7.png" /><figcaption>Fig. 5. The side (a) and top (b) views of force chains in the agglomerate composed of hexapods of aspect ratio $\alpha = 9$ and friction coefficient 0.6 under diametral compression after the onset of fracture. Line thickness is proportional to the total force acting between two neighboring hexapods. The white and red segments correspond to tensile and compressive contacts, respectively.</figcaption></figure>

### 3. Agglomerate microstructure

Before considering the mechanical behavior of the agglomerates under diametral compression, it is useful to characterize their microstructure in the unconfined initial state, i.e. after the removal of the ball and relaxation of the agglomerate under the action of adhesion forces with zero friction coefficient. All the features presented in this section are quite similar to those before the removal of the ball, which were characterized in detail in Ref. [25]. As stated previously, in the absence of friction forces, the hexapods within the agglomerate form a random close packing of highest packing fraction $\Phi$. The packing fraction is calculated by dividing the total volume $V_p$ of the particles by the total volume $V$ of the agglomerate. Fig. 6(a) shows $\Phi$ as a function of non-convexity parameter $\eta$. We observe an initial increase of $\Phi$ up to a maximum value 0.68 for $\eta = 0.5$ and a subsequent decrease to values as low as 0.25 at high values of $\eta$.

For $\eta = 0$ (spherical particles), we have $\Phi = \Phi_0 \approx 0.62$ which is very close to the known RCP density of spherical particle packings ($\approx 0.635$) with periodic boundary conditions ([49,54]). The lower packing fraction of the agglomerate, as compared to that with periodic boundary conditions, is due to the higher porosity at its boundary. This is a surface effect or a wall effect (since the agglomerate is prepared by a contracting ball), and the resulting perturbation

<figure><img src="image_8.png" /><figcaption>Fig. 6. (a) Packing Fraction $\Phi$ as a function of aspect ratio $\eta$ in frictionless aggregates. The continuous line represents the fitting form (17). (b) Coordination number $Z$ and contact coordination number $Z_c$ as a function of shape parameter $\eta$.</figcaption></figure>

propagates into the volume over a distance of nearly 10 particle diameters. The packing fraction is 0.635 in the core of the agglomerate, i.e., inside a spherical volume whose radius is about 10 particle diameters smaller than that of the agglomerate, and its measurement requires cutting the particles with the spherical control volume. We did not measure this core value of the packing fraction for other shapes since it requires cutting the hexapods by the control volume. In general, we may consider that the values of $\Phi$ as a function of $\eta$ in Fig. 6(a) are all slightly below their core values.

The non-monotonic variation of $\Phi$ as a function of $\eta$ has been observed for different particle shapes and seems to be a consequence of the competing effects of steric exclusions and decrease of particle volume as $\eta$ increases [55,49,14,56,57]. As shown in Fig. 6, the plot of $\Phi$ as a function of $\eta$ is well fit to a polynomial function of degree 3:

$$ \Phi = \Phi_0 + a\eta^2 - b\eta^3, \quad (17) $$

with $a \approx 0.85$ and $b \approx 1.5$. The positive term $a\eta^2$ accounts for increasing interlocking between hexapods whereas the negative term $-b\eta^3$ represents decrease of particle volume. The peak value occurs at $\eta = 2a/3b \approx 0.38$.

The lowest-order parameter for the description of the connectivity of particles is the coordination number [58]. There can be several contact points between two hexapods. For this reason, we distinguish the coordination number $Z$, defined as the average number of contact neighbors per particle (neighboring particles sharing at least one contact with a given particle), from the contact coordination number $Z_c$, defined as the average number of contacts per particle [25]. We will refer to a pair of particles sharing at least one contact as a 'link'. The number of contacts in a link is the 'vallence number'. Fig. 6(b) shows $Z$ and $Z_c$ as a function of $\eta$. The isostatic state of the agglomerate implies $Z_c = 6$ for spherical particles
and $Z_c = 12$ for hexapods in the limit of perfectly rigid particles [47]. We see that the values of $Z_c$ are indeed very close to their iso-static values except for $\eta \simeq 0.67$ ($\alpha = 3$) for which $Z_c$ is close to 14 for the confining pressure $p$ applied. Interestingly, the coordination number $Z$ has nearly a constant value ($\sim 7$) for $\alpha < 3$ (except for $\eta = 0$ where the particles are spherical and $Z = Z_c$), but increases to values above 8 for larger values of $\alpha$. This increase is a consequence of the ability of thin hexapods to touch their second neighbors (see below).

The higher value of $Z_c$ in the case $\alpha = 3$ reveals the higher sensitivity of the microstructure to the confining pressure $p$ for this aspect ratio. In fact, we find that when $p$ is reduced, $Z_c$ declines towards 12 while $Z$ keeps its value. In other words, the number of contacts of the hexapods declines with decreasing pressure without significant change of the number of neighbors. This higher sensitivity to pressure for $\alpha = 3$ reflects the trend of hexapods of this particular aspect ratio to form parallel ordered local structures of high valence, as discussed in more detail in Ref. [25].

It is also interesting to quantify the average valence number $\bar{\vartheta} = Z_c/Z$ (average number of contacts per link). Fig. 7(a) shows $\bar{\vartheta}$ as a function of $\eta$. The valence number $\bar{\vartheta}$ is, by definition, 1 for spheres and takes a value between 1.7 and 1.8 for thick hexapods. Its value reaches 2 for $\alpha = 3$ and declines to 1.5 for thin hexapods as $\alpha$ increases beyond 3. This latter decrease reflects the increase of coordination number $Z$ due to the increasing number of contacts with second neighbors. Therefore, the same number of contacts (12) is shared with more particles, leading to lower valence number.

The concept of “second neighbor” is intuitive in the case of spherical particles but must be defined for thin hexapods. Let us consider the distance $\ell$ between the centers of two hexapods (branch vector lengths). The values of $\ell$ are broadly distributed in a range varying from a mutual exclusion distance $\ell_{min}$ up to $d = L+h$. For hexapods with cylindrical arms ($\alpha > 2.5$), the mini-

mum distance $\ell_{min}$ between the centers of two hexapods occurs when their 6 arms are parallel and touch each other. In this limit, we have $\ell_{min} = \sqrt{3}h$ and thus $\ell_{min}/d = \sqrt{3}/\alpha$. We define the first neighbors as the hexapods whose centers are located within a distance [$\ell_{min}, 2\ell_{min}$] from the center of a hexapod. The second neighbors are the hexapods whose centers are located beyond the distance $2\ell_{min}$. With this definition, a hexapod can touch its second neighbors. This condition holds if the length $(L+h)/2$ of its arm is longer than $\ell_{min} = \sqrt{3}h$. Therefore, a hexapod may touch its second neighbors if $\alpha > 2\sqrt{3} \simeq 3.46$. This aspect ratio distinguishes ‘thick hexapods’ ($\alpha < 3.46$) from ‘thin hexapods ($\alpha > 3.46$). Only thin hexapods may touch their second neighbors. Fig. 8 displays an example of the first and second neighbors of a hexapod.

A simple way of demonstrating the contacts with second neighbors consists in considering the proportions $P_k$ of links (in the sense defined above) between hexapods with exactly $k$ contacts. Fig. 7(b) shows the distributions $P_k$ for different values of $\eta$. We can readily distinguish the group of thick hexapods from the group of thin hexapods. The proportion $P_1$ of links with only one contact is larger in the second group and increases with $\eta$. The significant increase of $P_1$ for larger values of $\alpha$ clearly shows the creation of contacts between hexapods and their second neighbors. Such contacts occur between the ends of the arms so that they cannot be linked with more than one contact. We also see that, for all values of $\alpha$, there is a significant number of links involving 2 contacts. Their proportion $P_2$ decreases from 0.4 for thickest hexapods down to 0.24 for thinnest hexapods. The same decreasing trend (from 0.2 to 0.1) is observed in the case of $P_3$. This means that a higher number of single contacts with second neighbors (increase of $P_1$) is sanctioned by a lower number of links involving 2 or 3 contacts as $\alpha$ increases. As we shall see, the increase of interlocking and $P_1$ has a dramatic effect on the strength of the agglomerates.

## 4. Effect of particle shape on the compressive strength

We first investigate the effect of particle shape during diametral compression with $\mu = 0.2$. Note that the bonding strength between

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** (a) Average valence number $\bar{\vartheta}$ (number of contacts per pair of touching hexapods) as a function of shape parameter $\eta$. (b) Proportions $P_k$ of valence numbers $k$ for different values of $\alpha$. The inset shows the variation of $P_1$, i.e. number of pairs of particles with exactly one contact, as a function of nonconvexity $\eta$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Three hexapods of aspect ratio 9. Based on the distances between particle centers, hexapod B is a first neighbor of hexapod A whereas hexapod C is a second neighbor of A; see text. In this example, hexapod A touches its second neighbor.</figcaption></figure>
particles is controlled by the cohesive stress $\sigma_c$ and the values of adhesion force are given by $f_c = \sigma_c L^2$. Dimensional analysis implies that the strength of the agglomerate, i.e. the largest value of the vertical stress during compression and, more generally, the level of stress at arbitrary instant of evolution of the agglomerate, are proportional to $\sigma_c$, which represents the only internal stress scale of the system. As we shall see below, the elastic deformations may make the strength of the agglomerate slightly depend also on the contact stiffness.

Fig. 9 shows the vertical stress $\sigma$ normalized by $\sigma_c$ as a function of strain during diametral compression. We see that, for all values of $\alpha$, $\sigma$ increases from zero to a peak value before decreasing and tending to a plateau. The peak is rather better pronounced for thick hexapods ($\alpha \le 3$) than for thin hexapods. In all cases, i.e. for both thick and thin hexapods, the first breaking event occurs at $\varepsilon \approx 0.01$. We refer to the corresponding stress as peak stress denoted by $\sigma_p$. For $\alpha \le 3$, the stress drops after this point whereas for $\alpha = 5$ and 7, the vertical stress remains nearly constant after this point and for $\alpha = 9$, $\sigma$ keeps increasing towards a higher peak, to which we will refer as secondary peak, denoted by $\sigma_p$, in Fig. 9. Considering also the evolution of $\sigma$ in the cases of $\alpha = 1.25$ and $\alpha = 3$, we see that the behavior after the first breakage event depends on the evolution of cracks inside the agglomerate and its contacts with the two plates. In fact, the agglomerate is never fully broken since its breakage is controlled by the slow motion of the top plate. The scenario would have been different if the top plate were controlled by stepwise increase of the stress, in which case the agglomerate would break apart instantly as a result of avalanching energy release as soon as cohesive bonds begin to break.

The above picture of primary versus secondary peak can be further highlighted by following the evolution of the connectivity of hexapods. Fig. 10 shows Z as a function of $\varepsilon$. For all values of $\alpha$, Z begins to decrease from the very beginning. At some value below $\varepsilon \approx 0.02$, Z declines rapidly and then keeps decreasing at a lower

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Vertical stress $\sigma$ normalized by the cohesive stress of particles as a function of vertical strain $\varepsilon$ during diametral compression of agglomerates composed of hexapods of different aspect ratios $\alpha$ and for friction coefficient $\mu = 0.2$. The case $\alpha = 1$ corresponds to spherical particles. The first vertical dashed line indicates the position of the first breakage event. The second dashed line indicates the position of the secondary peak in the case of $\alpha = 9$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Coordination number Z as a function of vertical strain $\varepsilon$ for all compression tests of Fig. 9.</figcaption></figure>

rate. We saw that the initial value of Z is a consequence of radial compaction with zero friction and increases with $\alpha$. But the amount of decline in Z is nearly the same in all cases. It is remarkable that, while the drop in the value of Z clearly preludes the first fracture event, the secondary peak observed in the evolution of $\sigma/\sigma_c$ in Fig. 9 for $\alpha = 9$ has no signature in the evolution of Z. The contact coordination number $Z_c$ shown in Fig. 11 reveals a similar behavior with the difference that, as we saw before, the initial values of $Z_c$ are nearly the same and equal to 12 except for spherical particles ($\alpha = 1$) and $\alpha = 3$.

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Contact coordination number $Z_c$ as a function of strain $\varepsilon$ for all compression tests of Fig. 9.</figcaption></figure>
The above observations suggest that the first peak value $\sigma_p$ of vertical stress occurring at $\varepsilon \approx 0.012$ represents the breakage threshold or the compressive strength of the agglomerates. Fig. 12 displays the evolution of $\sigma_p/\sigma_c$ as a function of $\eta$. This ratio increases from 0.1 at $\eta = 0$ (spheres) to 0.5 at $\eta = 0.67$ ($\alpha = 3$). Then, it continues to increase more rapidly to a value above 2 for $\alpha = 9$. In this way, the agglomerate strength grows by a factor above 20 as the aspect ratio varies by a factor 9. This enhanced strength of agglomerates at large aspect ratio is obviously a consequence of interlocking between thin hexapods and contacts with second neighbors. Interlocking has a strong hindering effect on the relative rotations between two hexapods, but the hexapods can still rotate by sliding along their arms. The contacts with second neighbors considerably reduce these degrees of freedom. These links are also responsible for the increase of Z with $\alpha$ in the group of thin hexapods.

Previously, we argued that, for perfectly rigid particles and for dimensional reasons, the compressive strength is proportional to $\sigma_c$ as the only internal stress scale of the system. To check this point, we repeated the simulations of diametral compression with larger values of $\sigma_c$ for five values of $\alpha$. Fig. 13 shows $\sigma_p/\sigma_c$ as a function of $\sigma_c/\sigma_e$. We see that $\sigma_p/\sigma_c$ is independent of $\sigma_c/\sigma_e$ only in the case of spherical particles. For hexapods, $\sigma_p/\sigma_c$ slightly increases with $\sigma_c/\sigma_e$ despite small values of the latter all the more that $\alpha$ is larger. This increase reflects the influence of the elastic deformations of hexapods on the microstructure. In particular, the contact coordination number $Z_c$ is expected to be sensitive to the increase of elastic overlaps between particles. Fig. 14 displays $Z_c$ both in the initial state and at the peak stress state as a function of $\sigma_c/\sigma_e$ for the simulated values of $\alpha$. We see that $Z_c$ slightly increases with $\sigma_c$ for all values of $\alpha$, but this increase is larger in the cases of $\alpha = 3$ and $\alpha = 9$. As previously stated, the hexapod packings of aspect ratio $\alpha = 3$ are more sensitive to the confining pressure due to their propensity for local ordering. This feature explains both the higher value of $Z_c$ and its faster increase in this case. The larger increase of $Z_c$ in the case $\alpha = 9$ reflects the higher level of interlocking in this case, which favors the closure of small gaps by reducing particle mobility.

Fig. 15 displays the probability density function $P_\sigma$ of the average particle stress $\sigma$ in the initial state and at the stress peak state for all values of $\alpha$. The particle stress is the stress tensor $\sigma_{ij}$ of a particle calculated from the forces exerted by its contact neighbors $c$: $\sigma_{ij} = (1/V_g) \sum_{c=1}^{N_g} f_i^c r_j^c$, where $f_i^c$ is the $i$-component of the force exerted by a neighboring particle at contact $c$, $r_j^c$ is the $j$-component of the vector joining the center of mass of the particle to the contact point $c$ (contact vector), $N_g$ is the number of contacts of the particles, and $V_g$ is its volume. The average particle stress is given by the trace $tr(\sigma_{ij})/3$. In Fig. 15(a) we see that $P_\sigma$ is centered

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Peak stress $\sigma_p$ normalised by the cohesive stress $\sigma_c$ of particles as a function of shape parameter $\eta$ for friction coefficient $\mu = 0.2$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Peak stress $\sigma_p$ normalized by cohesive stress $\sigma_c$ as a function of $\sigma_c/\sigma_e$ for different values of $\alpha$ and $\mu = 0.2$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Contact coordination number $Z_c$ as a function of cohesive stress $\sigma_c/\sigma_e$ (a) in the initial state and (b) in the peak stress state for several aspect ratios $\alpha$.</figcaption></figure>

on and symmetric around $\sigma = 0$ as expected for self-stresses in the absence of external forces. The shape of the distribution is exponential for both compressive and tensile pressures. The map of the percolating networks of particles of tensile and compressive pressures in the initial cohesive state were represented in Fig. 3. It is noteworthy that the largest compressive and tensile pressures in absolute value do not exceed $\sigma_c$. In contrast, at the peak state the probability density $P_\sigma$ shown in Fig. 15(b) is remarkably broader and asymmetric around $\sigma = 0$. The largest particle compressive pressures exceed $25\sigma_c$ for $\alpha = 9$. The high values of pressures specially for thin hexapods are consistent with the amplifying effect of interlocking. The largest forces occur in the central part of the agglomerate along the compression axis.

## 5. Effect of friction coefficient

We evidenced the effect of interlocking with increasing aspect ratio $\alpha$ for $\mu = 0.2$. In granular materials composed of spherical particles, the compressive strength is weakly dependent on friction
<figure><img src="image_17.png" /><figcaption>Fig. 15. Probability distribution function $P_σ$ of average particle stresses $σ$ normalized by the reference cohesive stress $σ_c$ for different values of the particle aspect ratio in the initial state (a) and at the peak stress state (b).</figcaption></figure>

coefficient. This is because the compressive strength essentially depends on the pull-off force between primary particles. In other words, debonding between particles takes place due to separation between particles (opening mode) rather than sliding. In contrast, the relative motion between the centers of two touching hexapods implies the mobilization of friction. For this reason, we expect that friction coefficient must amplify the effects of interlocking.

Fig. 16 shows stress-strain curves with four different values of $μ$ for $α = 3$ and $α = 9$. In the former case, we observe that the failure stress (located around the peak stress) is nearly the same for the three larger friction coefficients but is slightly higher as compared to $μ = 0.2$. In contrast, in the case of $α = 9$, the failure stress, whether the primary failure point at $ε ≃ 0.012$ or the secondary peak at larger strain, increase by orders of magnitude with $μ$. In particular, for $α = 9$, $σ'_{p}/σ_c$ is as high as ≃ 200 for $μ = 0.8$! In this case, and also in the case of $μ = 0.6$, we also observe a rapid drop of stress after failure despite our velocity-controlled setup. This elastic energy release is generally accompanied by the ejection of some hexapods from the surface of the agglomerate.

Fig. 17 shows the peak stress $σ_p$ as a function of shape parameter $η$ for the four values of $μ$. We see that, up to small fluctuations due to the finite size of the agglomerates, $σ_p$ increases almost linearly with $η$ for all values of $μ$ up to $α = 3$. Beyond this point a crossover occurs to a regime of dramatic increase of $σ_p$ with $α$. In the first regime, the peak stress is independent of $μ$ whereas in the second regime $μ$ has an amplifying effect on the peak strength.

To better understand this amplifying effect of interlocking and friction in the second regime, we consider the static equilibrium of the agglomerate as a whole. The compressive stress acting along the vertical axis of the agglomerate is balanced by the mobilization of tensile forces along the ortho-radial directions inside the agglomerate [59]. The peak stress $σ_p$ occurs when the tensile stress at some points inside the agglomerate reaches its maximum value and the first surface of fracture appears. In other words, the peak

<figure><img src="image_18.png" /><figcaption>Fig. 16. Vertical stress normalized by cohesive stress $σ_c$ as a function of axial strain $ε$ for (a) $α = 3$ and (b) $α = 9$, and for four different values of the friction coefficient $μ$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 17. Normalized compressive strength $σ_p/σ_c$ of agglomerates as a function of the shape parameter $η$ for four different values of the friction coefficient $μ$. The vertical dashed line indicates crossover from thick hexapods ($α ≤ 3.46$) to thin hexapods ($α > 3.46$).</figcaption></figure>

state corresponds to incipient failure. Assuming a Mohr-Coulomb behavior, the incipient failure requires that the stress state in a representative volume element around a given point under simple compression (vanishing lateral pressure) satisfies the following relation [60]:

$$ \frac{\sigma_1}{\sigma_t} = \frac{2 \sin \varphi}{1 - \sin \varphi}, \qquad (18) $$

where $σ_1$ is the principal compressive stress component, $σ_t$ is the tensile strength, and $φ$ is the internal (effective) friction coefficient.

Due to the spherical shape of the agglomerate and the action of the compressive stress over a small surface area between the plates and the agglomerate, the stress state is not uniform, as shown in Fig. 5. In particular, the peak compressive stress $σ_p$ is reached while the conditions of failure hold only at a small subset of points inside the agglomerate. We did not calculate the stress distribu-
tions since the total number of hexapods is not sufficiently large to allow for the calculation of statistically meaningful local stress tensors using sufficiently large volume elements placed at different points inside the agglomerate. However, the observation of crack surface development from the vertical axis of the agglomerate towards its periphery suggests that the stress state around the vertical axis is in a state of incipient failure when the peak stress is reached. Therefore, we may replace $\sigma_1$ in Eq. (18) by $\kappa\sigma_p$, where $\kappa$ is a geometrical factor depending only on the shape of the agglomerate. As a result, we have

$$ \frac{\sigma_p}{\sigma_c} = \frac{1}{\kappa} \frac{\sigma_t}{\sigma_c} \frac{2 \sin \varphi}{1 - \sin \varphi}. \qquad (19) $$

This relation involves three different contributions to the compression strength $\sigma_p$ of the agglomerates:

1. The factor $\kappa$ accounts for the inhomogeneous coarse-grained stress distribution inside the agglomerate and may be described as a stress concentration factor reflecting the shape of the agglomerate;

2. The tensile strength $\sigma_t$ of the material is expected to depend on the adhesion force $f_c$ and shape parameter $\alpha$ through the microstructure of the agglomerate;

3. The internal friction angle $\varphi$ is expected to increase with $\alpha$ and the friction coefficient $\mu$ between hexapods.

The values of $\sigma_p/\sigma_c$ obtained from simulations and shown in Fig. 17 can be predicted through Eq. (19) if the values of $\varphi, \sigma_t$ and $\kappa$ are independently measured from simulations or theoretically evaluated for different values of $\alpha$ and $\mu$. For example, using the expression of stress tensor (16), an upper bound of the tensile strength $\sigma_t$ is obtained by replacing the force component $\vec{f} \cdot \vec{n}$ along any direction $\vec{n}$ by $f_c$. Assuming that the structure of the agglomerate is isotropic and the branch vector length $\ell$ is uncorrelated with the value of the force, which is very well verified, we obtain the upper bound $\sigma_t = \frac{n_c}{3} \langle\ell\rangle f_c$. Normalizing by $\sigma_c = f_c/L^2$, yields the following expression [61,60]:

$$ \frac{\sigma_t}{\sigma_c} = \frac{n_c}{3} \langle\ell\rangle L^2. \qquad (20) $$

This expression indicates that the effect of particle shape and other characteristics of a granular packing are mediated only by the values of $n_c$ and $\langle\ell\rangle$. The contact density $n_c$ can be estimated as the average number of contacts per particle $Z_c/2$ divided by the average volume $V_p = \pi\langle\ell\rangle^3/6\Phi$ occupied by each particle. The normalized tensile strength in Eq. (20) can therefore be re-expressed as

$$ \frac{\sigma_t}{\sigma_c} = \frac{Z_c \Phi}{\pi} \left(\frac{L}{\langle\ell\rangle}\right)^2. \qquad (21) $$

It is instructive to evaluate the orders of magnitude of the tensile strength from Eq. (21) in special cases for which the information on the values of parameters is available. For example, for nearly spherical particles ($\alpha = 1.25$) and for $\mu = 0.4$, at failure we have $\Phi \approx 0.6, Z_c \approx 12, \varphi \approx 25^\circ$, and $\langle\ell\rangle \approx L$. Hence, Eq. (21) yields $\sigma_t/\sigma_c \approx 3.3/\kappa$. Equating this value with the measured value in Fig. 17 for $\alpha = 1.25$, implies $\kappa \approx 10$. This is a physically plausible value of $\kappa$ since $\sigma_p$ is defined for the largest cross section of the agglomerate while the volume element is located along the central axis which is the locus of high stress concentration.

For larger values of $\alpha$ in the first regime (thick hexapods), the product $Z_c\Phi$ at failure does not evolve much (see Figs. 6 and 11) and the contribution of $\sigma_t/\sigma_c$ depends mainly on the increase of $(L/\langle\ell\rangle)^2$. This ratio varies by a factor of 2 as $\alpha$ varies from 1.25 to 3.46. Furthermore, the internal friction angle varies from 25° to

$\approx 30^\circ$ in this range, implying an increase of $2 \sin \varphi/(1 - \sin \varphi)$ from 1.45 to 2. Putting together these two contributions due to the evolution of $\sigma_t$ and $\varphi$, and assuming that $\kappa$ does not change, Eq. (19) predicts an increase of $\sigma_p/\sigma_c$ by a factor of 2.75 in the first regime. This is consistent with the values obtained from numerical simulations in Fig. 17 for the first regime. In this range, $\langle\ell\rangle$ is independent of $\mu$ and $\varphi$ varies only slightly with $\mu$. This is the reason why the evolution of $\sigma_p/\sigma_c$ in the first regime is not sensitive to the value of $\mu$.

In the second regime (thin hexapods), the rapid increase of $\sigma_p/\sigma_c$ with both $\alpha$ and $\mu$ cannot be explained from the evolution of $\sigma_t/\sigma_c$. Indeed, as $\alpha$ varies from 3 to 9, both $Z_c$ and $\Phi$ decrease (in the approximate ratios 10/12 and 0.35/0.6, respectively) while $(L/\langle\ell\rangle)^2$ increases only by a factor of 2. Therefore, the enhanced compressive strength in the second regime can only be attributed to the increase of the internal friction angle $\varphi$ with $\alpha$ and $\mu$. In fact, Eq. (19) predicts that $\sigma_p/\sigma_c$ must diverge as $\varphi$ tends to $\pi/2$. This is in qualitative agreement with Fig. 17.

Large values of $\varphi$ (close to 90°) for thin cohesionless hexapods and other nonconvex shapes have been observed in experiments by measuring the angle of repose [21]. Such large angles are sometimes qualified as “geometric cohesion” in analogy with the high angles of repose of wet sand [62]. In our agglomerates, there is no geometric cohesion since the tensile strength $\sigma_t$ vanishes if $f_c$ is set to zero. A geometric cohesion would exist only in the case of ‘hooking’ between particles. However, even for a small nonzero value of $\sigma_t$, the cohesion is enhanced as a result of the large values of $\varphi$ due to interlocking and friction between hexapods. This giant scale-up of the local tensile strength into the macroscopic compressive strength of the agglomerate suggests that small amounts of adhesion may transform interlocked systems of nonconvex particles into strong cohesive solids under compressive load.

## 6. Conclusion

In this paper, we used DEM simulations to analyze the mechanical strength of cohesive agglomerates composed of hexpods of varying aspect ratio. The arms are rounded-cap cylinders and the aspect ratio is defined as the ratio of arm length to arm thickness, which is equal to the diameter of the caps. As a result, low aspect ratios (below 3) correspond to particles similar to 6 overlapped particles along three orthogonal axes (thick hexapods) whereas higher aspect ratios (above 3) have usual hexapod shapes with three mutually perpendicular arms (thin hexapods). We also defined a shape parameter $\eta$, which quantifies the degree of non-convexity of the hexapods. We used radial compaction to build dense agglomerates, which were then subjected to diametral compression between two platens.

The compressive strength of the agglomerates was measured at the peak stress, where the agglomerate begins to break. For thick hexapods the fracture is triggered at a small strain whereas thin hexapods deform up to larger strains before failure despite the irreversible nature of adhesion between hexapods. Our simulation data clearly evidence two regimes with a crossover at aspect ratio 3. In the first regime, the compressive strength (peak strength) increases linearly with the non-convexity parameter $\eta$ and is nearly independent of friction coefficient. In the second regime, a fast diverging increase of compressive strength occurs with increasing $\eta$ and friction coefficient between hexapods.

We discussed the origins of this giant scale-up of the adhesion force between hexapods into a high compressive strength in the light of a simple stress analysis based on the Mohr-Coulomb criterion. This analysis suggests that compressive strength of agglomerates depends on three factors: 1) stress distribution inside the agglomerate, which depends on the shape of the agglomerate, 2)
tensile strength, which can be evaluated from the microstructure, and 3) internal friction angle, which reflects interlocking and friction coefficient between hexapods. From the analytical expression of an upper bound on the tensile strength, we argued that the linear increase of the compressive strength in the first regime is mainly governed by the increase of the tensile strength as a result of the increase of contact density (or decrease of the average distance between particle centers) whereas in the second regime the fast diverging increase of the compressive strength or the Coulomb cohesion mainly reflects the increase of the internal friction angle as a function of nonconvexity and friction coefficient. In other words, in the presence of a minimum amount of adhesion between particles to endow the agglomerate with a minimum level of tensile strength, interlocking behaves as an amplification factor through the increase of the internal friction coefficient.

In pursuing this work, it would be interesting to consider agglomerates composed of hexapods of even larger aspect ratios (beyond 9) to see how the amplifying effect and the overall cohesion diverge as a function of aspect ratio. Another point in this regard is that the irreversible cohesion considered in this work is generally a property of solid bonds, which usually also involve finite rolling resistance. In our simulations, we set the rolling resistance to zero since simple adhesion at multiple contact points between two particles leads to an effective resistance to relative rolling between hexapods. But it is important to properly quantify the influence of rolling resistance at the contact points on the compressive strength.

In this work, we did not consider size effects, which can be characterized by comparing agglomerates of increasing diameter. Finite-size effects are expected to occur more strongly in the case of thin hexapods. We did not neither analyze the breakage of the agglomerates, which leads to the opening of a wedge parallel to the compression axis on one side of the agglomerate. In other words, the agglomerate is never divided into two or more fragments under quasi-static conditions. But we observe dynamic effects at failure such as hexapods detaching from the agglomerate. These phenomena can be studied in more detail by analyzing the evolution of cohesive bonds. Finally, with regard to Eq. (21), its limits and validity need a rigorous and systematic analysis by considering the evolution of the internal friction angle and tensile strength as a function of shape parameter and inter-particle friction coefficient. In the same way, the concentration factor can be evaluated by analyzing the stress fields inside the agglomerate for a larger number of particles than used in the present work.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

The authors thank the CEA SIFCO program for the funding of this project.

## References

[1] C. Nyström, G. Alderborn, M. Duberg, P.-G. Karehill, Bonding surface area and bonding mechanism—two important factors for the understanding of powder comparability, Drug Dev. Ind. Pharm. 19 (1993) 2143–2196.

[2] S.C. Endres, L.C. Ciacchi, L. Mädler, A review of contact force models between nanoparticles in agglomerates, aggregates, and films, J. Aerosol Sci. 153 (2021) 105719.

[3] J.-L. Song, G.-Y. Lee, J.-P. Choi, J.-S. Lee, Compaction behavior of bimodal iron nanopowder agglomerate, Powder Technol. 338 (2018) 333–341.

[4] T. Patey, A. Hintennach, F. La Mantia, P. Novák, Electrode engineering of nanoparticles for lithium-ion batteries—role of dispersion technique, J. Power Sources 189 (2009) 590–593.

[5] S. Palzer, Agglomeration of pharmaceutical, detergent, chemical and food powders—similarities and differences of materials and processes, Powder Technol. 206 (2011) 2–17.

[6] F. La Lumia, Granulation de suspensions concentrées UO2/PuO2: application à l'élaboration de compacts granulaires denses par pressage et à leu caractérisation structurale post frittage, Ph.D. thesis, Limoges, 2019.

[7] F. La Lumia, L. Ramond, C. Pagnoux, P. Coste, F. Lebreton, J.-R. Sevilla, G. Bernard-Granger, Dense and homogeneous mox fuel pellets manufactured using the freeze granulation route, J. Am. Ceram. Soc. 103 (2020) 3020–3029.

[8] S.H. Chien, G. Carmona, L.I. Prochnow, E.R. Austin, Cadmium availability from granulated and bulk-blended phosphate-potassium fertilizers, J. Environ. Qual. 32 (2003) 1911–1914.

[9] W. Stelte, A.R. Sanadi, L. Shang, J.K. Holm, J. Ahrenfeldt, U.B. Henriksen, Recent developments in biomass pelletization—a review, BioResources 7 (2012) 4451–4490.

[10] A. Nyembwe, R. Cromarty, A. Garbers-Craig, Prediction of the granule size distribution of iron ore sinter feeds that contain concentrate and micropellets, Powder Technol. 295 (2016) 7–15.

[11] A. Taboada, N. Estrada, F. Radjaí, Additive decomposition of shear strength in cohesive granular media from grain-scale interactions, PRL 97 (2006) 098302. https://doi.org/10.1103/PhysRevLett.97.098302, URL: https://link.aps.org/doi/10.1103/PhysRevLett.97.098302.

[12] V. Richefeu, M.S. El Youssoufi, F. Radjai, Shear strength properties of wet granular materials, Phys. Rev. E 73 (2006) 051304.

[13] J. Zhao, S. Zhao, S. Luding, The role of particle shape in computational modelling of granular matter, Nat. Rev. Phys. 5 (2023) 505–525.

[14] E. Azéma, F. Radjaí, Stress-strain behavior and geometrical properties of packings of elongated particles, Phys. Rev. E 81 (2010) 051304.

[15] B. Saint-Cyr, K. Szarf, C. Voivret, E. Azéma, V. Richefeu, J.-Y. Delenne, G. Combe, C. Nouguier-Lehon, P. Villard, P. Sornay, et al., Particle shape dependence in 2d granular media, EPL (Europhys. Lett.) 98 (2012) 44008.

[16] S.G. Bardenhagen, J.U. Brackbill, D. Sulsky, Numerical study of stress distribution in sheared granular material in two dimensions, Phys. Rev. E 62 (2000) 3882–3890.

[17] B.B. Dai, J. Yang, C.Y. Zhou, Observed effects of interparticle friction and particle size on shear behavior of granular materials, Int. J. Geomech. 16 (2016) 04015011.

[18] K.M. Salerno, D.S. Bolintineanu, G.S. Grest, J.B. Lechman, S.J. Plimpton, I. Srivastava, L.E. Silbert, Effect of shape and friction on the packing and flow of granular materials, Phys. Rev. E 98 (2018) 050901.

[19] B. Saint-Cyr, F. Radjai, J.-Y. Delenne, P. Sornay, Cohesive granular materials composed of nonconvex particles, Phys. Rev. E 87 (2013) 052207.

[20] H. Zheng, D. Wang, J. Barés, R. Behringer, Jamming by compressing a system of granular crosses, EPJ Web Conf. - Powders and Grains 140 (2017) (2017) 06014.

[21] H.M. Jaeger, Celebrating soft matter's 10th anniversary: Toward jamming by design, Soft matter 11 (2015) 12–27.

[22] N.A. Conzelmann, A. Penn, M.N. Partl, F.J. Clemens, L.D. Poulikakos, C.R. Müller, Link between packing morphology and the distribution of contact forces and stresses in packings of highly nonconvex particles, Phys. Rev. E 102 (2020) 062902.

[23] Y. Zhao, J. Barés, J.E. Socolar, Yielding, rigidity, and tensile stress in sheared columns of hexapod granules, Phys. Rev. E 101 (2020) 062903.

[24] A. Hafez, Q. Liu, T. Finkbeiner, R.A. Alouhali, T.E. Moellendick, J.C. Santamarina, The effect of particle shape on discharge and clogging, Scient. Rep. 11 (2021) 1–11.

[25] T.-D. Tran, S. Nezamabadi, J.-P. Bayle, L. Amarsid, F. Radjai, Contact networks and force transmission in aggregates of hexapod-shaped particles, Soft Matter 20 (2024) 3411–3424. https://doi.org/10.1039/D3SM01762A.

[26] D.-H. Nguyen, E. Azéma, P. Sornay, F. Radjai, Bonded-cell model for particle fracture, Phys. Rev. E 91 (2015) 022203.

[27] M. Dosta, S. Dale, S. Antonyuk, C. Wassgren, S. Heinrich, J.D. Litster, Numerical and experimental analysis of influence of granule microstructure on its compression breakage, Powder Technol. 299 (2016) 87–97.

[28] S.-Q. Yang, P.-F. Yin, Y.-H. Huang, Experiment and discrete element modelling on strength, deformation and failure behaviour of shale under brazilian compression, Rock Mech. Rock Eng. 52 (2019) 4339–4359.

[29] J. Horabik, J. Wiacek, P. Parafiniuk, M. Stasiak, M. Bańda, M. Molenda, Tensile strength of pressure-agglomerated potato starch determined via diametral compression test: Discrete element method simulations and experiments, Biosyst. Eng. 183 (2019) 95–109.

[30] G. Guillonneau, M. Mieszala, J. Wehrs, J. Schwiedrzik, S. Grop, D. Frey, L. Philippe, J.-M. Breguet, J. Michler, J.M. Wheeler, Nanomechanical testing at high strain rates: New instrumentation for nanoindentation and microcompression, Mater. Des. 148 (2018) 39–48.

[31] G.K. Reynolds, J.I. Campbell, R.J. Roberts, A compressibility based model for predicting the tensile strength of directly compressed pharmaceutical powder mixtures, International Journal of Pharmaceutics 531 (2017) 215–224. URL: https://www.sciencedirect.com/science/article/pii/S0378517317307822. doi:10.1016/j.ijpharm.2017.08.075.

[32] S. Antonyuk, S. Heinrich, J. Tomas, N.G. Deen, M.S. Van Buijtenen, J. Kuipers, Energy absorption during compression and impact of dry elastic-plastic spherical granules, Granular Matter 12 (2010) 15–47.
[33] S. Antonyuk, J. Tomas, S. Heinrich, L. Mörl, Breakage behaviour of spherical granulates by compression, Chem. Eng. Sci. 60 (2005) 4031–4044.

[34] V. Richefeu, P. Villard, Modeling Gravity Hazards from Rockfalls to Landslides, Elsevier, 2016.

[35] D.C. Vu, L. Amarsid, J.-Y. Delenne, V. Richefeu, F. Radjai, Macro-elasticity of granular materials composed of polyhedral particles, Granular Matter 26 (2024). https://doi.org/10.1007/s10035-023-01382-3.

[36] P. Cundall, O. Strack, The distinct element method as a tool for research in granular media, Part II, volume 20711, 1979.

[37] T. Pöschel, T. Schwager, Computational Granular Dynamics: Models and Algorithms, Springer Science & Business Media, 2005.

[38] F. Radjai, F. Dubois, Discrete-element modeling of granular materials, Wiley-Iste (2011).

[39] L. Liu, S. Ji, Bond and fracture model in dilated polyhedral dem and its application to simulate breakage of brittle materials, Granular Matter 21 (2019) 41. https://doi.org/10.1007/s10035-019-0896-4.

[40] K.L. Johnson, K. Kendall, A.D. Roberts, Surface energy and the contact of elastic solids, Proc. Royal Soc.A 324 (1971) 301–313.

[41] K. Ismail, W. Stronge, Impact of viscoplastic bodies: dissipation and restitution, J. Appl. Mech. 75 (2008) 061011.

[42] G. Hu, Z. Hu, B. Jian, L. Liu, H. Wan, On the determination of the damping coefficient of non-linear spring-dashpot system to model hertz contact for simulation by discrete element method, in: 2010 WASE International Conference on Information Engineering, volume 3, 2010, pp. 295–298.

[43] P. Müller, D. Krengel, T. Pöschel, Negative coefficient of normal restitution, Phys. Rev. E 85 (2012) 041306.

[44] N.V. Brilliantov, N. Albers, F. Spahn, T. Pöschel, Collision dynamics of granular particles with adhesion, Phys. Rev. E 76 (2007) 051302.

[45] A.K. Dubey, A. Bodrova, S. Puri, N. Brilliantov, Velocity distribution function and effective restitution coefficient for a granular gas of viscoelastic particles, Phys. Rev. E 87 (2013) 062202.

[46] D.C. Vu, L. Amarsid, J.-Y. Delenne, V. Richefeu, F. Radjai, Rheology and scaling behavior of polyhedral particle flows in rotating drums, Powder Technology 434 (2024) 119338. URL: https://www.sciencedirect.com/science/article/pii/S003259102301121X. doi:10.1016/j.powtec.2023.119338.

[47] J.-N. Roux, Geometric origin of mechanical properties of granular materials, Phys. Rev. E 61 (2000) 6802.

[48] H.P. Zhang, H.A. Makse, Jamming transition in emulsions and granular materials, Phys. Rev. E 72 (2005) 011301.

[49] P. Mutabaruka, M. Taiebat, R.J.-M. Pellenq, F. Radjai, Effects of size polydispersity on random close-packed configurations of spherical particles, Phys. Rev. E 100 (2019) 042906.

[50] V. Richefeu, F. Radjai, M.S. El Youssoufi, Stress transmission in wet granular materials, Eur. Phys. J. E 21 (2006) 359–369.

[51] L. Rothenburg, A. Selvadurai, A micromechanical definition of the cauchy stress tensor for particulate media, Mechanics of Structured Media (1981) 469–486.

[52] L. Staron, F. Radjai, J.-P. Vilotte, Multi-scale analysis of the stress state in a granular slope in transition to failure, The, Eur. Phys. J. E 18 (2005) 311–320.

[53] F. Radjai, J.-N. Roux, A. Daouadji, Modeling granular materials: century-long research across scales, J. Eng. Mech. 143 (2017) 04017002.

[54] J.G. Berryman, Random close packing of hard spheres and disks, Phys. Rev. A 27 (1983) 1053–1061. https://doi.org/10.1103/PhysRevA.27.1053, URL: https://link.aps.org/doi/10.1103/PhysRevA.27.1053.

[55] W. Xu, H. Chen, Z. Lv, A 2d elliptical model of random packing for aggregates in concrete, J. Wuhan Univ. Technol.-Mater Sci. Ed. 25 (2010) 717–720.

[56] B. Saint-Cyr, J.-Y. Delenne, C. Voivret, F. Radjai, P. Sornay, Rheology of granular materials composed of nonconvex particles, Phys. Rev. E 84 (2011) 041302.

[57] D.-H. Nguyen, É. Azéma, F. Radjai, P. Sornay, Effect of size polydispersity versus particle shape in dense granular media, Phys. Rev. E 90 (2014) 012202.

[58] R.M. German, Coordination number changes during powder densification, Powder Technol. 253 (2014) 368–376.

[59] T.-T. Vo, P. Mutabaruka, S. Nezamabadi, J.-Y. Delenne, E. Izard, R. Pellenq, F. Radjai, Mechanical strength of wet particle agglomerates, Mech. Res. Commun. 92 (2018) 1–7.

[60] J.-Y. Delenne, F. Soulié, M.S. El Youssoufi, F. Radjai, From liquid to solid bonding in cohesive granular media, Mechanics of Materials 43 (2011) 529–537. URL: https://www.sciencedirect.com/science/article/pii/S0167663611001037. doi:10.1016/j.mechmat.2011.06.008.

[61] F. Radjai, V. Richefeu, Bond anisotropy and cohesion of wet granular materials, Philos. Trans. Royal Soc. A: Math., Phys. Eng. Sci. 367 (2009) 5123–5138.

[62] S.V. Franklin, Extensional rheology of entangled granular materials, Europhys. Lett. 106 (2014) 58004. https://doi.org/10.1209/0295-5075/106/58004.