Anisotropy of force distributions in sheared soft-particle systems

The contact forces in disordered materials, such as colloidal suspensions, foams, emulsions, and granular media are remarkably organized into highly heterogeneous force networks [1]. A statistical mechanical description of stress transmission in disordered media should provide a way to understand and predict the contact force distributions. The tail behaviour of the normalized normal force distribution $\text{P}(f_{n} \equiv F_{n}/\langle F_{n} \rangle)$ has received much attention, and several theoretical models with different assumptions and approaches [2,3] predict an exponential as well as a Gaussian tail. While early experiments and numerical simulations [4–6] favoured the exponential decay, further studies revealed that the decay can also be faster than exponential [7–11]. A recent numerical study [12] of frictional soft-particle systems under pure compression showed that, independently of the distance from jamming, the tail behaviour can be described by a stretched exponential with an exponent around 1.8, which slightly depends on the choice of the contact force law, the friction coefficient, and the relative particle stiffness in tangential and normal directions.

In sheared systems, a slower decay of $\text{P}(f_{n})$ compared to isotropic packings has been observed [8,11,13,14], where increasing the shear stress enhances the broadening of $\text{P}(f_{n})$. This necessitates further efforts to provide a comprehensive description of the mechanisms underlying stress propagation in sheared systems. In this letter, the force distributions in periodic 2D granular systems under non-cyclic pure shear are studied. The shear-induced stress anisotropy is taken into account by categorizing the contacts in terms of their orientation. While the normal force distribution decays even faster than exponentially for the contacts oriented along the same direction, it is shown that averaging over all angle-resolved distributions leads to the broader shape of the overall distribution $\text{P}(f_{n})$ in the sheared system. Thus, a connection between the shear-induced stress anisotropy and the broadening of $\text{P}(f_{n})$, is established, which results in the saturation of broadening of $\text{P}(f_{n})$ at high shear deformations. In the asymptotic strain-independent regime, the distribution for any given direction nearly follows a Gaussian form. This enables one to integrate over all directions and obtain an approximate analytical expression for the invariant broad shape of $\text{P}(f_{n})$ at the limit of large shear deformations.

The distribution $\text{P}(f_{t})$ of tangential forces decreases monotonically (in contrast to $\text{P}(f_{n})$ that usually developsa peak) with a broad exponential-like tail [4,5,11,12,15,16]. Moreover, the collapse of $\text{P}(f_{t})$ curves has been reported for different values of inter-particle friction coefficients [5], and for different isotropic [12,15] or anisotropic [15] applied loads. Here, it is verified that the angle-resolved tangential force distributions nearly collapse onto a universal curve for different orientations. Therefore, the shape of $\text{P}(f_{t})$ remains approximately invariant with the applied load, in contrast to the distinctive shape of the normal force distribution $\text{P}(f_{n})$ under isotropic or shear strain.

The evolution of the contact forces during a quasi-static, pure shear deformation of an initially compressed packing of disks was studied numerically. The simulations were carried out by means of discrete element methods. The inter-particle forces are modelled by damped, linear springs, for both normal and tangential interactions, using the spring constant ratio $k_{t}/k_{n}=0.5$. Additionally, the tangential forces obey Coulomb's law with a friction coefficient set to $\mu=0.5$. The two-dimensional simulation cell with fully periodic boundary conditions contains nearly 20000 disks. The radii are uniformly distributed in the range $[0.8 \bar r,1.2 \bar r]$. The average particle radius $\bar r$ is the length unit, and $k_{n} \bar r$ is taken as the unit of force in the following.

The initial configuration is generated by placing the particles randomly into the simulation box, without accepting any overlap between them. This unjammed system is subsequently compressed quasi-statically by applying consecutive steps of incremental compression and relaxation. The compression is achieved by re-scaling the particle positions while keeping their radii fixed. The relaxation procedure ensures that the net force exerted on each particle is 8 orders of magnitude below the mean contact force. After the average normal overlap in the jammed state reaches a desired threshold, the system is sheared in an analogous manner with an applied pure shear deformation, so that the aspect ratio of the rectangular simulation box is changed, but the volume is kept constant. Each time the system is equilibrated, the force state of the packing is stored. Upon increasing shear strain, the fabric and the force network change and anisotropies develop. The texture $\text{P}(\alpha)$ and the average normal and tangential forces can be well approximated by a second-order Fourier expansion as [17]

$$\begin{eqnarray} \begin{array}{rcl} \text{P}(\alpha) &=& \frac{1}{2\pi}\left(1+a \cos(2\alpha)\right),\\ \bar f_n(\alpha) &=& 1+a_n \cos(2\alpha),\\ \bar f_t(\alpha) &=& a_t \sin(2\alpha), \end{array} \end{eqnarray} \tag{ 1 }$$

where a is the fabric anisotropy, and a_n and a_t represent the mechanical anisotropies in normal and tangential directions. Figure 1 shows how the anisotropies develop with the increase of the shear strain. The relation between the stress ratio and the anisotropies follows $\tau/p=(a+a_n+a_t)/2$ [17] with very small deviations throughout the shearing process (not shown).

Zoom In Zoom Out Reset image size

The components of the globally averaged stress tensor are measured using

$$\begin{equation} \sigma_{ij} = \frac{1}{A} \sum_{c=1}^{N_{c}} f_{i}^{c} r_{j}^{c}, \end{equation} \tag{ 2 }$$

where A is the area of the system, $f_{i}^{c}$ the i-th component of the force acting on contact c and $r_{j}^{c}$ the j-th component of the branch vector. The sum runs over all contacts N_c in the system. Denoting the eigenvalues of the stress tensor by $\sigma_{1}$ and $\sigma_{2}\ (\sigma_{1} \leq \sigma_{2})$, the pressure and shear stress are given by $p=\frac{1}{2} (\sigma_1 + \sigma_{2})$ and $\tau = \frac{1}{2}(\sigma_{2}-\sigma_{1})$, respectively. When the isotropic system with aspect ratio e = 1 is subject to a pure shear deformation, the engineering shear strain γ increases as the aspect ratio decreases since $\gamma = \frac {1-e^{2}}{2e}$. The principal axes of the stress rotate less than $1.13^{\circ}$ with respect to the biaxial deformation directions throughout the shearing process. Figure 1 shows that the stress ratio $\tau/p$ grows as γ increases, and eventually saturates for large shear strains $(0.01 \leq \gamma)$. Note that the shear strain at which $\tau/p$ and the shear stress saturate depends on the volumetric strain applied on the initial isotropic packing [18]. The fabric anisotropy in fig. 1 shows a similar saturation behavior as the normal force anisotropy a_n, while the pressure p did not saturate for the maximum applied shear strains $(\gamma = 0.015)$ in this study.

Upon increasing the shear deformation γ, the normalized normal force distribution $\text{P}(f_{n})$ broadens, as shown in fig. 2. Similar results were observed in numerical studies [7,8,10] as well as in experiments with photo-elastic particles [11]. A crucial question is how the shape of $\text{P}(f_{n})$ is influenced by the characteristics of the globally imposed stress, namely p and τ. In spite of the conserved volume during the pure shear deformation, the pressure can change and one may partially attribute the shape change of $\text{P}(f_{n})$ at different values of γ to the difference between their pressures. Moreover, it is known that shearing induces anisotropies, leading to spatial correlations between contact forces with direction-dependent correlation lengths [11]. While both p and τ seems to influence the shape of $\text{P}(f_{n})$, the evolution of $\text{P}(f_{n})$ in fig. 2 remarkably slows down at high γ and eventually saturates, which is reminiscent of the behaviour of stress ratio $\tau/p$ and anisotropy development in fig. 1. Note that when the stress anisotropy approaches an invariant state, the shape of $\text{P}(f_{n})$ does not vary any more.

Zoom In Zoom Out Reset image size

To elucidate the influence of shear-induced stress anisotropy, categorizing the contacts according to their orientation provides useful information about the angular dependence of force transmission. Therefore, the joint probability distribution $\text{P}(f_{n},\alpha)$ for the normal force and the contact angle is introduced. The contacts are divided into 12 angular bins of 15° each. Next, the angle-dependent conditional distribution $\text{P}(f_{n}|\alpha)=\frac{\text{P}(f_{n},\alpha)}{\text{P}(\alpha)}$ is calculated in each bin, where $\text{P}(\alpha)= \int_{0}^{\infty} \text{P}(f_{n},\alpha) \text{d} f_{n}$. Three examples of the resulting distributions are shown for different directions in fig. 3. Similar distributions have been recently reported in 3D packings under periodic uniaxial shear [19] as well as plane shear in a split-bottom Couette cell [20]. One finds that fits of the form

$$\begin{equation} \text{P}(f_{n}|\alpha)=\frac{1}{N(\alpha)} f_{n}^{\nu_{n}(\alpha)}\exp\left[-\left| \frac{f_{n}-b_{n}(\alpha)}{2 w_{n}(\alpha)}\right|^{\delta_{n}(\alpha)}\right], \end{equation} \tag{ 3 }$$

characterize the overall shape of $\text{P}(f_{n}|\alpha)$ along different directions and for different values of shear strain γ. Note that there are only three independent fit parameters in the above equation, due to the normalization constraint $\int_{0}^{\infty} \int_{0}^{2\pi} \text{P}(f_{n},\alpha) \text{d}\alpha \, \text{d}f_{n}=1$, and the constraint on the first moment of distribution $\int_{0}^{\infty} \int_{0}^{2\pi} f_{n} \text{P}(f_{n},\alpha) \text{d}\alpha \, \text{d} f_{n}=1$. The choice of force distribution in eq. (3) is inspired by the recent work by Tighe et al. [3], where a similar function (even though with slight differences) was proposed based on entropy maximization arguments with respect to the allowed force network ensemble [21] in isotropic systems. Note that other variants of the fit function have been also proposed, see, e.g., [22] for the force distributions in 3D isotropic packings. The dimensionality may play a crucial role in determining the shape of the contact force distributions, as, e.g., numerically verified for the tail behaviour of $\text{P}(f_{n})$ in [14].

Zoom In Zoom Out Reset image size

The angular dependence of the fit parameters in the initial isotropic packing $(\gamma=0)$ is compared with a highly sheared case $(\gamma=0.015)$ in fig. 4. One observes that the fit parameters are practically α-independent except for b_n which develops a pronounced angular dependence during shearing. This behaviour can be understood by comparing eq. (3) in the special case of $\delta=2$ with the derivation in [3]. There w and ν are related to the local force balance constraint on the grains, the friction coefficient, and the connectivity of the force network, thus, they are not expected to be angle dependent. On the other hand, b_n is set by a constraint to the pressure. Since the average normal force varies with α in the presence of stress anisotropy in the system, it becomes clear why b_n is α-dependent. However, note that pressure is not the only control parameter in determining the shape of the angle-resolved distributions. $\text{P}(f_{n}|\alpha)$ along a given direction in the sheared system notably differs from $\text{P}(f_{n})$ of an isotropic packing carrying the same average normal force [12]. Our generalized equation (3) captures the shape of sheared force distributions by allowing the exponent δ to act as an additional free parameter. Nevertheless, the global constraint on the applied shear stress/strain has to be taken into account to obtain an analytical expression for $\text{P}(f_{n})$ in sheared packings.

Zoom In Zoom Out Reset image size

The parameter b_n varies periodically with a peak in the direction of compression, which can be described by a second-order Fourier expansion of the form

$$\begin{equation} b_{n}(\alpha) = b_{n} \left[1- a_{b} \cos(2\alpha)\right]. \end{equation} \tag{ 4 }$$

This equation and the finding that the rest of fit parameters do not develop a clear angular dependence during shearing motivates us to propose a similar functional form for the joint probability distribution $\text{P}(f_{n},\alpha)$, assuming that only the shift parameter b_n has an angular dependence according to eq. (4), i.e. 

$$\begin{equation} \text{P}(f_{n},\alpha)=\frac{1}{N} f_{n}^{\nu_{n}} \exp\left[-\left|\frac{f_{n}- b_{n} \left[1- a_{b} \cos(2\alpha) \right]}{2w_{n}} \right|^{\delta_{n}}\right]. \end{equation} \tag{ 5 }$$

The results of the fits via eq. (5) are shown in fig. 5 by the evolution of the fit parameters of eq. (5) for increasing γ. Note that they approach an invariant state for large shear strains.

Zoom In Zoom Out Reset image size

By integrating the regulated form of $\text{P}(f_{n},\alpha)$ in eq. (5) over α, one obtains the overall distribution

$$\begin{equation} \text{P}(f_{n}) = \int_{0}^{2\pi} \text{P}(f_{n},\alpha)\, \text{d} \alpha, \end{equation} \tag{ 6 }$$

which can be compared to $\text{P}(f_{n})$ obtained from the simulations. Using numerical integration (due to the non-integer exponent δ), $\text{P}(f_{n})$ was obtained, e.g. for the packing with $\gamma=0.015$. The resulting curve, shown in fig. 2, matches the simulation results. It was also checked, whether the resulting $\text{P}(f_{n})$ reproduces the anisotropy a_n obtained directly from the simulation data. Figure 6 shows that both anisotropies are in good agreement, with small deviations for large anisotropies.

Zoom In Zoom Out Reset image size

Interestingly, the exponent of the stretched exponential $\delta_n$ increases during shearing and approaches two, i.e. it nearly follows a Gaussian tail at high γ (see fig. 5). This allows one to analytically integrate the angle-resolved distribution (i.e. by combining eqs. (5) and (6) using $\delta=2$) and obtain an approximate γ-invariant expression for the marginal distribution $\text{P}(f_{n})$ in the limit of large shear strains,

$$\begin{eqnarray} \lim_{\gamma \to \infty} \text{P}(f_{n}) &\approx \frac{2\pi}{N} f_{n}^{\nu} \exp \left[-\left(\frac{f_{n}-b_{n}}{2w}\right)^{2}\right] \nonumber \\ &\quad \times\,\text{I}_{0} \left(\frac{a_{b}b_{n} (b_{n}-f_{n})}{2w^{2}}\right), \quad \end{eqnarray} \tag{ 7 }$$

where $\text{I}_{0}$ is the modified Bessel function of the first kind. In the integration, all quadratic terms in $a _{b}\cos(\alpha)$ are neglected. The above expression is compared to the simulation data in fig. 2, which shows a satisfactory agreement. As expected, the decay is slightly faster than simulations, since a pure Gaussian exponent is used to obtain eq. (7) despite the fact that $\delta_{n}$ converges to an exponent slightly below 2.

Note that after extensive slow shearing, a system is expected to reach a critical state of flow, where the force state attains a statistically steady-state condition. However, we limit the application range of our results to quasi-static deformations because yielded systems may behave differently, as the flow properties are in general shear-rate dependent [23]. Investigation of granular flows is beyond the scope of this letter.

The marginal distribution of the tangential forces $\text{P}(\tilde f_{t})$, using $\tilde f_{t} = |f_{t}|/\left<|f_{t}|\right>$, decreases monotonically, as shown in fig. 7. There is no significant change of the distribution during shearing. This result, together with the fact that the anisotropy a_t of the average tangential forces is about an order of magnitude below the anisotropy a_n of the average normal forces in the sheared system (see fig. 1), shows that while friction stabilizes a packing at a lower coordination number, the main history dependence of the contact forces is observed in the normal forces.

Zoom In Zoom Out Reset image size

The normalized contact force distributions in sheared systems of soft frictional particles were studied numerically. The angle-resolved distribution $\text{P}(f_{n}|\alpha)$ along an arbitrary direction α decays faster than exponentially, similarly to the behaviour of the isotropic packings. However, different contact orientations carry different stresses, which affects the shape of $\text{P}(f_{n}|\alpha)$ at small forces as well as its tail behaviour. The resulting overall distribution $\text{P}(f_{n})$, obtained from averaging over different contact orientations, has a modified shape and an evident slower decay which does not necessarily follow the same form as the angle-resolved distributions. Therefore, a link between the broadening of the normalized normal force distribution $\text{P}(f_{n})$ and the shear-induced stress anisotropy was established.

The broadening is enhanced with increasing shear strain, as far as the stress anisotropy still develops. Eventually, the stress anisotropy saturates at high shear deformations, thus, $\text{P}(f_{n})$ reaches a strain-independent shape. Since the anisotropy development is considerably weaker in the tangential direction, $\text{P}(f_{t})$ remains approximately invariant through the deformation process. The fabric anisotropy remains small throughout the shearing process and has no major influence on the shape of the force distributions.

These findings show that the stress propagation in sheared systems can be better understood when angle-resolved distributions are considered. Analytical treatments [3] to obtain the force distributions in isotropic packings need to be reconsidered in sheared systems by taking the global constraint on the applied shear stress/strain into account, as a step forward towards fully describing the shape of force distributions in sheared packings. The results also help to better understand the mechanisms of deformation of granular materials at the microscopic level, which facilitates the development of stochastic approaches [24] for theoretical modelling of deformation and elastic behaviour of granular systems.

We would like to acknowledge the support by the German Research Foundation (DFG) via priority program SPP 1486 "Particles in Contact" and the Center for Computational Sciences and Simulation of the University of Duisburg-Essen.