1 Introduction
Suspensions, namely mixtures of solid particles and viscous liquids, can be considered as incompressible fluids as long as the volume fraction
$\unicode[STIX]{x1D719}$
of solid particles is less than a certain value, the jamming point, above which a solid-like behaviour is observed. The behaviour of suspensions is not usually captured by simple Newtonian models. As a primary example of a non-Newtonian effect, the viscosity can vary with the shear rate, exhibiting shear thinning and shear thickening (Laun Reference Laun1984; Barnes Reference Barnes1989; Bender & Wagner Reference Bender and Wagner1996; Guy, Hermes & Poon Reference Guy, Hermes and Poon2015). Moreover, non-vanishing normal stress differences
$N_{1}$
and
$N_{2}$
, further hallmarks of non-Newtonian behaviour, are often observed (Laun Reference Laun1994; Lootens et al.
Reference Lootens, van Damme, Hémar and Hébraud2005; Lee et al.
Reference Lee, Alcoutlabi, Magda, Dibble, Solomon, Shi and McKenna2006; Couturier et al.
Reference Couturier, Boyer, Pouliquen and Guazzelli2011; Dbouk, Lobry & Lemaire Reference Dbouk, Lobry and Lemaire2013; Cwalina & Wagner Reference Cwalina and Wagner2014). Discontinuous shear thickening is a particularly intriguing phenomenon of dense suspensions, and the underlying mechanism has raised significant debate (Brady & Bossis Reference Brady and Bossis1985; Hoffman Reference Hoffman1998; Melrose & Ball Reference Melrose and Ball2004; Fall et al.
Reference Fall, Huang, Bertrand, Ovarlez and Bonn2008; Brown & Jaeger Reference Brown and Jaeger2009). Analysis of the rheology of suspensions is a difficult task since forces of various types act among particles and the system lives mostly far from thermodynamic equilibrium. Particle simulations have been used to explore the microstructure emerging among particles in various flows and to estimate the importance of different interactions. Several particle simulations have recently succeeded in reproducing shear thickening by taking into account direct contact forces (Fernandez et al.
Reference Fernandez, Mani, Rinaldi, Kadau, Mosquet, Lombois-Burger, Cayer-Barrioz, Herrmann, Spencer and Isa2013; Heussinger Reference Heussinger2013; Seto et al.
Reference Seto, Mari, Morris and Denn2013). These works support the ‘stress-induced friction’ scenario (Mari et al.
Reference Mari, Seto, Morris and Denn2014; Wyart & Cates Reference Wyart and Cates2014), and the contribution of contact forces was also confirmed in experiments (Lin et al.
Reference Lin, Guy, Hermes, Ness, Sun, Poon and Cohen2015; Clavaud et al.
Reference Clavaud, Bérut, Metzger and Forterre2017). Thus, the particle-scale mechanism of shear thickening is, to a great extent, understood.
However, particle-scale simulations are not capable of reproducing engineering-scale flows of dense suspensions due to the practical limits imposed on the system size by computational tractability. For this reason, it is important to develop effective continuum models through the design of suitable non-Newtonian constitutive relations. Besides laboratory experiments, particle-scale simulations are an important source of indications for the development of such models. A complete model should describe the fluid response under any flow condition (Miller, Singh & Morris Reference Miller, Singh and Morris2009), not only in the simple shear flows in which most experimental and computational data are retrieved. Indeed, the response of non-Newtonian fluids can depend on the type of flow, as exemplified by the observations of shear thinning and extensional thickening in some viscoelastic fluids. Particularly important is the class of extensional flows of suspensions, for which few rheological characterisations are available (Dai & Tanner Reference Dai and Tanner2017) and the sole computational investigation of which the authors are aware was performed by Sami (Reference Sami1996), who studied semidilute Brownian suspensions. (We note that in his analysis, flow-type dependence was not evidenced.) A related computational method to treat hydrodynamic interactions in diluted suspensions was introduced by Ahamadi & Harlen (Reference Ahamadi and Harlen2008). For important developments regarding emulsions of deformable droplets, we refer the reader to the work of Zinchenko & Davis (Reference Zinchenko and Davis2015).
To study the material response, we simulate motions of particles in the bulk region under prescribed flow conditions. As usual, periodic boundary conditions are employed to minimise finite-size effects. The Lees–Edwards boundary conditions (Lees & Edwards Reference Lees and Edwards1972) are commonly used to impose simple shear flows in many contexts, including suspension rheology (Bossis & Brady Reference Bossis and Brady1984; Mari et al. Reference Mari, Seto, Morris and Denn2014). In this work, we also apply the Kraynik–Reinelt boundary conditions (Kraynik & Reinelt Reference Kraynik and Reinelt1992; Todd & Daivis Reference Todd and Daivis1998), originally devised to impose planar extensional flows in non-equilibrium molecular dynamics simulations. With these, we can provide a first assessment of the flow-type dependence of the response in dense suspensions.
In §§ 2.1 and 2.2, we describe our simulation technique, which operates in the inertialess approximation. To compare the results under different flow conditions consistently, we employ the rheometric framework introduced by Giusteri & Seto (Reference Giusteri and Seto2017) (summarised in § 2.3), which defines, for the case of planar flows, a dissipative response function,
$\unicode[STIX]{x1D705}$
, and two non-dissipative response functions,
$\unicode[STIX]{x1D706}_{0}$
and
$\unicode[STIX]{x1D706}_{3}$
. These are defined for any flow type (simple shear, extensional and mixed flows) and offer a unified description of the material response. The results of our analysis, discussed in § 3, highlight the presence of flow-type dependence in the microstructure and in the non-Newtonian effects observed for dense suspensions.
2 Methods
2.1 Bulk rheology with periodic boundary conditions
Non-Newtonian incompressible fluids obey the differential equations
$$\begin{eqnarray}\unicode[STIX]{x1D70C}\left[\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}\right]=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748},\quad \text{with }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$
where
$\boldsymbol{u}$
is the velocity field and
$\unicode[STIX]{x1D70C}$
is the density. To close the system of equations, the stress tensor
$\unicode[STIX]{x1D748}$
must be given in terms of the velocity gradient through a constitutive prescription. The local value of the stress tensor describes the material response and is determined by the local history of deformation. To investigate this response, we consider small volume elements in which the velocity gradient
$\unicode[STIX]{x1D735}\boldsymbol{u}$
is approximately uniform. By simulating motions of particles in the volume element with fixed
$\unicode[STIX]{x1D735}\boldsymbol{u}$
, we can find the typical stress for a certain deformation history.
Time-dependent periodic boundary conditions allow us to impose
$\unicode[STIX]{x1D735}\boldsymbol{u}$
and effectively simulate the bulk behaviour. Since we consider planar flows in a 3D geometry, we can describe our methods considering the 2D projections of the computational cells. The cell frame vectors
$\boldsymbol{l}_{1}(t)$
and
$\boldsymbol{l}_{2}(t)$
(see figure 1) are prescribed to follow the velocity field
$\boldsymbol{u}=\unicode[STIX]{x1D735}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{r}$
, and periodic images of a particle at
$\boldsymbol{r}$
are given by
$\boldsymbol{r}^{\prime }=\boldsymbol{r}+i\boldsymbol{l}_{1}(t)+j\boldsymbol{l}_{2}(t)$
with (
$i,j=\pm 1,\pm 2,\ldots$
). For simple shear flows (
$\unicode[STIX]{x1D735}\boldsymbol{u}=\dot{\unicode[STIX]{x1D6FE}}\boldsymbol{e}_{y}\boldsymbol{e}_{x}$
), this is equivalent to the Lees–Edwards boundary conditions. The initial periodic cells are rectangles in the flow plane (blue in figure 1
a). A simple shear flow deforms the cells to parallelogram shapes (red in figure 1
a). To avoid significantly deformed periodic cells, the initial rectangular cells can be recovered as shown in figure 1(a). To impose planar extensional flows (
$\unicode[STIX]{x1D735}\boldsymbol{u}=\dot{\unicode[STIX]{x1D700}}\boldsymbol{e}_{x}\boldsymbol{e}_{x}-\dot{\unicode[STIX]{x1D700}}\boldsymbol{e}_{y}\boldsymbol{e}_{y}$
) for long times, we employ the Kraynik–Reinelt periodic boundary conditions (Kraynik & Reinelt Reference Kraynik and Reinelt1992; Todd & Daivis Reference Todd and Daivis1998). If the initial master cell is a regular square oriented at a certain angle
$\unicode[STIX]{x1D703}^{\ast }$
from the extension axis (
$x$
axis), the deformed parallelogram cell after a certain strain
$\unicode[STIX]{x1D700}_{p}$
can be remapped to the initial regular shape, as shown in figure 1(b).
(a) The deforming periodic cells for simple shear flows. The initial rectangular shape can be recovered when the shear strain
$\unicode[STIX]{x1D6FE}$
equals
$1$
by removing a part
$\text{A}$
of the master cell and including the corresponding part
$\text{A}^{\prime }$
of a periodic image in the new master cell. It should be noted that the recovery can be performed at any value of the shear strain
$\unicode[STIX]{x1D6FE}$
if the periodic displacement of the rows is taken into account. (b) The deforming periodic cells for extensional flow (Kraynik–Reinelt boundary conditions). The initial rectangular cell, which is oriented at a certain angle
$\unicode[STIX]{x1D703}^{\ast }\approx 31.7^{\circ }$
, can be recovered when the strain
$\unicode[STIX]{x1D700}$
equals
$\unicode[STIX]{x1D700}_{p}\approx 0.962$
by removing parts
$\text{A}$
–
$\text{C}$
of the master cell and including the corresponding parts
$\text{A}^{\prime }$
–
$\text{C}^{\prime }$
of periodic images in the new master cell.
2.2 Inertialess particle dynamics for suspensions
We numerically evaluate the stress tensor
$\unicode[STIX]{x1D748}$
by using particle simulations with deforming periodic cells. Our simulation is analogous to rate-controlled rheological measurements in the sense that time-averaged stress responses
$\langle \unicode[STIX]{x1D748}\rangle$
are evaluated for imposed velocity gradients
$\unicode[STIX]{x1D735}\boldsymbol{u}$
.
We consider non-Brownian, density matched and dense suspensions. Suspended particles interact with each other in several ways. As discussed in Mari et al. (Reference Mari, Seto, Morris and Denn2014), we take into account contact forces
$\boldsymbol{F}_{C}$
(and torques
$\boldsymbol{T}_{C}$
) and stabilising repulsive forces
$\boldsymbol{F}_{R}$
, besides hydrodynamic interactions
$\boldsymbol{F}_{H}$
and
$\boldsymbol{T}_{H}$
. Since the inertia of sufficiently small particles is negligible in comparison with the hydrodynamic drag forces, the particles obey the quasi-static equations of motion
$$\begin{eqnarray}\boldsymbol{F}_{H}+\boldsymbol{F}_{C}+\boldsymbol{F}_{R}=\mathbf{0}\quad \text{and}\quad \boldsymbol{T}_{H}+\boldsymbol{T}_{C}=\mathbf{0}.\end{eqnarray}$$
Here, forces
$\boldsymbol{F}$
and torques
$\boldsymbol{T}$
represent the set of forces and torques for
$N$
particles. Flows around microscale particles are dominated by viscous dissipation, and the inertia of the fluid is negligible, so they are described by the Stokes equations. The imposed velocity gradient
$\unicode[STIX]{x1D735}\boldsymbol{u}$
gives the background flow via the velocity
$\boldsymbol{u}(\boldsymbol{r})$
, vorticity
$\unicode[STIX]{x1D74E}\equiv \unicode[STIX]{x1D735}\times \boldsymbol{u}$
and rate of deformation tensor
$\unicode[STIX]{x1D63F}$
such that
$\unicode[STIX]{x1D735}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{r}=\boldsymbol{u}(\boldsymbol{r})=\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\boldsymbol{r}+(\unicode[STIX]{x1D74E}/2)\times \boldsymbol{r}.$
In this case, the hydrodynamic interactions can be expressed as the sum of linear resistances to the relative velocities
$\unicode[STIX]{x0394}\boldsymbol{U}^{(i)}\equiv \boldsymbol{U}^{(i)}-\boldsymbol{u}(\boldsymbol{r}^{(i)})$
, angular velocities
$\unicode[STIX]{x0394}\unicode[STIX]{x1D734}^{(i)}\equiv \unicode[STIX]{x1D734}^{(i)}-\unicode[STIX]{x1D74E}/2$
, for
$i=1,\ldots ,N$
, and imposed deformation
$\unicode[STIX]{x1D63F}$
via
$$\begin{eqnarray}\left(\begin{array}{@{}c@{}}\boldsymbol{F}_{H}\\ \boldsymbol{ T}_{H}\end{array}\right)=-\unicode[STIX]{x1D64D}\boldsymbol{\cdot }\left(\begin{array}{@{}c@{}}\unicode[STIX]{x0394}\boldsymbol{U}\\ \unicode[STIX]{x0394}\unicode[STIX]{x1D734}\\ \end{array}\right)+\unicode[STIX]{x1D64D}^{\prime } :\unicode[STIX]{x1D63F}_{N},\end{eqnarray}$$
where
$\unicode[STIX]{x0394}\boldsymbol{U}$
and
$\unicode[STIX]{x0394}\unicode[STIX]{x1D734}$
represent the set of relative velocities for
$N$
particles,
$\unicode[STIX]{x1D63F}_{N}$
is block-diagonal with
$N$
copies of
$\unicode[STIX]{x1D63F}$
, and
$\unicode[STIX]{x1D64D}$
and
$\unicode[STIX]{x1D64D}^{\prime }$
are resistance matrices, which can, in principle, be derived from the Stokes equations once the particle configurations are given. In dense suspensions, the long-range hydrodynamic interactions are screened by crowds of particles. Therefore, we may approximately construct the resistance matrices by including only the contributions of Stokes drag and lubrication forces.
In real suspensions, the lubrication singularity in
$\boldsymbol{F}_{H}$
is absent due to factors such as the surface roughness of particles – direct contacts are not forbidden. Hence, we include the contact interactions
$\boldsymbol{F}_{C}$
and
$\boldsymbol{T}_{C}$
in (2.2). The contact forces between solid particles depend on the nature of the particle surfaces. This is effectively encoded in the friction coefficient
$\unicode[STIX]{x1D707}$
that enters a simple friction model. By denoting with
$F_{C}^{n}$
and
$F_{C}^{t}$
normal and tangential forces respectively, we prevent sliding if
$F_{C}^{t}\leqslant \unicode[STIX]{x1D707}F_{C}^{n}$
. The normal force depends on the overlap between particles through an effective elastic constant and the tangential force depends on the sliding displacement in a similar way. The details of the model employed are given in Mari et al. (Reference Mari, Seto, Morris and Denn2014).
The presence of the stabilising repulsive force
$\boldsymbol{F}_{R}$
in (2.2) generates the rate dependence of rheological properties in such suspensions. Indeed, while reaching the same strain,
$\boldsymbol{F}_{R}$
can work more to prevent particle contacts under lower deformation rates, but less under higher rates. As a result, the number of contacts depends on the rate of the imposed flow. In colloidal suspensions, Brownian forces may play a similar role, as discussed by Mari et al. (Reference Mari, Seto, Morris and Denn2015a
).
The bulk stress tensor is obtained as
$$\begin{eqnarray}\unicode[STIX]{x1D748}=-p_{0}\unicode[STIX]{x1D644}+2\unicode[STIX]{x1D702}_{0}\unicode[STIX]{x1D63F}+V^{-1}\left(\mathop{\sum }_{i}\unicode[STIX]{x1D64E}_{D}^{(i)}+\mathop{\sum }_{i>j}\unicode[STIX]{x1D64E}_{P}^{(i,j)}\right),\end{eqnarray}$$
where
$\unicode[STIX]{x1D702}_{0}$
is the viscosity of the solvent,
$\unicode[STIX]{x1D64E}_{D}^{(i)}$
is the stresslet on particle
$i$
due to
$\unicode[STIX]{x1D63F}$
, and
$\unicode[STIX]{x1D64E}_{P}^{(i,j)}$
is the stresslet due to non-hydrodynamic interparticle forces between particles
$i$
and
$j$
(Mari et al.
Reference Mari, Seto, Morris and Denn2015b
). It should be noted that, since the hydrostatic pressure
$p_{0}$
is arbitrary, we set
$p_{0}=0$
. However, the last term in (2.4) is not traceless and thus contributes to the total pressure
$p\equiv -(1/3)\text{Tr}\,\unicode[STIX]{x1D748}$
.
2.3 General response functions for steady flows of non-Newtonian fluids
The stress
$\unicode[STIX]{x1D748}$
is a tensorial quantity, and we need a procedure to extract the relevant information from it in terms of scalar quantities. We are interested in comparing the material response under different types of imposed flow conditions. For this reason, we need a framework in which it is possible to identify the dependence on the flow type of each independent non-Newtonian effect. To this end, we use the framework introduced by Giusteri & Seto (Reference Giusteri and Seto2017), in which the characteristic rate of the imposed flow is defined independently of the flow type and a complete set of response functions is given. These functions generalise to any flow-type standard quantities such as viscosity and normal stress differences.
The velocity gradient
$\unicode[STIX]{x1D735}\boldsymbol{u}$
is decomposed into symmetric and antisymmetric parts as
$\unicode[STIX]{x1D735}\boldsymbol{u}=\unicode[STIX]{x1D63F}+\unicode[STIX]{x1D652}$
. In the planar case, with
$\unicode[STIX]{x1D63F}\neq \unicode[STIX]{x1D7EC}$
, we denote by
$\dot{\unicode[STIX]{x1D700}}$
the largest eigenvalue of
$\unicode[STIX]{x1D63F}$
and express
$\hat{\unicode[STIX]{x1D63F}}\equiv \unicode[STIX]{x1D63F}/\dot{\unicode[STIX]{x1D700}}$
and
$\hat{\unicode[STIX]{x1D652}}\equiv \unicode[STIX]{x1D652}/\dot{\unicode[STIX]{x1D700}}$
on the basis of the eigenvectors
$\hat{\boldsymbol{d}}_{1}$
and
$\hat{\boldsymbol{d}}_{2}$
of
$\unicode[STIX]{x1D63F}$
(corresponding to the eigenvalues
$\dot{\unicode[STIX]{x1D700}}$
and
$-\dot{\unicode[STIX]{x1D700}}$
) as follows:
$$\begin{eqnarray}\hat{\unicode[STIX]{x1D63F}}=\hat{\boldsymbol{d}}_{1}\hat{\boldsymbol{d}}_{1}-\hat{\boldsymbol{d}}_{2}\hat{\boldsymbol{d}}_{2},\quad \hat{\unicode[STIX]{x1D652}}=\unicode[STIX]{x1D6FD}_{3}(\hat{\boldsymbol{d}}_{2}\hat{\boldsymbol{d}}_{1}-\hat{\boldsymbol{d}}_{1}\hat{\boldsymbol{d}}_{2}).\end{eqnarray}$$
The non-vanishing and positive rate
$\dot{\unicode[STIX]{x1D700}}>0$
is used to set the time scale of deformation in any flow type. With this definition, the standard rate
$\dot{\unicode[STIX]{x1D6FE}}$
for simple shear corresponds to the value
$2\dot{\unicode[STIX]{x1D700}}$
. The vorticity
$\unicode[STIX]{x1D714}_{z}$
is represented by the dimensionless parameter
$\unicode[STIX]{x1D6FD}_{3}$
through
$\unicode[STIX]{x1D714}_{z}=2\dot{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D6FD}_{3}$
. It should be noted that planar extensional flows are characterised by
$\unicode[STIX]{x1D6FD}_{3}=0$
and simple shear flows by
$\unicode[STIX]{x1D6FD}_{3}=1$
.
A general representation of the stress tensor in planar flows is then given by
$$\begin{eqnarray}\unicode[STIX]{x1D748}(\dot{\unicode[STIX]{x1D700}},\unicode[STIX]{x1D6FD}_{3})=-p(\dot{\unicode[STIX]{x1D700}},\unicode[STIX]{x1D6FD}_{3})\unicode[STIX]{x1D644}+\dot{\unicode[STIX]{x1D700}}[\unicode[STIX]{x1D705}(\dot{\unicode[STIX]{x1D700}},\unicode[STIX]{x1D6FD}_{3})\hat{\unicode[STIX]{x1D63F}}+\unicode[STIX]{x1D706}_{0}(\dot{\unicode[STIX]{x1D700}},\unicode[STIX]{x1D6FD}_{3})\hat{\unicode[STIX]{x1D640}}+\unicode[STIX]{x1D706}_{3}(\dot{\unicode[STIX]{x1D700}},\unicode[STIX]{x1D6FD}_{3})\hat{\unicode[STIX]{x1D642}}_{3}],\end{eqnarray}$$
where
$\hat{\unicode[STIX]{x1D640}}\equiv -(1/2)(\hat{\boldsymbol{d}}_{1}\hat{\boldsymbol{d}}_{1}+\hat{\boldsymbol{d}}_{2}\hat{\boldsymbol{d}}_{2})+\hat{\boldsymbol{d}}_{3}\hat{\boldsymbol{d}}_{3}$
,
$\hat{\boldsymbol{d}}_{3}$
is the eigenvector of
$\unicode[STIX]{x1D63F}$
orthogonal to the flow plane and
$\hat{\unicode[STIX]{x1D642}}_{3}\equiv \hat{\boldsymbol{d}}_{1}\hat{\boldsymbol{d}}_{2}+\hat{\boldsymbol{d}}_{2}\hat{\boldsymbol{d}}_{1}$
is introduced to complete an orthogonal basis for the space of symmetric tensors for planar flows. The functional dependence of
$\unicode[STIX]{x1D705}$
,
$\unicode[STIX]{x1D706}_{0}$
and
$\unicode[STIX]{x1D706}_{3}$
on the two kinematical parameters
$\dot{\unicode[STIX]{x1D700}}$
and
$\unicode[STIX]{x1D6FD}_{3}$
needs to be determined to characterise the response in generic flows. We remark that the response functions
$\unicode[STIX]{x1D705}$
,
$\unicode[STIX]{x1D706}_{0}$
and
$\unicode[STIX]{x1D706}_{3}$
can depend on any other quantity that characterises the system. For instance, in § 3, we will also show their dependence on the volume fraction
$\unicode[STIX]{x1D719}$
. The function
$\unicode[STIX]{x1D705}$
is the only one to carry information about dissipation. We therefore refer to it as the dissipative response function, a generalised viscosity. The functions
$\unicode[STIX]{x1D706}_{0}$
and
$\unicode[STIX]{x1D706}_{3}$
carry information about non-dissipative responses, and we call them non-dissipative response functions. The presence of a non-vanishing
$\unicode[STIX]{x1D706}_{0}$
leaves the eigenvectors of the stress
$\unicode[STIX]{x1D748}$
aligned with those of
$\unicode[STIX]{x1D63F}$
, as happens in Newtonian fluids, but gives a contribution to the stress in the form of a modified pressure that is isotropic in the flow plane but different in the direction normal to the flow plane. On the other hand, a non-vanishing
$\unicode[STIX]{x1D706}_{3}$
corresponds to a rotation of the eigenvectors of
$\unicode[STIX]{x1D748}$
in the flow plane with respect to those of
$\unicode[STIX]{x1D63F}$
, determining the reorientation angle
$$\begin{eqnarray}\unicode[STIX]{x1D711}\equiv \arctan \left[\unicode[STIX]{x1D706}_{3}\left/\left(\unicode[STIX]{x1D705}+\sqrt{\unicode[STIX]{x1D705}^{2}+\unicode[STIX]{x1D706}_{3}^{2}}\right)\right.\right].\end{eqnarray}$$
For the sake of comparison, the shear viscosity
$\unicode[STIX]{x1D702}$
and normal stress differences
$N_{1}$
and
$N_{2}$
, defined for simple shear flows with
$\unicode[STIX]{x1D6FD}_{3}=1$
as functions of
$\dot{\unicode[STIX]{x1D6FE}}=2\dot{\unicode[STIX]{x1D700}}$
, are given by
$$\begin{eqnarray}\unicode[STIX]{x1D702}(2\dot{\unicode[STIX]{x1D700}})=\unicode[STIX]{x1D705}(\dot{\unicode[STIX]{x1D700}},1)/2,\quad N_{1}(2\dot{\unicode[STIX]{x1D700}})=-2\dot{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D706}_{3}(\dot{\unicode[STIX]{x1D700}},1)\quad \text{and}\quad N_{2}(2\dot{\unicode[STIX]{x1D700}})=\dot{\unicode[STIX]{x1D700}}[\unicode[STIX]{x1D706}_{3}(\dot{\unicode[STIX]{x1D700}},1)-(3/2)\unicode[STIX]{x1D706}_{0}(\dot{\unicode[STIX]{x1D700}},1)].\end{eqnarray}$$
Moreover, the extensional viscosity, defined for extensional flows with
$\unicode[STIX]{x1D6FD}_{3}=0$
, is given by
$\unicode[STIX]{x1D702}_{E}(\dot{\unicode[STIX]{x1D700}})=2\unicode[STIX]{x1D705}(\dot{\unicode[STIX]{x1D700}},0)$
. Therefore, the Trouton ratio
$\unicode[STIX]{x1D702}_{E}(\dot{\unicode[STIX]{x1D700}})/\unicode[STIX]{x1D702}(2\dot{\unicode[STIX]{x1D700}})$
equals
$4$
only if
$\unicode[STIX]{x1D705}(\dot{\unicode[STIX]{x1D700}},0)=\unicode[STIX]{x1D705}(\dot{\unicode[STIX]{x1D700}},1)$
.
We want to stress that, to arrive at (2.6), no a priori assumption is made on the list of quantities on which the material functions can depend. Hence, the stress tensor for any non-Newtonian fluid model under steady flow conditions can be expressed in the form (2.6). For example, since
$\unicode[STIX]{x1D63F}^{2}$
in planar flows is a linear combination of
$\unicode[STIX]{x1D644}$
and
$\hat{\unicode[STIX]{x1D640}}$
, the class of Reiner–Rivlin fluids corresponds to choosing
$\unicode[STIX]{x1D706}_{3}=0$
, and assuming
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D706}_{0}$
to be independent of
$\unicode[STIX]{x1D6FD}_{3}$
. Similarly, second-order fluids under steady shear flows would produce constant values of
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D706}_{0}$
, and entail
$\unicode[STIX]{x1D706}_{3}\propto \unicode[STIX]{x1D6FD}_{3}\dot{\unicode[STIX]{x1D700}}$
, since, under such flows,
$\unicode[STIX]{x1D652}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}-\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\unicode[STIX]{x1D652}=\unicode[STIX]{x1D6FD}_{3}\dot{\unicode[STIX]{x1D700}}^{2}\hat{\unicode[STIX]{x1D642}}_{3}$
. A detailed discussion of the representation (2.6) and its relation to fluid models is given in Giusteri & Seto (Reference Giusteri and Seto2017).
3 Results
To obtain the numerical results, we mainly performed 50 independent 3D simulations with
$2000$
particles. The periodic cells were initially cuboids with ratio
$5:5:1$
. We also performed some simulations with 4000 particles using double-sized cells (
$5:5:2$
) to confirm the absence of significant finite-size effects (data are not shown). Regarding the friction coefficient, we set
$\unicode[STIX]{x1D707}=1$
, since it is the value that, in a previous paper (Mari et al.
Reference Mari, Seto, Morris and Denn2015a
), was found to give good agreement with the experimental data by Cwalina & Wagner (Reference Cwalina and Wagner2014). It is worth mentioning that Tanner & Dai (Reference Tanner and Dai2016) showed that
$\unicode[STIX]{x1D707}=0.5$
gives better quantitative agreement with different experimental data. Nevertheless, such a fine tuning of
$\unicode[STIX]{x1D707}$
is not necessary for our qualitative analysis.
The short-range repulsive force is given by
$|\boldsymbol{F}_{R}|=F_{R}(0)\exp [-(r-2a)/(0.02a)]$
, with
$a$
the particle radius. A reference rate is set as
$\dot{\unicode[STIX]{x1D700}}_{0}\equiv F_{R}(0)/(12\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{0}a^{2})$
. To estimate the importance of the inertial effects, we can use the Stokes number, given by
$$\begin{eqnarray}St\equiv \frac{2\unicode[STIX]{x1D70C}_{p}a^{2}\dot{\unicode[STIX]{x1D700}}}{\unicode[STIX]{x1D702}_{0}}=\frac{\unicode[STIX]{x1D70C}_{p}F_{R}(0)}{6\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{0}^{2}}\frac{\dot{\unicode[STIX]{x1D700}}}{\dot{\unicode[STIX]{x1D700}}_{0}},\end{eqnarray}$$
with
$\unicode[STIX]{x1D70C}_{p}$
the particle density. Hence, inertial effects can be neglected if
$St\ll 1$
; that is, when the ratio
$\dot{\unicode[STIX]{x1D700}}/\dot{\unicode[STIX]{x1D700}}_{0}$
is much smaller than
$6\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{0}^{2}/\unicode[STIX]{x1D70C}_{p}F_{R}(0)$
. The preceding threshold determines, for each specific system, the region in which the rheology curves obtained with our simulations can be expected to be in agreement with real data.
3.1 Dissipative response function
$\unicode[STIX]{x1D705}$
For the case of monodisperse suspensions, the dissipative response function
$\unicode[STIX]{x1D705}$
significantly increases with the rate
$\dot{\unicode[STIX]{x1D700}}$
in both simple shear and extensional flows (figure 2). Not only shear thickening but also extensional thickening occurs. However, below thickening, there is a clear flow-type dependence. The value of
$\unicode[STIX]{x1D705}$
in extensional flow is much higher than the value obtained in simple shear flow (the Trouton ratio is much larger than
$4$
). On the other hand, above thickening, the values of
$\unicode[STIX]{x1D705}$
in extensional and shear flows are almost indistinguishable (the Trouton ratio is very close to
$4$
) and the flow-type dependence is hindered.
The significant discrepancy observed below thickening is due to shear-induced ordering, which can occur only in simple shear flows, as we confirm by analysing the pair distribution functions in § 3.4. Since the streamlines of a simple shear flow are straight and parallel to each other, particles tend to be arranged in chain-like structures along the flow direction. We observe a gradual decrease of
$\unicode[STIX]{x1D705}$
over time (strain thinning) in simple shear flows, which indicates the growth of the ordered structure. It should be noted that the shear-induced ordering is enhanced by the periodic boundary conditions, since linear chains may connect with their own periodic images. By contrast, the streamlines of extensional flows are never parallel to each other. Therefore, there is no obvious ordered structure compatible with extensional flows. Indeed, we neither observe strain thinning nor any ordered microstructure in the extensional flow simulation.
Both shear thickening and extensional thickening are observed in the rate dependence of the dissipative material function
$\unicode[STIX]{x1D705}$
. The data are for monodisperse suspensions in extensional flows (filled symbols with solid lines) and in simple shear flows (open symbols with dashed lines). Below thickening,
$\unicode[STIX]{x1D705}$
in simple shear flows is much lower than
$\unicode[STIX]{x1D705}$
in extensional flows. In each simulation, the time average is taken over five strains after reaching the steady state. The error bars show the standard deviation for 50 independent simulations.
In the thickened regime, frictional contact forces are constantly activated. These contact forces are so strong that particles are easily prevented from following the background flow; thus, ordered structures cannot be developed. As long as the disordered structure is maintained under simple shear flows, the value of
$\unicode[STIX]{x1D705}$
remains very close to that observed in extensional flows.
The shear-induced ordering can be hindered by mixing particles with different sizes. To see this effect, we consider two types of bidisperse suspensions with different size ratios:
$a_{2}/a_{1}=1.2$
and
$a_{2}/a_{1}=1.4$
(named ‘weak’ and ‘strong’ respectively). The two populations occupy the same volume fractions, i.e.
$\unicode[STIX]{x1D719}_{1}=\unicode[STIX]{x1D719}_{2}=\unicode[STIX]{x1D719}/2$
. In the weakly bidisperse suspensions (figure 3
a), although the differences clearly become smaller, some flow-type dependence can still be seen, especially for
$\unicode[STIX]{x1D719}=0.5$
. In the strongly bidisperse suspensions (figure 3
b), we no longer see a noticeable flow-type dependence. Trouton ratios are always close to
$4$
.
3.2 Pressure and anisotropic response
The total stress tensor
$\unicode[STIX]{x1D748}$
is usually split into two parts: the isotropic pressure term and the traceless extra-stress term. Although only the extra-stress term determines the flows of incompressible fluids, the pressure
$p$
is also part of the material response. As seen in figure 4(a) for monodisperse suspensions in extensional flows, the pressure term
$p$
varies in a similar way to
$\dot{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D705}$
: the ratio
$\dot{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D705}/p$
remains of the order of unity even when
$\unicode[STIX]{x1D705}$
significantly increases by thickening. In our simulation, the volume of the periodic cells is fixed; therefore, the system can never dilate. However, such increase of
$p$
with
$\dot{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D705}$
suggests that extensional thickening (and shear thickening) of suspensions is a phenomenon related to that of dilatancy in granular materials.
The pressure term
$p$
contributes isotropically to
$\unicode[STIX]{x1D748}$
by definition. However, there is another contribution to the stress
$\unicode[STIX]{x1D748}$
sharing the same origin. The non-dissipative response associated with dilatancy can be anisotropic and activate the response function
$\unicode[STIX]{x1D706}_{0}$
. The dimensionless ratio
$\dot{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D706}_{0}/p$
represents such anisotropy. Its positive values reported in figure 4(b) indicate that the in-plane pressure is higher than the out-of-plane pressure. However, this anisotropy is not very strong, as it would be if the pressure dilatancy were only present in the flow plane.
Mixing particles of different sizes hinders the shear-induced ordering. Differences in
$\unicode[STIX]{x1D705}$
between extensional flow (filled symbols with solid lines) and simple shear (open symbols with dashed lines) are still present in the weakly bidisperse suspensions,
$a_{2}/a_{1}=1.2$
(a), but no longer significant in the strongly bidisperse suspensions,
$a_{2}/a_{1}=1.4$
(b).
(a) The ratio between
$\dot{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D705}$
and
$p$
for monodisperse suspensions in extensional flows remains of the order of unity even when the viscosity increases significantly. (b) The positive values of the ratio
$\dot{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D706}_{0}/p$
for monodisperse suspensions in extensional flows indicate some anisotropy in the pressure response; that is, the in-plane pressure is higher than the out-of-plane pressure.
3.3 Reorientation angle of stress eigenvectors
Besides the ordering in simple shear flow, we can see some flow-type dependence in the reorientation angle
$\unicode[STIX]{x1D711}$
, defined in (2.7). In extensional flows, the principal axes of the stress tensor
$\unicode[STIX]{x1D748}$
must be parallel to the eigenvectors of
$\unicode[STIX]{x1D63F}$
due to symmetry considerations. Indeed, the reorientation angle
$\unicode[STIX]{x1D711}$
fluctuates around zero in those simulations. In simple shear flows, the shear-induced ordering is accompanied by large negative values of
$\unicode[STIX]{x1D711}$
(figure 5
a). On the other hand, in the disordered states above thickening (figure 5
b) and with strong bidispersity (figure 5
c),
$\unicode[STIX]{x1D711}$
is always rather small but non-zero. In our inertialess simulation, this finite flow-type dependence indicates some characteristic microstructure (see § 3.4) due to the presence of vorticity in simple shear, which is absent in extensional flows.
It is worth commenting on the dependence of the angle
$\unicode[STIX]{x1D711}$
on the volume fraction
$\unicode[STIX]{x1D719}$
(figure 5). In the thickened regime, corresponding to higher values of
$\dot{\unicode[STIX]{x1D700}}$
, the angle
$\unicode[STIX]{x1D711}$
is always positive for
$\unicode[STIX]{x1D719}=0.5$
. The values of
$\unicode[STIX]{x1D711}$
become smaller and can take slightly negative values as
$\unicode[STIX]{x1D719}$
increases. This behaviour is consistent with some experimental measurements of
$N_{1}$
, which is proportional to
$-\unicode[STIX]{x1D706}_{3}$
. When the volume fraction is not very high, negative values have been observed for
$N_{1}$
(Lee et al.
Reference Lee, Alcoutlabi, Magda, Dibble, Solomon, Shi and McKenna2006; Cwalina & Wagner Reference Cwalina and Wagner2014), corresponding to positive
$\unicode[STIX]{x1D711}$
, while the sign of
$N_{1}$
turns positive (negative
$\unicode[STIX]{x1D711}$
) at higher volume fractions (Lootens et al.
Reference Lootens, van Damme, Hémar and Hébraud2005; Dbouk et al.
Reference Dbouk, Lobry and Lemaire2013).
The reorientation angle
$\unicode[STIX]{x1D711}$
is non-zero in simple shear flows for both monodisperse suspensions (a,b) and strongly bidisperse suspensions (c), while it fluctuates around zero in extensional flows (not shown). This angle is associated with the first normal stress difference
$N_{1}$
, and we have
$\unicode[STIX]{x1D711}\approx -N_{1}/(4\dot{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D705})$
when
$\unicode[STIX]{x1D706}_{3}\ll \unicode[STIX]{x1D705}$
. The large standard deviations present in the monodisperse case (a) below thickening are due to the existence of several types of stable ordered structures displaying rather different values of
$\unicode[STIX]{x1D711}$
. This is likely to be an effect originated by the finite size of the simulation cell. When the microstructure is disordered (b,c), the standard deviations are smaller and comparable.
3.4 Microstructure
As discussed in the modelling section § 2.2, it is reasonable to neglect particle and fluid inertia in the particle-scale dynamics. The response of such inertialess material elements to an imposed flow essentially depends on the microstructure built by the particles during the flow (Morris Reference Morris2009). To measure the correlation of particle positions, we evaluate the pair distribution function
$g(\boldsymbol{r})\equiv P_{1|1}(\boldsymbol{r}|\mathbf{0})/n$
, where
$n$
is the average number density of particles and
$P_{1|1}(\boldsymbol{r}|\mathbf{0})$
is the conditional probability of finding a particle at
$\boldsymbol{r}$
with the condition that another particle is at the origin
$\mathbf{0}$
. Figure 6(a,b) shows
$g(\boldsymbol{r})$
in the flow-plane slice
$|z|<0.1a$
. We also consider angular distributions
$g_{c}(\unicode[STIX]{x1D703})$
for contacting (and nearly contacting) particles such that
$|\boldsymbol{r}|<2.02a$
. The angle
$\unicode[STIX]{x1D703}$
is measured from the
$\hat{\boldsymbol{d}}_{1}$
axis.
(a) The pair distribution function
$g(\boldsymbol{r})$
highlights the presence of strong ordering below thickening (
$\dot{\unicode[STIX]{x1D700}}/\dot{\unicode[STIX]{x1D700}}_{0}=0.01$
) for monodisperse suspensions (
$\unicode[STIX]{x1D719}=0.54$
) under simple shear, and its absence above thickening (
$\dot{\unicode[STIX]{x1D700}}/\dot{\unicode[STIX]{x1D700}}_{0}=0.02$
). (b) No obvious ordered structure can be associated with
$g(\boldsymbol{r})$
for the case of extensional flow both below and above thickening. (c,d) The polar plot of the angular distribution
$g_{c}(\unicode[STIX]{x1D703})$
of contacting (and nearly contacting) particles such that
$|\boldsymbol{r}|<2.02a$
shows that the strong enhancement of contact interactions is mainly responsible for the thickening. It should be noted that, since practically no particles are in contact for
$\dot{\unicode[STIX]{x1D700}}/\dot{\unicode[STIX]{x1D700}}_{0}\leqslant 0.005$
in simple shear, the corresponding data in (c) are negligible. The difference between simple shear flow (c) and extensional flow (d) is discussed in the main text.
As seen in figure 6(a), a stripe-patterned correlation
$g(\boldsymbol{r})$
appears for the monodisperse suspensions in simple shear flow below thickening (
$\dot{\unicode[STIX]{x1D700}}/\dot{\unicode[STIX]{x1D700}}_{0}=0.01$
). The periodic peaks and striped correlation indicate the formation of chain-like structures by the particles. Once such chain-like structure is formed, particle interactions are rather weak, which leads to significantly low values of
$\unicode[STIX]{x1D705}$
, as seen in figure 2. The microstructure is totally different above thickening (
$\dot{\unicode[STIX]{x1D700}}/\dot{\unicode[STIX]{x1D700}}_{0}=0.02$
). The long-range correlation is no longer seen. The correlation pattern indicates some disordered anisotropic microstructure. In figure 6(c), the angular contact distribution
$g_{c}(\unicode[STIX]{x1D703})$
clearly shows that the number of contacting particles increases remarkably for
$\dot{\unicode[STIX]{x1D700}}/\dot{\unicode[STIX]{x1D700}}_{0}\geqslant 0.02$
. This observation is consistent with the idea that shear thickening is caused by the development of the contact network (Seto et al.
Reference Seto, Mari, Morris and Denn2013).
These results can be directly compared with those for the extensional flow simulation. As seen in figure 6(b), even below thickening (
$\dot{\unicode[STIX]{x1D700}}/\dot{\unicode[STIX]{x1D700}}_{0}=0.01$
), there is no long-range correlation in
$g(\boldsymbol{r})$
. The distribution pattern has horizontal and vertical mirror symmetries, and no vorticity skews the correlation in the extensional flow. The distribution pattern does not change much above thickening (
$\dot{\unicode[STIX]{x1D700}}/\dot{\unicode[STIX]{x1D700}}_{0}=0.02$
). However, a clear difference is present in the angular contact distribution
$g_{c}(\unicode[STIX]{x1D703})$
(figure 6
d). A flame-shaped distribution transforms into a fan-shaped distribution at the extensional thickening transition. Thus, just below the transition, we can find contacting particles only around the directions of the compression axis. Nevertheless, the width of the flame shape indicates that, differently from what we observed under shear, the contact chains do not correspond to stable ordered chains of particles. Rather, they are constantly rebuilt among new neighbouring particles. Such contact chains which are roughly parallel and oriented along the compression axis do not contribute to the viscosity significantly (figure 2). By contrast, above the thickening transition, contacting particles can be found in all directions, even in the directions of the extension axis (
$\unicode[STIX]{x1D703}=0$
and
$\unicode[STIX]{x03C0}$
). This distribution suggests an anisotropic network structure for the pattern of contacts, which enhances the viscosity. Thus, we can describe the essence of extensional thickening as a contact-chain to contact-network transition.
The fact that such a transition occurs in extensional flows without significant change of the long-range correlation (always absent) indicates that also in simple shear flows the contact-chain to contact-network transition is mainly responsible for the thickening. Indeed, while we observe a concurrent order–disorder structural transition in monodisperse suspensions under shear, this is not present in strongly bidisperse suspensions, which nevertheless display a strong thickening behaviour.
4 Conclusions
We numerically explored the non-Newtonian character of dense suspensions, which has a different origin from that of viscoelastic fluids. This character is manifested in three main aspects: rate dependence, non-dissipative responses and flow-type dependence. By analysing the thickening in both extensional and simple shear flows, we were able to confirm that the contact-chain to contact-network transition is its main cause. Non-dissipative responses, such as normal stress differences, are present in any flow regime. Flow-type dependence is evident in monodisperse suspensions below thickening, where ordering occurs under simple shear. The prevention of ordering through thickening or polydispersity hinders (but does not cancel) the flow-type dependence.
Acknowledgements
The authors acknowledge support from the Okinawa Institute of Science and Technology Graduate University. The research of R.S. is partially supported by JSPS KAKENHI Grant Number JP17K05618.
