Shear zones in granular mixtures of hard and soft particles with high and low friction

Abstract

Granular materials show inhomogeneous flows characterized by strain localization. When strain is localized in a sheared granular material, rigid regions of a nearly undeformed state are separated by shear bands, where the material yields and flows. The characteristics of the shear bands are determined by the geometry of the system, the micromechanical material properties, and the kinematics at the particle level. For a split-bottom shear cell, recent experimental work has shown that mixtures of hard, frictional and soft, nearly frictionless particles exhibit wider shear zones than samples with only one of the two components. To explain this finding, we investigate the shear zone properties and the stress response of granular mixtures using discrete element simulations. We show that both interparticle friction and elastic modulus determine the shear-band properties and packing density of granular mixtures of various mixing ratios, but their stress response depends strongly on the interparticle friction. Our study provides a fundamental understanding of the micromechanics of shear band formation in granular mixtures.

---

1 Introduction

Strain localization is a characteristic effect of sheared granular matter, which is not observed in the flow of other materials, such as liquids. Once strain localizes within a narrow shear zone (also referred to as the shear band), the stress, strain, and strain rate of the material assume high values within that zone. On the contrary, outside of it, the material remains nearly undeformed. The knowledge of the properties of shear zones and the factors that influence them is crucial for understanding the intrinsic properties of granular materials under shear.

In many systems where particles come into contact with boundaries, narrow shear zones appear at the interface because of particle–boundary interactions. These shear zones are triggered by boundary effects, which hinder the accurate characterization of bulk material properties based on the shear zone characteristics. In our work, we present a unique setup called a split-bottom shear cell, which is a modified version of a Couette cell. This setup involves splitting the cell's bottom at a specific radius (R_s), resulting in two parts consisting of the inner disk and outer cylinder components (see Fig. 1). After loading the cell with grains to a certain height (H), the inner disk is set to rotational motion. This design allows granular materials to be sheared to develop wide shear zones inside the bulk of the materials, far away from the rigid boundaries, which allows for a comprehensive characterization of their yielding properties and flow behavior at a macromechanical level. This is quite different from the usual narrow shear bands that we observe close to the boundary, which are formed due to particle–boundary interactions. This shear cell setup has been used in experiments¹ and simulations² to explore the bulk properties of granular systems.

Fig. 1 Schematic of a split bottom shear cell of radius R₀ with granular material at filling height H, illustrating the main system components. The bottom plate (dark-red color) rotates around the vertical axis at constant rotation frequency ω₀, while the outer wall (marked in grey) stays stationary. The differential motion triggers the formation of a shear band (marked in orange color) starting from the split radius R_s, which widens with increasing height. In contrast, the material outside the shear band (marked in yellow) stays undeformed.

The universal geometrical properties of shear bands in the split-bottom shear cell have been studied extensively in the literature.^3–5 The absolute value of the macroscopic normalized angular velocity of the granular material was measured as a function of the radial position, and it has been found to follow closely the shape of an error function, which is characterized by two parameters: the width W and the center position R_c of the shear zone.^4,5 The width of the shear band increases with increasing height within the bulk of the material, while the distance of its central position R_c from the center of rotation decreases (that is, the shear band moves inward with increasing height). The width of the shear band is affected by the size, shape, and interaction properties of the particles, such as friction between the particles, cohesion, and the shear rate.^6–9 The position of the shear band, R_c, at the surface of the granular material in the shear cell, depends on two geometrical parameters, the filling height H and split radius R_s.¹⁰ The bulk profile of the shear band as a function of height, R_c(z), can be linked to a variational problem. The solution of the function R_c(z) that minimizes the integral representing the dissipation rate, is related to the minimal mechanical torque and thus depends on the material properties. Therefore, the shape of the shear band within the bulk is not only influenced by geometric factors such as the split position (R_s) and the filling height (H) but is also closely connected to the material properties, e.g., the shear stress, and the packing fraction of the material, as highlighted in previous studies.^5,10

Interparticle friction and particle softness play a key role in influencing material bulk flows, which have been widely studied. Wang et al.¹¹ experimentally characterized the shear zones of binary mixtures made of hard, frictional and soft, slippery particles. Their study identified distinct shear band characteristics specific to the mixture, distinguishing them from those in pure-species systems. Luding⁸ explored the impact of interparticle friction on shear zones in split-bottom shear cells through numerical simulations. Increasing friction resulted in narrower shear bands that shifted towards the center. Ashour et al.¹² found that particle softness adds qualitatively new features to the dynamics and the packing characteristics of silo flows. Soft, deformable particles, such as hydrogel particles, can adapt to local stress conditions, leading to alterations in pressure profiles and flow fields. Hong et al.¹³ highlighted the importance of particle softness in hopper flows, where gravitational forces must be less than repulsive forces from particle stiffness to induce clogging. Harth et al.¹⁴ observed longer time scales for internal rearrangement in soft hydrogel particles during silo orifice flow compared to hard grain systems, where shorter time scales have been reported by Unac et al.¹⁵ This was attributed to mass flow rate differences, roundness, and increased viscous damping of hydrogel grains. The rate-dependent frictional regimes observed in hydrogel particle suspensions have been attributed to their high deformability and low friction.¹⁶ Wang et al.¹⁷ observed experimentally that, unlike granular materials with hard particles, soft slippery particles in a silo flow do not follow Beverloo's¹⁸ law, according to which the flow rate is independent of the overburden pressure. Adding just 5% hard particles to soft silo flows restored this behavior and significantly altered system flow properties. Based on the micromechanical observations, Saitoh et al.¹⁹ highlighted that the evolution of the probability distribution functions (PDFs) of forces in soft particle packings under quasi-static deformations extends beyond just particle contacts. Inter-particle gaps between nearest neighbors are also crucial and affect the stochastic evolution of the PDFs of forces. Recently, Luding et al.²⁰ studied the contact and force networks of frictionless soft granular materials, revealing that the evolution of the pressure in the system was correlated with the existence of loops formed in the force network during particle rearrangements. Götz et al.²¹ observed that soft particles with larger deformations and increased contact areas result in larger frictional forces and improved holding capacity of a granular gripper. In summary, both particle friction and softness are crucial factors influencing bulk flow behavior, and the intricate interplay between these properties and their combined effects on bulk behavior remains a subject of ongoing investigation.

The findings of Wang et al.¹¹ established that mixtures of hard and soft particles exhibit wider shear zones than assemblies made of only one of the individual materials (only hard or only soft particles). However, the underlying reason for this behavior remains unclear. Our objective here is to gain insights into the combined effects of microscopic parameters, namely particle stiffness and inter-particle friction, on the macroscopic shear zone characteristics for ensembles of soft, slippery particles and hard, frictional particles sheared in a split-bottom shear cell.

2 Numerical setup and materials

We use an open-source discrete element method (DEM) code MercuryDPM²² to simulate sheared systems of hard and soft particles. The numerical setup used for this study is a split-bottom shear cell, as shown in Fig. 1, where the inner bottom plate rotates around a central axis at a constant angular velocity, ω₀ = 0.52 rad s⁻¹, while the outer boundary of the cell remains stationary. The granular material located in the outer parts remains static, while the material near the center of the cell rotates with the bottom disk. The shear in this system is localized in an intermediate region, defining a wide shear band as indicated by the dark orange colour in Fig. 1. The position and the width of the shear zone depend on the ratio H/R_s, where H is the filling height and R_s is the split radius of the rotating bottom plate.^4,23 For small and moderate filling height (H < 0.6R_s), the shear zone starts from the bottom split position, evolves through the entire height of the granular layer, and ends up at the surface,³ as shown in Fig. 1, while the inner part of the granular material co-rotates synchronously with the bottom plate. In our numerical setup, the filling height is H = 20 mm and H/R_s = 0.42, which is considered a moderate filling height. The inner radius of the shear cell is R₀ = 75 mm. This numerical setup is similar to that employed in the experimental work by Wang et al.,¹¹ with the difference that the experimental configuration incorporates an additional small rod in the center of the cell, where the shear load was imposed.

2.1 Material parameters

We consider a mixture of soft low-frictional particles and hard frictional particles. In this section, we describe the generation of binary granular mixtures of such particles. We take the properties of hydrogel spheres as the soft, low-frictional material component, and mustard seeds represent hard, frictional particles. A reference sample of pure soft hydrogel particles consists of N⁰_A = 26 795 polydisperse particles of species A with mean diameter, d_A = 2.5 mm. Starting from this baseline sample volume, we create mixtures by substituting N⁰_A − N^T_A particles of species A with particles of species B of diameter d_B = 3 mm, such that the same volume of material A is replaced by B, i.e. V_B = (N⁰_A − N^T_A)(π/6)d_A³ = N^T_B(π/6)d_B³. For both species A and B, we use particles of homogeneous size distribution of mean particle diameter d_A and d_B, respectively, with polydispersity of ±8%. To characterize the mixtures with different ratios of species B, we define the mixing ratioWe analyze systems with X ∈ {0, 0.05, 0.25, 0.50, 0.75, 0.85, 1.0}. Table 1 presents the number of particles of the two species in each sample. Fig. 2 shows some examples of different mixtures.

Table 1 Number of particles of the species in each sample

X	0	0.05	0.25	0.50	0.75	0.85	1
Species A	26 975	24 675	17 104	9757	4244	2497	0
Species B	0	1361	5695	9917	13 160	14 133	15 465

---

Fig. 2 Top view of the shear cell for granular mixtures of different fractions of hard particles X, where the blue and red particles represent species A (soft) and B (hard), respectively.

The properties of soft particles of species A correspond to that of hydrogel spheres with friction coefficient μ_A = 0.005 and elastic modulus E_A = 50 kPa. The frictional properties of hard particles correspond to mustard seeds with friction coefficient μ_B = 0.20 and the elastic modulus E_B = 5000 kPa. It is noteworthy that the true density of mustard seeds and hydrogel spheres is approximately the same.^24–27 For simplicity, we considered the density of both species to be equal, set at 1017 kg m⁻³. All relevant material and interaction properties of the two species are given in Table 2.

Table 2 Microscopic material parameters for the model

Properties	Species A	Species B
Mean particle diameter	2.5 mm	3 mm
Density	1017 kg m^−3	1017 kg m^−3
Elastic modulus	50 kPa	5000 kPa
Poisson's ratio	0.40	0.40
Friction coefficient	0.005	0.20

---

2.2 Contact model

In the soft-particle DEM, the evolution of particle positions and velocities is calculated by integrating Newton's equation of motion. When two particles i and j interact, they can slightly overlap. Only viscoelastic interactions are considered, i.e., after separating the contacting particles, they recover their initial shape. The compression of contacting particles ξ is written as:where R_i, R_j are the particle radii and r⃑_i, r⃑_j are the particle positions. The normal force is calculated as in ref. 28, where an elastic Hertzian normal force is considered in conjunction with a viscous damper of velocity-dependent coefficient of restitution.withwhere E* is the effective elastic modulus, R* = R_iR_j/(R_i + R_j) is the effective radius, is the normal unit vector, and A_n = 4.1 × 10⁻⁶ s is the normal dissipative parameter, calculated as in ref. 29, considering a coefficient of restitution of 0.70 for the tangential translational velocity of the rotating bottom plate at the split radius. The effective elastic modulus E* is given bywhere E_i, E_j and ν_i, ν_j are the elastic moduli and Poisson ratios of particles i, j, respectively.

We model the tangential viscoelastic forces considering a path-dependent calculation, following the no-slip solution of Mindlin³⁰ for the elastic part and Parteli and Pöschel³¹ for the tangential dissipative parameter A_t ≈ 2A_nE*, which are capped by the static friction force between two particles. The tangential force is thereby given by

where μ is the friction coefficient, e⃑_t is the unit vector in the tangential contact direction and G* the effective shear modulus, given byThe shear modulus is calculated assuming isotropic material. For the interaction between different species, the effective friction coefficient is considered as . It is typical for many DEM formulations to consider the minimum coefficient of friction when two particles of different materials come into contact. Since the harmonic average of two values favors the smaller one, considering this calculation of the friction coefficient for contacts between high and low friction coefficients gives values close to considering the minimum value. This assumption is chosen as a milder approach to simply assuming the minimum friction coefficient of the two, and more experimental evidence is required to understand the friction between hard and soft particles. Considering the harmonic average here, we have three scenarios of possible contacts between frictional–frictional, slippery–slippery and frictional–slippery (with friction very close to slippery–slippery) particles.

3 Segregation effects

When granular mixtures with particles differing in size and/or density are subjected to external forces such as shaking, stirring, or shearing, they often segregate to form complex patterns. These patterns of segregation can lead to different and inhomogeneous macroscopic behavior. In many applications where a homogeneous packing is sought, this is undesirable, and significant effort is often invested in attempts to avoid or control segregation during processing and handling of granular materials. Shear-induced percolation segregation that occurs due to kinetic sieving and squeeze expulsion^32,33 is the dominant mechanism in dense granular flows. In this study, we study a bi-disperse mixture of two species that differ slightly in terms of particle size. Therefore, in this section, we examine the extent to which particles segregate in this system of hard frictional and soft slippery particles, sheared in the split-bottom shear cell.

First, we analyze the probability distribution of the particles of each species in the radial r and vertical z directions for a given snapshot after λ = 104 radians of shear cell rotation (200 s of simulation time), as shown in Fig. 3(a) and (b), respectively, for X = 0.50. The distribution increases in the radial direction for both species due to the increased measuring volume as the probability is calculated using co-centric cylinders of increasing radius, aligned with the center of the cell. Note that the distribution of the two species overlaps with each other in the radial r direction, indicating that there is no segregation effect of the two species horizontally. The distribution in the vertical z direction for species B is nearly constant. However, the distribution of species A decreases in vertical z direction from 0.065 to 0.035, indicating a higher probability of species A near the bottom of the shear cell than near the free surface. This is a signature of the segregation of species A in the vertical direction due to kinetic sieving, and is attributed to the smaller particle size compared to species B. To further analyze this effect, Fig. 4(a) and (b) illustrate the mean radial positions based on the weighted surface area and perpendicular positions z_mean of the two species as a function of the shear displacement λ for X = 0.50. Here, it is to be noted that the mean position of Species A is about 0.40 mm less in the r direction and 1.20 mm less in the z direction than that of Species B. However, the differences are negligibly small and of less than one particle diameter.

Fig. 3 Probability of number of particles, P(N), as a function of (a) radial distance and (b) distance from the bottom surface for species A and B after λ = 104 radians of shear cell rotation (200 s of simulation time) for granular mixtures of X = 0.50.

Fig. 4 Average radial and bottom surface distances (a) and (b) for species A and B as a function of shear displacement λ = ωt_rot in granular mixtures with X = 0.50. Insets show top view (a) and cross-sectional front view in yz plane (b) at λ = 0 and λ = 104 radians for left and right inset figures, respectively.

4 Micro–macro transition

The current work aims to derive appropriate macroscopic fields commensurate with the shear band features and macromechanical stress analysis, e.g., strain rate, velocity profile, and stress, based on the given micromechanical properties. The flow that rapidly relaxed to a steady state is purely azimuthal and is proportional to the driving rate of the shear cell, ω₀. The averaging is thus performed in toroidal volumes over many snapshots of time, leading to generic continuum fields Q(r,z) as a function of the radial and vertical positions. We choose the representative elementary volume according to the ref. 34 as equal to the minimum particle diameter d_p = 1.37 mm.

Shearing of a granular system leads to volume dilation and build-up of shear stress and anisotropy in fabric. Since we are interested in the steady-state flow profile, we probe for the minimum strain required to reach the steady-state flow regime. We define the shear displacement λ = ωt_rot through which the system rotates (where t_rot is the simulation time) as the shear control parameter. As reported by Singh et al.,³⁵ both the global quantities of kinetic energy and average number of contacts reach a steady state very fast (λ ∼ 5 radians). We analyze the relaxation of local quantities, specifically the local velocity profiles, and conclude that λ ∼ 20–30 radians is the amount of shear displacement required to form a stable shear band, which is in agreement with Ries et al.³⁶ The data depicting the initial, transient evolution of local shear band properties is not presented here, as our focus is solely on the steady-state characteristics of the shear band. Consequently, we perform the local averaging over almost 600 snapshots distributed over λ ∼ 80 to 100 radians.

The consistency of the local averaged quantities also depends on the accumulated local shear strain during the averaging time. We concentrate our interest in the region where the system can be considered to reach a critical state. The critical state is a unique steady state reached after extensive shear, where the material deforms with applied strain without any further change in normal stress, shear stress, and volume fraction. This is a state where the system forgets its loading, packing, or other sample-preparation history.³⁷ In accordance with previous experimental and numerical results in the same setup,^38,39 the criterion of large strain rate identifies the shear band region, e.g., higher than a critical strain rate of [small gamma, Greek, dot above]_c = 0.08 s⁻¹. The shear band center at a given height is defined by the region where strain rates higher than 80% of the maximum strain rate at that height in the shear cell. Therefore, the critical strain rate here is a function of height in the shear cell and is denoted as [small gamma, Greek, dot above]_c(z) = 0.8[small gamma, Greek, dot above]_max(z).

5 Flow profile of granular mixtures

Once the flow of the mixtures relaxes to steady state shearing, the azimuthal velocity is proportional to the driving rate ω₀. Fig. 5(a) shows the scaled angular velocity profile ω(r, z)/ω₀. Within the limits of assumptions in our simulations, we do not see a clear difference in the flow profiles of the three cases of pure particles and mixtures for X = 0, 0.50, and 1. However, we observe that the shear zone of the mixture (X = 0.50) is slightly more diffused compared to the pure case simulations involving only soft (X = 0) or only hard particles (X = 1). This indicates a wider shear band for the mixture, although it is not clearly concluded from this figure. A closer and more precise introspection on the shear band properties is therefore essential to conclude the shear zone properties of the mixture. On this ground, we analyze the dimensionless ratio of the measured angular velocity ω/ω₀ as a function of the radial coordinate r. The dimensionless velocity profile in a split bottom shear cell is well approximated by an error function^5,10where R_c and W are the position and width of the shear band developed above the split at R_s, respectively. Both the position and the width of the shear band depend on the height z in the system. Therefore, we extract the position and the width of the shear band at layers of varying height along the filling height in the shear cell by fitting the scaled ω(r) against the radial position r data. With the goal to extract quantitative data for the shear band features, Fig. 5(b) shows values of non-dimensional angular velocity ω(r)/ω₀ at different heights in the bulk of granular materials against the radial distance.

Fig. 5 (a) Shear band profiles for mixtures of soft and hard particles in different mixing ratios. The color scale illustrates the normalized ratio of macroscopic angular velocity ω(r, z) scaled by the rotation velocity of the bottom disk ω₀. R_s marks the split radius, which serves as the origin of the shear band in space. (b) Graphical representation of macroscopic angular velocity ω(r) normalized by the bottom disk angular velocity ω₀, as a function of radial position r for a fraction of hard particles X = 0.50 at different heights. The inset shows the scaled angular velocity ω(r)/ω₀vs. radial position r fitted by eqn (8), also showing the features of the shear band.

A quantitative evaluation of the strain rate [small gamma, Greek, dot above], expressed in terms of the velocity gradient components ∂v_θ/∂r and ∂v_θ/∂z, is given as:⁴⁰

Fig. 6(a) shows the strain rate [small gamma, Greek, dot above](r, z) evaluated from eqn (9) and plotted as a function of radial and vertical position for the sample with X = 0.50. The dashed lines overlaid on top of the contour graph indicate the center R_c and the edges R_c + W and R_c − W of the shear band, indicated in the figure, as obtained from the fit function in eqn (8).

Fig. 6 (a) Contour plot of strain rate [small gamma, Greek, dot above](r, z) for X = 0.50, with the overlaid lines indicating shear band center R_c and edges R_c ± W from the fit function in eqn (8). (b) Shear band center R_c and edges R_c ± W for varying fractions of hard particles X at z_top = 16 mm, (c) scaled shear bandwidth W/W_topvs. scaled height z/H and (d) relative shear band position R_c/R_svs. scaled height z/R_s for X = 0,0.50 and 1. (e) Shear bandwidth W vs. fraction of hard particles X and (f) shear band position R_cvs. fraction of hard particles X at z = z_top. Error bars represent fitting errors, determined by fitting ω/ω₀vs. r to eqn (8). Markers denote X = 0 (), X = 0.50 (◊) and X = 1 ().

An important objective of this present study is to characterize the shear zone profiles of granular materials for ensembles of soft slippery and hard frictional materials. Fig. 6(b) shows the shear band profiles for samples with X = 0, 0.50 and 1. As observed from the figure, the width of shear band corresponding to a mixture sample X = 0.50 is wider than that of the pure mixtures, corresponding to X = 0 and 1. Also, the position of the shear band corresponding to a mixture sample X = 0.50 is shifted closer to the center than the pure mixtures. Note that the profile for the scaled width W/W_top as a function of z/H shown in Fig. 6(c) does not follow a quadrant shape, as described by the analysis of Ries et al.,³⁶ for any of the samples. The reason behind this deviation needs further exploration. However, a possible reason for this deviation can be due to the low W/d and z/d ratios, where d is the mean particle diameter. The normalized position of the shear band for the different mixtures as a function of bulk height is shown in Fig. 6(d), and the trends agree with that shown by Unger et al.¹⁰ From both Fig. 6(c) and (d), we observe a marked deviation in the trend for X = 0.50 compared to the pure samples with X = 0 and 1. Next, we analyze the width and position of the shear band for a horizontal slice through the axis of the shear cell at z = z_top = 16 mm, as represented by the dashed line in Fig. 6(b), for all samples of different fraction of hard particles X and show them in Fig. 6(e) and (f), respectively. Here, we find that the shear band gets wider as well as moves inwards in a mixture of soft and hard spheres up to a fraction X = 0.50, which corresponds to an intermediate mixture of 50% soft −50% hard particles. With a further increase of X, the shear band narrows and moves outwards in the shear cell. The non-monotonic trend of the width of the shear band with increasing X was also observed by Wang et al.¹¹

From the above observations, we conclude that the ratio of hard spheres in a mixture of hard frictional and soft slippery particles significantly affects the flow profiles at the bulk scale. For up to intermediate ratios of hard particles, the shear band gets wider by up to 8%, and also the center shifts inwards by up to 2%, compared to the values for pure soft particles. Both of these features are influenced by the effect of inter-particle friction and stiffness, which are dominant for frictional hard particles.⁸

Wang et al.¹¹ observed experimentally a notable difference in the shear band properties between samples made of pure materials and mixtures. It was observed that the shear zones of the mixtures are wider than those of the pure materials. Unlike the experimental observations, which show a significant widening of the shear zone for mixtures, our simulations show a rather small effect of widening the shear zone. However, the distinct non-monotonic trend of the shear bandwidth as a function of the mixing ratio is the feature observed in experiments, which is also captured in our simulations. In addition, we observe the shift in shear band position for the mixtures, which were not reported in the experimental work. It is important to note that the capillary bridges resulting from small amounts of liquid between hydrogel particles, and adhesive forces in general, potentially have an impact on experiments. Singh et al.⁷ showed that interparticle cohesion strongly affects the flow profiles, specifically shear banding where the width of the shear band increases with cohesion, and the internal force structures of granular materials. Typical velocity profiles emerging during shear are strongly affected by capillary bridge cohesion.⁴¹ Yang et al.⁴² experimentally studied the effects of cohesion on normal force distribution, showing a wider force distribution for higher cohesion between particles, which, in turn, influences the flow behavior. However, the effects of cohesive interactions between particles are not taken into consideration in this numerical study. For larger quantities of liquid in the mixtures, beyond the capillary regime, the materials exhibit a more dispersed shear behavior, which is also not considered in this numerical study.

6 Simplified cases to explore the role of inter-particle friction and elastic modulus

Following the findings of this study and the qualitative agreement of our results with the observations noted by Wang et al.,¹¹ it becomes evident that mixtures of soft, slippery and hard frictional grains demonstrate shear bands of variable width for increasing mixing ratio, which evolves non-monotonically with increasing values of the ratio. Since the materials of interest have two independent material parameters differing distinctly, namely friction and elastic modulus, it is interesting to look into their individual effects on the shear band characteristics. Hence, to explore the role of inter-particle friction and elastic modulus on the flow profile of the granular mixture systematically and independently of each other, we vary one parameter at a time in additional discrete element simulations of a parametric nature, while keeping the other parameter constant. Therefore, we conducted four case studies in addition to the base simulations, as outlined in Table 3. In these studies, our primary focus is on analyzing the shear band characteristics at a specific height, denoted as z = z_top = 16 mm in Fig. 6(b) as a function of mixing ratio X. We then proceed to compare the obtained results with those from the real simulations.

Table 3 Base simulations and simplified cases to explore the effects of friction and elastic modulus on shear band characteristics separately

Cases	Elastic modulus (E)	Friction coefficient (μ)
Base simulations	E soft = 50 kPa|Ehard = 5000 kPa	μ soft = 0.005|μhard = 0.20
E soft μ base	E = 50 kPa for both species	μ soft = 0.005|μhard = 0.20
E hard μ base	E = 5000 kPa for both species	μ soft = 0.005|μhard = 0.20
E base μ soft	E soft = 50 kPa|Ehard = 5000 kPa	μ = 0.005 for both species
E base μ hard	E soft = 50 kPa|Ehard = 5000 kPa	μ = 0.20 for both species

---

6.1 Shear band characteristics of granular mixtures

Luding⁸ examined the effects of friction on shear bands via simulations of a split-bottom shear cell. He showed that the presence of inter-particle friction does not affect the qualitative behavior of the shear band. However, the shear band moves inwards and gets narrower in the presence of friction. This is justified by the fact that with increasing friction, there is less relative displacement between the particles, particles tend to move more as a rigid block, and hence, the shear bandwidth decreases. Here, we vary the inter-particle friction of the mixture by varying the fraction of hard particles X, keeping the elastic modulus constant in cases E_softμ_base and E_hardμ_base as implied by the red and blue curves, respectively, in Fig. 7. For both the cases, we observe that the shear bandwidth decreases with increasing X in Fig. 7(a). Also, the shear band position moves inwards with increasing X as shown in Fig. 7(b). Our results agree with the findings of Luding on the effects of friction on shear band characteristics.

Fig. 7 Shear band (a) width W and (b) position R_c for fixed elastic modulus and real inter-particle friction, represented by cases E_softμ_base and E_hardμ_base, as compared with the real simulations, and shear band (c) width W and (d) position R_c for real elastic modulus and fixed inter-particle friction, represented by cases E_baseμ_soft and E_baseμ_hard, as compared with the real simulations, where all the features are plotted as a function of fraction of hard particles X. The width and position of the shear band are obtained at z_top = 16 mm. The error bars represent the fitting error, determined by fitting the numerical data for angular velocity to eqn (8).

On the contrary, varying the elastic modulus and considering same friction coefficient for the two materials leads to a wider shear band formation with increasing fraction of hard particles X, as demonstrated by cases E_baseμ_soft and E_baseμ_hard. This is implied by the green and orange curves in Fig. 7(c). However, for the same friction between the two species, varying only the elastic modulus of the species has only a small effect on the position of the shear band as a function of X as shown by the green and orange trends in Fig. 7(d).

Overall, there is a dual effect of friction and elastic modulus that influences the shear band properties. The width of the shear band decreases with increasing X for E_softμ_base simulations. In contrast, the width of the shear band increases with increasing X for E_baseμ_soft and E_baseμ_hard. In case of the base simulations where friction and elastic modulus of both the species are different, for smaller values of X, the friction of soft particles is smaller and hence is less dominant, i.e., it is easier for inter-particle contacts to slide. With increasing X, there are more contacts between hard frictional particles, and hence friction dominates over the effect of elastic modulus.

6.2 Macroscopic friction coefficient of granular mixtures

The local shear stress τ(r, z) plotted against the local normal stress P(r, z) is shown in Fig. 8(a) for different X. Since the system is inhomogeneous in nature, for a given pressure, we observe a wide range of local strain rate values [small gamma, Greek, dot above], and we find that shear stress τ is higher for increasing strain rate [small gamma, Greek, dot above]. However, for values above the critical strain rate, i.e., [small gamma, Greek, dot above] > [small gamma, Greek, dot above]_c, τ becomes almost independent of the local strain. This means that τ/p is nearly constant for all data points with strain rate [small gamma, Greek, dot above] > [small gamma, Greek, dot above]_c. A linear trend is observed for the shear stress as a function of the normal stress, which is fitted well by:

τ = μP, (10)

where μ is the macroscopic friction coefficient obtained from linear fitting of the local shear stress and normal stress values.

Fig. 8 (a) Shear stress τ plotted against normal stress P for granular mixtures of different fractions of hard particles X for the base simulations. The different opacity and size of markers are according to the magnitude of local strain rate, with opacity changing from high to low and size changing from small to large for increasing magnitude of local strain rate. The lines represent linear regression fits given by τ = μP. (b) Macroscopic friction coefficient μ as a function of the fraction of hard particles X for the base simulations and the simplified cases. The solid lines indicate linear and quadratic regression fits for the different cases. The error bars in the graph illustrate the fitting error, which was determined by fitting the numerical data of μ vs. P to the equation presented in eqn (10).

We further analyze the macroscopic friction coefficient μ for different mixtures as a function of X for the real simulations and the simplified cases and show the results in Fig. 8(b). The macroscopic friction coefficient μ increases as a quadratic function with X for the simplified cases of E_softμ_base and E_hardμ_base, i.e., same elastic modulus for both materials and varying inter-particle friction coefficient. On the contrary, μ increases linearly with X for cases E_baseμ_soft and E_baseμ_hard, i.e. varying elastic modulus and same inter-particle friction coefficient for both materials. Furthermore, the macroscopic friction coefficient μ for the real simulations also increases as a quadratic function with X, which suggests that μ for the real simulations is strongly dominated by the varying inter-particle friction coefficient between the two species. The elastic modulus does not play a significant role in influencing the μ of the bulk behavior. Focusing on the functional form of the quadratic function μ vs. X for the real simulations, the black solid line represents the quadratic fitting equation μ = μ_soft + aX², where μ_soft ≈ 0.14 is the macroscopic friction coefficient corresponding to the pure soft particle mixtures and a = 0.18 here depends on the inter-particle friction coefficient values of pure soft and hard species.

6.3 Packing density of granular mixtures

The local packing density ϕ plotted against the stress ratio τ/P is shown in Fig. 9 for different X. The data points correspond to observations within the shear band, where the strain rate [small gamma, Greek, dot above] exceeds the critical value [small gamma, Greek, dot above]_c. We observe a dilation effect, i.e., a decrease in the packing density with increasing τ/P from approximately ϕ = 0.70 for X = 0 to ϕ = 0.58 for X = 1 and all the data for different X collapse to a singular trend. We further analyze the packing density ϕ for different mixtures as a function of the stress ratio τ/P for the simplified cases and show the results in Fig. 10(a)–(d). For all the cases, we observe that the packing density decreases as a function of the stress ratio τ/P, indicating that both the friction coefficient and elastic modulus of particles contribute to dilatancy, although a weak dilatancy effect is observed for the cases E_baseμ_soft and E_baseμ_hard as compared to that of E_softμ_base and E_hardμ_base. Thus, the decrease in packing density with increasing stress ratio τ/P observed in Fig. 9 for base simulations is a collective contribution from the increasing friction coefficient and elastic modulus of particles.

Fig. 9 Packing density ϕ plotted against stress ratio τ/P for granular mixtures of different fractions of hard particles X for the base simulations.

Fig. 10 Packing density ϕ as a function of the fraction of hard particles X for the simplified cases (a) E_softμ_base (b) E_hardμ_base (c) E_baseμ_soft and (d) E_baseμ_hard.

7 Conclusions

This work investigated the shear zone characteristics in binary mixtures of soft low frictional and hard frictional particles, via discrete element simulations of mixtures inside a split-bottom shear cell. The ratio of the volume of hard particles to the soft particles was varied parametrically to observe its effects on measures of interest, including the shear bandwidth W, shear band position R_c and the macroscopic friction coefficient μ.

By comparing the width and position of the shear zone for samples of pure soft low frictional and for samples of pure hard frictional particles, it becomes evident that both systems show somewhat similar shear band features. However, mixtures of the two materials exhibit broader shear bands and a shifting of the shear band position towards the center of the shear cell. Both these effects demonstrate non-monotonicity for increasing fraction of hard particles. Although the changes in shear band characteristics for different mixtures are small, the non-monotonic trend of the shear band features as a function of mixing ratio is an interesting effect of the shear behavior of binary mixtures of this nature. We attribute the effects of shear band widening in these mixtures to an interplay of the role of inter-particle friction and particle stiffness.

Contrary to the shear band characteristics, the macroscopic friction coefficient shows a monotonically increasing trend as a quadratic function with an increasing fraction of hard particles in the mixtures. The macroscopic friction coefficient is strongly influenced by the inter-particle friction coefficient of the species present in the mixture, and the elastic modulus seems to play a less significant role.

The packing density of the granular mixtures exhibits a consistent decrease as the stress ratio τ/P increases. Notably, the data from different X values collapse into a singular trend. Both the particle friction coefficient and the elastic modulus collectively contribute to the dilatancy effect of the granular mixtures.

This work captures critical features of shear band formation in soft-hard particle mixtures of low-high inter-particle friction, such as the non-monotonic evolution of the shear bandwidth and position for increasing mixing ratios. However, the sensitivity of these features in our simulations is minor as compared to that observed in the experiments by Wang et al.¹¹ Unlike the simulation results showing a curved inward shear zone in the velocity profile in the rz-plane, the experimental velocity profile exhibits a predominantly vertical shear zone. Additionally, the velocity profiles from simulations reveal no significant differences among X = 0, 0.50, and 1. From the experimental findings, it is evident in the velocity profile that the shear zone significantly widens by over two times near the free surface for 50 : 50 particle mixtures compared to both 100 : 0 (soft : hard) and 0 : 100 (soft : hard) particle mixtures. However, a quantitative analysis is lacking in the experiments, making it challenging to provide an accurate estimation of the shear zone width from the velocity profile in the rz-plane. Some of the factors that might play a role in the experiments but are not captured in our simulations include the following: (i) our simulations only consider spherical particles and do not account for non-spherical shapes or the actual particle size distribution which can significantly affect the packing density and flow behaviour of granular materials,^43,44 (ii) the real stiffness of the mustard seeds is not taken into account; a sufficiently stiffer material is chosen for the mustard seeds instead, (iii) in experiments, shearing is controlled through angular displacements, whereas in simulations, we achieve shearing by continuously rotating the disk at a constant angular velocity, and (iv) capillary bridges due to small amounts of liquid between hydrogel particles and in general adhesive forces might play a role in experiments which are not considered here. The material would exhibit a more dispersed shear behavior for larger quantities of water in the mixtures beyond the pendular regime. Even with these discrepancies, our work captures the shear band behavior observed experimentally and the evolution of its characteristics for varying mixing ratios of soft and hard particles.

Author contributions

SR and VA designed the research and planned the parametric DEM simulations which APS carried out. SR and VA wrote the first draft of the manuscript, and TP reviewed and wrote with SR and VA the second draft. All authors reviewed the results and approved the final version of the manuscript. SR secured the research funding.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We acknowledge Stefan Luding for providing the code for micro–macro data transition and carefully reviewing this work. We thank Joshua Dijksman for stimulating discussion on this work. Financial support through the Emerging Talents Initiative (ETI Project no. ETI-Antrag_2022-1_Tech_13_Roy) and the DFG grant for Granular Weissenberg Effect (Project no. PO472/40-1) are acknowledged by SR and TP, respectively.

Notes and references

1. G. Wortel, T. Börzsönyi, E. Somfai, S. Wegner, B. Szabó, R. Stannarius and M. van Hecke, Soft Matter, 2015, 11, 2570–2576.
2. D. Fischer, T. Börzsönyi, D. S. Nasato, T. Pöschel and R. Stannarius, New J. Phys., 2016, 18, 113006.
3. D. Fenistein and M. Van Hecke, Nature, 2003, 425, 256.
4. D. Fenistein, J. W. van de Meent and M. van Hecke, Phys. Rev. Lett., 2004, 92, 094301.
5. J. A. Dijksman and M. Van Hecke, Soft Matter, 2010, 6, 2901–2907.
6. D. L. Henann and K. Kamrin, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 6730–6735.
7. A. Singh, V. Magnanimo, K. Saitoh and S. Luding, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 022202.
8. S. Luding, Part. Sci. Technol., 2007, 26, 33–42.
9. D. Faroux, K. Washino, T. Tsuji and T. Tanaka, Phys. Rev. Fluids, 2022, 7, 084306.
10. T. Unger, J. Török, J. Kertész and D. E. Wolf, Phys. Rev. Lett., 2004, 92, 214301.
11. J. Wang, Z. Farmani, J. A. Dijksman, C. Lübeck, O. Speck and R. Stannarius, Granular Matter, 2022, 24, 103.
12. A. Ashour, T. Trittel, T. Börzsönyi and R. Stannarius, Phys. Rev. Fluids, 2017, 2, 123302.
13. X. Hong, M. Kohne, M. Morrell, H. Wang and E. R. Weeks, Phys. Rev. E, 2017, 96, 062605.
14. K. Harth, J. Wang, T. Börzsönyi and R. Stannarius, Soft Matter, 2020, 16, 8013–8023.
15. R. O. Uñac, A. M. Vidales, O. Benegas and I. Ippolito, Powder Technol., 2012, 225, 214–220.
16. R. E. Rudge, J. P. Van De Sande, J. A. Dijksman and E. Scholten, Soft Matter, 2020, 16, 3821–3831.
17. J. Wang, B. Fan, T. Pongó, K. Harth, T. Trittel, R. Stannarius, M. Illig, T. Börzsönyi and R. C. Hidalgo, Soft Matter, 2021, 17, 4282–4295.
18. W. A. Beverloo, H. A. Leniger and J. Van de Velde, Chem. Eng. Sci., 1961, 15, 260–269.
19. K. Saitoh, V. Magnanimo and S. Luding, Soft Matter, 2015, 11, 1253–1258.
20. S. Luding, K. Taghizadeh, C. Cheng and L. Kondic, Soft Matter, 2022, 18, 1868–1884.
21. H. Götz, A. Santarossa, A. Sack, T. Pöschel and P. Müller, Granular Matter, 2022, 24, 1–9.
22. T. Weinhart, L. Orefice, M. Post, M. P. van Schrojenstein Lantman, I. F. Denissen, D. R. Tunuguntla, J. Tsang, H. Cheng, M. Y. Shaheen and H. Shi, et al. , Comput. Phys. Commun., 2020, 249, 107129.
23. K. Sakaie, D. Fenistein, T. J. Carroll, M. van Hecke and P. Umbanhowar, Europhys. Lett., 2008, 84, 38001.
24. S. S. Rohokale, Y. D. Dhanorkar, V. Pahuja and G. T. Kulkarni, J. Chronother. Drug Delivery, 2012, 3, 41–54.
25. P. Grewal and A. Singh, Int. J. Eng. Dev. Res., 2016, 4, 464–470.
26. V. Singh, S. C. Moses, R. N. Aalam and P. M. D'souza, Int. J. Environ. Clim. Change, 2023, 13, 1392–1399.
27. E. Ropelewska, K. J. Jankowski, P. Zapotoczny and B. Bogucka, Zemdirbyste-Agriculture, 2018, 105, 257–264.
28. N. V. Brilliantov, F. Spahn, J.-M. Hertzsch and T. Pöschel, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 53, 5382.
29. P. Müller and T. Pöschel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 021302.
30. R. D. Mindlin, J. Appl. Mech., 1949, 16, 259–268.
31. E. J. Parteli and T. Pöschel, Powder Technol., 2016, 288, 96–102.
32. J. M. N. T. Gray, P. Gajjar and P. Kokelaar, C. R. Phys., 2015, 16, 73–85.
33. J. W. Vallance and S. B. Savage, IUTAM Symposium on Segregation in Granular Flows: Proceedings of the IUTAM Symposium held in Cape May, NJ, USA June 5-10, 1999, 2000, pp. 31–51.
34. M. Lätzel, S. Luding and H. J. Herrmann, Granular Matter, 2000, 2, 123–135.
35. A. Singh, V. Magnanimo and S. Luding, New J. Phys., 2015, 17, 043028.
36. A. Ries, D. E. Wolf and T. Unger, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 051301.
37. D. M. Wood, Soil behaviour and critical state soil mechanics, Cambridge University Press, 1990.
38. S. Luding, Particuology, 2008, 6, 501–505.
39. B. Szabó, J. Török, E. Somfai, S. Wegner, R. Stannarius, A. Böse, G. Rose, F. Angenstein and T. Börzsönyi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 032205.
40. M. Depken, W. van Saarloos and M. van Hecke, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 73, 031302.
41. R. Mani, D. Kadau and H. J. Herrmann, Granular Matter, 2013, 15, 447–454.
42. R. Yang, R. Zou, A. Yu and S. Choi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 031302.
43. D. B. Nagy, P. Claudin, T. Börzsönyi and E. Somfai, Phys. Rev. E, 2017, 96, 062903.
44. D. B. Nagy, P. Claudin, T. Börzsönyi and E. Somfai, New J. Phys., 2020, 22, 073008.

---