Consititutive modeling of time-dependent flows in frictional suspensions

This paper summarizes recent joint work towards a constitutive modelling framework for dense granular suspensions. The aim is to create a time-dependent, tensorial theory that can implement the physics described in steady state by the Wyart-Cates model. This model of shear thickening suspensions supposes that lubrication films break above a characteristic normal force so that frictional contact forces come into play: the resulting non-sliding constraints can be enough to rigidify a system that would flow freely at lower stresses [1]. Implementing this idea for time-dependent flows requires the introduction of new concepts including a configuration-dependent ‘jamming coordinate’, alongside a decomposition of the velocity gradient tensor into compressive and extensional components which then enter the evolution equation for particle contacts in distinct ways. The resulting approach [2, 3] is qualitatively successful in addressing (i) the collapse of stress during flow reversal in shear flow, and (ii) the ability of transverse oscillatory flows to unjam the system. However there is much work required to refine this approach towards quantitative accuracy, by incorporating more of the physics of contact evolution under flow as determined by close interrogation of particle-based simulations.


Introduction
Recent years have seen a shift in our understanding of the rheology of dense suspension of hard, non-Brownian particles in a Newtonian solvent. Naive application of Stokesian physics would imply that such particles, if spherical and governed by a no-slip boundary condition, never come into direct contact because the force needed to create a finite interparticle normal velocity diverges as the fluid film between them becomes thin [4]. However, this is an approximation too far: real particles have aspherities and a finite (nanometre scale) slip length, and therefore suspensions do contain direct interparticle contacts. In a Stokesian system, sudden flow reversal changes the sign but not the magnitude of the shear stress. In reality, the magnitude drops discontinuously before slowly recovering to steady state, because the direct contact forces that are present carry compressive but not tensile forces (unlike a lubrication film). This was demonstrated long ago in the experiments of Gadala-Maria and Acrivos [5] on flow reversal, and recently confirmed in forensic detail by Lin et. al. [6]. Although it is possible to recover some aspects of the Stokesian picture by replacing the word 'contact' in the above discussion with a more general concept involving intact lubrication films and explicit aspherities [7], we stick with the simpler picture of direct particle contacts.
Once these contacts are present, one has essentially a wet granular system and friction enters the stage with a central role. To people familiar with granular systems this may appear obvious, yet the idea that friction matters has only recently achieved consensus in the dense suspension community through a combination of numerical (e.g. [8]), experimental (e.g. [6,9]) and theoretical (e.g. [1]) studies that are by now already too numerous for a full survey.

Wyart-Cates Theory
In the absence of friction, or if the interparticle friction is linear (i.e., Coulombic), dense hard-core suspensions without Brownian motion or inertia are 'quasi-Newtonian'. In simple shear flow all stresses are then linear in the strain-rateγ ≡ dγ/dt. This includes normal stresses, hence the 'quasi': such stresses vanish only in the pure Stokes limit to give a truly Newtonian fluid. The quasi-Newtonian outcome is readily established by dimensional analysis given that the physics of the problem does not admit an intrinsic time scale (no Brownian motion or inertia) or stress scale (hard-core forces are scale-free). Now suppose that there is a soft-core repulsion which is overcome at large normal interparticle forces. One can then consider two quasi-Newtonian branches, a frictionless one occupied at low stress σ ≤ b 1 σ * , say, and a frictional one at high stress, σ ≥ b 2 σ * where b 1 < b 2 and both are of order unity. Here σ * is stress scale set by the softcore forces, the particle volume fraction φ, and the mean particle radius. Crucially, as the number density increases, the viscosity η 2 (φ) of the frictional branch diverges at a jamming point φ J 2 that is lower than the jamming point φ J 1 of the frictionless branch whose viscosity is η 1 (φ). Hence for φ approaching φ J 2 from below, η 2 η 1 . Accordingly, the underlying steady-state flow curve has an unavoidably sigmoidal shape as shown in figure 1 [1]. The region of negative slope is unstable so that at some stress at or below b 1 σ * there will be a discontinuous jump in the observed stress σ(γ) from low to high values. This is discontinuous shear thickening (DST). Backing off the density so the flow curve is no longer sigmoidal, one still expects a very large, but now smooth increase in viscosity on traversing the crossover, giving continuous shear thickening (CST).  [1]. For densities just below the frictional jamming point φ J 1 , the slope of the upper (frictional) branch of the flow curve is vastly larger than that of the lower (frictionless) branch. Thus any reasonable stress-dependent friction law will create a sigmoidal curve.
A simple quantification of these ideas is to choose η = η s ν(φ J − φ) −2 with η s the solvent viscosity and ν some constant, and then set Here f (Π) is the fraction of contacts at which friction is mobilized enough to enforce a rolling rather than a sliding motion. This is written here as a function of the particle pressure Π = −TrΣ/3, with Σ the particle stress tensor and Π * of order σ * , although a similar dependence on the shear stress σ = Σ xy would be sufficient for pure shear flows. The Wyart-Cates theory captures much of the phenomenology of shear thickening in dense suspensions [10] and it has also been extended to allow for Brownian motion [11], interparticle adhesion [12] and other effects. However a fundamental limitation is that it is a steadystate theory only.

Gillissen-Wilson Theory
Let us set aside the physics of shear thickening, and consider a quasi-Newtonian dense suspension belonging to (say) the frictional branch -effectively setting σ * = 0 in figure 1. How can one construct a constitutive model of the time-dependent rheology of such a material? If this question is answered, we can hope to include the concept of stress-dependent friction afterwards. And, of course, there are many dense suspensions that are indeed quasi-Newtonian rather than shear thickening, for which such a model can be directly useful [9].
The problem is that a direct particle contact, unlike a lubrication film, is a highly nonlinear object: it supports compressive stress but not tensile ones. This issue is of course familiar to those working on granular matter. A key insight, however, is that this nonlinearity is so extreme it can be treated like a switch [13]. At the level of macroscopic constitutive equations, this idea is encoded by a decomposition of the macroscopic rate of strain tensor E = (L + L T )/2, where L i j = ∂ j v i (with v the fluid velocity) is the velocity gradient tensor, into its compressive and extensile components, E = E c + E e . These have distinct effects on the evolution of the particle contact distri-bution, which we describe by a microstructure (fabric) tensor nn created from the vectors n between neighbours. These are taken to be interparticle vectors that lie within some coarse-graining shell, thin compared to the particle radius a and thick compared to the range of any direct interparticle force F(h) with h the interparticle separation.
Specifically, the compressive strain rate E c advects, into the coarse-graining shell, previously non-contacting particle pairs whose separation-vector distribution can be modelled as isotropic. It imports preferentially contacts along the compressive axis (or axes). Conversely the extensional strain rate E e advects existing particle pairs, whose separation vectors are anisotropically distributed, out of the coarse graining shell. It exports preferentially those oriented along the extensional axis (or axes).
Combining these ideas with a standard advective model for the n vectors, Gillissen and Wilson were thereby led to the following equation for the fabric tensor [13,14]: with δ the unit tensor. (The precise prefactor of the φ term is inessential and can be absorbed by a rescaling [2].) To complete their model, Gillissen and Wilson adopted a standard Hinch-Leal closure [15] to express nnnn in terms of nn -a choice inspired by the fact that the fabric should be almost isotropic close to jamming [13,14]. Finally they proposed an instantaneous dependence of the particle stress on microstructure and strain rate as In this equation the α term encodes lubrication forces and is therefore linear in E. The χ term represents direct interparticle forces comprising any soft repulsion F(h), hardcore terms, and friction; this neglects tangential forces which are known to be subdominant when calculating the stress [16]. (Note that this choice does not preclude a crucial role for tangential forces, such as friction, in controlling the jamming point.) Because E c and E e interchange on flow reversal, this contact term can capture the sudden loss in stress when compressive contacts suddenly open up. In contrast the lubrication term changes sign at fixed magnitude as required for purely Stokesian physics.

Shear Thickening Suspensions
A marriage of the Wyart-Cates and the Gillissen-Wilson theory was proposed in [2]. To achieve this one should allow the stress parameter χ in (3) to depend on f , the fraction of contacts that are frictional. In this context, we can retain the choice of f (φ) in (1) but can no longer pretend that the viscosity depends only on φ − φ J ( f ) as assumed there. Such a form makes sense for steady states but in unsteady flow the structure undergoes changes at fixed volume fraction and the viscosity clearly depends on these changes. We thus need a scalar function of nn (which can also involve the flow) to describe the distance from jamming in the way φ − φ J ( f ) does for steady states. We call this the 'jamming coordinate'. Proximity to jamming is linked to the number of coarse-grained contacts that are close to becoming frictional. Such contacts are the subset that lie within the range of the soft-core interactions; most of these are oriented along the compression axis or axes. In [2] we proposed on this basis the following jamming coordinate: and showed that, in particle simulations, ξ evolves similarly to a coordination number Z that counts direct (h < ) contacts only. (The latter might be an equally good choice for the jamming coordinate, but it cannot help us close the description at coarse-grained level until approximated by something like ξ.) Replacing (1) with the extremal jamming coordinates ξ J 1,2 can be found from φ J 1,2 by solution of the fabric evolution equations [2].
To complete the picture, we replace the viscosity law with a similar divergence, on the approach to jamming, of the non-hydrodynamic stress coefficient in (3): Here the exponent −2 is retained for simplicity, with the empirical support of particle-based simulations [2].

Results
For reversing shear flows (L i j = ±γδ iy δ jx ) we define a reduced viscosity η r = Σ xy /(γη s ). Sample results from the model are shown in figure 2 alongside corresponding results for particle-based simulations at various reduced shear ratesγ r ≡γη s /Π * . These results use α = 120 and χ 0 = 2.4 which are fitted using steady-state flow curve data. Those fits are done with β = 50 which is in turn fitted to the simulated reversal data in figure 2 itself, to match the observed strain scale for recovery of the stress after reversal. Other parameters are φ, Π * , and φ J 1,2 , with the latter fixing χ J 1,2 via matching of the steady-state solutions (see above). These are already present in the Wyart-Cates theory and since our χ 0 replaces its ν parameter, the extension from steady shear to arbitrary time-dependent flow requires only two new parameters, α and β. Given this, the agreement between the constitutive model and the simulation data is good. However, various discrepancies for other quantities, particularly normal stress differences, are found , as described in [2].
A topical aspect of time-dependent flow in shearthickening suspensions concerns the imposition of a transverse oscillatory flow on top of a steady shear. It is known that the oscillatory component can mitigate or prevent shear thickening, essentially because the compressive contacts leading to large frictional forces are opened by the transverse oscillation, effectively moving the system further from the jamming point than it otherwise  [2]. Lower panel: corresponding data from particle-based simulations. See [2] for further details.
would be [17,18]. Interestingly a transverse oscillation is much more effective than superposing an oscillation in the same direction as the main flow. The unjamming effect is strongly dependent on friction, and indeed barely detectable when non-frictional materials (σ * → ∞) are simulated; but there is no requirement that the material be shear-thickening [18]. The transverse vibration, by opening contacts, is broadly similar in its effects to a reduction in the volume fraction φ.
The model of [2], as outlined above, is used in [3] to address such steady shear flows with superposed transverse oscillations. That paper also provides a more microscopic derivation of the governing equations and introduces slight technical modifications, without changing the conceptual basis of the model. For instance α, which can be treated as constant in the neighbourhood of φ J 2 , is replaced with α 0 (1 − φ/φ J 1 ) 2 . This captures the growth of lubrication forces near φ J 1 which, with transverse shear applied, becomes a more accessible region of parameter space (since jamming at φ 2 is now avoided) even in systems that are frictional from the outset (σ * → 0). This setting offers a stringent test of the new modelling framework because there is no further parameter fitting available and one is testing a model developed to describe flows of fixed orientation but variable rate (reversal of simple shear) by applying to flow where the extensional and compression axes are fully time dependent. The results are somewhat mixed. Figure 3 shows the reduced viscosity η r as a function ofγ ⊥ = ωγ/γ (where ω and γ refer to the transverse flow andγ to the steady shear) alongside its contact and hydrodynamic contributions for both the constitutive model and the particle-based simulations.  Figure 3. Reduced viscosity as a function of nondimensional transverse strain rate from the constitutive model of [3]. Full curve: total value; dashed curve, contact contribution; dotted curve, hydrodynamic contribution. Circles, squares and inverted triangles: the same quantities measured in particle-based simulations. See [3] for further details.

Discussion
The results shown in figure 2 for flow reversal are most encouraging, but those for steady normal stress differences (see [2]), like the data from [3] in figure 3 for transverse oscillations, show qualitative agreement at best when confronted with particle-based simulation data. (Such data, as emphasized in [19], offers stringent tests for constitutive modelling approaches to dense suspensions in general.) Discrepancies reported in [2] include poor prediction of normal stress differences which were attributed to a systematic over-prediction of microstructural anisotropy for any given flow conditions. As reported in [3] the model accounts well for the decrease in contacts on increasing oscillation frequency ω, and for both the build-up of transverse shear stress under transient conditions and its phase relation to the transverse shear rate. On the other hand, figure 3 shows a qualitatively incorrect prediction that the contact contribution to the viscosity collapses to negligible levels at highγ ⊥ . In practice it does fall substantially but continues to dominate the lubrication contribution. An encouraging element of the work summarized above is the way that close interrogation of microstructural data from simulations can inform and constrain the continuum modelling. By looking more closely at the birth and death statistics of contacts, among other quantities, one can hope to improve upon the modelling assumptions laid out above without either altering their basic structure or introducing a proliferation of uncontrolled parameters. Among obvious targets for improvement are the absence of any frictional term in the microstructural evolution equation (2); the simplified treatment of contact birth and death processes reflected in that equation; the ansatz of a near-isotropic (Hinch-Leal [15]) closure; and the specific definition of the jamming coordinate ξ (4) in terms of the fabric tensor and the compressional flow.