Microstructural description of shear-thickening suspensions

Dynamic particle-scale numerical simulations are used to study the variation of microstructure with shear stress during shear thickening in dense non-Brownian suspensions. The microscale information is used to characterize the differences between the shear thickened (frictional) and non-thickened (lubricated, frictionless) states. Here, we focus on the force and contact networks and study the evolution of associated anisotropies with increase in shear stress. The force and contact networks are both more isotropic in the shear-thickened state than in non-thickened state. We also find that both force and structural anisotropies are rate independent for both low and high stress, while they are rate (or stress) dependent for the intermediate stress range where the shear thickening occurs. This behavior is similar to the evolution of viscosity with increasing stress, showing a clear correlation between the microstructure and the macroscopic rheology.


Introduction
Suspensions are dispersions of rigid particles in fluids.The fluid allows numerous possible interactions between the particles, but here we consider only hydrodynamic interactions and short-range or contact forces.Even with this simple system, a rich rheology is found, including shear thinning or shear thickening, normal stress differences, particle migration, shear jamming and yield stress behavior [1][2][3].For moderate volume fractions, the viscosity increases with increase in shear rate (or shear stress) and is termed continuous shear thickening (CST).However, at high volume fraction the viscosity can increase by orders of magnitude at a critical shear rate in what is known as discontinuous shear thickening (DST) [1,3].
The mechanism for DST has been of interest for a number of years, with an order-disorder transition [4] and hydrodynamic clustering [5] postulated.Recent numerical and experimental studies have linked DST to a transition between 'lubricated' low-viscosity and 'frictional' high viscosity states as the shear stress increases [1,2,6,7].The friction requires contact and this scenario thus requires an added repulsive force F 0 to stabilize the material in the lubricated regime.F 0 sets the stress scale F 0 /a 2 for particles of size a for the transition from the lubricated to frictional behavior.We explore features of this lubricatedfrictional transition through numerical simulations.
Suspension rheology is often discussed in the framework of shear-induced microstructure [2,8].For smoothly shearing suspensions, representative of the lubricated state, the pair distribution is strongly anisotropic, with accumulation of pairs near contact in the compressional quadrant [2,8].Other particulate systems, such as dry granular materials, show similar behavior; isotropically generated assemblies develop anisotropic structure under shear, which results in the finite shear strength of the material [9,10].Recent studies [1,2] have shown that, when contact with friction is simulated, the contact network in the shear-thickened state is more isotropic than the one in the non-thickened state; this is seen especially in the normal stress response, as the normal stress differences grow through the DST transition, but less so than the isotropic particle normal stress (particle pressure) [1].Here, we report the evolution of both force and contact network anisotropies, toward the goal of bridging between the microscopic and macroscopic scales in dense suspension flows.

Model system
We simulate a two-dimensional or monolayer suspension of non-Brownian spherical particles.We consider the motion to be inertialess, so that the equation of motion reduces to a force balance between hydrodynamic (F H ) and contact (F C ) forces: where u and r are many-body position and velocity vectors (u ≡ ṙ).Rather than employing a repulsive force to set the stress scale for the transition to the frictional interactions, we embed this scale in the contact force as noted below.The hydrodynamic forces generally consist of two terms: Stokes drag on the particle due to its motion relative to the surrounding fluid, −R FU •(u − u ∞ ), and resistance associated with the bulk straining motion, γR FE : Ê∞ , where γ Ê∞ is the rate of strain (details are elsewhere [1,2]).The hydrodynamic forces scale linearly with shear rate γ, and thus their contribution to the relative viscosity, given by η 0 η r = σ/γ, with σ the shear stress and η 0 the fluid viscosity, has no rate dependence for purely hydrodynamic interactions.Similarly, frictional interactions do not introduce a rate dependence.To incorporate rate dependence we employ a Critical Load Model (CLM) with an inter-particle threshold normal force F CL used to activate friction between particles [1,2].This force scale sets the stress scale F CL /a 2 for the transition between 'frictionless' interactions that take place via a liquid film and those making unlubricated 'frictional' contacts at high stress.CLM can be thought of as a small Debye length limit of the electrostatic repulsion model [2].The interparticle contacts are defined by a linear spring model in both normal and tangential directions, as commonly used in granular physics [1,2,11].The tangential force between two particles is subjected to Coulomb's friction law In this study we use μ p = 1.0.The strain rate and stress are nondimensionalized by γ0 ≡ F CL /6πη 0 a 2 and σ 0 ≡ η 0 γ0 .

Results and discussion
In Fig. 1a, the average steady state shear rate is shown as a function of the applied shear stress for an area fraction φ A = 0.78.The simulations are performed at controlled shear stress σ, and the shear rate is thus variable [12].The error bars represent the standard deviation taken over many samples in a simulation of typically γt = 20 strain units.At small and large shear stresses (σ/σ 0 < 1 and σ/σ 0 ≥ 20), the shear rate is linear with σ, i.e. the suspension viscosity is rate independent, while it shows a non-monotonic S-shaped flow curve for the intermediate values of shear stress.A rate-controlled simulation at this volume fraction would display discontinuous shear thickening, because such S-shaped curves are not accessible under a rate-controlled simulation [6,12].Figure 1b shows the evolution of the direct particle-particle coordination number or number of contacts per particle Z = N C /N as a function of stress.The number of frictional contacts increases with an increase in the shear stress, whereas the number of frictionless contacts decreases.For low stress all contacts are frictionless, but at high stress contacts are predominantly frictional, which is consistent with the evolution of frictional contacts previously reported [2].The number of contacts in the shear-thickened state exceeds the frictionless interactions that were present in the nonthickened state.
We can gain more insight by examining the contact network anisotropy.There are several methods that have been used to quantify the structural anisotropy of a granular packing [9,10].One common method is to consider the probability distribution P(n i j ) of the normal vectors n i j oriented from particle i to particle j.In a two dimensional setting, the unit vector n can be described by its orientation angle θ, measured here relative to the flow direction.Now the anisotropy can be characterized by a probability density function P(θ) of finding a vector aligned along the direction θ.This probability density function is π-periodic, as the contact normal could be defined in either direction.A Fourier expansion with only even terms provides an exact representation of the geometrical anisotropy where a i (with i = 2, 4, 6..) are the anisotropy descriptors along principal direction θ c .When a i = 0, the distribution is isotropic and P(θ) is a constant equal to 1/2π (as indicated by the dashed line in Fig. 2a). Figure 2 displays a polar representation of P(θ) for various values of stress (only a few are shown here for the sake of clarity).Each of these probabilities is normalized by the respective number of interactions in each class.Before we discuss our results, we introduce the nomenclature used.We use 'contact' when particles touch each other (distance between their centers is less than the sum of their radii), and 'interaction' is used otherwise (hydrodynamic forces would fall here).
The generalized network including both contacts and interactions (non-contact hydrodynamic lubrication and direct particle-particle contacts) is shown in Fig. 2a.The network is highly non-uniform, with no preferred orientation at low stress and tending to become isotropic with an increase in stress (for high stress the data almost collapse on the dashed circle).We observe at low stress a large accumulation of pairs aligned along the flow direction, consistent with examination of structures seen in sheared Brownian suspensions [13].
To gain better insight the general network is split into various classes, i.e. lubrication and direct particle-particle contacts (which are further divided into frictional and frictionless contacts).For low stress values, both hydrodynamic and particle-particle contacts have peaks in P(θ) along the flow direction (Figs.2b and 2c), and as the stress increases the magnitude of these peaks decreases.With increasing stress, the lubrication interactions develop the elongation axis as the dominant direction of anisotropy, while the direct particle-particle contacts are found to align along the compressional axis; i.e., the directions of anisotropy in the two types of interactions are orthogonal to each other.Eventually, at high enough values of the imposed stress, the contact network tends to become closer to isotropic.We further categorize the direct particle-particle contacts as frictional and frictionless, as shown in Figs.2d and  2e, respectively.As these polar plots are normalized by the total number of interactions in each class, Figs.2d and 2e should be analyzed together with Fig. 2f, where we show the fraction of frictional contacts f = N f C /N C (this is the ratio between contacts for which the normal force exceeds the critical load force F CL and total number of particleparticle contacts).The frictional contacts are rare at low stress and are concentrated in small force chains along the compression direction, whereas the frictionless contacts behave similarly to the lubrication interactions (being aligned along the flow direction).As the stress increases, the number of frictional contacts increases (Fig. 2f), and at the same time these contacts tend to become more isotropically distributed.Figure 3 shows the contact anisotropy descriptors a 2 and a 4 as functions of shear stress for all, frictional and frictionless contacts.Opposite signs of a 2 and a 4 signify the orthogonal anisotropic directions (a 2 captures P(θ) along the compression direction, a 4 along the elongational axis).Interestingly, a 4 is large for low stress and diminishes with increasing stress, signifying it is required to capture the nonuniform and highly anisotropic distributions at low stress as seen in Figs.2c, 2e, and 2f.We observe that the contact anisotropy (for all and frictionless contacts) is almost a constant for low stress σ < 1 (as the number of frictional contacts for σ = 0.1 is almost zero, we do not have P(θ) for that case).For all the cases shown, the contact anisotropy decreases with increasing stress and once again saturates for high stress σ > 20 the magnitudes of both a 2 and a 4 decrease and tend to saturate to constant values with increase in stress).θ c for all types of contacts nearly coincide and show a slight decrease θ c ∈ [135, 130] (in degrees) with increasing stress.
Next we examine the anisotropy of the force network in the non-thickened (frictionless) and thickened (frictional) states.The angular distribution of normal contact forces between particles can be represented by average force f n (n) as a function of contact normal direction n.This has been previously defined by [9] as f c n , where f c n represents the normal contact force acting at the contact c and N c (θ) represents the number of contacts in S (θ) (set of contacts with direction θ).Similar to the contact normal vector, the normal force could also be defined in either direction, hence the probability density is π-periodic.In the steady state, f n (θ) can be described by a similar function as used for contact orientation distribution in Eq. ( 2): where f 0 is the mean of the magnitude of normal contact forces, b i (with i = 2, 4, 6..) are the anisotropy descriptors of the normal force network, and θ f denotes the direction of the maximum average normal force.In Fig. 4, f n (θ) (for only direct particle-particle contact forces) is displayed in polar coordinates.Once again, these distributions are normalized by the mean force for each class of contacts.We observe that all particle-particle contacts display pronounced force anisotropy, unlike P(θ) (as shown in Figs 2c).The force polar plots do not show any alignment of force chains along the flow direction as previously seen in the case of contact normal polar representation.We also observe that the overall force network tends to become less anisotropic with increasing stress.The polar representations of frictional and frictionless normal forces exhibit highly contrasting features.The anisotropy in the frictional force network increases with increase in the applied shear stress.On the other hand, the force network of frictionless contacts is initially anisotropic and the anisotropy decreases with increase in stress.Figure 5 shows the evolution of force anisotropy parameters b 2 and b 4 with the applied shear stress σ.The overall normal force network is initially anisotropic and its anisotropy decreases and tends toward constant values of b 2 and b 4 with increasing stress.The force anisotropy of frictionless normal forces decreases to zero with σ.At the same time, the anisotropy of frictional normal forces increases sharply initially and then saturates as a function of σ. θ f for the three types of contacts are nearly identical (θ f 3π/4) and are independent of stress.

Conclusions
We have studied shear thickening in dense suspensions by considering the microstructure and its correlation to the bulk rheology.Estimation of anisotropy parameters revealed that the anisotropy of both contact and force networks plateau for low and high stress values and show a rate dependence for intermediate stress values, much as the viscosity does.This shows that the rate dependence in the viscosity can be attributed to both increase in the number of frictional contacts as well as the local rearrangements in contact and force networks.

Figure 1 .
Figure 1.Macroscopic and microscopic characterization.(a) The macroscopic shear rate as a function of stress for φ A = 0.78 with N = 2000 particles.(b) Number of contacts as functions of stress.Squares show all contacts, triangles represent frictionless contacts , while circles represent frictional contacts.

Figure 2 .
Figure 2. Polar representation of the probability density function P(θ) of the contact normal vectors for (a) generalized network, (b) only lubrication interaction, (c) all direct particle-particle contacts, (d) frictional contacts, and (e) frictionless contacts.Different symbols represent simulations with σ/σ 0 as shown in the legend.All the data shown are sampled in the steady state.(f) Fraction of frictional contacts f (σ) as a function of shear stress.As can be seen from Fig. 1, σ = 0.5 corresponds to nonthickened state, while σ ≥ 20 correspond to the thickened state, and σ = 1, 5 belong to the transitional state.In Figs.(a-e) the direction of shear is from right to left at the top and left to right at the bottom.

4 Figure 3 .
Figure 3. Contact anisotropy descriptors a 2 and a 4 as functions of stress for (a) all direct particle-particle contacts, (b) frictional contacts, and (c) frictionless contacts.All the data shown are sampled in the steady state.

Figure 4 . 4 Figure 5 .
Figure 4. Polar representation of the angle-averaged normal forces f n (θ) for (a) all, (b) frictional, and (c) frictionless direct particle-particle contacts for various stress values.All the data shown are sampled in the steady state.