Discontinuous shear thickening of dense suspensions under confining pressure

TL;DR

This study uses 2D simulations to analyze how confining pressure influences discontinuous shear thickening in dense suspensions, highlighting the role of frictional contacts and shear protocols.

Contribution

It introduces a simple model linking frictional contact fraction to viscosity and jamming, revealing pressure-dependent shear thickening mechanisms.

Findings

Viscosity diverges at a packing fraction controlled by frictional contacts.

A relationship between frictional contacts, viscous number, and pressure ratio is established.

Viscosity curves depend on shear protocol, with potential negative compressibility.

Abstract

We use 2D numerical simulations to study dense suspensions of non-Brownian hard particles using the Critical Load Model (CLM) under constant confining pressures. This simple model shows discontinuous shear thickening (DST) as the tangential forces get activated upon increased shear stresses. By parameterizing a simple binary system of frictional and non-frictional particles of different proportions we show that the jamming packing fraction, at which the viscosity diverges, is controlled by the fraction of frictional contacts. The viscosity of dense suspensions can thereby be expressed as a function of the fraction of frictional contacts as well as the packing fraction of solid particles. In addition, we show that there exists a simple relationship between the fraction of frictional contacts and the two control parameters (under confining pressure): the viscous number J and the ratio between the repulsive barrier force and confining pressure. Under confining pressures the viscosity curves are found to depend on the shear protocol, with the possibility of yielding negative dynamic compressibility.

Figures

b) so must $\phi_{\mathrm{c}}\to\phi_{\mathrm{c}}^{\mathrm{f}}$, yielding a negative *dynamic* compressibility.

Equations

$${\mu }_{\mathrm{p}}^{ij}=\{\begin{array}{ll}1,& f_{\mathrm{n}}^{ij}>f_{\mathrm{n}}^{\mathrm{cl}};\\ 0,& f_{\mathrm{n}}^{ij}<f_{\mathrm{n}}^{\mathrm{cl}}.\end{array}$$

$$\eta /{\eta }_0\simeq k({\varphi }_{\mathrm{c}}-\varphi {)}^{-m},$$

$$\phi_{\mathrm{c}}(\chi_{\mathrm{f}})=\phi_{\mathrm{c}}^{\mathrm{f}}g(\chi_{\mathrm{f}})+\phi_{\mathrm{c}}^{\mathrm{nf}}\big{(}1-g(\chi_{\mathrm{f}})\big{)},$$

$$\eta/\eta_{0}=\eta(\phi,\chi_{\mathrm{f}})=k(\chi_{\mathrm{f}})\big{(}\phi_{\mathrm{c}}(\chi_{\mathrm{f}})-\phi\big{)}^{-m(\chi_{\mathrm{f}})}.$$

$$\eta =\eta (f_n^{cl}/Pd,J).$$

$$D=\frac{f_{\mathrm{n}}^{\mathrm{cl}}}{Pd}\Big{(}\frac{J}{J_{0}}\Big{)}^{\alpha},$$

$${\chi }_f=k_1+(1-k_1)\cdot \frac{J^{k_2}}{k_3+J^{k_2}},$$

$$\varphi ={\varphi }_{\mathrm{c}}+k^{\prime }{J}^{0.45},$$

Full text

Dense suspensions of rigid particles are of great importance both from a geological and an industrial point of view, including materials such as mud, quicksand, slurry, tooth-paste, paint, etc. Such suspensions a variety of rheological behaviors and a substantial amount of research have been done trying to facilitate our understanding of these behaviors Mari et al. (2014).

In the simplest model describing non-Brownian suspensions, where particles are assumed to have no inertia or Brownian motion and being immersed in a highly viscous fluid, the suspension’s viscosity will be a sole function of packing fraction $\phi$ and the fluid viscosity Boyer et al. (2011). While the viscosity of dilute suspensions is given by the hydrodynamic stresses on single particles in a shear flow and distant hydrodynamic interactions between pairs Batchelor and Green (1972) the viscosity in dense suspensions is governed by viscous dissipation due to particle departures from the affine flow due to geometrical constraints Andreotii et al. (2012); Lerner et al. (2012); Trulsson et al. (2012). The path of departures will increase as the packing fraction increases and finally diverge at the jamming transition. While this is a good model for most of suspensions, it does neither include the possibility of shear thinning nor shear thickening. The latter phenomenon has lately acquired an increasing interest in basic research. In these suspensions the viscosity increases as a function of shear rate Barnes (1989). This increase is labelled as continuous or discontinuous depending on how significant and how sharp the increase is. Transition with significant and sharp viscosity increase is defined as discontinuous shear thickening (DST). DST is usually observed when packing fraction of the suspension exceeds a certain threshold while lower packing fraction gives a more continuous shear thickening Seto et al. (2013). Several scenarios have been proposed to explain shear thickening behavior. It is known that the inertia of the particles Bagnold (1954) and/or fluid Picano et al. (2013) gives a continuous shear thickening above certain shear rate, usually described by the Stokes and/or Reynolds number(s). However, other scenarios are also proposed to explain shear thickening with inertia being subdominant.

Hydroclustering has been a dominating explanation for a couple of decades, where shear thickening is driven by clustering of particles upon increasing shear rate leading to effectively larger quasi-particles. These transient hydroclusters have little internal re-arrangement but a large collective rotational and translational movement which leads to a high dissipation, hence increased viscosity Brady and Bossis (1988); Farr et al. (1997). Although hydro-clustering works well in describing weak shear thickening, simulations based on this phenomenon have not been able to produce the dramatic viscosity increase observed in discontinuous shear thickening Brown and Jaeger (2014). Another alternative explanation relies on the onset of tangential force with increasing shear rate or, equivalently increasing stresses Fernandez et al. (2013); Seto et al. (2013). It is now well-studied that the rheology of frictional and non-frictional dense suspension differ in the jamming packing fraction and the viscosity at the same shear-rate, where frictional suspensions diverge at lower packing fractions and have higher viscosities at the same packing fraction compared to non-frictional suspensions Mari et al. (2014). The former results from a counting argument where extra frictional forces can mechanically stabilise the packings with fewer contacts per particle, hence lower packing fractions. These extra tangential forces also increase the dissipation and hence the viscosity.

In the frictional explanation of DST one goes from a non-frictional to a frictional suspension as the frictional contacts are ”activated”. This activation is initiated once the particles have enough energy to overcome a certain repulsive barrier that protect the particles from being in contact and feeling the frictional forces. Such repulsive and lubricative forces could be both electrostatic and/or steric, e.g. polymer brushes. Several recent papers have reported simulation results based on this explanation/model which produce DST behavior. Although there has been a lot of attention following this explanation/model, most works are done under constant packing fraction Wyart and Cates (2014); Seto et al. (2013); Mari et al. (2015); Ness and Sun (2012) and very little is known about how these suspensions behave under constant confining pressure, which might be a more relevant boundary condition for geological processes or non-planar shear cells A.Fall et al. (2015).

In our work, we explore behavior of dense suspensions under confining pressure based on the scenario of ”activation” of tangential force. We focus on non-Brownian suspensions in highly viscous fluids. Besides of finding relations between viscosity, packing fraction and fraction of frictional contact, we also explore different shear protocols yielding different viscosity versus packing fraction curves.

Simulations are run in 2D with a constant number of disks, $N_{\mathrm{p}}=948$, with average diameter $d$ and polydispersity $\pm 50\%$ to avoid crystallisation. Periodic boundary condition is applied along the $x-$direction. The cell size is approximately $47d$ in the $x$-direction (along the shear) but is allowed to vary in the other direction (typically $\sim 20d$). The disks are confined between two rough walls and immersed in a highly viscous fluid with viscosity $\eta_{0}$. The disk dynamics is hence strictly overdamped without Brownian noise which gives us a set of equations of force and torque balance for each disk $i$. The force balance is given as $\vec{F}^{\mathrm{ext}}_{i}+\vec{F}^{\mathrm{v}}_{i}-\sum_{j}{\vec{f}_{ij}}=0$, where $F^{\mathrm{ext}}$ is the external force from walls, $\vec{F}^{\mathrm{v}}$ is the viscous drag force, and $f_{ij}$ is the contact force exerted on disk $j$ by disk $i$. The drag force is given by Stokes drag $\vec{F}^{\mathrm{v}}_{i}=-\frac{3\pi\eta_{0}}{1-\phi_{0}}(\vec{V}_{i}-\vec{V}_{i}^{\mathrm{a}})$, where $\vec{V}_{i}^{\mathrm{a}}=\dot{\gamma}y\hat{x}$ is the affine fluid velocity at a shear-rate $\dot{\gamma}$ and $\phi_{0}=0.76$. The torque balance is $\tau^{\mathrm{ext}}_{i}+\tau^{\mathrm{v}}_{i}-\sum_{j}{\tau_{ij}}=0$. The force between two disks consist of two components: normal force $\vec{f}_{\mathrm{n}}$ and tangential force $\vec{f}_{\mathrm{t}}$. Here, we assume elastic stiff disks in our simulation so the normal force is given by a linear force, $f_{\mathrm{n}}^{ij}=k_{\mathrm{n}}\delta_{\mathrm{n}}^{ij}$, where $\delta_{\mathrm{n}}^{ij}$ is the overlap between two disks and $k_{\mathrm{n}}$ is the normal spring constant. The tangential force is obtained in a similar way but using tangential spring and a tangential spring constant $k_{\mathrm{t}}=0.5k_{\mathrm{n}}$. The relation between normal component and tangential component is constrained by Coulomb friction, $|{f_{\mathrm{t}}}|\leq\mu_{\mathrm{p}}|{f_{\mathrm{n}}}|$, where $\mu_{p}$ is the friction coefficient (for more information see e.g. Trulsson et al. (2017)).

We use the Critical Load Model (CLM) Mari et al. (2014); Seto et al. (2013); Ness and Sun (2012) which generates DST behavior. This model describes transition between lubrication (non-frictional) and frictional contacts. We first define a critical normal force $f_{\mathrm{n}}^{\mathrm{cl}}$ which represents the magnitude of the repulsive force between two disks at which frictional contacts set in. The friction coefficient between disk $i$ and $j$ is then determined by comparing the normal force between two disks with $f_{\mathrm{n}}^{\mathrm{cl}}$,

[TABLE]

Accordingly, a contact is defined as frictional when $\mu_{\mathrm{p}}^{ij}=1$; otherwise, it is frictionless. We define the fraction of frictional contact $\chi_{\mathrm{f}}$ as the ratio between number of frictional contacts and the number of total contact. In order to study the influence of $\chi_{\mathrm{f}}$ on viscosity of the suspension, we also run simulation with constant $\chi_{\mathrm{f}}$, where disks have a binary distribution of friction coefficient. On this condition, we set certain fraction of disks, $N_{\mathrm{f}}$, as frictional, i.e. $\mu_{\mathrm{p}}^{i}=1$, while the others are frictionless, i.e. $\mu_{\mathrm{p}}^{i}=0$. The friction coefficient between disk $i$ and $j$ is defined as $\mu_{\mathrm{p}}^{ij}=\sqrt{\mu_{\mathrm{p}}^{i}\mu_{\mathrm{p}}^{j}}$ (i.e. only between two frictional disks can one have frictional contacts). In this way, it is easy to control the fraction of frictional contacts $\chi_{\mathrm{f}}\simeq{N_{\mathrm{f}}^{2}}$.

The simulations are shear rate-controlled, by imposing constant velocities of the walls, either at constant packing fraction or constant pressure. The latter is fulfilled by imposing a constant pressure $P$ to both walls, allowing the shear cell to dilate or compress during the shearing. The disk velocities follow well the linear profile, with slope $\dot{\gamma}$, set by the viscous fluid. All the measurements are done considering only the central part of the shear cell, excluding five layers of disks close to each wall. In our simulation, we vary the viscous number $J=\eta_{0}\dot{\gamma}/P$ by changing either the confining pressure $P$ or the shear rate $\dot{\gamma}$ and measure how the viscosity $\eta=\sigma/\dot{\gamma}$, where $\sigma$ is shear stress, and $\chi_{\mathrm{f}}$ changes with $J$. The stiffness of disks is maintained by keeping $k_{\mathrm{n}}/P\approx 2\cdot 10^{3}$ (giving in principal a hard disk behavior). To get more systematic understanding, we run simulations with either $f_{\mathrm{n}}^{\mathrm{cl}}/(Pd)$ or $f_{\mathrm{n}}^{\mathrm{cl}}/(\eta_{0}\dot{\gamma}d)$ constant corresponding to two different shear protocols; keeping either $P$ constant and varying $\dot{\gamma}$ or keeping $\dot{\gamma}$ constant and varying $P$. We also run two reference simulations; in one case, all contacts are frictional (i.e. $f_{\mathrm{n}}^{\mathrm{cl}}=0$) while in the other case, all contacts are frictionless (i.e. $f_{\mathrm{n}}^{\mathrm{cl}}=\infty$).

To validate that the model used is able to describe DST behavior, we first run the simulation with four constant packing fraction $\phi=0.43,0.54,0.63,0.76$. We plot viscosity as a function of shear rate in normalised units, $\eta/\eta_{0}$ against $\dot{\gamma}\eta_{0}/(f_{\mathrm{n}}^{\mathrm{cl}}d)$ (see Supplementary Information SI ). All the curves clearly show a discontinuous transition in viscosity when the shear rate reaches a critical value $\dot{\gamma_{0}}$. At $\phi=0.43$, the increase is quite small. As $\phi$ increases, the transition becomes more pronounced. This result is consistent with previous reported works Seto et al. (2013); Wyart and Cates (2014).

We then run simulations with a constant confining pressure $P$ and study how the viscosity $\eta$ changes with packing fraction $\phi$. Divergence of dense suspension close to jamming can be estimated as

[TABLE]

where $\phi_{\mathrm{c}}$ is critical packing fraction where jamming transition happens. For a constant $\chi_{\mathrm{f}}$, $\phi$ is a function of $J$ Boyer et al. (2011); $\phi\to\phi_{\mathrm{c}}$ when $J\to 0$. Since frictional contact will influence packing of disks, it is reasonable to think that $\phi_{\mathrm{c}}=\phi_{\mathrm{c}}(\chi_{\mathrm{f}})$ and bound by the two limits $\phi_{\mathrm{c}}(0)=\phi_{\mathrm{c}}^{\mathrm{nf}}$ and $\phi_{\mathrm{c}}(1)=\phi_{\mathrm{c}}^{\mathrm{f}}$, where the subscripts denote either the fully frictional or non-frictional references. All intermediate $\chi_{\mathrm{f}}$ cases can then be expressed as

[TABLE]

where $g(\chi_{f})\in[0,1]$ is a combination function conjectured by Wyart and Cates Wyart and Cates (2014). Furthermore, we assume that $k=k(\chi_{\mathrm{f}})$ and $m=m(\chi_{\mathrm{f}})$ so that

[TABLE]

To test our assumption, we first run simulations with constant $\chi_{\mathrm{f}}$ (binary system). The results are plotted in Fig. 1. We fit the data using Eq. 2 and get a set of $\phi_{\mathrm{c}}$ corresponding to different $\chi_{\mathrm{f}}$. We now fit these data to Eq. 3 with $g(\chi_{\mathrm{f}})=\chi_{\mathrm{f}}^{b_{\phi_{\mathrm{c}}}}$ and $b_{\phi_{\mathrm{c}}}\simeq 2$. We apply same strategy to $k$ and $m$ and get $b_{k}\simeq 0.4$ and $b_{m}\simeq 2.2$ respectively (fits are presented in Supplementary Information SI ). Since we have expressions for $k(\chi_{\mathrm{f}})$, $\phi_{\mathrm{c}}(\chi_{\mathrm{f}})$ and $m(\chi_{\mathrm{f}})$, we now plot Eq. 4 with corresponding $\chi_{\mathrm{f}}$. The plots are presented in Fig. 1 as dashed lines, which shows that Eq. 4 fit well with simulation data. This indicates that our assumption works well to predict viscosity for given $\phi$ and $\chi_{\mathrm{f}}$.

Now we apply the same assumptions to the CLM model. We run simulations by varying $J$ with constant $f_{\mathrm{n}}^{\mathrm{cl}}/(Pd)$ or constant $f_{\mathrm{n}}^{\mathrm{cl}}/(\dot{\gamma}\eta_{0}d)$. Fig. 2 shows both simulation results and plots using Eq. 4 with parameters from the binary system with $\phi$ and $\chi_{\mathrm{f}}$ measured from simulations with astonishing agreement and shows a robustness of these reduced parameters irrespectively of the detailed microscopic structure (see Supplementary Information SI for two snapshots with the same $\phi$ and $\chi_{\mathrm{f}}$). For both shear protocols (keeping constant either $f_{\mathrm{n}}^{\mathrm{cl}}/(Pd)$ or $f_{\mathrm{n}}^{\mathrm{cl}}/(\dot{\gamma}\eta_{0}d)$), curves show a tendency of crossing over from frictional to frictionless behavior as the control parameters are changed. Both $\eta$ and $\phi$ increase monotonically for constant $f_{\mathrm{n}}^{\mathrm{cl}}/(Pd)$. Similar is found controlling $f_{\mathrm{n}}^{\mathrm{cl}}/(\dot{\gamma}\eta_{0}d)$ if the parameter is low enough (typically $\lesssim 300$). On the other hand, for high and constant $f_{\mathrm{n}}^{\mathrm{cl}}/(\dot{\gamma}\eta_{0}d)$, one finds cusp-shape behaviors with an increasing cusp as $f_{\mathrm{n}}^{\mathrm{cl}}/(\dot{\gamma}\eta_{0}d)$ increases. For these cups one finds two (or more) viscosities for each packing fraction (except possibly at the maximum $\phi$). These two states do however have different pressures (assuming fixed $f_{\mathrm{n}}^{\mathrm{cl}}$), with a higher pressure for the high viscosity case. This indicates a range of pressures where one has a negative dynamic compressibility, i.e. $\phi$ decreases with increasing $P$. This appears when $\phi$ is above $\phi_{\mathrm{c}}^{\mathrm{f}}$ since then $\phi_{\mathrm{c}}^{\mathrm{f}}<\phi<\phi_{\mathrm{c}}\leq\phi_{\mathrm{c}}^{\mathrm{nf}}$; as $\chi_{\mathrm{f}}$ increases with $P$ (see Fig. 2b) so must $\phi_{\mathrm{c}}\to\phi_{\mathrm{c}}^{\mathrm{f}}$, yielding a negative dynamic compressibility.

Unlike in our reference binary system, $\chi_{\mathrm{f}}$ can no longer be considered as control parameter in the CLM model. Instead it will vary with $J$ and the control parameter as determined by the shear protocol. It is reasonable to presume that $\chi_{\mathrm{f}}$ is related to $f_{\mathrm{n}}^{\mathrm{cl}}/\langle f_{\mathrm{n}}\rangle$ according to the presumption of the CLM model, where $\langle f_{\mathrm{n}}\rangle$ is the average of normal force between two disks. Given that $P\sim\langle f_{\mathrm{n}}\rangle$, we assume $\chi_{\mathrm{f}}=\chi_{\mathrm{f}}(f_{\mathrm{n}}^{\mathrm{cl}}/(Pd),J)$. As $\lim_{J\to 0}\phi=\phi_{\mathrm{c}}$ and given that $\phi_{\mathrm{c}}=\phi_{\mathrm{c}}(\chi_{\mathrm{f}})$, we can assume $\phi=\phi(f_{\mathrm{n}}^{\mathrm{cl}}/(Pd),J)$. We obtain expression for $\chi_{\mathrm{f}}$ and $\phi$ by fitting corresponding simulation data (in Supplementary Information SI ). Now we can rewrite Eq. 4 as

[TABLE]

Remembering that $J=\dot{\gamma}\eta_{0}/P$ one can reformulate both control parameters in functions of $J$ and $f_{\mathrm{n}}^{\mathrm{cl}}/(Pd)$: $f_{\mathrm{n}}^{\mathrm{cl}}/(Pd)=(f_{\mathrm{n}}^{\mathrm{cl}}/(Pd)){J^{0}}$ and $f_{\mathrm{n}}^{\mathrm{cl}}/(\dot{\gamma}\eta_{0}d)=(f_{\mathrm{n}}^{\mathrm{cl}}/(Pd)){J^{-1}}$. This shows the generality of Eq. 5. Following the above reasoning we actually find a whole family of shear protocols that can be encoded in the new parameter:

[TABLE]

where $J_{0}$ is a reference point. With a given $\alpha$, we can obtain relation between $J$ and $P$ which can then be put into Eq. 5. Different $\alpha$’s can be divided into three parts around the two previous reference points $-1$ and [math], which correspond to varying $\dot{\gamma}$ or varying $P$. The other $\alpha$’s indicate that both $\dot{\gamma}$ and $P$ are varied simultaneously. From Eq. 6 we see that pressure should be varied as $P\sim\dot{\gamma}^{\alpha/(1+{\alpha})}$ and noticing that $J\sim\dot{\gamma}/P$ we can obtain $J\sim\dot{\gamma}^{1/(1+\alpha)}$. To observe divergence in viscosity (i.e. $J\to 0$ and consequently, $\phi\to\phi_{c}$), there are three possibilities to realize this. In the first case, both $\dot{\gamma}$ and $P$ increases with $P$ increasing more rapidly, which corresponds to $\alpha<-1$. In the second case, $\dot{\gamma}$ decreases while $P$ increases, which corresponds to $-1<\alpha<0$. In the third cases, both $\dot{\gamma}$ and $P$ decreases with $P$ decreasing faster, which corresponds to $\alpha>0$. These combination of $\dot{\gamma}$ and $P$ can be easily mapped to different shear protocol in experiments. Fig. 3 illustrates plots of Eq. 5 with $J_{0}=0.001$ and $D=1$ for various $\alpha$’s. Frictional and frictionless curves are also produced by setting $D=0$ and $\infty$ respectively. The intersection point is controlled by both $J_{0}$ and $D$. Negative $\alpha$ causes $\phi$ to behave non-monotonically while positive $\alpha$ causes $\eta$ to become non-monotonic. Cusped curves become more significant with larger $|\alpha|$.

In this work, we simulate shear thickening of dense suspensions under confining pressure. By plotting $\eta(\phi)$ and fitting the data to the empirical equation $\eta/\eta_{0}=\eta(\phi,\chi_{\mathrm{f}})$, we illustrate that both $\chi_{\mathrm{f}}$ and $\phi$ determine the viscosity. In particular, the jamming point $\phi_{\mathrm{c}}$ is determined by $\chi_{\mathrm{f}}$ and hence the divergence of viscosity. We find different transition paths between the fully frictional and frictionless curves depending on the shear protocol. Negative dynamic compressibility is observed when varying $P$ at $\phi>\phi_{\mathrm{c}}^{\mathrm{f}}$, resulting from a decreasing $\phi_{\mathrm{c}}$ due to increasing $\chi_{\mathrm{f}}$. We further find that both $\chi_{\mathrm{f}}$ and $\phi$ are functions of $f_{\mathrm{n}}^{\mathrm{cl}}/(Pd)$ which indicates that a relation between repulsive forces and confining pressure is a key factor of shear thickening behavior and could be used in continuum modelling. We also introduce a new dimensionless constant $D=\frac{f_{\mathrm{n}}^{\mathrm{cl}}}{Pd}(\frac{J}{J_{0}})^{\alpha}$. Encoded in $\alpha$ is a whole set of different shear protocols, each which generating its own flow curve. This observation could be tested and verified using current experiment protocols for shear-thickening suspensions Boyer et al. (2011).