1. Introduction
Suspensions are one of the simplest complex fluids to make, consisting solely of solid particles suspended in a viscous fluid. Despite their simplicity, driven suspensions can show a multitude of fascinating and complex physics (Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018; Ness, Seto & Mari Reference Ness, Seto and Mari2022). For example, under deformation such as shear, they can jam, thicken or thin, depending on the particle–particle interactions and the solid content, i.e. volume fraction (Seto et al. Reference Seto, Singh, Chakraborty, Denn and Morris2019). Understanding when and how the suspension mechanical properties, such as viscosity and normal stress differences, change upon changing conditions is crucial for an effective processing of these materials, not least as they are widely used in industry (e.g. in food, cosmetic, personal care and drug processing).
For dense suspensions composed of rigid particles under excluded-volume interactions, with or without frictional forces, the steady-state viscosity is only a function of the (solid) volume fraction and shows a power-law divergence around the shear jamming fraction (Zarraga et al. Reference Zarraga, Hill and Leighton2000; Ovarlez, Bertrand & Rodts Reference Ovarlez, Bertrand and Rodts2006; Olsson & Teitel Reference Olsson and Teitel2007; Bonnoit et al. Reference Bonnoit, Darnige, Clement and Lindner2010; Andreotti, Barrat & Heussinger Reference Andreotti, Barrat and Heussinger2012; Trulsson, Andreotti & Claudin Reference Trulsson, Andreotti and Claudin2012). However, in the case of unsteady flow conditions, for example in oscillatory shear (OS) with only excluded-volume interactions, the viscosity is found to depend not just on the packing fraction but also on the strain amplitude (Breedveld et al. Reference Breedveld, van den Ende, Jongschaap and Mellema2001; Bricker & Butler Reference Bricker and Butler2006, Reference Bricker and Butler2007; Lin, Phan-Thien & Khoo Reference Lin, Phan-Thien and Khoo2013; Ness, Xing & Eiser Reference Ness, Xing and Eiser2017; Martone, Carotenuto & Minale Reference Martone, Carotenuto and Minale2020; Ge & Elfring Reference Ge and Elfring2022). The strain amplitude dependence shows a non-monotonic behaviour of the viscosity, with minima at an intermediate strain amplitude which depends on the volume fraction.
The particle dynamics also depends on the strain amplitude in OS. In
$1966$
, G. I. Taylor illustrated the reversibility feature of Stokes flow in a now-classic experiment: by reversing the direction of the flow, the dye particles suspended in a viscous Newtonian fluid returned to their initial positions. However, later studies have shown that neutrally buoyant non-Brownian particles in fact diffuse in a sheared Stokesian fluid (Eckstein, Bailey & Shapiro Reference Eckstein, Bailey and Shapiro1977; Leighton & Acrivos Reference Leighton and Acrivos1987; Davis Reference Davis1996; Marchioro & Acrivos Reference Marchioro and Acrivos2001; Sierou & Brady Reference Sierou and Brady2004). Such a paradox could be reconciled by the chaotic nature of the hydrodynamic interactions between the particles (Drazer et al. Reference Drazer, Joel, Boris and Andreas2002), or the presence of non-hydrodynamic forces that dominate at particle contacts (Metzger, Pham & Butler Reference Metzger, Pham and Butler2013; Pham, Metzger & Butler Reference Pham, Metzger and Butler2015). More recently, Pine et al. (Reference Pine, Gollub, Brady and Leshansky2005) and others (Corte et al. Reference Corté, Chaikin, Gollub and Pine2008; Menon & Ramaswamy Reference Menon and Ramaswamy2009; Metzger & Butler Reference Metzger and Butler2010; Metzger et al. Reference Metzger, Pham and Butler2013; Nagasawa, Miyazaki & Kawasaki Reference Nagasawa, Miyazaki and Kawasaki2019; Ness & Cates Reference Ness and Cates2020; Ge et al. Reference Ge, Martone, Brandt and Minale2021; Mari, Bertin & Nardini Reference Mari, Bertin and Nardini2022) demonstrated that non-Brownian suspensions subject to periodic shear flow can undergo a dynamical phase transition, called ‘reversible–irreversible transition’ (RIT). Remarkably, when the amplitude is low, particles self-organise at the end of each period to avoid further collisions. The phenomenon, termed random organisation, leads to reversible particle trajectories in idealised situations and was first discovered for dilute suspensions, but recently was also seen in the semi-dilute and the dense regimes. In many real cases, the dynamics does not become perfectly reversible but rather sub-diffusive, due to thermal or other noise sources in the system. Nevertheless, both perfect reversibility and sub-diffusivity are a signature of the onset of random organisation (Pine et al. Reference Pine, Gollub, Brady and Leshansky2005). On the other hand, above a critical strain amplitude, the particle motion becomes diffusive and irreversible, consistently with the chaotic dynamics observed in steady shear (SS). Furthermore, it was found that this critical strain amplitude (
$\gamma _{0,c}$
) decreases with the volume fraction (
$\phi$
) as
$\gamma _{0,c} \sim \phi ^{-2}$
(Pine et al. Reference Pine, Gollub, Brady and Leshansky2005), suggesting that diffusion results from interactions of any particle with at least two neighbours, including both hydrodynamic (lubrication) and non-hydrodynamic (contact) contributions.
The above studies clearly show that suspensions in OS display intricate rheological and dynamical behaviours, both of which depend on the volume fraction and driving strain amplitude. As these properties are ultimately determined by the suspension microstructure (assuming no friction or time-dependent interactions), it is reasonable to expect a relation between the suspension rheology and the dynamics. Indeed, Ge & Elfring (Reference Ge and Elfring2022) have recently shown that at the strain amplitude where the suspension viscosity is minimal, the suspension dynamics also transitions from reversible to diffusive. Furthermore, these authors found other rheological signatures, including an enhanced intracycle shear thinning and a finite second normal stress difference, at the onset of RIT. These observations are intriguing and suggestive of a unique relationship between the dynamics and rheology; however, we note that they were made in periodically sheared suspensions where the underlying flow is reversible. To determine whether there is indeed such a unique relationship requires exploring new flow protocols with a more complex dynamics.
Although most of the previous works have been devoted to OS flow, with inherent periodicity, time reversibility and sudden shear reversal, the characterisation of the mechanical response to changes in shear axes has seldom been considered. The only exceptions are perhaps the recent efforts of Blanc et al. (Reference Blanc, Peters, Gillissen, Cates, Bosio, Benarroche and Mari2023) and Acharya & Trulsson (Reference Acharya and Trulsson2024), who investigated the rheology (but not the dynamics) of dense suspensions after a sudden rotation of the shear axes. These authors observed a drop in shear viscosity after shear rotation on strain scales of the order of unity or less. The magnitude of the viscosity drop was found to depend on the rotation angle, from no drop for 0-degree rotation, since this case corresponds to no change, to a maximum drop at 180-degree rotation around the gradient axis, corresponding to a classical shear reversal experiment as previously studied by e.g. Peters et al. (Reference Peters, Ghigliotti, Gallier, Blanc, Lemaire and Lobry2016), with a gradual and smooth transition between these two extremes. Alternating shear rotations also reduce the average viscosity and dissipation, if a small strain is chosen between repeated rotations (Acharya & Trulsson Reference Acharya and Trulsson2024). Despite the novelty of the shear rotation protocol, the change in the shear direction always occurs suddenly. This raises an interesting question: What if, instead, the gradient shear axis is rotated smoothly and gradually?
In this work, we aim to fill this gap by investigating both the rheology and dynamics of dense suspensions beyond the sudden shear reversal and sudden shear rotation. To this end, we propose a novel shear protocol called rotary shear (RS), where we slowly change the shear direction in a continuous way. This new RS flow, which is a linear combination of two out-of-phase OSs perpendicular to each other, has the benefit of not relying on any shear reversals (i.e. a sudden change in sign/direction), has a constant shear magnitude and it is still periodic. Rotary shear finds its origin in many real-world flows that undergo multidirectional shear deformations. One everyday example is the rubbing motion while applying an ointment or cream to the skin: when applying skincare products people usually do not use linear strokes or oscillatory motion but a circular rubbing movement, which resembles our RS protocol. Another conceivable example, which may be experimentally realised in a Taylor–Couette rheometer, is to impose oscillations both in the circumferential and axial directions of the inner or outer cylinder. If the two oscillations are out of phase, the net motion will also resemble a RS flow.
We note that our new RS protocol is inherently irreversible in contrast to OS, as the fluid flow moves continuously circularly in a directional manner (i.e. clockwise or anticlockwise) without reversing its path. However, when viewed stroboscopically, both RS and OS appear to be reversible because both are periodic, prompting one to ask whether the suspension dynamics and rheology would share any similarity.
Therefore, we investigate the RS protocol in terms of both mechanical properties and microstructural changes, numerically comparing it with the extensively studied OS protocol to find the correspondence between the rheology and the dynamics. Ultimately, our findings will help to link the rheological and dynamical behaviours of a dense suspension through its microstructure.
The outline of the paper is as follows: we start by discussing our numerical model in § 2. In § 3 we show the RS’s suspension rheology, its dynamics and the microstructure in comparison with OS and try to understand the difference and similarity between the two protocols using a reversible combination of the two: reversible-RS (RRS). We summarise our main findings in § 4.
2. Modelling and method
Figure 1(a) shows a schematic view of the dense suspension under consideration. As a reference, we will expose the suspension to the OS protocol, which has been widely used in the literature to investigate the time-dependent rheology of particulate suspensions. In our coordinate system, we consider particles subject to a time-periodic shear rate in the
$y$
direction with
$z$
being the gradient direction,
$\dot {\gamma }_{yz} = \gamma _0 \omega \cos (\omega t)$
, as shown in figure 1(b). Where
$\gamma_0$
is strain amplitude,
$\omega$
is angular frequency of oscillation, and
$t$
is time. To investigate a more complex flow, we propose to add to the
$yz$
shear a second shear in the
$xz$
plane with a non-zero phase difference to the first shear.
If the second shear is lagging the first one by
$\pi /2$
in phase,
$\dot {\gamma }_{xz} = \gamma _0 \omega \sin (\omega t)$
, the direction of the effective applied shear would rotate in the clockwise direction (see figure 1
c), which we call RS. In practice, this might be realised, for example, in a cylindrical rheometer by oscillating the inner cylinder in both axial and circumferential directions. The external rate-of-strain tensors corresponding to the OS and RS protocols are
\begin{equation} {\mathbb{E}}^\infty _{\textit{OS}} = \left [ \begin{array}{ccc} 0 & \quad 0 & \ 0 \\ 0 & \quad 0 & \ {\dot {\gamma }_{yz}}/{2} \\ 0 & \quad {\dot {\gamma }_{zy}}/{2} & \ 0 \\ \end{array} \right ] , \quad {\mathbb{E}}^\infty _{\textit{RS}} = \left [ \begin{array}{ccc} 0 & \quad 0 & \ {\dot {\gamma }_{xz}}/{2} \\ 0 & \quad 0 & \ {\dot {\gamma }_{yz}}/{2} \\ {\dot {\gamma }_{zx}}/{2} & \quad {\dot {\gamma }_{zy}}/{2} & \ 0 \\ \end{array} \right ]\! . \end{equation}
Figure 1. (a) Schematic view of the dense suspension. Top view showing the directions of shear rate for (b) OS, (c) RS, and (d) RRS, the dashed red line shows shear direction just before the reversal and the solid red line shows the instantaneous shear direction. Here, we denote the vorticity, streamwise and velocity-gradient directions as
$x$
,
$y$
and
$z$
for OS, but the vorticity and streamwise directions rotate in the case of RS and RRS, whereas the gradient direction remains along
$z$
.
In the RS protocol, the suspension is acted upon by shear of constant magnitude but varying direction. A fluid parcel would return to its original position after a shear cycle in RS but it does not reverse its path like that which occurs in OS. This leads to a flow that is not time reversible. To nonetheless keep reversibility and also examine the effect of sudden shear reversal, we consider a rotary deformation where we reverse the flow direction and the direction of rotation of the shear, which we call RRS, see figure 1(d). Specifically, we apply
$\dot {\gamma }_{xz} \to - \dot {\gamma }_{xz}$
,
$\dot {\gamma }_{yz} \to - \dot {\gamma }_{yz}$
and
$\omega \to - \omega$
after a rotation of angle
$\theta _0$
or after a time
$t=\theta _0/\omega$
. For consistency and ease of comparison, we keep the same time period (
$T = 2\pi /\omega$
) in OS, RS and RRS such that the shear rate amplitude
$\dot {\gamma } = \gamma _0\omega$
remains a small constant. Thus, the time-dependent shear rate in the RRS protocol is expressed as
\begin{equation} \dot {\gamma }_{xz} = \gamma _0 \hat \omega \sin (\hat \omega t), \quad \dot {\gamma }_{yz} = \gamma _0 \hat \omega \cos (\hat \omega t), \end{equation}
where
$\hat \omega =\omega$
if
$nT \leqslant t\lt (n+1/2)T$
, and
$\hat \omega =-\omega$
if
$(n+1/2)T \leqslant t\lt (n+1)T$
, for
$n \in {\mathbb{N}}$
.
The motion of a rigid particle in a particle suspension is governed by the Newton–Euler equations
\begin{equation} \begin{array}{ll} \boldsymbol {F}_{\!i} = m_i \dfrac {{\textrm d} \boldsymbol {u}_i}{{\textrm d}t},\\[6pt] \boldsymbol {T}_{\!i} = \boldsymbol {I}_{\!i} \dfrac {{\textrm d}\boldsymbol {w}_i}{{\textrm d}t} + \boldsymbol {w}_i \times (\boldsymbol {I}_{\!i} \boldsymbol {w}_i),\end{array} \end{equation}
where
$\boldsymbol {F}_{\!i}$
and
$\boldsymbol {T}_{\!i}$
are, respectively, the total force and torque exerted on particle
$i$
of mass
$m_i$
and moment of inertia tensor
$\boldsymbol {I}_{\!i}$
, and
$\boldsymbol {u}_i$
and
$\boldsymbol {w}_i$
its translational and rotational velocities, respectively. For spheres, rotation has no effect on the translational dynamics of the suspension, thus we omit the associated degrees of freedom from now on. When the suspension is dense, long-range hydrodynamic interactions may also be negelcted (Seto et al. Reference Seto, Mari, Morris and Denn2013; Mari et al. Reference Mari, Seto, Morris and Denn2014), and the total force on a particle is simply a sum of hydrodynamic forces, such as the Stokes drag
$(\boldsymbol {F}_{\!i}^S)$
and pairwise, short-range lubrication
$(\boldsymbol {F}_{\!i}^L)$
, and any non-hydrodynamic forces, including contact forces
$(\boldsymbol {F}_{\!i}^C)$
and friction. Specifically, we write
\begin{equation} \boldsymbol {F}_{\!i} = \boldsymbol {F}_{\!i}^S + \sum _j^{N_l} \boldsymbol {F}^L_{\textit{ij}} + \sum _j^{N_c} \boldsymbol {F}^C_{\textit{ij}}, \end{equation}
where
$N_l$
and
$N_c$
are the numbers of lubricating and contact pairs. The expressions for the individual terms on the right-hand side of the above equation are taken from the hybrid lubrication/granular dynamics model derived for dense particle suspensions in low Reynolds number flow (Cheal & Ness Reference Cheal and Ness2018; Ge & Brandt Reference Ge and Brandt2020). Their detailed expressions are given in Appendix A.
Finally, we note that there are several model parameters in our system, including the particle radius
$a$
, density
$\rho$
, stiffness
$k_n$
, dynamic viscosity of the fluid
$\eta _0$
and shear rate
$\dot \gamma$
. For the model output to correspond physically to a dense suspension of inertialess, rigid particle, we require the Stokes number
$\textit{St} =\rho \dot {\gamma } a^2/\eta _0 \ll 1$
and a stiffness-scaled shear rate
$\hat {\dot {\gamma }} = \dot {\gamma }a_1/\sqrt {k_n/(\rho a_1)} \ll 1$
. Throughout this work, we have
$\textit{St} \approx 10^{-2}$
and
$\hat {\dot {\gamma }} \approx 10^{-4}$
. See Ge & Brandt (Reference Ge and Brandt2020) for more details.
2.1. Quantities of interest
To examine the rheology of the suspension, we consider the bulk stress tensor, defined as
\begin{equation} \sigma = 2\eta _0{\mathbb{E}}^\infty + \frac {1}{V} \left ( \sum _i^N {\mathbb{S}}_i^S + \sum _{j\gt i}^{N_l} {\mathbb{S}}_{\textit{ij}}^L + \sum _{j\gt i}^{N_c} {\mathbb{S}}_{\textit{ij}}^C \right )\!, \end{equation}
where
$V = L_x L_y L_z$
is the volume of the simulation box (cf. Figure 1
a). The single-body hydrodynamic stresslet (
${\mathbb{S}}^S$
), the pairwise lubrication stresslet (
${\mathbb{S}}^L$
) and the stress tensor due to contact force (
${\mathbb{S}}^C$
) are given as
\begin{align} {\mathbb{S}}_i^S &= \big (20\pi \eta _0 a_i^3/3 \big ){\mathbb{E}}^\infty , \end{align}
\begin{align} {\mathbb{S}}_{\textit{ij}}^L &= \boldsymbol {F}_{\textit{ij}}^L(\boldsymbol {r}_j- \boldsymbol {r}_i), \end{align}
\begin{align} {\mathbb{S}}_{\textit{ij}}^C &= \boldsymbol {F}_{\textit{ij}}^C(\boldsymbol {r}_j- \boldsymbol {r}_i), \end{align}
where,
$i$
and
$j$
are particle indices. In SS flow, the viscosity is simply
$\eta = \sigma _{yz}/\dot {\gamma }$
, where
$\sigma _{yz}$
is the shear stress component of the output stress tensor, and
$\dot {\gamma }$
is the shear rate. In the unsteady shear case, the concept of complex viscosity is adopted, which includes both viscous and elastic contributions. It is customary to express the complex viscosity
$\eta ^*$
as
\begin{equation} \eta ^* = \eta ' - i\eta '', \end{equation}
where
$\eta '$
denotes the dynamic viscosity, and
$\eta ''$
the non-dissipating elastic contribution (Mewis & Wagner Reference Mewis and Wagner2012). The complex viscosity,
$\eta ^*$
can be calculated using Fourier analysis for the OS protocol. The shear stress along the flow direction within a cycle can be decomposed in a Fourier series
\begin{equation} \sigma _{yz}(t) = \gamma _0 \sum _{n=1}^{\infty }(G^{\prime}_n \sin (n\omega t) + \omega \eta ^{\prime}_n \cos (n\omega t)), \end{equation}
where
$G^{\prime}_n$
and
$\eta ^{\prime}_n$
denote the
$n$
th elastic and viscous coefficients, respectively.
In the case of RS and RRS, the flow-vorticity plane rotates with respect to the gradient direction. Therefore, we need to use tensor rotation to calculate the stress tensor in the frame of reference rotating with the imposed shear
\begin{equation} \sigma ' = R \sigma\! R^T, \end{equation}
where
$R$
is the rotation matrix (cf. C2), and
$\sigma$
and
$\sigma '$
are the stress tensor in the laboratory and rotating frames, respectively. The normal and shear stress components are given as
\begin{align} \sigma^{\prime}_{\textit{ff}} &= \sigma _{{xx}}\sin ^2(\theta ) + ( \sigma _{{xy}} + \sigma _{{yx}})\sin (\theta )\cos (\theta ) + \sigma _{{yy}}\cos ^2(\theta ),\nonumber \\\sigma^{\prime}_{\textit{gg}} &= \sigma _{{zz}},\nonumber \\\sigma^{\prime}_{vv} &= \sigma _{{xx}}\cos ^2(\theta ) - (\sigma _{{yx}}+\sigma _{{xy}})\sin (\theta )\cos (\theta ) + \sigma _{{yy}}\sin ^2(\theta ),\nonumber \\\sigma^{\prime}_{\textit{fg}} &= \sigma _{{xz}}\sin (\theta )+\sigma _{{yz}}\cos (\theta ),\nonumber \\\sigma^{\prime}_{v\textit{g}} &= \sigma _{{xz}}\cos (\theta )-\sigma _{{yz}}\sin (\theta ), \end{align}
where
$\theta = \omega t$
is the angle of rotation of the flow-vorticity plane,
$\sigma^{\prime}_{\textit{ff}}$
,
$\sigma^{\prime}_{\textit{gg}}$
and
$\sigma^{\prime}_{vv}$
are, respectively, the normal stresses along the instantaneous flow, gradient and vorticity directions and
$\sigma^{\prime}_{\textit{fg}}$
, and
$\sigma^{\prime}_{v\textit{g}}$
are the viscous and elastic components of the shear stress in the rotating frame. Note that the normal stresses in the laboratory frame do not contribute to dissipation, while the shear stress in the laboratory frame appears in the normal stresses in the rotating frame. The suspension viscosity under RS and RRS can be calculated as
\begin{align} \eta ' &= \sigma^{\prime}_{\textit{fg}}/\gamma _0\omega , \end{align}
\begin{align} \eta '' &= -\sigma^{\prime}_{v\textit{g}}/\gamma _0\omega , \end{align}
\begin{align} \eta ^* &= \sqrt {\eta '^{2}+\eta ''^{2}}. \\[9pt] \nonumber \end{align}
In this work, we report the relative complex viscosities,
$\eta _r^* = \eta ^*/ \eta _0$
. In our particulate system, the stress response is closely in phase with the shear rate irrespective of the strain amplitude, as shown in figure 2. Therefore, the complex viscosity has a dominant viscous contribution, and the complex viscosity and dynamic viscosity are typically not distinguishable in the plots.
Figure 2. Stress–strain evolution for OS (a)
$\gamma _0=0.05$
, (b)
$\gamma _0=10.0$
. The signals are normalised for comparison. Here,
$T$
represents the period of a shear cycle.
In order to understand the non-Newtonian nature of the flow, we also examine the first,
$N_1$
, and second,
$N_2$
, normal stress differences, defined as
\begin{align} N_1 &= \sigma^{\prime}_{\textit{ff}} - \sigma^{\prime}_{\textit{gg}},\nonumber\\ N_2 &= \sigma^{\prime}_{\textit{gg}} - \sigma^{\prime}_{vv}. \end{align}
For the OS protocol, the stress terms in the above equation are simply
$\sigma^{\prime}_{\textit{ff}} = \sigma _{{yy}}$
,
$\sigma^{\prime}_{\textit{gg}} = \sigma _{{zz}}$
,
$\sigma^{\prime}_{vv} = \sigma _{xx}$
. For the RS and RRS protocol, on the other hand, the normal stresses are obtained from (2.10).
To investigate the suspension dynamics we calculate the mean squared displacement (MSD) of the particle ensemble stroboscopically, after an integral number of shear cycles
\begin{equation} \langle [\Delta r (\gamma _t)/d]^2 \rangle = 6 D_{\textit{eff}} \gamma _t, \end{equation}
where
$\Delta r$
is the particle displacement,
$d=2a$
the particle diameter,
$\langle \boldsymbol{\cdot }\rangle$
an average over all particles and the time,
$D_{\textit{eff}}$
the effective diffusivity and
$\gamma _t$
the total strain:
$\gamma _t = |\dot {\gamma }|t$
, for SS;
$\gamma _t = 4 \gamma _0 n$
, for OS;
$\gamma _t = 2 \pi \gamma _0 n$
, for RS and RRS, for
$n$
number of shear cycles.
To understand the evolution of the suspension microstructure we calculate the coordination number
$(Z)$
, i.e. the average number of particles in contact with any particle in the suspension, as well as the pairwise distribution of particles
$g(h,\theta )$
which provides information on the relative distribution of particle pairs in the suspension in terms of their radial surface separation
$(h)$
and relative orientation
$(\theta )$
when projected on the shear plane. Here,
$h\leqslant 0$
shows particle pairs in contact, and
$h \gt 0$
shows that particles are separated. The angle
$\theta$
is calculated with respect to the shear direction, where
$\theta \in [0,\pi /2]$
and
$[\pi , 3\pi /2]$
represent the two extensional quadrants, and
$\theta \in [\pi /2, \pi ]$
and
$[3\pi /2, 2\pi ]$
the two compressional quadrants.
2.2. Problem set-up
In our simulations, we consider a bidisperse suspension of
$500$
particles with a size ratio of
$1.4$
and particle number such that the suspension has an equal volume fraction of large and small particles. The computational domain is a cubic box with Lees–Edwards boundary conditions (Lees & Edwards Reference Lees and Edwards1972). The initial state is obtained for each simulation by preshearing a random configuration with a constant shear rate for a total strain of
$\dot {\gamma }t = 40$
. Note that this results in an initially anisotropic microstructure with many particles in contact. The suspension is then subjected to the different shear protocols under consideration (OS, RS and RRS) and data are collected for at least 200 accumulated strains, after an initial pre-conditioning for
$10$
shear cycles using the same protocol and strain amplitude.
The control parameters are the volume fraction of the solid particles
$\phi$
, the friction coefficient for particle–particle contact
$\mu _c$
and the strain amplitude
$\gamma _0$
. In this work, we present results for
$\phi =$
0.40, 0.50 and 0.55, and
$\mu _c=$
0.0, 0.2 and 0.5 to focus on the dense regime, where frictional contact is most prominent. The strain
$\gamma _0$
is varied from 0.05 to 10, while keeping the maximal shear rate
$\gamma _0\omega$
constant, to explore both small- and large-amplitude oscillatory/rotary shear. For RRS, the shear direction is reversed every half a cycle i.e. after a rotation of
$\pi$
radians, if not specified otherwise.
3. Results
3.1. Rheology
First, we report the main rheological observables: the complex viscosity and normal stress differences. Figure 3 shows the suspension relative complex viscosity,
$\eta _r^*$
, as a function of the strain amplitude
$\gamma _0$
at different volume fractions
$\phi$
and friction coefficients
$\mu _c$
, subject to OS (circle), RS (star) or SS (dotted). The main qualitative observation is a non-monotonic dependence of
$\eta _r^*$
on
$\gamma _0$
: specifically,
$\eta _r^*$
first decreases with the applied strain until it reaches a minimum at an intermediate strain
$\gamma _{0,m}$
(denoted hereafter as the minimum strain amplitude), then increases and converges to the steady shear viscosity at large
$\gamma _0$
(reported by dotted lines in the figure). In OS, a non-monotonic complex viscosity with respect to strain amplitude has long been observed (Breedveld et al. Reference Breedveld, van den Ende, Jongschaap and Mellema2001; Bricker & Butler Reference Bricker and Butler2006, Reference Bricker and Butler2007), and can be explained by the shear-induced microstructure and its evolution under repeated flow reversal (Ge & Elfring Reference Ge and Elfring2022). The observation of a similar viscous response in RS suggests that rotating the shear while keeping a constant shear rate has, on average, a similar effect on the suspension microstructure as oscillating the shear rate in one direction, since the suspension rheology ultimately depends on its microstructure. Despite the qualitative similarity of the viscous responses, the specific value of the minimum viscosity is higher in RS as compared with OS. As we will discuss later, this can be attributed to the absence of microstructure breakage in RS, resulting in the continuous resistance to the flow unlike that which occurs within a shear cycle in OS.
Figure 3. Complex viscosity for OS (circles), RS (stars), SS (dotted lines). Black lines
$\mu _c=0.0$
, red lines
$\mu _c=0.2$
, blue lines
$\mu _c=0.5$
.
The specific value of the minimum strain amplitude
$\gamma _{0,m}$
depends on the volume fraction
$\phi$
, the shearing protocol and the friction coefficient
$\mu _c$
in a complex way. Generally,
$\gamma _{0,m}$
is smaller under the RS protocol than under the OS protocol. At
$\phi =0.40$
, the minimum viscosity is attained at
$\gamma _{0,m} \approx 1$
for all three values of the friction coefficient
$\mu _c$
in OS (as also observed in the experiments by Bricker & Butler Reference Bricker and Butler2006, Reference Bricker and Butler2007), whereas the minimum moves to smaller
$\gamma _0$
values in RS. Varying
$\mu _c$
has only a weak effect on
$\gamma _{0,m}$
, slightly shifting
$\gamma _{0,m}$
to smaller values as it increases. With the increase in volume fraction
$\phi$
,
$\gamma _{0,m}$
decreases and the dependence of
$\gamma _{0,m}$
on the shearing protocol reduces. Since
$\gamma _{0,m}$
is a purely geometric parameter in the absence of attractive interparticle interactions (Corte et al. Reference Corté, Chaikin, Gollub and Pine2008; Ge et al. Reference Ge, Martone, Brandt and Minale2021), its dependence on
$\mu _c$
and the shearing protocol indicates a subtle interplay of the surface roughness (modelled here by
$\mu _c$
) and the microstructure formed under shear (Lemaire et al. Reference Lemaire, Blanc, Claudet, Gallier, Lobry and Peters2023).
Figure 4. Viscosity budget, total stress (solid lines), contribution from contact (circles), contribution from hydrodynamic (stars) for OS (a,c,e) and RS (b,d,f). Friction coefficient,
$\mu _c=0.0$
(black lines),
$\mu _c=0.2$
(red lines),
$\mu _c=0.5$
(blue lines). (a,b)
$\phi =0.40$
, (c,d)
$\phi =0.50$
, and (e,f)
$\phi =0.55$
.
To better understand the non-monotonic variation of the viscosity, we show the contributions from the hydrodynamic and non-hydrodynamic (contact) stresses to
$\eta _r^*$
, obtained from the total stress, for all values of
$\phi$
and
$\mu _c$
under OS or RS in figure 4. The decomposition reveals two rheological regimes. First, at low
$\gamma _0$
, the stress is always dominated by the hydrodynamic contribution (mainly lubrication), which is independent of
$\mu _c$
and decreases with increasing
$\gamma _0$
, leading to the reduction of the total viscosity at moderate
$\gamma _0$
. Second, from moderate to high
$\gamma _0$
, the contact stress also contributes to the rheology (especially at large
$\phi$
or
$\mu _c$
), and its role increases with
$\gamma _0$
. Note that, new to the RS protocol, the contact contribution is already non-zero for
$\gamma _0 \lt \gamma _{0,m}$
. This indicates the presence of particle contact at an earlier strain amplitude and explains why the minimal
$\eta _r^*$
is higher in magnitude and occurs at smaller
$\gamma _{0}$
in RS than in OS. For the low strain amplitudes,
$\mu _c$
has an almost negligible effect in either protocol as friction does not affect the hydrodynamic forces and, at smaller amplitudes, particles should not be in contact. On the other hand,
$\mu _c$
has a more prominent effect on
$\eta _r^*$
at large strain amplitudes, when the viscosity increases dramatically with
$\mu _c$
for the same
$\phi$
. This is expected as the presence of friction tends to enhance the importance of contacts. On the other hand, both the contact and hydrodynamic stresses increase with
$\phi$
, increasing the overall complex viscosity.
Figure 5. The first (a,c,e) and second (b,d,f) normal stress differences under RS at three volume fractions: (a,b)
$\phi =0.40$
, (c,d)
$\phi =0.50$
and (e,f)
$\phi =0.55$
. The inset figures show the individual contribution from the contact and hydrodynamic forces at
$\phi =0.50$
and
$\mu _c=0.5$
. The yellow dashed line shows the zero line for reference.
In addition to the complex viscosity, the first (
$N_1$
) and second (
$N_2$
) normal stress differences (NSDs) are typically used to characterise the rheology of suspensions. The cycle-averaged values of
$N_1$
and
$N_2$
are shown in figures 5 and 6 as a function of
$\gamma _0$
for different values of
$\phi$
and
$\mu _c$
for RS and OS, respectively. While both
$N_1$
and
$N_2$
are negligible at small
$\gamma _0$
, they become finite and negative when
$\gamma _0 \gt \gamma _{0,m}$
. The negligible
$N_1$
at smaller
$\gamma _0$
suggests that the pair distribution function characterising the suspension microstructure is nearly fore–aft symmetric (Morris Reference Morris2009); its finite values at larger
$\gamma _0$
indicates the breaking of such symmetry. Another observation is that
$N_1$
is slightly positive (more so at larger
$\mu _c$
) at
$\gamma _0 \approx \gamma _{0,m}$
. As discussed in Ge & Elfring (Reference Ge and Elfring2022), this is due to the more frequent particle collisions (resulting in compressive normal stresses) in the velocity-gradient direction than in the flow direction; cf. figure 6 inset and figure 7 in Lemaire et al. (Reference Lemaire, Blanc, Claudet, Gallier, Lobry and Peters2023). Here, our results further suggest that the effect can be enhanced by surface roughness (i.e. larger
$\mu _c$
). As for
$N_2$
, it is approximately zero at small strain amplitudes, but increases (in magnitude) sharply around
$\gamma _{0,m}$
, until it saturates to its corresponding SS value. A negative value of
$N_2$
is widely observed in the literature (Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018), which is attributed to the dominant contribution from contact stresses that occur mostly in the shear plane, i.e. there are more collisions in the velocity-gradient direction than in the vorticity direction. The NSDs for RS are consistent with the behaviour observed in our OS protocol, shown in figure 6, and also in the previous studies on OS.
Figure 6. The first (a,c,e) and second (b,d,f) normal stress differences under OS at three volume fractions: (a,b)
$\phi =0.40$
, (c,d)
$\phi =0.50$
and (e,f)
$\phi =0.55$
. The inset figures show the individual contribution from the contact and hydrodynamic forces at
$\phi =0.50$
and
$\mu _c=0.5$
. The yellow dashed line shows the zero line for reference.
Moreover, we note that both
$|{N_1}|$
and
$|{N_2}|$
increase with the friction coefficient
$\mu _c$
at higher applied strain amplitude
$\gamma _0$
for both OS and RS (with an exception at 40 %). Because the NSDs are mainly determined by contact forces at higher
$\gamma _0$
, the presence of friction enhances the effect of contacts. However, the strain amplitude at which
$N_2$
becomes finite is only weakly sensitive to the value of
$\mu _c$
and
$\phi$
, as the (absolute) slopes of the curves barely increase with
$\mu _c$
or
$\phi$
near the onset of
$N_2$
(cf. the relatively larger difference in OS, figure 6
b,d,f). Finally, we note that recent work by Ge & Elfring (Reference Ge and Elfring2022) has shown that certain rheological features – such as minimum complex viscosity and the onset of the second NSD – coincide with the critical strain amplitude associated with the RIT under OS. While the RS protocol is intrinsically irreversible due to the absence of sudden shear reversal, we observe analogous rheological signatures at an intermediate strain amplitude. This motivates our examination of the suspension’s dynamical response under RS in the next section.
3.2. Dynamics
In this section, we report the stroboscopic investigation of the particle dynamics in RS, comparing it with the standard OS protocol. In the case of OS, previous works have shown that periodically sheared suspensions can display a dynamical phase transition called RIT or absorbing-phase transition (Pine et al. Reference Pine, Gollub, Brady and Leshansky2005; Corte et al. Reference Corté, Chaikin, Gollub and Pine2008; Menon & Ramaswamy Reference Menon and Ramaswamy2009). In particular, there is a
$\phi$
-dependent critical strain amplitude
$\gamma _{0,c}$
below which the suspension can fall into a reversible ‘absorbing’ state, i.e. all particles return to their original positions after each strain cycle. On the other hand, for strain amplitudes larger than
$\gamma _{0,c}$
, the suspension may remain in an irreversible fluctuating state, where the particles do not return to their original positions and display a biased random walk-like dynamics if observed at the end of each strain cycle (Pine et al. Reference Pine, Gollub, Brady and Leshansky2005).
We therefore quantify the MSDs of suspended particles in three dimensions, as shown in figure 7. For OS, our results show that at large strain amplitudes the MSD scales linearly with the total strain at large strains, which indicates a diffusive dynamics. On the other hand, at lower strain amplitudes, the MSD increases faster than linearly; however, its value is much smaller and the dynamics is practically negligible. As for RS, the MSD scales linearly for all simulated strain amplitudes, see e.g. data for
$\gamma _0=0.05$
in figure 7.
Figure 7. Mean square displacements for (a) OS and (b) RS at different
$\gamma _0$
for
$\phi =0.55$
and
$\mu _c= 0.5$
. The black lines are for reference to show the slopes.
We calculate the effective diffusivity
$D_{\textit{eff}}$
for both the OS and RS protocols from the linear fit of the MSD data similar to Pine et al. (Reference Pine, Gollub, Brady and Leshansky2005); see figure 8. Here,
$D_{\textit{eff}}$
is calculated after the suspension has been sheared for
$300$
strain units. For OS, the dynamics clearly displays two regimes and a transition in between, consistent with the RIT; however, we cannot rigorously determine the critical strain amplitude
$\gamma _{0,c}$
at which the transition occurs. As an alternative, we plot the experimentally fitted
$\gamma _{0,c}$
by Pine et al. (Reference Pine, Gollub, Brady and Leshansky2005),
$\gamma _{0,c} = C\phi ^{-\alpha }$
, where
$C = 0.14 \pm 0.03$
, and
$\alpha = 1.93 \pm 0.14$
(shaded regions in figure 8). The fitted values (based on the experimental data at
$\phi \leqslant 0.4$
) appear to be slightly higher than our simulations, particularly at
$\phi \gt 0.4$
and
$\mu _c \gt 0$
. This could be due to differences between the experimental system and our numerical model or uncertainties in the fitting procedures. Nevertheless, the general trend of an earlier transition at higher
$\phi$
is consistently observed.
For suspension undergoing RS,
$D_{\textit{eff}}$
increases gradually with
$\gamma _0$
and is not negligible even at the smallest strain amplitudes, in stark contrast to the dynamics in OS. As an example, note that
$D_{\textit{eff}}$
can be approximately 2 orders of magnitude larger than for the corresponding OS case at the same strain below
$\gamma _{0,m}$
. Although the rate of increase of
$D_{\textit{eff}}$
with respect to
$\gamma _{0}$
seems to be higher at around
$\gamma _{0,m}$
(cf. the
$\phi =0.55$
case in the figure), the dynamics is always diffusive even at the smallest
$\gamma _0$
(cf. figure 7
b). Therefore, we show that a non-monotonic rheological response under periodic shear, with a minimal viscosity at an intermediate
$\gamma _0$
, does not necessarily imply the existence of RIT. Another interesting observation is that, although
$D_{\textit{eff}}$
depends on
$\mu _c$
in the diffusive regime in OS (and to a less extent over the same range of
$\gamma _0$
in RS), the dependence diminishes at larger
$\gamma _0$
in the irreversible regime. As irreversible dynamics arises from non-hydrodynamic particle contacts, we examine their effect on the microstructure next.
Figure 8. Effective diffusivity (
$D_{\textit{eff}}$
) vs
$\gamma _0$
for OS and RS. The shaded region shows the region of prediction of the critical amplitude from the empirical scaling,
$\gamma _{0,c} = C\phi ^{-\alpha }$
, where
$C = 0.14 \pm 0.03$
, and
$\alpha = 1.93 \pm 0.14$
(Pine et al. Reference Pine, Gollub, Brady and Leshansky2005).
3.3. Microstructure
Non-Newtonian stresses are generally caused by an anisotropic microstructure. Because hard spheres do not deform and are isotropic in shape, the microstructure is entirely dictated by the spatial arrangement in the suspension. To find the origin of the differences in the suspension dynamics between the sudden shear reversal (OS) and the gradual shear reversal (RS) protocols, we first examine the particle coordination number (
$Z$
), defined as the average number of particles in contact with one particle in suspension.
Figure 9. Coordination number (Z). (a) OS and (b) RS. The coordination number is averaged over roughly
$100$
strain units in the steady state.
Figure 10. Pairwise particle distribution
$g(h,\theta )$
. The angles
$\theta = 0^{\circ }$
and
$90^{\circ }$
represent the flow and velocity-gradient directions, respectively,
$\theta \in [0,\pi /2]$
or
$[\pi , 3\pi /2]$
represents extensional quadrants and
$\theta \in [\pi /2, \pi ]$
or
$[3\pi /2, 2\pi ]$
represents the compressional quadrants. The contact point
$h = 0$
is marked in red,
$h\leqslant 0$
shows particle pairs in contact and
$h \gt 0$
shows particles are separated. Figures are shown for suspension at
$\phi =0.40, \mu _c=0.0$
.
The coordination number
$Z$
as a function of
$\gamma _{0}$
is displayed in figure 9 for all the values of volume fraction and friction coefficient under investigation. For OS,
$Z \approx 0$
at smaller
$\gamma _0$
, suggesting that the particles self-organise into contact-free absorbing states. In contrast,
$Z\gt 0$
even at low strain amplitudes in RS. This reveals why RIT is not observed in RS: in this case, we observe non-hydrodynamic particle contacts, leading to a diffusive dynamics, at all strain amplitudes. The particle contacts also affect the suspension viscosity. As shown in figure 4, the contribution from the contact stress remains negligible for
$\gamma _0\lt \gamma _{0,c}$
in OS, but it is non-negligible for
$\gamma _0\lt \gamma _{0,m}$
in RS. Since the total viscosity is mainly due to the hydrodynamic interactions at small amplitudes, the rheology is nevertheless similar in OS and RS.
Next, we investigate the pairwise particle distribution
$g(h,\theta )$
, in the vicinity of the particle contacts,
$h \in [\pm 0.0002]$
(see figure 10). Here,
$g(h,\theta )$
provides information about the relative position of particle pairs in the suspension, projected onto the shear plane, as described in § 2.1. The figure shows that particles are isotropically distributed and well separated (thus the sparsity of data within a narrow gap) at the smallest strain amplitude under OS. Increasing the strain amplitude above
$\gamma _{0,m}$
particles come into closer contact
$(h\leqslant 0)$
and form an anisotropic structure with an accumulation of particles around the compression axes (
$\theta =3\pi /4$
and
$7\pi /4$
). These observations are in line with the previous findings for suspensions under simple shear (Morris Reference Morris2009). Moreover, we note that, within the compression quadrants, neighbouring particles tend to accumulate more towards
$\theta \approx \pi /2$
and
$3\pi /2$
more so than for
$\theta \approx \pi$
or
$\theta \approx 2\pi$
, leading to more frequent collision along the velocity-gradient direction than in the flow direction. As discussed by Ge & Elfring (Reference Ge and Elfring2022) (see their figure 3, inset), such anisotropic collision might explain the slightly positive first NSD
$N_1$
near
$\gamma _{0,m}$
seen in figure 6(a,c,e).
Although collisions are enough to cause the suspension dynamics to transition from reversible to irreversible, our analysis provides a more detailed understanding of how the dynamics is related to the rheology as the strain amplitude is increased. For RS, in contrast, the figure shows the development of an anisotropic microstructure and particle contact even at the smaller strain amplitudes. Both the number of contacting pairs and the extent of anisotropy increase with
$\gamma _0$
, suggesting an analogy between RS and SS in terms of microstructural reorganisation.
So far, we have examined and compared the rheology, dynamics and microstructure of OS and RS. These are the two extremes: the suspension undergoes sudden shear reversal in OS, whereas it undergoes gradual shear rotation and, ultimately, reversal in RS. Moreover, OS is a reversible protocol, while RS is not. To bridge the gap between the two, we therefore consider a combination of both, i.e. RRS, and present the main results in the following section.
Figure 11. (a,c) Viscosity budget: total stresses (solid lines), contact stresses (circles), hydrodynamic stresses (stars) for OS (black lines), RS (red lines) and RRS (blue lines) at
$\phi =0.55$
. (b,d) Normal stress differences comparing OS, RS and RRS for
$\phi =0.55$
and
$\mu _c=0.0, 0.5$
. (b) first NSD (
$N_1$
) and (d) second NSD (
$N_2$
); (c) uses the same legend as for (a), and (d) the same as for (b).
3.4. Reversible rotary shear
Reversible rotary shear, as already introduced in § 2, is a combination of OS and RS, where the directions of both the flow and the shear rotation are reversed at half-cycles. For comparison, we investigate the suspension under RRS protocol for
$\phi =0.55$
and two friction coefficients
$\mu _c = 0.0$
and
$0.5$
.
As shown in figure 11(a,c), suspensions undergoing RRS have a non-monotonic dependence of the complex viscosity on the strain amplitude similar to OS and RS: the complex viscosity decreases at smaller strain amplitudes, attains a minimum value at an intermediate strain amplitude and increases afterwards. As for the NSDs, we observe that both
$N_1$
and
$N_2$
are negative at large strain amplitudes, with
$|N_1| \ll |N_2|$
, as shown in figure 11(b,d). From these observations we conclude that the suspension rheology is not sensitive to the different shearing protocols considered.
Next, we investigate the dynamics of suspensions undergoing RRS and compare it with the dynamics of the same suspension undergoing either OS or RS. The MSDs versus strain are therefore shown in figure 12(a); the data reveal a similar trend as that observed for OS: the MSD does not scale linearly with the strain at smaller strain amplitudes, cf. the curve pertaining the case
$\gamma _0 = 0.05$
; it is, however, linear for larger strain amplitudes. Interestingly, the suspension under RRS experiences RIT, as shown in figure 12(b): the suspension is in the reversible ‘absorbing’ state before a critical strain amplitude,
$\gamma _0 \lt \gamma _{0,c}$
, and in the irreversible ‘diffusive’ state at
$\gamma _0 \gt \gamma _{0,c}$
.
Figure 12. (a) Particle MSD versus strain for suspensions undergoing RRS and the different values of the strain amplitude in the legend. (b) Comparison of the effective particle diffusivity for the three protocols under investigation: OS, RS and RRS.
The above results show that suspensions driven by RRS can undergo the RIT, as observed in OS but not in RS, despite the fact that they all share essentially the same rheology. The different dynamical behaviour of the suspension may be ascribed to the fundamental difference between the protocols: shear reversibility. The OS and the RRS have shear reversibility, whereas the RS does not. This suggests that absorbing states can only be reached via protocols that involve shear reversal. Hence, we can conclude that the difference that the shear reversal brings is in the suspension microstructure. The coordination number
$(Z)$
versus the strain amplitude
$(\gamma _0)$
is shown in figure 13 for the three protocols. The RRS protocol is, from all perspectives, similar to OS. At smaller strain amplitude the microstructure needs a long time to evolve but is broken before it develops by reversing the shear, which causes an immediate drop in the coordination number i.e. particle contacts; this mechanism is essential for the irreversible dynamics. When applying the RS protocol, the microstructure does not break, and the particle contacts grow continuously until saturation, similarly to the behaviour in simple shear flows, which do not include shear reversal. This suggests that the differences between OS and RS are primarily due to the presence of sudden shear reversals in the former but not the latter, affecting both the dynamics at low strain amplitudes and the saturation to SS at higher strain amplitudes.
Figure 13. Coordination number (Z) in OS, RS and RRS for
$\phi =0.55$
,
$\mu _c=0.0$
(a) and
$\mu _c=0.5$
(b). The coordination number is averaged over roughly
$100$
strain units in the steady state.
4. Conclusion
In this work, we performed numerical simulations of dense, bidisperse suspensions of rigid, spherical particles at volume fractions ranging from 0.40 to 0.55 under three periodic flow protocols: the classical OS, the novel RS and the RRS, a combination of the first two.
Summarising the main findings, we observe the suspension rheology, both the viscosity and the NSDs, to be only weakly sensitive to the flow protocol, but strongly dependent on the applied strain amplitude. The viscosity shows a non-monotonic behaviour with the strain amplitude, attaining a minimum viscosity at an intermediate strain amplitude (
$\gamma _{0,m}$
) that depends on the volume fraction. At the same strain amplitude
$\gamma _{0,m}$
, we observed the onset of second NSDs regardless of the shear protocol. However, the suspension dynamics changes with the flow protocol and the strain amplitude. Specifically, under OS and RRS, the suspensions show reversibility at low strain amplitudes and irreversibility at large strain amplitudes, thus undergoing a RIT at an intermediate strain amplitude (
$\gamma _{0,c}$
). In contrast, the suspension dynamics under the RS protocol is inherently irreversible at all strain amplitudes.
The reason for this difference lies in the flow-induced suspension microstructure: for shear reversal protocols (OS and RRS), the particle pairs are more separated and isotropically distributed at smaller strain amplitudes, only to become anisotropic and in closer contact as the strain amplitude increases beyond a critical value. On the other hand, under the gradually varying RS protocol, particle pairs are anisotropically distributed and at contact at all strain amplitudes, even if the number of contacts is lower at low strain amplitudes. The presence of particle contacts at even the smallest strain amplitudes reflects the origin of the inherent irreversibility of RS, as also suggested in the literature based on simpler descriptions of particle interactions (Corte et al. Reference Corté, Chaikin, Gollub and Pine2008; Menon & Ramaswamy Reference Menon and Ramaswamy2009; Pham et al. Reference Pham, Metzger and Butler2015).
For all protocols considered, we observe the build-up of particle contacts and, consequently, contact stresses, which are closely linked to the emergence of irreversibility in the stroboscopic dynamics. In shear reversal protocols such as OS and RRS, the contact network breaks upon shear reversal and then gradually rebuilds. At low strain amplitudes, only the nearest particles can collide, allowing the system to reach a self-avoiding and reversible (or absorbing) state after multiple shear reversals. This state is characterised by a microstructure with almost no contacts. At higher strain amplitudes, however, the contact network has sufficient time to form and persist before being disrupted, where particles collide with multiple particles before any shear reversal, leading to an irreversible dynamics. In contrast, under the RS protocol, which lacks any shear reversal, the contact network evolves continuously and saturates into a steady-state structure without undergoing periodic breaking. As a result, irreversibility is observed at all strain amplitudes. This suggests that the formation of reversible or absorbing states is intrinsically tied to the presence of shear reversal, which periodically resets the microstructure and enables the system to avoid sustained contact-driven irreversibility.
An important finding of the current work is that the presence of rheological characteristics such as minimum viscosity and the onset of second NSDs is not a sufficient condition for the suspensions to show dynamical characteristics such as RIT. As demonstrated by the RS protocol, these rheological features can coexist with the irreversible behaviour. Finally, the current understanding in the literature suggests that a time-reversible flow protocol can lead to either a stroboscopically reversible or irreversible particle dynamics, depending on parameters such as strain amplitude and volume fraction. Our results extend this understanding by showing that, in the absence of time reversibility in the flow protocol, as in the RS case, stroboscopic reversibility does not emerge, even at low strain amplitudes. Based on this, we argue that time reversibility of the flow protocol is a necessary (although not sufficient) condition for achieving an absorbing state in dense suspensions. For future works, it would be interesting to examine how the absence or presence of an actual shear reversal step (as in OS and RRS) controls the emergence of an absorbing state. For example, tuning the number of rotations after which shear is reversed in the RRS protocol could offer a mechanism to interpolate between the fully reversible and inherently irreversible flow regimes.
Funding
This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 955605, YIELDGAP ITN network. Z.G. is supported by an international postdoc grant from the Swedish Research Council (Grant No. 2021-06669VR). M.T. acknowledges funding from the Swedish Research Council under grant No. 2021-04997. O.T. acknowledges funding from European Research Council through the Starting Grant MUCUS (ERC-StG2019-852529).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Forces
The individual force terms in the right-hand side of (2.4) acting on the
$i^{}$
th particle of radius
$a_i$
, translational velocity
$\boldsymbol {u}_i$
and angular velocity
$\boldsymbol {w}_i$
are described here. The
$j$
index is used for the interacting particle. For details, please refer to the previous studies (Seto et al. Reference Seto, Mari, Morris and Denn2013; Mari et al. Reference Mari, Seto, Morris and Denn2014; Cheal & Ness Reference Cheal and Ness2018).
-
(i) Stokes’ drag: the Stokes drag force acts on each particle
$i$
due to the surrounding fluid. It is given by(A1)
\begin{equation} \boldsymbol {F}_{\!i}^S = - 6 \pi \eta _0 a_i (\boldsymbol {u}_i - \boldsymbol {U}_{\!\infty} ),\end{equation}
where
$\eta _0$
is the dynamic viscosity of the fluid,
$\boldsymbol {U}_{\!\infty}$
is the translational velocity of the underlying shear flow at the location of the particle centre of mass and
$a_i$
is particle radius. -
(ii) Lubrication force: the lubrication force acting on each lubricating particle pair
$i,j$
is given as(A2)
\begin{align} \boldsymbol {F}^L_{\textit{ij}} &= - \big(X_{\textit{ii}}^A {\mathbb{P}}_n + Y_{\textit{ii}}^A {\mathbb{P}}_t \big)(\boldsymbol {u}_i - \boldsymbol {u}_{\!j}) + Y^B_{\textit{ii}}(\boldsymbol {w}_i \times \boldsymbol {n}_{\textit{ij}}) + Y^B_{\textit{ji}}(\boldsymbol {w}_j \times \boldsymbol {n}_{\textit{ij}}),\nonumber\\[4pt] \boldsymbol {F}^{\boldsymbol {L}}_{\textit{ji}} &= -\boldsymbol {F}^{\boldsymbol {L}}_{\textit{ij}}, \end{align}
where
$\boldsymbol {n}_{\textit{ij}}$
denotes the unit normal vector pointing from particle
$i$
to particle
$j$
,
${\mathbb{P}}_n = \boldsymbol {n}_{\textit{ij}}\boldsymbol {n}_{\textit{ij}}$
and
${\mathbb{P}}_t = {\mathbb{I}} - \boldsymbol {n}_{\textit{ij}}\boldsymbol {n}_{\textit{ij}}$
represent the normal and tangential projection matrices. The
$X$
and
$Y$
are scalar resistances depending on
$\eta _0, a_i, a_j$
and the gap between the two particles, whose complete expressions are given in Ge & Brandt (Reference Ge and Brandt2020). -
(iii) Contact force: in case of small overlap between particles
$i$
and
$j$
, equal and opposite contact forces are exerted on the pair, following Coulomb’s friction law(A3a)
\begin{align} \boldsymbol {F}^C_{\textit{ij}} &= -k_n\boldsymbol {h}_{\textit{ij}} - k_t \boldsymbol{\xi }_{\textit{ij}}, \end{align}
(A3b)
\begin{align} \boldsymbol {F}^C_{\textit{ji}} &= -\boldsymbol {F}^C_{\textit{ij}}, \end{align}
(A3c)
\begin{align} |k_t \boldsymbol{\xi }_{\textit{ij}}| &\leqslant \mu _c |k_n\boldsymbol {h}_{\textit{ij}}| . \end{align}
In the expressions above,
$\boldsymbol {h}_{\textit{ij}} = h_{\textit{ij}}\boldsymbol {n}_{\textit{ij}}$
is the normal surface gap between particle i and j,
$k_n$
the normal spring constant,
$k_t$
the tangential spring constant and
$\mu _c$
the friction coefficient. The friction force (
$k_t {\boldsymbol{\xi}}_{\textit{ij}}$
) comes as a consequence of particle roughness. Here, we adopt the stick/slide model for friction force. The tangential stretch vector is calculated as(A4)
\begin{equation} \boldsymbol{\xi }_{\textit{ij}} = \left \{\!\! \begin{array}{ll} \int _{t_0}^{t} -{\mathbb{P}}_t [(\boldsymbol {u}_i-\boldsymbol {u}_{\!j})+(a_i \boldsymbol {w}_i+a_j \boldsymbol {w}_j) \times \boldsymbol {n}_{\textit{ij}}]{\textrm d}t', & \mbox{if} \ |\boldsymbol{\xi }_{\textit{ij}}| \lt |\boldsymbol{\xi }_{\textit{max}}|,\\[4pt] \boldsymbol{\xi }_{\textit{max}} & \mbox{otherwise}, \end{array} \right . \end{equation}
where
$\boldsymbol{\xi }_{\textit{max}}$
comes from the Coulomb law of friction in (A3c).
Appendix B. Numerical implementations
The governing equations described in (2.4) are marched in time using the modified velocity-Verlet algorithm (Groot & Warren Reference Groot and Warren1997)
\begin{align} \boldsymbol {x}_i^{n+1} &= \boldsymbol {x}_i^{n} + \Delta t \boldsymbol {u}_i^{n} + \frac {\Delta t^2}{2} \boldsymbol{ \alpha }_i^{n} \end{align}
\begin{align} \boldsymbol {u}_i^{n+\frac {1}{2}} &= \boldsymbol {u}_i^{n} + \frac {\Delta t}{2} \boldsymbol{\alpha }_i^{n}, \end{align}
\begin{align} \boldsymbol{\alpha }_i^{n+1} &= {\mathcal{F}} \{\boldsymbol {x}_i^{n+1},\boldsymbol {u}_i^{n+\frac {1}{2}}\}, \end{align}
\begin{align} \boldsymbol {u}_i^{n+1} &= \boldsymbol {u}_i^{n} + \frac {\Delta t}{2} \big(\boldsymbol{\alpha }_i^{n} + \boldsymbol{\alpha }_i^{n+1} \big). \end{align}
Here,
$\boldsymbol {x}_i,\boldsymbol {u}_i, \boldsymbol{\alpha _i}$
denote the position, velocity and acceleration vector of the
$i^{}$
th particle. The time steps are denoted by:
$n, n+({1}/{2}), n+1$
are the current, half-time step forward and one time step forward, respectively, and
${\mathcal{F}}$
denote the force functional as in (2.4). To simulate the shear flow and remove the wall effect we imposed a Lees–Edwards (Lees & Edwards Reference Lees and Edwards1972) boundary condition on particle positions and velocity components as follows.
For OS
\begin{align} y &= \left \{\!\! \begin{array}{ll} \textit{mod}\{(y + L_y - y'), L_y\}, &\mbox{if}\ z\gt L_z,\\ \textit{mod}\{(y + L_y + y'), L_y\}, &\mbox{if}\ z\lt 0,\\ \textit{mod}\{(y + L_y), L_y\}, &\mbox{otherwise},\\ \end{array} \right . \end{align}
\begin{align} z &= \textit{mod}\{(z + L_z), L_z\}, \end{align}
\begin{align} x &= \textit{mod}\{(x + L_x), L_x\}, \end{align}
\begin{align} u_y &= \left \{\!\! \begin{array}{ll} u_y - u^{\prime}_y, &\mbox{if}\ z\gt L_z,\\ u_y + u^{\prime}_y, &\mbox{if}\ z\lt 0.\\ \end{array}\right. \\[9pt] \nonumber \end{align}
For RS and RRS, (B2a
), (B2b
), (B2d
) are used in addition to the shear periodicity in
$x$
(secondary shear) direction
\begin{align} x &= \left \{ \begin{array}{ll} \textit{mod}\{(x + L_x - x'), L_x\}, &\mbox{if}\ z\gt L_z,\\ \textit{mod}\{(x + L_x + x'), L_x\}, &\mbox{if}\ z\lt 0,\\ \textit{mod}\{(x + L_x), L_x\}, &\mbox{otherwise},\\ \end{array} \right . \end{align}
\begin{align} u_x &= \left \{ \begin{array}{ll} u_x - u^{\prime}_x, &\mbox{if}\ z\gt L_z,\\ u_x + u^{\prime}_x, &\mbox{if}\ z\lt 0,\\ \end{array}\right . \end{align}
where
$x'$
and
$y'$
are the position shifts in the
$x$
, and
$y$
directions, respectively, and
$u^{\prime}_x$
and
$u_y^{\prime}$
are the velocity shifts when crossing the periodic boundaries, given by
\begin{equation} \begin{array}{llll} x' = \textit{mod}\{\dot {\gamma }_{xz} L_zt, L_x \}, &u_x^{\prime} = \dot {\gamma }_{xz}L_z, \\ y' = \textit{mod}\{\dot {\gamma }_{yz} L_zt, L_y \}, &u_y^{\prime} = \dot {\gamma }_{yz}L_z. \\ \end{array} \end{equation}
Appendix C. Stress rotation
The stress tensor in continously rotating frame of flow-vorticity directions is given as
\begin{equation} \sigma ' = R \sigma \!R^T, \end{equation}
where, the rotation matrix
$R$
is given as
\begin{equation} R = \begin{bmatrix} \cos (-\theta ) & \sin (-\theta ) & 0 \\ -\sin (-\theta ) & \cos (-\theta ) & 0 \\ 0 & 0 & 1 \end{bmatrix}\!, \end{equation}
where
$\theta = \omega t$
in clockwise direction (therefore negative) for RS. Therefore, the stress tensor in the new basis is
\begin{align} \sigma ' &\!=\! \begin{bmatrix} \sigma^{\prime}_{vv} & \sigma^{\prime}_{v\textit{f}} & \sigma^{\prime}_{v\textit{g}}\\ \sigma^{\prime}_{\textit{f}v} & \sigma^{\prime}_{\textit{ff}} & \sigma^{\prime}_{\textit{fg}}\\ \sigma^{\prime}_{\textit{g}v} & \sigma^{\prime}_{\textit{g}\kern1pt\textit{f}} & \sigma^{\prime}_{\textit{gg}} \end{bmatrix}\!, \end{align}
\begin{align} \sigma ' &\!=\! \begin{bmatrix} \begin{array}{l} \sigma _{xx}\cos ^2(\theta ) + \sigma _{yy}\sin ^2(\theta )\\ \!\!\!\!\!\!\!\! - (\sigma _{xy} + \sigma _{yx})\sin (\theta )\cos (\theta ) \end{array} & \begin{array}{l} \sigma _{xx}\sin (\theta )\cos (\theta ) - \sigma _{yx}\sin ^2(\theta ) \\ \!\!\! + \sigma _{xy}\cos ^2(\theta ) - \sigma _{yy}\sin (\theta )\cos (\theta ) \end{array} & \begin{array}{l} \sigma _{xz}\cos (\theta ) \\ \!\!- \sigma _{yz}\sin (\theta ) \end{array}\\ \begin{array}{l} \sigma _{xx}\sin (\theta )\cos (\theta ) + \sigma _{yx} \cos ^2(\theta ) \\ \!\!\! - \sigma _{xy}\sin ^2(\theta ) - \sigma _{yy}\sin (\theta )\cos (\theta ) \end{array} & \begin{array}{l} \sigma _{xx}\sin ^2(\theta ) + \sigma _{yy}\cos ^2(\theta ) \\ \!\!\!\! +(\sigma _{yx}+\sigma _{xy})\sin (\theta )\cos (\theta ) \end{array} &\begin{array}{l} \!\!\!\!\!\sigma _{xz}\sin (\theta ) \\ \!\!\!\!\!\!\!\!+ \sigma _{yz}\cos (\theta ) \end{array}\\ \begin{array}{l} \sigma _{zx}\cos (\theta ) - \sigma _{zy}\sin (\theta ) \end{array} & \begin{array}{l} \sigma _{zx}\sin (\theta ) + \sigma _{zy}\cos (\theta ) \end{array} & \sigma _{zz}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \end{bmatrix}\!. \end{align}
For RRS,
$\theta$
takes the form
\begin{equation} \theta = \left \{ \begin{array}{ll} \omega t & \mbox{if}\ 0 \leqslant t \lt T/2, \\ \pi + \omega (T-t) & \mbox{if}\ T/2 \leqslant t \lt T. \end{array} \right . \end{equation}
Appendix D. Validations
To validate our simulations we compared the simulation results with those in literature. The rheology; relative complex viscosity ratio is compared with the experimental results by Bricker & Butler (Reference Bricker and Butler2006) and the dynamics; effective diffusivity is compared with the experimental and numerical results by Pine et al. (Reference Pine, Gollub, Brady and Leshansky2005).
Figure 14. Validation at 40 % volume fraction in OS. (a) Complex viscosity ratio against the experiments by Bricker and Butler (Bricker & Butler Reference Bricker and Butler2006) and dynamics against experiments by (Pine et al. Reference Pine, Gollub, Brady and Leshansky2005), effective diffusion in the (b) flow direction, (c) gradient direction for volume fraction
$\phi =0.40$
.
Figure 14(a) compares the relative complex viscosity obtained from our simulations and experimental result from Bricker & Butler (Reference Bricker and Butler2006) for different strain amplitudes for a volume fraction of
$\phi =0.40$
. The simulation results are shown for three friction coefficients
$\mu _c = 0.0, 0.2$
and
$0.5$
. The complex viscosity of both works shows qualitative agreement. The relative complex viscosity shows non-monotonic behaviour with respect to strain amplitude: it decreases at smaller
$\gamma _0$
, attains a minimum value at an intermediate
$\gamma _0$
and then increases at the larger
$\gamma _0$
. The precise value of the viscosity not only depends upon
$\gamma _0$
but on the physiochemical interactions between the particles as well.
Figures 14(b) and 14(c) show effective diffusivity in flow and gradient directions, respectively, comparing our simulation results with the experimental and numerical results from Pine et al. (Reference Pine, Gollub, Brady and Leshansky2005). It can be seen that our simulation results closely match the numerical results from Pine et al. (Reference Pine, Gollub, Brady and Leshansky2005). However the quantitative difference could arise from the difference in shear protocols: they considered square wave shear protocol whereas we considered the OS protocol, the particle dispersity: we have bi-dispersed particles with a size ratio of
$1.4$
, whereas they considered a monodispersed system. However, their own numerical and experimental results have quantitative discrepancies owning to the machine precession, etc. Overall, the rheology and dyanamics from our simulations closely follow the behaviour observed in the literature.
We checked the initial condition dependency of mechanical response. Figure 15 shows the relative complex viscosity with respect to strain amplitude for two different initial conditions obtained by
$(i)$
a random distribution of particles and
$(ii)$
after
$40$
strain units of pre-shearing for both OS and RS. For both OS and RS, we see that the steady-state results are independent of the initial conditions i.e. the suspension’s mechanical response is independent of the particles’ initial positions.
Figure 15. Effect of initial condition on suspension rheology. Two initial conditions are compared for packing fraction (
$\phi = 0.50$
), and friction coefficient (
$\mu _c=0.5$
) for both (a) OS and (b) RS.
Open access
