Rheology in dense assemblies of spherocylinders: frictional vs. frictionless

TL;DR

This study uses molecular dynamics simulations to compare the rheological behavior of dense assemblies of frictional and frictionless spherocylinders, revealing the impact of friction on shear-thickening, particle alignment, and jamming density.

Contribution

It provides new insights into how friction influences the rheology and jamming of anisotropic particles, highlighting differences from frictionless systems.

Findings

Frictional particles exhibit shear-thickening behavior.

Friction suppresses particle alignment along flow.

Jamming density in frictional systems is insensitive to aspect ratio.

Abstract

Using molecular dynamics simulations, we study the steady shear flow of dense assemblies of anisotropic spherocylindrical particles of varying aspect ratios. Comparing frictionless and frictional particles we discuss the specific role of frictional inter-particle forces for the rheological properties of the system. In the frictional system we evidence a shear-thickening regime, similar to that for spherical particles. Furthermore, friction suppresses alignment of the spherocylinders along the flow direction. Finally, the jamming density in frictional systems is rather insensitive to variations in aspect-ratio, quite contrary to what is known from frictionless systems.

Figures

Equations

$${\mathbf{f}}_{ij}^n$$

$${\mathbf{f}}_{ij}^t$$

$$I_i{\dot{\omega }}_i=\sum _j{\mathbf{a}}_{ij}\times {\mathbf{f}}_{ij}$$

$$\sigma ={\sigma }_y+c{\dot{\gamma }}^x,$$

$$\sigma (\varphi ,\dot{\gamma })=|\delta \varphi {|}^{a}F(\dot{\gamma }\tau )$$

$$T_{\mu\nu}=\Bigl{\langle}\frac{3}{2N}\sum_{i=1}^{N}\Bigl{[}\hat{\text{l}}_{i\mu}\hat{\text{l}}_{i\nu}-\frac{\delta_{\mu\nu}}{3}]\Bigr{\rangle}$$

$$s_i=\frac{3{\hat{\text{l}}}_{ix}{\hat{\text{l}}}_{ix}-1}{2}$$

Full text

I Introduction

Flow and arrest in dense dispersions and granular systems is an important current theme in both materials science and fundamental research. Computer simulations have been developed to better understand the very broad range of observed rheological phenomena. One such phenomenon is discontinuous shear thickening, which signals the discontinuous increase of flow resistance (viscosity) when the forcing is increased by only small amounts. Such a behavior is highly relevant for many industrial operations, like the mixing or pumping of cementitious paste or concrete Feys et al. (2009). Recent work hints at a mechanism that relies on the direct particle frictional interactions upon contact Brown and Jaeger (2014); Heussinger (2013); Seto et al. (2013); Wyart and Cates (2014); Comtet et al. (2017); Clavaud et al. (2017); Hsiao et al. (2017).

Given the relevance of particle rotations for the frictional forces and the coupling to translational motion Maiti et al. (2016) we expect particle shape to be a determining factor for the rheological behavior, in particular for shear-thickening. Thus, in this work we study the flow of dense assemblies of frictional particles with non-spherical shape. Particle shape has been investigated in connection with the jamming transition in several publications Williams and Philipse (2003); Hoy (2017); Marschall and Teitel (2018); VanderWerf et al. (2018).

In granular systems numerous parameters such as aspect ratio, convexity, angularity play important roles to design transport, conduction, diffusive properties of a material Börzsönyi and Stannarius (2013). One of the fundamental observation in systems of non-spherical elongated particles is the orientation in shear flow, the angle between flow direction and average orientation of the particles being nonzero Campbell (2011); Guo et al. (2012, 2013); Börzsönyi et al. (2012a); Boton et al. (2013). Alignment decreases with aspect ratio of the particle, but has zero or no dependence on strain rate Börzsönyi et al. (2012a, b). The nematic order parameter is either monotonic Reddy et al. (2009); Nagy et al. (2017), or nonmonotonic function of density depending on the particle shape or aspect ratio Marschall et al. (2018). Here, we study the interplay of rheological properties and shear-induced alignment for model spherocylindrical particles. With an interest in the shear-thickening response we aim at comparing frictionless and frictional particles.

II Model

We perform computer simlations on a 3D system of spherocylindrical particles using interaction forces developed e.g. in Ref. Pournin et al. (2005). A spherocylinder in 3D consists of a cylindrical part of length $L$ and two hemispherical caps of diameter $D$ at two ends (see Fig. 2). The spherocylinders are parametrized by their length-to-diameter (aspect) ratio $\mathrm{AR}=L/D$. Spheres correspond to $\mathrm{AR}=0$.

The system has equal numbers of big and small particles with the same aspect ratio $\mathrm{AR}$ and a diameter ratio of 1.4. We have simulated systems of particles with $\mathrm{AR}$=0.001, 0.1, 0.5, 1.0, 1.5 and 2.0, where the particles with $\mathrm{AR}$=0.001 could be considered as spherical particles.

Homogeneous shear flow is implemented with periodic boundary conditions in the flow-direction ($\hat{x}$) and in the direction of average vorticity ($\hat{z}$), and Lees-Edwards boundary conditions in the direction of the velocity-gradient ($\hat{y}$) Lees and Edwards (1972). The volume fraction $\phi$, hence the box size is kept constant throughout the simulation.

The particles interact via repulsive, finite-range contact forces $\mathbf{f}$. A contact is made between particles $i$ and $j$, when the shortest distance between the spines of the particles, $r_{ij}$ is less than $D_{ij}$, where $D_{ij}=(D_{i}+D_{j})/2$. The normal $\mathbf{f}^{n}_{ij}$ and tangential $\mathbf{f}^{t}_{ij}$ components of $\mathbf{f}_{ij}$, force on particle $i$ due to contact between $i$ and $j$, are calculated by a force-overlap relation as in a linear-viscoelastic model,

[TABLE]

Here, the normal direction $\mathbf{n}_{ij}$ is a unit vector pointing from particle $j$ to $i$, $k_{n}$ and $k_{t}$ are spring constants, $c_{n}$ and $c_{t}$ are viscous damping constants, $\boldsymbol{\xi}^{t}_{ij}$ is elastic shear displacement pointing tangential to the contacting particles, and $\mathbf{v}^{n}_{ij}$ and $\mathbf{v}^{t}_{ij}$ are relative velocity comnonents in normal and tangential direction. The relative velocities include translational and angular components at the point of contact. Particle $j$ experiences the same but opposite force $-\mathbf{f}_{ij}$. We have set $k_{n}$ to unity, $k_{t}/k_{n}$ as 2/7, $c_{t}$ as zero and $c_{n}$ is such that the coefficient of restitution equals 0.7. All particles have the same mass set at 1.0.

Solid sliding friction is taken into account replacing $\mathbf{f}^{t}$ by $\mathbf{f}^{t}=\mu(\frac{|\mathbf{f}^{n}|}{|\mathbf{f}^{t}|})\mathbf{f}^{t}$, whenever the Coulomb inequality $|\mathbf{f}^{t}|>\mu|\mathbf{f}^{t}|$ holds. We have considered $\mu$=0.0, 0.5 and 1.0.

The motion of each particle is goverened by the collective force gathered from all particles in contact. The elastic and dissipative forces give rise to torques on the particles, and rotation of particle $i$ is dictated by

[TABLE]

where $I_{i}$ is the moment of inertia of particle $i$ and $\mathbf{a}_{ij}$ is moment of arm from the center of particle $i$ to the point of contact with particle $j$.

We integrate the equations of motion on a GPU using a velocity Verlet algorithm for the translational degrees of freedom, and a Richardson-like iteration for the rotational degrees of freedom, which are represented as quaternions. Normalization of the quaternions is ensured by rescaling at each time step. The time-step of the integration is 0.01 and number of particles $N$ is 1000.

III Results

III.1 Rheology for aspect ratio $\mathrm{AR}=1.0$

We start the discussion with the rheological properties of spherocylinders with a given aspect-ratio $\mathrm{AR}=1.0$. In Fig. 3 we present the flowcurves, i.e. the relation between strainrate $\dot{\gamma}$ of the shear flow and the stress $\sigma$ necessary to drive the flow. In the figure, we compare the rheology of a frictional system ($\mu=1.0$, left) with that of a frictionless system (i.e. $\mu=0$, center).

These flowcurves look rather similar to the much studied system of spherical particles Grob et al. (2016, 2014); Otsuki and Hayakawa (2011). The frictional system has three distinctive regimes. At small volume-fractions $\phi$ and at small strainrates $\dot{\gamma}$ the stress shows Bagnold scaling, $\sigma\propto\dot{\gamma}^{2}$, characteristic of inertial dynamics of hard particles. At high volume-fractions, beyond a jamming threshold, the flow follows the typical Herschel-Bulkley form (HB),

[TABLE]

with a yield stress $\sigma_{y}$ and an HB exponent $x$. Delimited by these two asymptotic regimes we find an additional shear-thickening regime (ST), which is a consequence of the frictional interactions. ST implies a stronger increase of the stress with strainrate than given by Bagnold scaling. In agreement with previous work Heussinger (2013); Seto et al. (2013); Wyart and Cates (2014) it seems that the onset of ST is not given by a characteristic strainrate but by a stress $\sigma_{0}\approx 10^{-2}\sigma_{y}$. The physical origin of this value is unclear at present.

The frictionless system does not have such a ST regime. Instead it displays a continuous crossover-scenario from Bagnold to shear-thinning behavior below a critical volume fraction $\phi_{c}$, and from yield-stress to the same shear-thinning behavior above $\phi_{c}$. Such a scenario corresponds to a continuous transition and is well known from spherical-particle systems Olsson and Teitel (2007). On the other hand, the properties of the frictional system are characteristic of a discontinuous transition like, for example, a liquid-gas coexistence. For the continuous transition, we can use a scaling Ansatz for the flowcurves

[TABLE]

with a time-scale $\tau=|\delta\phi|^{-b}$, $\delta\phi=\phi-\phi_{c}$, and approximately determine the exponents $a$, $b$ as well as $\phi_{c}$ (see Fig. 3 right). The obtained values imply, e.g. for the viscosity $\eta=\sigma/\dot{\gamma}^{2}\sim|\delta\phi|^{-\mu}$ with $-\mu=a-2b\approx-3$. Likewise, the HB exponent defined in Eq. (4) is given by $x=a/b\approx 0.5$. However, care must be taken as the values of the exponents sensitively depend on the value for $\phi_{c}$. Furthermore, our systems are rather small, and finite-size effects are likely to strongly influence these values.

III.2 Shear-induced alignment

It is to be expected that the spherocylindrical particles align with the flow. To quantify this effect we compute the average orientational ordering tensor Börzsönyi et al. (2012c),

[TABLE]

where $\mathbf{\hat{\text{l}}}_{i}$ is the unit vector pointing along the long axis of particle $i$. The quantity $T$ is a $3\times 3$ matrix, the largest eigenvalue of which is called the global nematic order parameter $\mathbf{S2}$. The local nematic order parameter is given by,

[TABLE]

which is the $(x,x)^{th}$ entry of the local orientational tensor of particle $i$, and represents a measure for the alignment with the flow direction. Snaphshots of the system highlighting the local ordering are given in Fig. 1.

The global nematic order parameter is given in Fig. 4, again comparing frictional and frictionless systems. For the frictional system, in the fluid regime (small $\phi$ and $\dot{\gamma}$) the nematic order is rather high, at least at the $\phi$ values close to jamming that we consider. In the HB regime (high $\phi$), on the other hand, $\mathbf{S2}$ is small, possibly because of plastic events leading to efficient randomization. In the ST regime $\mathbf{S2}$ rapidly interpolates between the high value in the fluid and the small value in HB. The strainrate at which this transition happens corresponds to the onset of ST in Fig. 3 and represents the characteristic stress $\sigma_{0}$.

If we now turn to the frictionless system, then the first observation is that the scale of nematic order is strongly enhanced as compared to the frictional system. Apparently, friction suppresses or frustrates flow alignment. The transition from fluid to jammed state is also different than for frictional systems: as in the flow-curves it features a continuous crossover scenario from the fluid or the HB branch into a critical branch at a volume-fraction dependent strainrate. The strainrate-independent value of S2 in the fluid branch therefore continuously decreases upon approaching the jamming limit. In contrast, the frictional scenario is discontinuous. Two branches are connected via a rapid crossover at the onset stress for ST. Similar data for frictionless ellipsoids in 3d has been presented in Ref.  Marschall et al. (2018).

III.3 Distribution of alignment angles

For ellipsoidal and cylindrical particles with small friction the typical angle between the ellipsoid major axis and the flow has been shown to be non-zero Campbell (2011); Guo et al. (2012, 2013). Similar results have been obtained for dumbbell-shaped particles Reddy et al. (2009) and in experiments of shear flows of elongated particles in a cylindrical split bottom shear cell Börzsönyi et al. (2012a).

We display the distribution of orientation angles in the fluid regime in Fig. 5. We define $\beta_{\rm xy}$ as the angle between the flow direction and the projection of the particle orientation into the shear-gradient (xy) plane. For the frictionless system we observe a clear peak at a finite angle which, for the frictional system, is strongly degraded.

The origin of the peak can be explained with the thermal tumbling dynamics of a single particle, which spends a rather long time at small positive angles Hinch and Leal (1972); Lee and Kim (2009); Gerashchenko and Steinberg (2006). A negative angle, on the other hand, which may be achieved by a thermal fluctuation, is unstable and immediately leads to full rotation of the particle back into a stable positive angle. The smaller the temperature (the larger the Peclet number) the more pronounced is this effect. In our system, the temperature may be due to the interaction with the surrounding particles. The frictional system then providing higher effective temperature as the frictionless system. A similar trend is observed when going into the plastic flow HB regime (not shown). There the peak height is reduced in agreement with the presumed higher randomization due to plastic events as mentioned above.

III.4 Contacts for short spherocylinders, $\mathrm{AR}=0.1$

Two recent publications Nagy et al. (2017); Marschall et al. (2018) describe an interesting phenomenon relating to the number of contacts between (frictionless) spherocylinders/ellipsoids during flow. Apparently, nearly-spherical particles tend to establish an excessive number of contacts at their sides. Notably, the area of the side is rather small as compared to the spherical cap region.

In Fig. 6 we plot the ratio $R_{\rm con}$ of side-to-end contacts comparing the frictionless with the frictional system. From the range of values it is readily observed that the number of side contacts in the frictionless simulations are increased as compared to the frictional case. The fraction of side contacts rises to about $25\%$ without friction and $~{}8\%$ with friction. The latter value can be understood from the ratio $r$ of surface areas on the side ($\pi DL$) and on the caps ($\pi D^{2}$), giving $r=\mathrm{AR}=0.1$. Thus, the number of end contacts in the frictional scenario is to be considered normal, while that in the frictionless scenario is markedly enhanced.

We also set up additional simulations for the case of isotropic jamming without shear. To this end we run a compression-decompression cycle with $\phi_{\rm min}=0.64$ and $\phi_{\rm max}=0.74$. From the decompression branch the jamming transition can be located at $\phi_{c}\approx 0.68\ldots 0.685$ (the error stems from the step-size). At jamming we find roughly $5$ end contacts per particle, as well as $1.6$ side contacts. For all volume fractions above $\phi_{c}$ the ratio of side-to-end contacts is constantly about $32\%$ and therefore even higher than in our shearing simulations. With the strong increase towards lower strainrates, the latter high value may eventually be reached, however.

This strong increase visible in Fig. 6 is indeed remarkable, as it indicates the presence of a long (but necessarily finite) time scale, at which the contact numbers cross over into their zero-strainrate limit. This time-scale cannot be the same as the $\tau$ extracted from the scaling of the flowcurve, Eq. (5), as that diverges at the jamming transition. One may speculate that this time-scale is connected to small-scale rotational motion. In packings of nearly-spherical ellipsoids rotational modes occur at the lower end of the frequency spectrum Zeravcic et al. (2009); Mailman et al. (2009), thus possibly leading to a small time-scale visible under shear. Longer simulations at smaller strainrates are necessary to answer this question in detail.

III.5 Shear thickening

If we compare data from different aspect ratios $\mathrm{AR}$ we expect the amount of shear alignment to increase with increasing $\mathrm{AR}$ Börzsönyi et al. (2012a); Marschall et al. (2018). This may already be apparent from the snapshots in Fig. 1. What is the consequence for the rheological properties? At first sight stronger alignment should make flow easier. This is what we have observed in the fluid regime in Fig. 4, increasing alignment corresponds to decreasing stress. Note, however, that the opposite is true for the HB regime, where plastic events supposedly counteract alignment and lead to small $\mathbf{S2}$.

The rheology for aspect ratios $\mathrm{AR}\neq 1$ is not much different from Fig. 3. The relevant range of volume fractions changes with the aspect ratio. We chose to use the viscosity in the fluid regime to compare these data. In Fig. 7 we therefore compare effective viscosities $\eta\equiv\sigma/\dot{\gamma}^{2}$ for different $\mathrm{AR}$ but with similar low-strainrate limit $\eta(\dot{\gamma}\to 0)=\eta_{0}$. What is most apparent for these data is the effect of aspect ratio on thickening. While short particles display only a mild viscosity increase in the ST regime, long particles have a much stronger increase – all starting from the same small-strainrate viscosity. At the same time the ST regime spans to smaller stresses. Apparently, the onset stress for ST decreases with the length of the particles. While there is no general understanding of the parameter dependencies of this onset stress, a negative correlation with particle length at least is plausible, when one considers the effects of particle rotations. In Ref. Maiti et al. (2016) we have argued that increasing stress leads to the build-up of rotationally frustrated structures that may bear higher loads. As to the anisotropic shape of longer particles, such an effect is likely to be reinforced. Thus, smaller stresses are sufficient to lead to the same “strength” of ST.

In this context we also want to mention the two experiments Brown et al. (2011); Egres and Wagner (2005), which deal with the rheology of cylindrical particles in the range $\mathrm{AR}=2\ldots 9$. While in Ref. Brown et al. (2011) the onset stress is masked by a shear-thinning regime at small stress, Ref.Egres and Wagner (2005) seem to find the same onset stress for the different particle lengths (which are longer than ours).

III.6 Aspect-ratio dependence of $\phi_{C}$

It is well known Williams and Philipse (2003) that the jamming density $\phi_{c}$ of isotropic packings of spherocylinders first increases with $\mathrm{AR}$, then presents a maximum at $\mathrm{AR}\approx 0.5$ before declining asymptotically as $\phi_{c}\propto\mathrm{AR}^{-1}$ (see Fig. 8, data marked as “Williams”).

The latter can be understood with an excluded volume argument Semenov (1986); Philipse (1996) that relates the excluded volume of a section $L_{e}$ of a long cylinder $v_{\rm excl}\sim(\pi/2)L_{e}^{2}D$ to the available volume $v_{\rm avail}\sim 1/\rho z$, where $\rho$ is the number density of cylinders and $z=L/L_{e}$ is the number of sections of length $L_{e}$ per particle. In terms of volume fraction $\phi\sim\rho(\pi/4)LD^{2}$ one arrives at $\phi\sim z/2\mathrm{AR}$. The number of sections $z$ is equal to the number of contacts, which –at jamming– should be $z\approx 8$ for long frictionless cylinders. The jamming volume fraction thus is given by $\phi_{c}\approx 4/\mathrm{AR}$ Rodney et al. (2005). We do not probe this asymptotic regime here, but remain in the range of the maximum.

As can be seen in the figure, frictionless particles in isotropic packings (Williams et al. Williams and Philipse (2003)) display a rather strong increase of the jamming density when increasing the length of the particles from the sphere limit $\mathrm{AR}=0$. It has been argued Donev et al. (2007) that nearly-spherical particles can make more efficient use of the available space than spheres, because orientations may be optimized within a given local surrounding. In fact, the slope of an initial linear increase in $\phi_{c}(\mathrm{AR})$ can be well approximated by satisfying torque-balance constraints for each particle, given the local surrounding of the equivalent sphere packing Donev et al. (2007). This sets an “optimal” orientation of the particle which allows for calculating the jamming density. On the other hand, chosing random orientations for the particles one would obtain a vanishing slope, i.e. insensitivity of $\phi_{c}$ towards aspect-ratio Donev et al. (2007).

Under shear the local surrounding is constantly changing, thus additional rotations are necessary to fulfil the torque-balance constraints. Recent pressure-controlled shear simulations with frictionless particles Nagy et al. (2017) indicate that the strong increase for small $\mathrm{AR}$ is preserved, albeit at slightly elevated volume fractions $\phi$ (see Fig. 8 data marked as “Nagy”). We confirm these data within our framework of $\phi$-controlled simulations. Note, that a similar difference between the jamming density in isotropic systems and that under shear is well known from systems of spherical particles O’Hern et al. (2002); Heussinger and Barrat (2009). At larger $\mathrm{AR}$ the jamming density does not decrease as in the isotropic packings. As discussed above these systems display alignment. This induces correlations between particles such that the excluded volume argument is invalidated Philipse (1996).

With friction, the relevant volume fractions are strongly reduced, which is to be expected. What may be surprising, however, is that the non-monotonic behavior is much less pronounced and the initial slope is strongly decreased. At large enough friction coefficient $\mu$ the maximum even seems to disappear. A similar effect is observed in a two-dimensional flow of ellipses in Ref. Trulsson (2018). Given the role of particle rotations for the increasing $\phi_{c}$ for small $\mathrm{AR}$ as discussed above, it is tempting to ascribe this effect to the tangential frictional forces acting parallel to the surface of the particle. These very effectively resist particle rotation. Thus, additional packing optimization due to particle rotation, which seems very prominent in nearly-spherical frictionless particles, is absent in the presence of friction. As remarked in Ref. Donev et al. (2007), insensitivity of the jamming density towards aspect-ratio could result from randomly orienting the particles irrespective of their local surrounding.

IV Summary

We have studied, by molecular dynamics simulations, the shear rheology of dense packings of soft spherocylinders, highlighting the differences between frictional and frictionless systems. We concentrate on the regime of short spherocylinders with aspect ratios $\mathrm{AR}=0.001\ldots 2.0$, close to the sphere limit. The flowcurves we obtain are rather similar to those of spherical particles. In particular, we find a shear-thickening regime (ST) in the frictional system, just as with spheres. The frictionless system, on the other hand, only shows shear-thinning behavior and no shear thickening. When comparing different particle aspect ratios we find that longer particles have a smaller onset stress for ST. Starting with the same small-stress viscosity, the maximal increase in viscosity can thus be modulated by particle shape. Friction also affects the strength of shear-induced alignment, with frictional systems showing a much smaller alignment than frictionless systems. Furthermore, frictionless systems display a rather large percentage of side-contacts between spherocylinders. It is much larger than in frictional systems, where the percentage roughly reflects the surface area of the side of the spherocylinder relative to its total area.

In order to rationalize these phenomena we highlight the importance of particle rotations, which follows from the nature of frictional forces to act tangential to particle surfaces. Frictional torques from the fluctuating dense environment may act as additional noise on the particle orientation, thus reducing overall alignment. We have discussed this question in terms of the probability distribution of orientation angle of the particles. This is much flatter for frictional particles than for frictionless particles. In addition frictional forces may be very effective in suppressing small-scale rotations that could optimize local packing arrangements. This could also be responsible for the observation that in frictional systems the jamming threshold is quite insensitive to aspect ratio, in contrast to what is known from frictionless systems.

Acknowledgements.

We acknowledge discussions with L. Pournin.