# Independence of shear strength with particle size dispersity still valid in polyhedral particle assemblies

Emilien Azéma, David Cantor, Itthichai Preechawuttipong

► To cite this version:

Emilien Azéma, David Cantor, Itthichai Preechawuttipong. Independence of shear strength with particle size dispersity still valid in polyhedral particle assemblies. Powders & Grains, Jul 2021, Buenos Aires, Argentina. pp.06009, 10.1051/epjconf/202124906009. hal-03259804

HAL Id: hal-03259804

https://hal.science/hal-03259804

Submitted on 14 Jun 2021

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
# Independence of shear strength with particle size dispersity still valid in polyhedral particle assemblies

Emilien Azéma<fnref n="1" />,<fnref n="2" />,*, David Cantor<fnref n="3" />,**, and Itthichai Preechawuttipong<fnref n="4" />,***

**Abstract.** A very staggering result that has been constantly highlighted in granular media is that the shear strength of granular assemblies is independent of the particle size dispersity. In other words, a packing composed of monodisperse particles has similar strength properties to those of polydisperse systems. This has been shown numerically for the simplified case of disc and polygon assemblies in 2D and spheres in 3D. In this paper, we use three-dimensional contact dynamics simulations to revisit these results for the more complex case of assemblies composed of highly polydisperse rigid polyhedra. Although non-spherical shapes induce more intricate spatial correlations than spherical shapes because of the multiple contact types (i.e., vertex-face, edge-edge, edge-face, face-face), our numerical data provide evidence that the shear strength independence as the particle size dispersity increases still holds up for assemblies of polyhedra. We explain this finding from compensation mechanisms at the micro-scale between geometrical and mechanical anisotropies developed within the assemblies.

## 1 Introduction

Most granular materials in their natural state or in industrial processes are characterized by a very wide range of grain shapes and sizes. Polydispersity can result from various natural processes by which grains are formed (by fragmentation or mineral precipitation). Polydispersity may also be a target in the production, transformation or use of grains. Although the space-filling properties of highly polydisperse systems and their elastic properties under uniaxial compression tests have widely been studied in the past [1-4], it is remarkable to note that very few studies are devoted to their shear behavior.

In this case, the few existing works have constantly highlighted that the shear strength of dry granular assemblies is independent of the particle size polydispersity. Such surprising result has been evidenced in 2D simulations (via discrete-element modeling) using disks [5-7], using regular and irregular polygons [8], and in 3D simulations with spheres [9-11]. The present paper is not intended to review on the practical or industrial implications of phenomenon of the shear strength independence on grain size dispersion. It is instead a generalization of such counterintuitive observation to general angular three-dimensional grains.

In the following, we present the numerical procedures to build and test the shear strength of samples composed

of polyhedra (Sec. 2). In Sec. 3, we focus on the evolution of shear strength and packing fraction with particle size polydispersity. The microstructure, described in terms of force and fabric anisotropies, is discussed in Sec. 4. We conclude, in Sec. 5 with a summary and perspectives of this work.

## 2 Numerical procedures

We use the Contact Dynamics (CD) [12, 13] to simulate assemblies of angular rigid grains in interaction. The CD method is based on implicit time integration of the equations of motion and a nonsmooth formulation of mutual exclusion and dry friction between particles. The contact forces and body velocities are found by means of an iterative nonlinear Gauss-Seidel algorithm. This method requires no elastic repulsive potential and no smoothing of the Coulomb friction law to determine interaction forces. For this reason, simulations can be performed with larger time steps compared to molecular dynamics simulations. We used LMGC90 [14], a multipurpose software developed in Montpellier, capable of modeling a collection of deformable or rigid particles of various shapes by different algorithms.

We consider grains of octahedral shape whose size, defined by the diameter *d* of the circumscribed sphere, is varied in the range of [*d*<sub>min</sub>, *d*<sub>max</sub>] with a uniform distribution of particle volume fractions. We define the size span *S* of

<footnote n="1">¹LMGC, Université de Montpellier, CNRS, Montpellier, France</footnote>

<footnote n="2">²Institut Universitaire de France (IUF), Paris, France</footnote>

<footnote n="3">³Department of Civil, Geological and Mining Engineering, Polytechnique Montréal, Québec, Canada</footnote>

<footnote n="4">⁴Department of Mechanical Engineering, Chiang Mai University, Chiang Mai, Thailand</footnote>

<footnote>*e-mail: emilien.azema@umontpellier.fr</footnote>

<footnote>**e-mail: david.cantor@polymtl.ca</footnote>

<footnote>***e-mail: itthichai_p20@hotmail.com</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1. Screenshots of samples with a grain size dispersion S = 0.8.</figcaption></figure>

the distribution as

$$S = \frac{d_{\max} - d_{\min}}{d_{\max} + d_{\min}}, \qquad (1)$$

In this work, we systematically varied S in the range [0, 0.8] in steps of 0.1. For each value of S, 19000 frictionless grains where geometrically deposited in a cubic box and then isotropically compressed to reach a dense state. Figure 1 shows a screenshots for S = 0.2 and 0.8 after the isotropic compression.

Then, we removed the four lateral walls, we set periodic boundary conditions and the friction coefficient between particles is set to 0.4. The shear tests were performed by fixing the bottom wall and applying on the upper wall a pressure P in the vertical direction and a constant velocity $v_x$ in the x-direction so the inertial number I is below $10^{-3}$ for all the simulations. All the particles in contact either with the upper or the lower wall are bonded and follow the boundary constrains. Gravity was set to 0 to avoid a force gradient within the samples and to prevent particle size segregation. The shearing was undertaken up to reaching the so-called steady flow state in which the volume and shear strength remain constant or present only small fluctuations.

## 3 Shear strength and packing fraction

The shear strength is computed from the well known formula of the granular stress tensor $\sigma$, defined as $\sigma_{ij} = (\sum_c f_i^c \ell_j^c)/V$ where, f is the contact force, $\ell$ is the branch vector (i.e., the vector joining the centers of particles interacting at contacts c) and the sum run over all contacts in a control volume V. Due to the geometry of the simple shear test, the stress state is considered invariant along the y direction. Thus, we may only consider the stresses of the tensor $\sigma$ on the shear plane $xz$.

According to the Mohr-Coulomb model, the effective macroscopic friction coefficient $\varphi$ is given by $\sin \varphi = q/p$,

where $p = (\sigma_1 + \sigma_2)/2$ is the mean stress, $q = (\sigma_1 - \sigma_2)/2$ is the deviatoric stress, and $\sigma_1$ and $\sigma_2$ are the principal values of the stress tensor on the shear plane [15]. The sample's density is quantified by the packing fraction $v = \sum_p V_p/V$, where $V_p$ is the volume of a particle p.

As a general observation (see Fig.2a), $\sin \varphi$ increases rapidly at the beginning of the shearing, and reaches a maximum value before decreasing towards a constant value at large strains (steady state). Along the same line, starting with an initially dense state, all packings dilate during shear and, hence, the volume increases (i.e., $v$ declines) and tend to a constant value presenting only minor fluctuations. Figure 2 shows $\sin \varphi^*$ and $v^*$, the averaged steady state values of $\sin \varphi$ and $v$, as a function of S. For the set of samples tested, we see that the average shear strength remains almost constant and around $\sin \varphi^* \approx 0.5$. In contrast the packing fraction $v^*$ increases with S. This finding thus extends to polyhedral particles assembly a similar behavior previously observed in the case of circular/polygonal particles in 2D and spherical particles in 3D.

The fact that the shear strength remains constant while the particle size polydispersity is increased reflects the organization of the microstructure and the features of force transmission that we discuss below.

<figure><img src="image_3.png" /><figcaption>Figure 2. (a) Evolution of the shear strength $q/p$ as a function of the shear deformation $\gamma$. In the inset: The solid fraction $v$ as function of $\gamma$. (b) In black squares: the average shear strength in the steady state as $\sin \varphi^*$ as a function of S. With red circles: the microstructural prediction of strength strength using Eq. (3). In the inset: The solid fraction at the steady state as a function of S.</figcaption></figure>
## 4 Micro-scale compensation mechanisms

At the contact scale, the local geometry and the mechanical equilibrium associated with two touching particles can be statistically characterized via the analysis of the angular distributions of the branch vector $\langle l \rangle(\Omega)$ and the contact force $\langle f \rangle(\Omega)$, where $\Omega$ is the solid angle. For coherence with the previous section, the branch and force vectors are projected on the shear plane 'xy' and can be thus split into their normal contributions $\langle l_n \rangle(\theta)$, $\langle f_n \rangle(\theta)$, respectively, and tangential contributions $\langle l_t \rangle(\theta)$, $\langle f_t \rangle(\theta)$, respectively, where $\theta$ is the orientation of $n'$, the projection of the real normal contact vector on the shear plane. Along with these distributions, we also define the angular distribution of contacts along the shear plane as $P(\theta)$. As shown in a number of previous studies, the analysis of these distributions should reveal the key microscopic elements in play as the particle size polydispersity increases.

<caption>**Figure 3.** Angular distribution of (a) contact orientations, (b) normal force and (c) tangential force intensities, and (d) normal branch and (e) tangential branch components for different values of $S$ at the steady state.</caption>

Under simple shear, these distribution take relatively simple shapes. Figure 3 displays polar representations of the five angular distributions $P_c(\theta)$, $\langle l_n \rangle(\theta)$, $\langle l_t \rangle(\theta)$, $\langle f_n \rangle(\theta)$ and $\langle f_t \rangle(\theta)$ as a function of the angle $\theta$, and averaged in the steady state for three different values of $S$. For all $S$ values, we observe an anisotropic behavior in $P_c(\theta)$ and $\langle f_n \rangle(\theta)$, where the peak value occurs, in average, along the major principal stress direction $\theta_{cr} \approx 3\pi/4$. Although less pronounced, an anisotropic orientation along the major principal stress is also seen for $\langle l_n \rangle(\theta)$. In contrast, the peak values for $\langle l_t \rangle(\theta)$ and $\langle f_t \rangle(\theta)$ occurs at $\theta_{ft} \approx \pi/2$. These distributions are well approximated by their lowest-order Fourier expansions [16, 17]

$$ \left\{ \begin{aligned} P_c(\theta) &= \frac{1}{2\pi}\{1 + a_c \cos 2(\theta - \theta_c)\} & (a) \\ \langle l_n \rangle(\theta) &= \langle l_n \rangle\{1 + a_{ln} \cos 2(\theta - \theta_{ln})\} & (b) \\ \langle l_t \rangle(\theta) &= \langle l_n \rangle\{-a_{lt} \sin 2(\theta - \theta_{lt})\} & (c) & (2) \\ \langle f_n \rangle(\theta) &= \langle f_n \rangle\{1 + a_{fn} \cos 2(\theta - \theta_{fn})\} & (d) \\ \langle f_t \rangle(\theta) &= \langle f_n \rangle\{-a_{ft} \sin 2(\theta - \theta_{ft})\} & (e), \end{aligned} \right. $$

where $a_c$ is the contact orientation anisotropy, $a_{ln}$ is the normal branch anisotropy, $a_{lt}$ is the tangential branch anisotropy, $a_{fn}$ is the normal force anisotropy, and $a_{ft}$ is

tangential force anisotropy. From Fig. 3, we get the corresponding privileged directions $\theta_c$, $\theta_{ln}$, $\theta_{lt}$, $\theta_{fn}$, and $\theta_{ft}$, that are nearly equal and aligned respect to the shear direction.

Figure 4 displays the averaged contact and branch anisotropies in the steady state, as function of $S$. We see that $a_c$ declines from 0.55 to 0.4 as $S$ increases. The decrease of $a_c$ traduce the fact that the mean connectivity of the particles increases. In fact, a detailed analysis of the connectivity as a function of particle size evidence that large particles are increasingly surrounded by a number of small particles as the size span is increased [9, 18]. In same time, the increases of normal branch anisotropies, from 0 to $\approx 0.09$, while the tangential branch anisotropy remains close to 0, evidence that longer branch appears preferentially in the direction of shear. This suggest that biggest particles self-organise along the strong force chains.

<figure><img src="image_4.png" /><figcaption>**Figure 4.** (a) Contact orientation anisotropy as a function of $S$ at the steady state. (b) The normal and tangential branch anisotropies ($a_{ln}$, $a_{lt}$) as a function of the particle size span $S$. Bars present the standard deviation of the data.</figcaption></figure>

The variations of the normal and tangential force anisotropies with $S$ are shown in Fig. 5. Basically, $a_{fn}$ slightly increases with $S$, from 0.4 to 0.5. In contrast the tangential force anisotropy remains independent with $S$ and close to 0.2. In fact, the above anisotropic parameters are the key elements behind the macroscopic friction coefficient $\varphi$. Indeed, it can be shown that the general expression of the stress tensor leads to the following simple relation [17, 19]:

$$ \sin \varphi^* \approx \frac{1}{2}(a_c + a_{ln} + a_{lt} + a_{fn} + a_{ft}). \quad (3) $$

The predicted values of $\sin \varphi^*$ by Eq. (3) are shown in Fig. 2 together with the measured values as a function of $S$. We see that the approximation of our data by Eq.
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Evolution of the normal and tangential force anisotropies ($a_{fn}, a_{ft}$) as a function of the particle size span $S$. Bars represent the standard deviation of the data.</figcaption></figure>

(3) is outstanding for all values of $S$. By virtue of this equation, the independence of $\sin \varphi^*$ with respect to $S$ results from a compensation mechanism between geometrical and mechanical source of anisotropy. More precisely, the decreases of $a_c$ is compensated by the increases of $a_{fn} + a_{ft}$ together with the fact that $a_{ft}$ is independent of $S$. Such compensation mechanism was also observed in previous studies with assemblies of discs, polygons and spheres but was unexpected for polyhedra particles assemblies. Indeed, when dealing with three-dimensional regular grains, a rich contact configuration appears given the different contact types that can occur: vertex-face, edge-face, edge-face, face-face and punctually, vertex-vertex. In particular, it is now well known that these contact contribute differently to the shear strength [17]. For example, face-face contacts are more likely to mobilize friction while vertex-face contacts, as more unstable interactions, are supporting relatively lighter stresses. But a surprising observation, detailed in [18], is that, in fact, the contributions of each family of contacts to the shear strength and anisotropies remain independent of the grain size span too.

## 5 Conclusion

In this paper, we developed a set of numerical experiments using the three dimensional Contact Dynamics method to generalize - or the opposite - the shear strength independence on particle size using very angular polyhedral grains. We built a set of samples composed of octahedra, where the grain size span varied from mono-size particles, to highly polydisperse particles with a size ratio (maximum/minimum) equivalent to ~ 10. The systems were sheared until a steady state, under quasi-static conditions, was reached.

Our numerical data clearly allow us to generalize to three-dimensional angular grains the observation that grain size span does not affect the shear strength although the solid fraction increases as a broader grain size distribution is considered. This unexpected result, consistently evidenced with more simple particle geometry, appear thus

to be a robust observation. By means of an additive decomposition of the shear strength in terms of anisotropic parameters related to particle connectivity, force transmission and friction mobilization, we show that subtle compensation mechanisms between anisotropies explains the invariance of shear strength as the particle size polydispersity is increased.

This long series of numerical experiments that include more and more realistic properties (from 2D discs assembly to 3D polyhedra systems), suggest that it is now time to develop controlled and systematic physical experiments to definitively conclude on this phenomenon.

## References

[1] T. Aste, Phys. Rev. E **53**, 2571 (1996)

[2] R.M. Baram, H.J. Herrmann, N. Rivier, Phys. Rev. Lett. **92**, 044301 (2004)

[3] J. Wiącek, M. Molenda, Int. J. Solids Struct. **51**, 4189 (2014)

[4] J. Wiącek, M. Stasiak, Granul. Matter **20**, 1 (2018)

[5] C. Voivret, F. Radjai, J.Y. Delenne, M.S. El Youssefi, Phys. Rev. Lett. **102**, 178001 (2009)

[6] N. Estrada, Phys. Rev. E **94**, 062903 (2016)

[7] E. Azéma, S. Linero, N. Estrada, A. Lizcano, Phys. Rev. E **96**, 022902 (2017)

[8] D.H. Nguyen, E. Azéma, P. Sornay, F. Radjai, Phys. Rev. E **91**, 032203 (2015)

[9] D. Cantor, E. Azéma, P. Sornay, F. Radjai, Phys. Rev. E **98**, 052910 (2018)

[10] S. Zhao, J. Zhao, N. Guo, Phys. Rev. E **101**, 012906 (2020)

[11] J. Hao, Y. Li, Y. Guo, H. Jin, J.S. Curtis, Powder Technol. **361** (2020)

[12] M. Jean, J. Moreau, Dynamics in the presence of unilateral contacts and dry friction: A numerical approach, in Unilateral problems in structural analysis, edited by G. Del Piero, F. Maceri (1987), pp. 151–196

[13] M. Jean, Comput. Methods in Appl. Mech. Eng. **177**, 235 (1999)

[14] F. Dubois, M. Jean, et al, LMGC90 wiki page, https://git-xen.lmgc.univ-montp2.fr/lmgc90/lmgc90_user/wikis/home (2019), [Online; accessed 28-Jul-2019]

[15] J.K. Mitchell, K. Soga, Fundamentals of Soil Behavior (Wiley, New-York, 2005), ISBN 978-0-471-46302-3

[16] L. Rothenburg, R.J. Bathurst, Geotechnique **39**, 601 (1989)

[17] E. Azéma, F. Radjai, F. Dubois, Phys. Rev. E **87**, 062203 (2013)

[18] D. Cantor, E. Azéma, I. Preechawuttipong, Phys. Rev. E **101**, 062901 (2020)

[19] H. Ouadfel, L. Rothenburg, Mech. Mater. **33**, 201 (2001)