Two-scale evolution during shear reversal in dense suspensions

We use shear reversal simulations to explore the rheology of dense, non-Brownian suspensions, resolving lubrication forces between neighbouring particles and modelling particle surface contacts. The transient stress response to an abrupt reversal of the direction of shear shows rate-independent, nonmonotonic behaviour, capturing the salient features of the corresponding classical experiments. Based on analyses of the hydrodynamic and particle contact stresses and related contact networks, we demonstrate distinct responses at small and large strains, associated with contact breakage and structural re-orientation, respectively, emphasising the importance of particle contacts. Consequently, the hydrodynamic and contact stresses evolve over disparate strain scales and with opposite trends, resulting in nonmonotonic behaviour when combined. We further elucidate the roles of particle roughness and repulsion in determining the microstructure and hence the stress response at each scale.

We use shear reversal simulations to explore the rheology of dense, non-Brownian suspensions, resolving lubrication forces between neighbouring particles and modelling particle surface contacts. The transient stress response to an abrupt reversal of the direction of shear shows rate-independent, nonmonotonic behaviour, capturing the salient features of the corresponding classical experiments. Based on analyses of the hydrodynamic and particle contact stresses and related contact networks, we demonstrate distinct responses at small and large strains, associated with contact breakage and structural re-orientation, respectively, emphasising the importance of particle contacts. Consequently, the hydrodynamic and contact stresses evolve over disparate strain scales and with opposite trends, resulting in nonmonotonic behaviour when combined. We further elucidate the roles of particle roughness and repulsion in determining the microstructure and hence the stress response at each scale.

I. INTRODUCTION
The flow behaviour of dense suspensions is strongly affected by details of the microstructure and interparticle forces [1]. Recent theoretical [2], experimental [3,4] and computational [5,6] work suggests that particle surface contacts make a major contribution to suspension rheology, though their precise role, and importance relative to hydrodynamic interactions, is still debated. To this end, shear reversal experiments, in which the flow direction is suddenly reversed, prove to be an elegant means of estimating the contact stress (CS) contribution, while probing microstructural anisotropy. For example, it was shown that following a flow cessation period, the shear [7] (or normal [8]) stress in a suspension of ∼ 40 µm polystyrene spheres reaches an "immediate" peak upon reversal, then evolves nonmonotonically over a strain of ≈ 3 to its steady state. The large-strain evolution was attributed to microstructural realignment [7]; the initial stress peak was hypothesized to represent the hydrodynamic stress (HS), leading to a suggestion of the larger role played by particle contacts in denser suspensions [8]. It remains difficult to isolate the evolution of the CS and HS contributions and link the microstructural effect to the puzzling nonmonotonic behaviour. In the present paper, we reveal responses at two disparate strain scales, an elegent manifestation of the microfragile vs. macrofragile distinction proposed by Cates et al. [9]. The small strain stress peak is shown to be a hydrodynamic response to surface-contact breakage, but is distinct from the steady state HS. The large strain scale is determined by microstructural reorientation, as predicted [7]. The nonmonotonic behaviour is a combined effect of the CS and HS evolutions. We show that different surface characteristics control the stress response at different strain scales, meaning our two-scale explanation, and hence the micro-versus macro-fragility paradigm, can be applied usefully to a wide range of suspended systems.

II. SIMULATION MODEL
We solve the equations of motion numerically [10] for neutrally buoyant suspended non-Brownian particles, subject to forces and torques arising due to hydrodynamics and particle surface contact [11]. For dense suspensions, in which the average surface separation between neighbouring particles becomes very small, the full hydrodynamic resistance matrix [12] can be suitably approximated by resolving pairwise, frame-invariant lubrication forces [13], which diverge at contact and significantly exceed the long-range force components. Such a simplification has been proven to be effective in capturing the behaviour of dense suspensions [6,14,15]. For an interaction between particles i and j, (with particle and fluid density ρ) the force and torque on particle i due to hydrodynamic lubrication can be expressed as for particle diameter d i , fluid viscosity η f , particle translational and rotational velocity vectors v i and ω i respectively, centre-to-centre unit vector n ij pointing from particle j to i and identity tensor I, with the squeeze a sq , shear a sh and pump a pu resistance terms as derived by arXiv:1509.01530v1 [cond-mat.soft] 4 Sep 2015 [16], for β = d j /d i , as: The separation between particles i and j is calculated for centre-to-centre vector r. We calculate the lubrication force when the interparticle gap h is smaller than h max = 0.05d (where d is the harmonic average particle diameter). An increasing body of evidence [3,4] shows that direct particle-particle surface contacts can play a major role in suspension rheology; indeed, simulations that strictly resolve lubrication forces (treating particles as ideally hard and ideally smooth) [17] have proven to be inadequate for capturing dense suspension rheology for cases where particleparticle contacts are presumed to be important. Therefore we truncate the lubrication divergence and regularize the contact singularity at a typical asperity length scale h min (= 0.001d unless specified otherwise), i.e., setting h = h min in the force calculation, when h < h min . We use a value of η f = 0.1 [viscosity unit: ρd 2 /t].
Mechanical contact occurs at h ≤ 0, giving normal repulsive and tangential forces described by a linear spring model and related through a Coulomb friction coefficient µ p (= 0.2 unless specified otherwise) [11]. A linear (as opposed to Hertzian) spring is chosen for convenience, though we expect Hertzian results to lead to identical conclusions reagarding the respective roles of contacts and lubrication. The normal (F c,n ) and tangential (F c,t ) contact force and torque Γ c are given by for a collision between particles i and j with normal and tangential spring stiffnesses k n and k t respectively [k n = 20000, unit: ρd 3 /t 2 and k t = (2/7)k n ], particle overlap δ and tangential displacement u ij .
The bulk stress tensor is calculated from the particle force and velocity data. It is decomposed into contributions due to the hydrodynamic interaction and the particle-particle interaction, given by Eqs. 4a and 4b, respectively, In the following discussion, we consider the shear components of the above stress tensors corresponding to the direction of the applied deformation, σ l and σ c , as well as the mean of the diagonal components, namely the "pressures" P l and P c . The hydrodynamic stress σ l is further decomposed in two ways. In the first, we isolate the contributions from normal forces (the squeeze a sq terms) and tangential forces (the shear a sh and pump a pu terms) [13] as σ l normal and σ l tangential respectively. In the second, we isolate contributions from opening and closing particle pairs (pairs for which dh/dt > 0 and dh/dt < 0 respectively), presented as σ l opening and σ l closing . It is noted that σ l opening + σ l closing = σ l and σ l normal + σ l tangential = σ l . Assemblies of 5000 spheres are sufficiently large to achieve system size independence, and bidispersity with diameter ratio 1 : 1.4 prevents crystallisation [18]. Simulation results are ensemble averaged over 20 realizations with different initial particle configurations. We note that although the overlap is exceedingly small, typically of order 10 −7 d in the Stokesian regime [19], it can lead to qualitatively different rheology from that produced using the "ideal" hard-sphere model, as demonstrated later. The present technique produces results at the dense limit (solid volume fraction φ 0.45) closely approximating those that would be obtained by fully resolving the hydrodynamics (e.g. [20]), but assuming particles co-move with fluid at the mean flow level [13], valid for shear flows. We verify this by incorporating an additional drag force, similar to [14], which leads to a negligible increase in the calculated suspension viscosity. The particle assemblies are subjected to rate (γ) controlled simple shear flow in a 3-dimensional periodic domain at constant φ (= 0.54) and Stokes numbers St (= ργd 2 /η f ) < 10 −2 , inhibiting particle inertia. The suspension is first sheared froṁ γt = −8 → −2, reaching steady flow. No shear is applied forγt = −2 → 0. Fromγt = 0, the suspension is sheared in the opposite direction until a new steady state is obtained.

III. STRESS AND MICROSTRUCTURE EVOLUTION
The total stress (σ = σ c + σ l ) evolution, Fig 1(a), is strikingly reminiscent of classical experiments [7,8,21]. Rate-independence is demonstrated by collapsing stress components with the respective steady state total showing total (σ, P ), contact (σ c , P c ) and fluid (σ l , P l ) contributions, each scaled by the steady state total stress (σ,P ). Multiple lines for σ illustrateγ independence. Inset: Same data with logarithmic x-axis; (b) σ l components arising from (i) normal and tangential forces; (ii) opening and closing interparticle gaps. Inset: Evolution of the mean fluid film thickness h and scaled mean neighbouring-particle normal Vn / V n and tangential Vt / V n velocity magnitudes. Embedded: particle-pair configurations corresponding to different times. (c) Evolution of coordination number Zc and shear fabric component A12. Inset: Coordination Zc and surface coordination Z h min evolution, with logarithmic x-axis.
stressσ, for multipleγ. Fig. 1(c) shows the microstructural evolution, characterised by a mechanical coordination number Z c , the mean number of per particle contacts which support a contact stress greater than 10 −6 of the mean steady-state stressP , a surface coordination number Z hmin counting all pairs with h < h min and a fabric tensor [19,22], A = 2/(Z hmax N ) h<hmax n ij n ij − 1 3 I. Under shear flow, particles preferentially align along the compressive axis at 45 • with the shear component of A, |A 12 | = 0.5 representing perfect alignment of all contacts and A 12 = 0 representing perfect isotropy. The evolution of additional microstructure variables is shown in Appendix A, Fig 4.

A. Steady flow and cessation
In the steady state (γt < −2, Fig 1(a)), the contact stress contribution is surprisingly large given the small contact overlaps, representing 60% of the total shear stress. The relative contribution is φ-dependent, e.g., at φ = 0.47, we found σ c ≈ 0.3σ [19]. We find Z c ≈ 1.5 and the shear component of A, A 12 ≈ −0.01, indicating persistent mechanical contacts and an anisotropic network of lubrication films. In this condition, relative particle motions are as illustrated in particle-pair diagram A, Fig 1(b), a configuration that results in the mean relative normal velocity of neighbouring particles V n being smaller than the mean relative tangential velocity V t , highlighted in Fig 1(b) Inset, which gives these quantities scaled by the steady state value of V n , V n . This leads to comparable normal and tangential lubrication forces (and corresponding stress contributions σ l normal and σ l tangential as decomposed in Fig. 1(b)), in spite of the order of magnitude difference expected from their respective 1/h versus ln(1/h) dependence. The average number of particle pairs moving together or apart is equal at steady state, as required to satisfy the constant volume constraint, resulting in constant mean lubrication film thickness h , Fig. 1(b) Inset, equal stress contributions σ l opening and σ l closing , Fig. 1(b), andP dominated by P c , Fig 1(a).
Upon flow cessation (γt = −2), the contact stress relaxes together with the hydrodynamics stress, suggesting that caution should be exercised when interpreting the "instantaneous" stress loss as entirely hydrodynamic in such experiments [1,23]. Correspondingly, Z c drops to zero in the relaxation period, though a small portion of weak contacts relax more slowly due to confinement and fluid overdamping. The shear-induced anisotropic microstructure pertaining to hydrodynamics, however, remains intact throughout the relaxation period, evidenced by constant A 12 and Z hmin (Fig 1(c) and Inset, Z hmin (γt → 0) = Z hmin (γt = 4)), implying the steady state HS can be recovered instantaneously (with opposite sign) upon shear reversal.

B. Shear reversal: micro-and macro-strain responses
Indeed, σ l does resume it's steady state magnitude upon reversal for strains ≤ 10 −4 , Fig. 1(a) Inset. It then surges to a significant peak, around 50% greater than the steady value, at strain 10 −3 , sustaining until about 10 −2 where it starts to subside. Resumption of the steady value followed by a demonstrable peak is also observed for P l over the same strain scale. We attribute this smallstrain surge, the manifestation of a microfragile response [9], to the pulling apart of particle surfaces at the h min (= 10 −3 d) scale due to the new (reversed) load being incompatible with the present microstructural alignment. This is clearly demonstrated by the coincident decrease of Z hmin , which reaches a minimum near 10 −2 . The mechanism is further evidenced by the significantly greater V n than under steady flow, the dominance of σ l normal and σ l opening , and the tensile nature of P l . This is further illustrated by the evolution of the h distribution, given in Appendix B. Relative particle motions during this time are illustrated in particle-pair diagram B, Fig. 1(b). Such microfragile events in, e.g., a dry granular system, would be subtle to detect or difficult to distinguish from the macroscopic process. These events in dense suspensions, nonhydrodynamic in nature, however, lead to the spectacular hydodynamic responses, which have been measured robustly in experiments [7,8]. We note that a microfragile hydrodynamic response is absent in Stokesian Dynamics simulations of shear reversal that strictly inhibit fluid films smaller than 0.01d [24], strengthening the argument for direct surface contacts in addition to hydrodynamics, as a crucial contributor to the rheology observed by [7,8].
The subsequent building up of Z hmin after a strain of 0.01 is coupled to re-orientation of the microstructural anisotropy A 12 , corresponding to macrofragile evolution at a larger strain scale of order unity. The initial subsidence of σ l from its peak untilγt ≈ 0.5 (while A 12 < 0), corresponds to a net opening of lubrication films (see h and σ opening , Fig 1(b)), consistent with the leading 1/h dependence of the lubrication forces, combined with a reduction in V n . At larger strains, a new contact network establishes in the now-compatible, oppositely aligned, compressive direction (evidenced by A 12 > 0) with net repulsive lubrication forces during 0.5 <γt < 2, restoring h to its steady value thereby producing positive P l and a marginally dominant σ l closing . The consequent mean relative particle motion is highlighted in paritcle-pair diagram C , Fig 1(b). Although σ l evolves continuously during this large-scale period, the responsible mechanism therefore switches as the anisotropy changes sign. The stress presented by [24] has a comparable macrofragile evolution, but exhibits nonmonotonic behaviour due to the absence of a microfragile response.
The contact stresses (σ c , P c ) follow a similar macrofragile evolution, their associated microfragile contact breakage having occurred at flow cessation as dis-cussed. The stress evolution is closely correlated with the building up of the mechanical coordination number Z c , which occurs on a similar strain scale as the fabric reorientation described above, as illustrated in Fig 1(c). The separation of scales in the evolution of the contact stress and the hydrodynamic peak ensures dominance of the hydrodynamic stress at small strains after reversal, an assertion made in [7,8], though overlooking the microfragile hydrodynamic response. Combining the increasing σ c with the decreasing σ l atγt > 0.01 gives rise to the nonmonotonic total stress, meaning microfragility in the hydrodynamic response is crucial in capturing the experimental behaviour.
The above analysis sheds light on the two-scale nature of the stress evolution, linked to configurational change at small strains and anisotropy re-orientation at large strains. The importance of particle contacts in achieving the nonmonotonic stress response naturally leads to the question of the sensitivity of the evolution at each scale to particle interactions and surface properties.

IV. ROLE OF PARTICLE PROPERTIES
In order to test the applicability of the above described mechanism to a wide range of particle systems, we address two factors pertaining to well studied suspensions, namely surface roughness and stabilising repulsion. For suspensions of large particles (e.g., d > 10 µm), such as the 40 − 50µm polystyrene spheres suspended in density matched silicon oils studied in [7], surface roughness is perhaps the more relevant factor; for those of small particles (e.g., d < 10 µm) steric or electrostatic repulsion may give well-defined repulsive forces.

A. Surface roughness
Surface roughness is represented numerically by an asperity length scale (by means of h min ) and the friction coefficient µ p . h min contributes to the strain scale of the microfragile HS response and should also affect the HS magnitude. We explore such effects by varying h min between 10 −4 d and 10 −2 d, considering the physical size of surface asperities and bounded numerically by the singularity and the overdamping requirement at the lower and upper limits, respectively. The resulting σ l scaled by the steady state total stressσ, following a reversal atγt = 0, is plotted against strain on a log-linear scale in Fig. 2(a). The strain scale of the microfragile peak decreases rather linearly with decreasing h min in the 0.01d to 0.001d range, but saturates approaching 10 −4 d. We verified that the saturation is not due to inertial effects, but is perhaps due to the nonlinear coupling between particle configuration and dynamics. Decreasing h min also significantly increases the peak magnitude, but only weakly affects the macrofragile HS response, the evolution of which interestingly collapses relative toσ. The surface roughness effect of h min thus controls the microfragile, but not the macrofragile HS response.
In a manner often employed in dry granular studies (see e.g. [25]), we vary the friction coefficient µ p incrementally between 0 and 1, exploring particle surfaces from ideally frictionless to very frictional. The total shear stress evolution is given on a linear x-scale in Fig. 2(b). On the contrary to the effects of varying h min , the microfragile response is largely insensitive to µ p , which is unsurprising given the small scale stress is dominated by attractive lubrication forces. The invariance of the hydrodynamic stress gives strong support to the central role of particle contacts in achieving the very different viscosities observed in such systems. The stresses differ hugely at larger strains, however, indicating that the µ p -dependent contact stress is important during 0.3 <γt < 2, coinciding with recovery of Z c to its steady value. The increase of the contact stress with increasing friction can be understood from the increase of tangential contact forces and the decreased departure from the jamming volume fraction φ c [26], which is known to decrease as friction increases [22,27]. The latter effect is also consistent with the experimental observation that the peak immediately after reversal becomes lower relative to the steady state stress when increasing volume fraction [8]. The interparticle friction thus mainly affects the large scale microstructure and contact stress and hence the macrofragile response. In reality, µ p and h min are probably simultaneously coupled to the surface roughness variation, though the combined effect may be deduced from the present separate analyses, exploiting the marked separation of scales associated with our two-scale description.

B. Surface stabilisation
We next probe the effect of a generic stabilising repulsive force, extending the above analysis to consider particles in the size range d < 10µm. It is assumed, based on previous simulation results [5,6], that a static, short range, normal repulsive potential is sufficient to capture the essence of a stabilising mechanism such as electrostatic repulsion or a grafted polymer hair coating. Enhanced dissipation in the lubrication forces, a phenomena described by [28], is neglected for simplicity. A generic form of the repulsive force model derived by Fredrickson et al. [29] is used, F r = k 1 h 5/4 n ij , where k is some constant that encapsulates (among other things) the chemical properties of the hairs and their density on the surface, essentially quantifying the "strength" of the static repulsion. We apply the same singularity regularisation as in the lubrication model, and the same values h min , h max . Coupling to the mechanical contact model is as before. The total shear stress response to reversal is given in Fig. 2(c), for k spanning 2 orders of magnitude (quantified by the relative magnitude of the steady-state repulsiveσ r and contactσ c stresses). As expected for small k, the additional static repulsion is insufficient to separate particles, so the stress response closely resembles that for the base case in Fig. 1(a). As k (orσ r /σ c ) is increased, we note that while the large strain scale for the evolution appears to be unchanged, the steady value of σ decreases. This is attributed to increasing inhibition to mechanical contacts (for which h < 0) as the repulsion becomes stronger. We note that this trend is valid whenσ c is of comparable magnitude toσ r . For very largeσ r /σ c , an opposite trend is observed [6] due to a shear thinning mechanism-the polymer hair length can begin to contribute to an effectively larger total particle diameter, leading to a higher effective volume fraction and therefore a higher shear stress, as explained in detail by [6] and references therein. This then leads to shear thinning behaviour with reducing k, rather than with increasing k as we observe here. For small strains after reversal, we observe a marked loss of the microfragile stress peak as k is increased.
To further understand this loss, we present the full evolution of shear stress contributions for large k and a flow cessation period sufficient to relax to steady state, Fig. 3(a), with h increasing from around 0.06d to 0.08d and Z hmin decreasing modestly (middle Inset). To further characterise the relaxation period, we provide the associated h distributions in Appendix B, under the action of the repulsive potential. Upon reversal, some remaining h min contacts are opened, resulting in a Z hmin decrease and a HS response over a 10 −3 strain scale (middle and right Insets, respectively), consistent with the microfragile response in Fig. 1. In this sense, a microfragile HS response still occurs; though it starts from a "loosened" microstructure, producing a HS lower than its steady value, rather than the surge noted previously. Following this reasoning, a HS peak would be recovered if the relaxation period were shortened sufficiently to disallow any increase in h . We verify this in Fig 3(b) using a very short relaxation periodγt = −0.01 → 0. A HS surge of about 100% of its steady value is observed, although it does not result in an appreciable peak in the total stress since the HS contribution is small. The repulsion and contact stresses similarly follow a macrofragile evolution. This again creates a strain window for the HS to be measured separately from other components. In short, the repulsive force magnitude together with the associated scales provides extra control over particle configurations and hence the stress response. The two-scale evolution concept is, however, still robustly helpful in understanding this more complicated behaviour.

V. CONCLUDING REMARKS
We provide a robust characterisation of a two-scale response to shear reversal in dense suspensions, that is highly reminiscent of the micro-versus macro-fragility proposed by Cates et al. [9]. Links are established between stress responses at small and large strains with microfragile contact breakage and macrofragile microstructural (re-)building respectively, resolving the hitherto unexplained nonmonotonic stress evolution following shear reversal. This substantiates the emerging understanding about the importance of particle contacts in suspension rheology -they not only provide a significant contact stress at steady state, but also give rise to a pronounced small strain transient hydrodynamic response. This understanding provides a sound theoretical framework from which to formulate constitutive models with appropriate two-scale characteristics, and previous attempts at such models [30] might be revised to correctly link the stress and microstructure at each scale. The evidence that different microstructural features control the contact and hydrodynamic stresses respectively and in an analogous way to that in dense granular flows [22], supports further unification of dense suspension and granular rheology extending from steady [31] to unsteady state. The findings on surface features and interactions also open doors to either devising new experiments and protocols, e.g., varying relaxation time, to characterise particle surface properties and stress contributions; or designing new particles, e.g., with different grafted polymer hairs, to realise certain desired rheological properties. Additional microscopic quantities to back up the findings reported in Figure 1(a) of the main article are presented in Figure 4 of this document. Z is a lubrication contact number, counting all pairs with h < h max . All such pairs contribute to the fabric component A 12 . A c12 is the shear component of the mechanical fabric tensor, which we define as A c = 2/(Z c N ) h<0 n ij n ij − 1 3 I, and omitting those pairs which support a contact stress less than 10 −6 of the mean steady-state stressP , consistent with the definition of Z c given in the main article.

Appendix B: Distributions of h
We plot the distribution of the particle-particle separation length h for the simulations in Fig. 1(a) and Fig. 3 of the main article. It is noted that in the case with significant polymer hair repulsion, there remains a peak in P DF (h) at very small h. This is consistent with the corresponding stress evolution, which demonstrates that there is still a non-negligible contribution from direct particle-particle contacts, σ c . In addition, we find that the peak in P DF (h) at small h remains even after the relaxation period. We attribute this somewhat counter-intuitive finding to the repulsive force magnitude and cut-off scales and confinement effects. It is the subsequent opening of these remaining small h particle pairs that is responsible for the very rapid evolution of σ l reported in Fig. 3 of the main article.