Discontinuous shear thickening of dense suspensions under confining pressure
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 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.
pacs:
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 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, , with average diameter and polydispersity to avoid crystallisation. Periodic boundary condition is applied along the direction. The cell size is approximately in the -direction (along the shear) but is allowed to vary in the other direction (typically ). The disks are confined between two rough walls and immersed in a highly viscous fluid with viscosity . 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 . The force balance is given as , where is the external force from walls, is the viscous drag force, and is the contact force exerted on disk by disk . The drag force is given by Stokes drag , where is the affine fluid velocity at a shear-rate and . The torque balance is . The force between two disks consist of two components: normal force and tangential force . Here, we assume elastic stiff disks in our simulation so the normal force is given by a linear force, , where is the overlap between two disks and is the normal spring constant. The tangential force is obtained in a similar way but using tangential spring and a tangential spring constant . The relation between normal component and tangential component is constrained by Coulomb friction, , where 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 which represents the magnitude of the repulsive force between two disks at which frictional contacts set in. The friction coefficient between disk and is then determined by comparing the normal force between two disks with ,
(1) |
Accordingly, a contact is defined as frictional when ; otherwise, it is frictionless. We define the fraction of frictional contact as the ratio between number of frictional contacts and the number of total contact. In order to study the influence of on viscosity of the suspension, we also run simulation with constant , where disks have a binary distribution of friction coefficient. On this condition, we set certain fraction of disks, , as frictional, i.e. , while the others are frictionless, i.e. . The friction coefficient between disk and is defined as (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 .
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 to both walls, allowing the shear cell to dilate or compress during the shearing. The disk velocities follow well the linear profile, with slope , 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 by changing either the confining pressure or the shear rate and measure how the viscosity , where is shear stress, and changes with . The stiffness of disks is maintained by keeping (giving in principal a hard disk behavior). To get more systematic understanding, we run simulations with either or constant corresponding to two different shear protocols; keeping either constant and varying or keeping constant and varying . We also run two reference simulations; in one case, all contacts are frictional (i.e. ) while in the other case, all contacts are frictionless (i.e. ).
To validate that the model used is able to describe DST behavior, we first run the simulation with four constant packing fraction . We plot viscosity as a function of shear rate in normalised units, against (see Supplementary Information SI ). All the curves clearly show a discontinuous transition in viscosity when the shear rate reaches a critical value . At , the increase is quite small. As 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 and study how the viscosity changes with packing fraction . Divergence of dense suspension close to jamming can be estimated as
(2) |
where is critical packing fraction where jamming transition happens. For a constant , is a function of Boyer et al. (2011); when . Since frictional contact will influence packing of disks, it is reasonable to think that and bound by the two limits and , where the subscripts denote either the fully frictional or non-frictional references. All intermediate cases can then be expressed as
(3) |
where is a combination function conjectured by Wyart and Cates Wyart and Cates (2014). Furthermore, we assume that and so that
(4) |
To test our assumption, we first run simulations with constant (binary system). The results are plotted in Fig. 1. We fit the data using Eq. 2 and get a set of corresponding to different . We now fit these data to Eq. 3 with and . We apply same strategy to and and get and respectively (fits are presented in Supplementary Information SI ). Since we have expressions for , and , we now plot Eq. 4 with corresponding . 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 and .
Now we apply the same assumptions to the CLM model. We run simulations by varying with constant or constant . Fig. 2 shows both simulation results and plots using Eq. 4 with parameters from the binary system with and 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 and ). For both shear protocols (keeping constant either or ), curves show a tendency of crossing over from frictional to frictionless behavior as the control parameters are changed. Both and increase monotonically for constant . Similar is found controlling if the parameter is low enough (typically ). On the other hand, for high and constant , one finds cusp-shape behaviors with an increasing cusp as increases. For these cups one finds two (or more) viscosities for each packing fraction (except possibly at the maximum ). These two states do however have different pressures (assuming fixed ), with a higher pressure for the high viscosity case. This indicates a range of pressures where one has a negative dynamic compressibility, i.e. decreases with increasing . This appears when is above since then ; as increases with (see Fig. 2b) so must , yielding a negative dynamic compressibility.
Unlike in our reference binary system, can no longer be considered as control parameter in the CLM model. Instead it will vary with and the control parameter as determined by the shear protocol. It is reasonable to presume that is related to according to the presumption of the CLM model, where is the average of normal force between two disks. Given that , we assume . As and given that , we can assume . We obtain expression for and by fitting corresponding simulation data (in Supplementary Information SI ). Now we can rewrite Eq. 4 as
(5) |
Remembering that one can reformulate both control parameters in functions of and : and . 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:
(6) |
where is a reference point. With a given , we can obtain relation between and which can then be put into Eq. 5. Different ’s can be divided into three parts around the two previous reference points and , which correspond to varying or varying . The other ’s indicate that both and are varied simultaneously. From Eq. 6 we see that pressure should be varied as and noticing that we can obtain . To observe divergence in viscosity (i.e. and consequently, ), there are three possibilities to realize this. In the first case, both and increases with increasing more rapidly, which corresponds to . In the second case, decreases while increases, which corresponds to . In the third cases, both and decreases with decreasing faster, which corresponds to . These combination of and can be easily mapped to different shear protocol in experiments. Fig. 3 illustrates plots of Eq. 5 with and for various ’s. Frictional and frictionless curves are also produced by setting and respectively. The intersection point is controlled by both and . Negative causes to behave non-monotonically while positive causes to become non-monotonic. Cusped curves become more significant with larger .
In this work, we simulate shear thickening of dense suspensions under confining pressure. By plotting and fitting the data to the empirical equation , we illustrate that both and determine the viscosity.
In particular, the jamming point is determined by 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 at , resulting from a decreasing due to increasing . We further find that both and are functions of 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 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).
. Encoded in
Acknowledgements.
M.T. acknowledges financial support by the Swedish Research Council (621-2014-4387) and Crafoord Foundation (20160568). The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the center for scientific and technical computing at Lund University (LUNARC). J.D. acknowledges Jasenko Gavran for initial help and instruction running the binary systems.References
- Mari et al. (2014) R. Mari, R. Seto, J. F. Morris, and M. M. Denn, Journal of Rheology 58 (2014).
- Boyer et al. (2011) F. Boyer, Élisabeth Guazzelli, and O. Poulique, Phys. Rev. Lett. 107, 188301 (2011).
- Batchelor and Green (1972) G. Batchelor and J. Green, J. Fluid Mechanics 56, 401 (1972).
- Andreotii et al. (2012) B. Andreotii, J.-L. Barrat, and C. Heussinger, Phys. Rev. Lett. 109, 105901 (2012).
- Lerner et al. (2012) E. Lerner, G. Düring, and M. Wyart, PNAS 109, 4798 (2012).
- Trulsson et al. (2012) M. Trulsson, B. Andreotti, and P. Claudin, Phys. Rev. Lett. 109, 118305 (2012).
- Barnes (1989) H. Barnes, Journal of Rheology 33 (1989).
- Seto et al. (2013) R. Seto, R. Mari, J. F. Morris, and M. M. Denn, Phys. Rev. Lett. 111 (2013).
- Bagnold (1954) R. Bagnold, Progress of the Royal Society A 225 (1954).
- Picano et al. (2013) F. Picano, W.-P. Breugem, D. Mitra, and L. Brandt, Phys. Rev. Lett. 111, 098302 (2013).
- Brady and Bossis (1988) J. F. Brady and G. Bossis, Annual Review of Fluid Mechanics 20 (1988).
- Farr et al. (1997) R. Farr, J. Melrose, and R. Ball, Physical Review E 55 (1997).
- Brown and Jaeger (2014) E. Brown and H. M. Jaeger, Review on Progress in Physics 77 (2014).
- Fernandez et al. (2013) N. Fernandez, R. Mani, D. Rinaldi, D. Kadau, M. Mosquet, H. Lombois-Burger, J. Cayer-Barrioz, H. J. Herrmann, N. D. Spencer, and L. Isa, Phys. Rev. Lett. 111 (2013).
- Wyart and Cates (2014) M. Wyart and M. Cates, Phys. Rev. Lett. 112, 098302 (2014).
- Mari et al. (2015) R. Mari, R. Seto, J. F. Morris, and M. M. Denn, Phys. Rev. E 91 (2015).
- Ness and Sun (2012) C. Ness and J. Sun, Soft Matter 12 (2012).
- A.Fall et al. (2015) A.Fall, F. Bertrand, D. Hautemayou, C. Mezi re, P. Moucheront, A. Lemaître, and G. Ovarlez, Phys. Rev. Lett. 114, 098301 (2015).
- Trulsson et al. (2017) M. Trulsson, E. DeGiuli, and M. Wyart, Phys. Rev. E 95, 012605 (2017).
- (20) Supplementary Information: available online.
Fig. 4 illustrates the dependence between viscosity and shear rate at various constant ’s.
Fig. 5 shows curves for , , and . The fitting function has been described before. Fig. 6 shows how varies with increasing using different shear protocols. The data in Fig. 6(a) are fit into
(7) |
where , and are all functions of . After getting expressions for , and , we obtain expression for . Plots of for both protocols are illustrated in Fig. 6 as dashed lines. The expression for is obtained with same strategy. The fitting function is
(8) |
where is a function of . The results are shown in Fig. 7 where dashed lines are plots of Eq. 8 and symbols are simulation data. Fig. 8 shows how macroscopic friction varies with using two shear protocols.
Fig. 9 shows microscopic structure of suspensions under shearing using CLM model and binary system (two animations are also attached). Corresponding and are , ; , . Although and for two different systems are quite close to each other, the microscopic structures are rather different as illustrated in Fig. 9. For system using CLM model, disks with frictional contacts form chains or even networks while for binary system, these disks are distributed more separately. This observation shoes the subdominant effect of the chains for viscosity.