2. Micro-mechanics

A steady, planar extensional flow of a dense suspension of identical spheres with radius $a$ is characterized by an average rate of deformation tensor $\boldsymbol{\mathsf{D}}$ with non-zero components ${\mathsf{D}}_{11}=-{\mathsf{D}}_{22}=\dot {\gamma }$, where $x_{1}$ and $x_{2}$ are the axes in the directions of greatest extension and compression, respectively, and $\dot {\gamma }$ is the constant shear rate. We focus on a typical pair of spheres and their near-contacting neighbours, and take $\skew4\hat {\boldsymbol {d}}^{(BA)}$ to be the unit vector directed from the centre of sphere $A$ to that of sphere $B$, with $\skew4\hat {\boldsymbol {d}}^{(AB)} = -\skew4\hat {\boldsymbol {d}}^{(BA)}$ (see figure 1). Then, with $\theta ^{(BA)}$ the time-dependent angle between $\skew4\hat {\boldsymbol {d}}^{(BA)}$ and the $x_{2}$ axis,

(2.1)\begin{equation} \skew4\hat{d}_{\alpha}^{(BA)}=(\sin \theta^{(BA)},\cos \theta^{(BA)}) \end{equation}

and the components of the unit tangent vector, $\boldsymbol {\hat {t}}^{(BA)}=-\boldsymbol {\hat {t}}^{(AB)}$, perpendicular to it, are

(2.2)\begin{equation} t_{\alpha }^{(BA) }=( \cos \theta^{(BA)},-\sin \theta^{(BA)}), \end{equation}

or $t_{\alpha }^{(BA)}=\varepsilon _{\alpha \beta } d_{\beta }^{(BA)}$, where $\varepsilon _{12} = -\varepsilon _{21} = 1$ and $\varepsilon _{11} = \varepsilon _{22} =0$. The unit vectors $\skew4\hat {\boldsymbol {d}}^{(BA)}$ and $\boldsymbol {\hat {t}}^{(BA)}$ are indicated in figure 1.

Figure 1.

The pair $AB$ and their near-contacting neighbours, with the angle $\theta ^{(BA)}$ and the unit vectors $\skew4\hat {\boldsymbol {d}}^{(BA)}$ and $\boldsymbol {\hat {t}}^{(BA)}$ along and perpendicular to the line of centres, respectively.

2.1. Kinematics

In planar extensional flow, the relative motion of the centre of particle $B$ with respect to the centre of particle $A$ is

(2.3)\begin{equation} v_{\alpha }^{(BA)}=\frac{{\textrm{d}s}^{(BA)}}{\textrm{d}t}\skew4\hat{d}_{\alpha }^{(BA)} +2a\frac{\textrm{d}{\theta}^{(BA)}}{\textrm{d}t}\hat{t}_{\alpha }^{(BA)}, \end{equation}

where $s$ is the separation of the edges along the line of centres. The relative velocity of their points of near contact is then,

(2.4)\begin{equation} v_{\alpha }^{(BA)}+a(\omega ^{(A)}+\omega ^{(B)})\hat{t}_{\alpha }^{(BA)} \equiv v_{\alpha }^{(BA)}+aS\hat{t}_{\alpha }^{(BA)}, \end{equation}

where $\omega$ is the angular velocity of the sphere and $S$ is their sum.

The interaction of $A$ with a near-contacting neighbours $n$, other than $B$, is treated differently; the sphere $n$ is assumed to move relative to $A$ with the average flow. Then, neglecting fluctuations in translational velocity, the relative velocity of centres of pair $nA$ is

(2.5)\begin{equation} v_{\alpha }^{(nA)}=2a{\mathsf{D}}_{\alpha \beta }\skew4\hat{d}_{\beta }^{(nA)}, \end{equation}

and the relative velocity of the points of near contact $nA$ is

(2.6)\begin{equation} v_{\alpha }^{(nA)}+a\omega ^{(A)}\hat{t}_{\alpha }^{(nA)}. \end{equation}

2.2. Force

The force of interaction between a typical pair $AB$ of particles is related to the relative velocity and distance between their points of near contact. According to Jeffrey & Onishi (Reference Jeffrey and Onishi1984) and Jeffrey (Reference Jeffrey1992), the force $\boldsymbol {F}^{(BA)}$ exerted by sphere $B$ on sphere $A$ through a fluid with viscosity $\mu$, is

(2.7)\begin{align} F_{\alpha }^{(BA) } &= 6{\rm \pi} \mu a{\mathsf{K}}_{\alpha \beta }^{(BA)} v_{\beta }^{(BA) } -\frac{F_{0}}{s^{(BA)}}\skew4\hat{d}_{\alpha}^{( BA)} -9.54{\rm \pi} \mu a^{2} (\hat{t}_{\beta} {\mathsf{D}}_{\beta \xi } \skew4\hat{d}_{\xi }) \hat{t}_{\alpha }^{(BA) }\nonumber\\ &\quad +{\rm \pi} \mu a^{2}\left[ \ln \left( \frac{a}{s^{(BA) }}\right) -0.96\right] \omega^{(A)} \hat{t}_{\alpha}^{(BA)}+{\rm \pi} \mu a^{2}\ln \left( \frac{a}{s^{(BA)}} \right) \omega^{(B)} \hat{t}_{\alpha}^{(BA)} , \end{align}

where

(2.8)\begin{equation} {\mathsf{K}}_{\alpha \beta }^{(BA)}=\frac{1}{4}\frac{a}{s^{(BA)}} \skew4\hat{d}_{\alpha}^{(BA)} \skew4\hat{d}_{\beta}^{(BA)}+\left[ \frac{1}{6}\ln \left( \frac{a}{s^{(BA)}}\right) +0.64 \right] \hat{t}_{a}^{(BA)} \hat{t}_{\beta}^{(BA)} \end{equation}

and the constant terms have been retained because they are of a similar size, unless $s$ is extremely small. The interaction force also includes a short-range repulsion of strength $F_{0}$ force times length (e.g. Singh & Nott Reference Singh and Nott2000).

We take the near-contacting neighbours, $m\neq B$, to be those that most influence equilibrium and make the greatest contribution to the stress. There are $k-1$ of these per sphere and we assume that the separation between their edges is $\bar {s}$. The number, $k$, of near-contacting neighbours is expected to be less, perhaps far less, than the number of nearest neighbours and to depend upon the area fraction, or concentration, $c$.

For the near-contacting neighbours, $m$, the corresponding force is based on the average motion and the separation $\bar {s}$

(2.9)\begin{align} F_{\alpha }^{(mA) } &= \frac{3}{\bar{s}}a^{3}{\rm \pi} \mu \left({\mathsf{D}}_{\beta \xi} \skew4\hat{d}_{\xi}^{(mA)} \skew4\hat{d}_{\beta}^{(mA)}\right) \skew4\hat{d}_{\alpha }^{(mA)}+{\rm \pi} \mu a^{2}\left[\ln \left( \frac{a}{\bar{s}}\right) -0.96\right] \omega^{(A)} \hat{t}_{\alpha}^{(mA)} \nonumber\\ &\quad + 2 a^{2}{\rm \pi} \mu \left[ \ln \left( \frac{a}{\bar{s}}\right) -0.96\right] \left({\mathsf{D}}_{\beta \xi}\hat{t}_{\xi}^{(mA)} \skew4\hat{d}_{\beta}^{(mA)}\right) \hat{t}_{\alpha }^{(mA)}-\frac{F_{0}}{\bar{s}}\skew4\hat{d}_{\alpha}^{(mA)}. \end{align}

2.3. Force and torque balances

In the more complicated two-dimensional simple shear flow, Jenkins & La Ragione (Reference Jenkins and La Ragione2015) require the balance of forces for a typical pair of spheres under the action of their near-contacting neighbours. Here, for the less complicated planar extensional flow, we consider the balance of both force and torque. The focus on a flow in which there is no average rotation makes this easier to do; and the possibility of solving for both the translational and rotational degrees of freedom of a typical pair should increase the accuracy of the modelling.

The balance of forces for particle $A$ is

(2.10)\begin{equation} F_{\alpha}^{(BA)} + \sum_{m{\neq} B}^{N^{(A)}}F_{\alpha}^{(mA)}=0; \end{equation}

and that for particle $B$ is

(2.11)\begin{equation} F_{\alpha}^{(AB)} + \sum_{m{\neq} A}^{N^{(B)}}F_{\alpha}^{(mB)}=0 , \end{equation}

with $F_{\alpha }^{(BA) }=-F_{\alpha }^{(AB)}$. The difference in the force balances projected along $\skew4\hat {\boldsymbol {d}}^{(BA) }$ is

(2.12)\begin{equation} 3{\rm \pi} \mu a\frac{a}{s^{(BA) }} \frac{\textrm{d}s^{(BA) }}{\textrm{d}t} -2\frac{ F_{0}}{s^{(BA) }}+6{\rm \pi} \mu a^{2}\frac{a}{\bar{s}}\skew4\hat{d}_{\alpha }^{(BA) }{\mathsf{J}}_{\alpha \beta \gamma }D_{\beta \gamma }-2\frac{F_{0}}{ \bar{s}}{\mathsf{Y}}_{\alpha }\skew4\hat{d}_{\alpha }^{(BA) }=0; \end{equation}

while along $\hat {t}_{\alpha }^{(BA) }$ it is

(2.13)\begin{align} 0 &=4\left[ \ln \left( \frac{a}{s^{(BA)}}\right) +3.84\right] \frac{\textrm{d}\theta^{(BA)}}{\textrm{d}t}-19\hat{t}_{\beta }^{(BA) }{\mathsf{D}}_{\beta \xi }\skew4\hat{d}_{\xi }^{(BA) } \nonumber\\ &\quad +\left[ 2\ln \left( \frac{a}{s^{(BA)}}\right) -0.96\right] S+\left[ \ln \left( \frac{a}{\bar{s}}\right) -0.96\right] S\varepsilon _{\alpha \beta }{\mathsf{Y}}_{\beta }^{(BA) }t_{\alpha }^{(BA) } \nonumber\\ &\quad +6\frac{a}{\bar{s}}{\mathsf{D}}_{\beta \xi }{\mathsf{J}}_{\alpha \xi \beta }^{(BA)}t_{\alpha }^{(BA) }+2\left[ 2\ln \left( \frac{a}{\bar{s}}\right) -1.92 \right] {\mathsf{D}}_{\beta \xi }{\mathsf{J}}_{\alpha \xi \beta }^{(BA)}t_{\alpha }^{(BA) } , \end{align}

with

(2.14)\begin{equation} {\mathsf{J}}_{\alpha \xi \beta }^{(BA)}=\sum_{m{\neq} B}^{N^{(A) }}\skew4\hat{d} _{\alpha }^{(mA) }\skew4\hat{d}_{\beta }^{(mA) }\skew4\hat{d} _{\xi }^{(mA) } \end{equation}

and

(2.15)\begin{equation} {\mathsf{Y}}_{\alpha }^{(BA) }=\sum_{m{\neq} B}^{N^{(A) }}\skew4\hat{d} _{\alpha }^{(mA) }. \end{equation}

In writing (2.12) and (2.13), we assume that ${\mathsf{J}}_{\alpha \xi \beta }^{(BA)}=-{\mathsf{J}}_{\alpha \xi \beta }^{(AB)}$ and ${\mathsf{Y}}_{\alpha }^{(BA) }=-{\mathsf{Y}}_{\alpha }^{(AB) };$ that is, the arrangement of near-contacting neighbours of $B$ is the reflection of that of $A$. The terms proportional to $S$ incorporate the influence of the rotations on the force balance.

The balance of torques for particle $A$ is

(2.16)\begin{equation} \varepsilon _{\alpha \beta}d^{(BA)}_{\alpha}F_{\beta}^{(BA) }+\varepsilon _{\alpha \beta}\sum_{m{\neq} B}^{N^{(A) }}d^{(mA)}_{\alpha}F_{\beta }^{(mA)}=0, \end{equation}

and that for particle $B$ is

(2.17)\begin{equation} \varepsilon _{\alpha \beta }d^{(AB)}_{\alpha}F_{\beta }^{(AB) }+\varepsilon _{\alpha \beta }\sum_{m{\neq} A}^{N^{\left( B\right) }}d^{(mB)}_{\alpha}F_{\beta }^{(mB)}=0; \end{equation}

so their sum is

(2.18)\begin{align} 0 &=4\left[ \ln \left( \frac{a}{s^{(BA)}}\right) +3.84\right]\frac{\textrm{d}\theta^{(BA)}}{\textrm{d}t}-19\hat{t}_{\mu }^{(BA) }D_{\mu \xi }\skew4\hat{d}_{\xi}^{(BA)} \nonumber\\ &\quad +\left[ 2\ln \left( \frac{a}{s^{(BA)}}\right) -0.96\right] S+\left[ \ln \left( \frac{a}{\bar{s}}\right) -0.96\right] S\left( k-1\right) \nonumber\\ &\quad +2\left[ 2\ln \left( \frac{a}{\bar{s}}\right) -1.92\right] \varepsilon _{\xi \nu }{\mathsf{A}}_{\nu \mu }^{(BA) }{\mathsf{D}}_{\mu \xi } , \end{align}

with

(2.19)\begin{equation} {\mathsf{A}}_{\nu \mu }^{(BA) }=\sum_{m{\neq} B}^{N^{(A) }}\skew4\hat{d} _{\nu }^{(mA) }\skew4\hat{d}_{\mu }^{(mA)} \end{equation}

and, again, ${\mathsf{A}}_{\nu \mu }^{(BA) }={\mathsf{A}}_{\nu \mu }^{(AB)}$.

The tensors $\boldsymbol {A}$, $\boldsymbol {J}$ and $\boldsymbol {Y}$ provide information on the distribution of spheres in the plane about a typical pair $AB$, as shown in figure 1. We assume here that the distributions about a pair at a given orientation is the average over all pairs at that orientation. These average distributions should depend on both $\skew4\hat {\boldsymbol {d}}^{(AB)}$ and $\boldsymbol {D}$. As do Jenkins & La Ragione (Reference Jenkins and La Ragione2015), we treat the local equilibrium with the approximation that $\boldsymbol {A}$, $\boldsymbol {J}$ and $\boldsymbol {Y}$ are independent of $\boldsymbol {D}$. Then,

(2.20)\begin{gather} {\mathsf{A}}_{\nu \mu }^{(BA) }=b_{1}\delta _{\nu \mu }+b_{2}\skew4\hat{d}_{\mu }^{(BA) }\skew4\hat{d}_{\nu }^{(BA) }, \end{gather}

(2.21)\begin{gather}{\mathsf{J}}_{\alpha \xi \beta }^{(BA)}=b_{3}\skew4\hat{d}_{\alpha }^{(BA)} \skew4\hat{d}_{\xi}^{(BA) }\skew4\hat{d}_{\beta }^{(BA) }+b_{4}\left(\skew4\hat{d}_{\alpha}^{(BA) } \delta_{\xi \beta}+\skew4\hat{d}_{\xi}^{(BA)} \delta_{\alpha \beta} + \skew4\hat{d}_{\beta}^{(BA)} \delta _{\xi \alpha}\right), \end{gather}

and

(2.22)\begin{equation} {\mathsf{Y}}_{\alpha }^{(BA) }=b_{5}\skew4\hat{d}_{\alpha }^{(BA)}. \end{equation}

To calculate the coefficients, Jenkins et al. (Reference Jenkins, La Ragione, Johnson and Makse2005) assume that given sphere $B$, the remaining near-contacting neighbours of $A$ are distributed uniformly around itscircumference. The results are given as a function of coordination number $k$ through

(2.23)\begin{equation} b=-\frac{3\sqrt{3}\left(k-1\right) }{16{\rm \pi} }, \end{equation}

by

(2.24a,b)\begin{equation} b_{1}=\frac{k-1}{2}-b,\ b_{2}=2b,\ b_{3}=0,\quad \text{and}\quad b_{4}=b,\ b_{5}=4b. \end{equation}

In the planar extensional flow of interest,

(2.25a,b)\begin{equation} \hat{t}_{\mu }^{(BA)} {\mathsf{D}}_{\mu \xi}\skew4\hat{d}_{\xi}^{(BA)} = \dot{\gamma}\sin 2\theta \quad\text{and}\quad \skew4\hat{d}_{\mu }^{(BA)} {\mathsf{D}}_{\mu \xi}\skew4\hat{d}_{\xi}^{(BA)}=-\dot{\gamma}\cos 2\theta . \end{equation}

We use these in the differences of the components of the force balances, make lengths dimensionless by the sphere radius $a$, time by the inverse of the shear rate, forces by ${\rm \pi} a^{2}\mu \dot {\gamma }$, write the dimensionless strength of the repulsion as $\hat {F}=F_{0}/({\rm \pi} a^{3}\mu \dot {\gamma })$ and remove the superscript $(BA)$. Then, the normal component becomes

(2.26)\begin{equation} \frac{1}{s}\frac{\textrm{d}s}{\textrm{d}\gamma }=\frac{2}{3}\hat{F}\left( \frac{1}{s}+\frac{4b}{\bar{s}}\right) +\frac{4b}{\bar{s}}\cos 2\theta; \end{equation}

and the tangential component is

(2.27)\begin{equation} \left[ \ln \left( \frac{1}{s}\right) +3.84\right] \frac{\textrm{d}\theta }{\textrm{d}\gamma }= \left[ c_{1}+c_{2}\ln \left( \frac{1}{s}\right) \right] \sin 2\theta , \end{equation}

with

(2.28a,b)\begin{equation} c_{1}=4.77-3\frac{b}{\bar{s}}\quad\text{and}\quad c_{2}=\frac{6b}{(4b-k+1) } \frac{1/\bar{s}}{\ln \left( 1/\bar{s}\right) -0.96}, \end{equation}

and we have employed the difference in the force balances and the sum of the torque balances to write

(2.29)\begin{equation} S=-2c_{2}\sin 2\theta . \end{equation}

The balances of force and torque, (2.26), (2.27) and (2.29), employed in (2.7), provide an expression for $\boldsymbol {F}^{(BA)}$ in terms of average quantities

(2.30)\begin{align} F_{\alpha }^{(BA) } &=4b\frac{F_{0}}{\bar{s}}\skew4\hat{d}_{\alpha}^{(BA)} +{\rm \pi} \mu a^{3}\frac{6b}{\bar{s}}\cos 2\theta \dot{\gamma} \skew4\hat{d}_{\alpha}^{(BA) } +{\rm \pi} \mu a^{2}\left(2c_{1}+0.96c_{2}\right) \sin 2\theta \dot{\gamma}\hat{t}_{a}^{(BA)} \nonumber\\ &\quad -9.54{\rm \pi} \mu a^{2}\sin 2\theta \dot{\gamma}\hat{t}_{\alpha}^{(BA)}. \end{align}

This is later used in the calculation of the stress.

2.4. Representative trajectory

The representative trajectory is a single trajectory that incorporates the influence of those that contribute most to the stress. Along the representative trajectory particle $B$ moves with respect to particle $A$ in a succession of states in which force and torque are balanced. The other near-contacting particles, $m$, of the pair are assumed to move with the average flow, at the constant distance $\bar {s}$ from the pair. The equation that determines this trajectory results from the balances of force and torque and is a function of two parameters: the average number of near-contacting particles, $k$, and the distance, $\bar {s}$. Upon combining (2.26) and (2.27), it is

(2.31)\begin{equation} \frac{\textrm{d}s}{\textrm{d}\theta }=\frac{2}{3\bar{s}}\frac{\hat{F}\left( \bar{s}+4bs\right) +6bs\cos 2\theta }{\left[ c_{1}+c_{2}\ln \left( 1/s\right) \right] \sin 2\theta }\left[ \ln \left( \frac{1}{s}\right) +3.84\right] . \end{equation}

Within the $\theta$ interval 0 to ${\rm \pi} /2$, the trajectory begins at $\theta _{0}$ and ends at $\theta _{1}$, and both angles must be determined. Because of the presence of $\hat {F}$, the trajectory is asymmetric about ${\rm \pi} /4$, and $\theta _{0}$ differs from ${\rm \pi} /2-\theta _{1}$.

The amount of total strain, $\hat {\gamma }$, necessary to complete the trajectory may be calculated from the pair interaction in the average flow. From (2.27)

(2.32)\begin{equation} \frac{\textrm{d}\gamma }{\textrm{d}\theta }=\frac{\ln \left( 1/s\right) +3.84}{[c_{1}+c_{2}\ln ( 1/s)] \sin 2\theta }; \end{equation}

so,

(2.33)\begin{equation} \hat{\gamma}=\int_{\theta _{0}}^{\theta _{1}}\frac{\ln \left( 1/s\right) +3.84}{[c_{1}+c_{2}\ln ( 1/s)] \sin 2\theta }\textrm{d}\theta. \end{equation}

2.5. Particle distribution

We next introduce the distribution of near-contacting neighbours along the trajectory, $A(\theta )$, defined so that $A( \theta ) \,\textrm {d}\theta$ is the average number of such particles within the element $\textrm {d}\theta$. At steady state, the flux, $A( \theta ) \,\textrm {d}\theta /\textrm {d}\gamma$, of these equilibrated particles along the trajectory is constant. That is, particles are more likely to be where the velocity along the trajectory is least. Because the repulsive force breaks the symmetry of approach and departure, the distribution is anticipated to be asymmetric about ${\rm \pi} /4$. In computations, we implement the flux condition as a differential equation

(2.34)\begin{equation} \frac{\textrm{d}A}{\textrm{d}\theta }=-\frac{A}{\dot{\theta}}\frac{\textrm{d}\dot{\theta}}{\textrm{d}\theta }, \end{equation}

with

(2.35)\begin{equation} \frac{\textrm{d}\dot{\theta}}{\textrm{d}\theta }=\frac{\partial \dot{\theta}}{\partial \theta } +\frac{\partial \dot{\theta}}{\partial s}\frac{\textrm{d}s}{\textrm{d}\theta }. \end{equation}

The distribution $A( \theta )$ is related to the average number near-contacting neighbours per particle by

(2.36)\begin{equation} 4\int_{\theta_{0}}^{\theta_{1}} A\left( \theta \right) \textrm{d}\theta = k. \end{equation}

We implement this as a differential equation for the partial number of near-contacting neighbours

(2.37)\begin{equation} I\left( \theta \right) \equiv \int_{\theta _{0}}^{\theta }A(\theta^{\prime }) \textrm{d}\theta^{\prime}, \end{equation}

as

(2.38)\begin{equation} \frac{\textrm{d}I}{\textrm{d}\theta }= A \left( \theta \right) , \end{equation}

with boundary conditions $I( \theta _{0}) =0$ and $I(\theta _{1})=k/4$.

Given that the beginning and ending angles of the trajectory differ, we take the beginning and ending values of the particle separation to be the same. There are three first-order differential equations, (2.31), (2.34) and (2.38), for $s$, $A$ and $I$ as functions of $\theta$, and four boundary conditions: one for each of $s_{0}$ and $s_{1}$ that introduce a single parameter, and two for $I$. Consequently, $\theta _{1}$ may be determined as part of the solution. The inputs are $\theta _0$, $s_0 = s_1$, $\bar {s}$ and $k$. In appendix B, we provide the Matlab code that is employed in the solver. We generate solutions and compare them with the results of Stokesian dynamics simulations in a later section.

3. Particle stress

Knowledge of the distribution of near-contacting neighbours $A( \theta )$ and the contact forces along the trajectory permits the calculation of the macroscopic particle stress in the suspension. The stress tensor is, according to Cauchy (Love (Reference Love1944), appendix, Note B),

(3.1)\begin{equation} {\mathsf{T}}_{\alpha \beta }=na\int_{0}^{2{\rm \pi} }A (\theta) F_{\alpha }\skew4\hat{d}_{\beta}\,\textrm{d}\theta, \end{equation}

where $n$ is the number of particles per unit area and $F_{\alpha }$ is given by (2.30). Then, the two-dimensional viscosity is $2a\mu$. The dimensionless form, $t_{\alpha \beta }=T_{\alpha \beta }/(2a\mu \dot {\gamma })$, with $n=c/({\rm \pi} a^{2})$, where c is the concentration, is

(3.2)\begin{align} {\mathsf{t}}_{a\beta } &={c} \frac{b}{\bar{s}}\int_{0}^{2{\rm \pi} }A\left( \theta \right) ( 2\hat{F}+3\cos 2\theta ) \skew4\hat{d}_{\alpha }\skew4\hat{d}_{\beta} \,\textrm{d}\theta \nonumber\\ &\quad +{c} \left( c_{1}+0.48c_{2}-4.77\right) \int_{0}^{2{\rm \pi} }A\left( \theta \right) \sin 2\theta \hat{t}_{a}\skew4\hat{d}_{\beta }\,\textrm{d}\theta. \end{align}

The particle shear stress,

(3.3)\begin{equation} \tau\equiv \tfrac{1}{2}\left({\mathsf{t}}_{11}-{\mathsf{t}}_{22}\right), \end{equation}

is

(3.4)\begin{align} \tau&=-2c\frac{b}{\bar{s}}\int_{0}^{{\rm \pi}/2 }A( \theta )( 2\hat{F}+3\cos 2\theta ) \cos 2\theta \,\textrm{d}\theta \nonumber\\ &\quad +c\left( 2c_{1}+0.96c_{2}-9.54\right) \int_{0}^{{\rm \pi}/2 }A\left( \theta \right) \sin^{2} 2\theta \,\textrm{d}\theta, \end{align}

where $b$, $c_1$ and $c_2$ are given in terms of $k$ in (2.24a,b) and (2.28a,b), respectively. The shear stress depends on the separation, $\bar {s}$, of near-contacting neighbours other than $B$, and on the area fraction, explicitly and through the coordination number, $k$. Because the direct contribution of the repulsive force to the integral is very small and the trigonometric factors associated with the other contributions are even about ${\rm \pi} /4$, the shear stress is independent of the asymmetry of the particle distribution about ${\rm \pi} /4$. In contrast, this asymmetry is crucial to the determination of the particle pressure.

The particle pressure,

(3.5)\begin{equation} p\equiv -\tfrac{1}{2}\left({\mathsf{t}}_{11}+{\mathsf{t}}_{22}\right) , \end{equation}

is

(3.6)\begin{equation} p=-2c\frac{b}{\bar{s}}\int_{0}^{{\rm \pi}/2 }A\left( \theta \right) ( 2\hat{F}+3\cos 2\theta ) \textrm{d}\theta . \end{equation}

This pressure also depends on $\bar {s}$ and $c$ and its existence is due to the asymmetry of $A$ about ${\rm \pi} /4$. This asymmetry is due to that of the separation along the representative trajectory created by $\hat {F}$ and the influence of the asymmetry of the separation on the angular velocity, $\dot {\theta }.$ The particle pressure and the mechanisms responsible for it are a focus of this paper; a particle shear stress may be calculated based on the average flow, although that determined here is several times less than this, because of the approximate satisfaction of equilibrium.

Particle stresses associated with motion along the representative trajectories are compared with those measured in Stokesian dynamics simulations after a discussion of the simulations.

4. Stokesian dynamics

We determine the trajectories of spherical particles in the flows by performing simulation with the same conditions as the theory (a monolayer with no inertia). We impose a planar extensional flow with shear rate $\dot {\gamma }$ (note that this is equivalent to the extensional rate $\dot {\varepsilon }$ in Seto, Giusteri & Martiniello (Reference Seto, Giusteri and Martiniello2017)),

(4.1a,b)\begin{equation} \boldsymbol{u}^{\infty} (\boldsymbol{r}) = \boldsymbol{\mathsf{D}} \boldsymbol{\cdot} \boldsymbol{r},\quad \boldsymbol{\mathsf{D}}= \begin{pmatrix} \dot{\gamma} & 0 \\ 0 & - \dot{\gamma} \end{pmatrix}. \end{equation}

A simulation box with periodic boundary conditions constantly deforms according to this velocity gradient $\boldsymbol{\mathsf{D}}$. Significant deformations of the simulation box can be avoided by using the Kraynik–Reinelt periodic boundary conditions, which rearrange the deformed box to the original square box after a constant strain interval (Kraynik & Reinelt Reference Kraynik and Reinelt1992; Todd & Daivis Reference Todd and Daivis1998; Seto et al. Reference Seto, Giusteri and Martiniello2017). Thus, the flow can be applied for a sufficiently long time to evaluate its steady states.

Due to the negligible inertia of the particles, translational and angular velocities can be determined by solving the force and torque balance equations for the respective particles ($i = 1, \dotsc , N$)

(4.2)\begin{equation} \begin{pmatrix} \boldsymbol{0} \\ \boldsymbol{0} \end{pmatrix} = \begin{pmatrix} \boldsymbol{F}_{{H}} \\ \boldsymbol{T}_{{H}} \end{pmatrix} + \begin{pmatrix} \boldsymbol{F}_{{R}} \\ \boldsymbol{0} \end{pmatrix}.\end{equation}

Here, a vector, such as $\boldsymbol {F}_{{H}}$, represents all $N$ particles, $\boldsymbol {F}_{{H}} \equiv (\boldsymbol {F}_{{H}}^{(1)}, \dotsc , \boldsymbol {F}_{{H}}^{(N)} )$.

The hydrodynamic interactions in the Stokes, zero Reynolds number, regime are linear in the velocities

(4.3)\begin{equation} \begin{pmatrix} \boldsymbol{F}_{{H}} \\ \boldsymbol{T}_{{H}} \end{pmatrix} = - \boldsymbol{\mathsf{R}}_{{FU}} \boldsymbol{\cdot} \begin{pmatrix} \boldsymbol{U}- \boldsymbol{u}^{\infty} \\ \boldsymbol{\varOmega} \end{pmatrix} + \boldsymbol{\mathsf{R}}_{{FE}}:\boldsymbol{\mathsf{D}}_{{N}},\end{equation}

where $\boldsymbol{\mathsf{D}}_{{N}}$ is block diagonal of $N$ copies of $\boldsymbol{\mathsf{D}}$. There exist several levels of approximations to construct the resistance matrices $\boldsymbol{\mathsf{R}}_{{FU}}$ and $\boldsymbol{\mathsf{R}}_{{FD}}$. Brady & Bossis (Reference Brady and Bossis1988) constructed them using truncated multipole expansions for the far-field interactions and a pairwise solution for lubrication interactions. In this work, we focus on a special situation in which repulsive forces are very weak in comparison with viscous drag forces. Under such conditions, particles tend to approach their neighbours very closely.

Because the resistance coefficients diverge at contacts ($s = 0$), the nearly touching hydrodynamic interactions dominate the dynamics. This is why we construct the approximate resistance matrices with the leading $1/s$ term in the normal component and the logarithmic term $\log (1/s)$ and following constants in the tangential component, using the solution for two nearly touching rigid spheres (Jeffrey & Onishi Reference Jeffrey and Onishi1984; Jeffrey Reference Jeffrey1992). (A detailed description can be seen elsewhere, cf. Mari et al. (Reference Mari, Seto, Morris and Denn2014).) The hydrodynamic interaction is effective only when $s < 0.10$; thus, the resistance coefficients remain positive in this range.

The repulsive force employed in this work is the same as that used in Nott & Brady (Reference Nott and Brady1994)

(4.4)\begin{equation} \boldsymbol{F}_{{R}} = F_0 \frac{\lambda^{-1} \,\textrm{e}^{- s/\lambda}}{1-\textrm{e}^{- s/\lambda}} \boldsymbol{n}, \end{equation}

where the range of repulsive force is set by a parameter $\lambda$. Because the repulsive force diverges as $F_0/s$ in the limit of contact, $s \to 0$, some force balance can occur at a finite gap. Because the repulsive force diverges as $F_0/s$ in the limit of contact, $s \to 0$, some force balance can occur at a finite gap. Thus, the gap $s$ remains positive, and contact forces do not appear in the current system. Note that the divergence in the lubrication coefficient does not guarantee the presence of a minimum $s >0$, thus it leads to a pathologic singularity in theoretical models (Ball & Melrose Reference Ball and Melrose1995).

By solving the force and torque balance equations (4.2) with the hydrodynamic interaction (4.3) and repulsive force (4.4), the linear and angular velocities $(\boldsymbol {U}, \boldsymbol {\varOmega })$ can be determined at each time step. Integrating these velocities $\boldsymbol {U}$ with a discretized time step, we obtain trajectories of particles. The particle stress tensor $\hat {\boldsymbol{\mathsf{T}}}$ is given by the symmetrized first moment

(4.5)\begin{equation} \hat{\boldsymbol{\mathsf{T}}} =\frac{1}{V} \sum_{j} \frac{ \boldsymbol{r}^{ij} \boldsymbol{F}^{ij} + \boldsymbol{r}^{ji} \boldsymbol{F}^{ji}}{2}, \end{equation}

with the pairwise forces $\boldsymbol {F}^{ij} \equiv \boldsymbol {F}^{ij}_{{Lub}}+\boldsymbol {F}^{ij}_{{R}}$ and relative positions $\boldsymbol {r}^{ij} \equiv \boldsymbol {r}^{i} -\boldsymbol {r}^{j}$ of all interacting particle pairs. Here, $V = 2a L^2$ is the volume of the monolayer system. Normalizing the symmetrized first moment with the shear stress of the suspending fluid $2 \mu \dot {\gamma }$, gives the dimensionless stress ${\mathsf{t}}_{\alpha \beta } \equiv \hat {{\mathsf{T}}}_{\alpha \beta }/2 \mu \dot {\gamma }$. Thus, we have the dimensionless particle pressure $p \equiv -({\mathsf{t}}_{11}+{\mathsf{t}}_{22})/2$ and the dimensionless particle shear stress $\tau \equiv ({\mathsf{t}}_{11}-{\mathsf{t}}_{22})/2$, respectively.

6. Tensorial formulation

As elaborated upon by Onat & Leckie (Reference Onat and Leckie1988) and Advani & Tucker (Reference Advani and Tucker1987, Reference Advani and Tucker1990), the distribution of near-contacting neighbours can be represented by an infinite series with respect to basis functions, such as

(6.1)\begin{equation} {\mathsf{f}}_{\alpha \beta }=\skew4\hat{d}_{\alpha }\skew4\hat{d}_{\beta }-\tfrac{1}{2}\delta _{\alpha \beta } \end{equation}

and

(6.2)\begin{align} {\mathsf{g}}_{\xi \eta \rho \beta } &=\skew4\hat{d}_{\xi }\skew4\hat{d}_{\eta }\skew4\hat{d}_{\rho }\widehat{ d}_{\beta }-\tfrac{1}{6}( \delta _{\xi \eta }\skew4\hat{d}_{\rho }\skew4\hat{d} _{\beta }+\skew4\hat{d}_{\xi }\skew4\hat{d}_{\eta }\delta _{\rho \beta }+\skew4\hat{d}_{\xi } \skew4\hat{d}_{\rho }\delta _{\eta \beta }+\skew4\hat{d}_{\eta }\skew4\hat{d}_{\beta }\delta _{\xi \rho }+\delta _{\xi \beta }\skew4\hat{d}_{\eta }\skew4\hat{d}_{\rho }\nonumber\\ &\quad +\skew4\hat{d}_{\xi }\skew4\hat{d}_{\beta }\delta _{\eta \rho }) +\tfrac{1}{24}( \delta _{\xi \eta }\delta _{\rho \beta }+\delta _{\xi \rho }\delta _{\eta \beta }+\delta _{\xi \beta }\delta _{\eta \rho }): \end{align}

(6.3)\begin{equation} A(\theta )=\frac{k}{2{\rm \pi} }(1+4\mathcal{Z}_{\alpha \beta }f_{\alpha \beta }+16\mathcal{B}_{\xi \eta \rho \beta }g_{\xi \eta \rho \beta }+\cdots). \end{equation}

The coefficients $\boldsymbol {\mathcal {Z}}$ and $\boldsymbol {\mathcal {B}}$ are completely traceless and completely symmetric tensors, related to the distribution through

(6.4)\begin{equation} \mathcal{Z}_{\alpha \beta }=\int_{0}^{2{\rm \pi} }A(\theta){\mathsf{f}}_{\alpha \beta }\,\textrm{d}\theta \end{equation}

and

(6.5)\begin{equation} \mathcal{B}_{\xi \eta \rho \beta }=\int_{0}^{2{\rm \pi} }A(\theta ){\mathsf{g}}_{\xi \eta \rho \beta }\,\textrm{d}\theta . \end{equation}

These are the second and fourth moments of the distribution, respectively.

The stress of (3.1) may be written in terms of these as

(6.6)\begin{align} {\mathsf{T}}_{\alpha \beta } &=4nak\frac{F_{0}}{\bar{s}}b\delta _{a\beta }-n{\rm \pi} \mu a^{3}\left( M+N\right) \mathcal{B} _{\alpha \beta \gamma \eta } {\mathsf{D}}_{\gamma \eta } -n{\rm \pi} \mu a^{3}\left( M-N\right) \frac{k}{4}\delta_{\alpha \eta} \delta_{\beta \gamma}{\mathsf{D}}_{\gamma \eta } \nonumber\\ &\quad -n{\rm \pi} \mu a^{3}\left(M+N\right) \frac{1}{6}\delta _{\alpha \beta }\mathcal{Z}_{\eta \gamma }{\mathsf{D}}_{\gamma \eta } -n{\rm \pi} \mu a^{3}\frac{2M-N}{6} (\delta_{\alpha \eta}\mathcal{Z}_{\gamma \beta } +\delta_{\beta \eta}\mathcal{Z}_{\gamma \alpha}) {\mathsf{D}}_{\eta \gamma }, \end{align}

where

(6.7a,b)\begin{equation} M=\frac{6b}{\bar{s}/a} \quad \text{and} \quad N=2c_{1}+0.96c_{2}-9.54; \end{equation}

or, more compactly, as

(6.8)\begin{equation} {\mathsf{T}}_{\alpha \beta }=na\left( 4kb\frac{F_{0}}{\bar{s}}\delta _{a\beta }+{\mathsf{G}}_{\alpha \beta \gamma \eta}{\mathsf{D}}_{\gamma \eta }\right), \end{equation}

where

(6.9)\begin{align} {\mathsf{G}}_{\alpha \beta \gamma \eta} &=-n{\rm \pi} \mu a^{3} \left[\frac{M-N}{4} k\delta_{\alpha \eta} \delta_{\beta \gamma} +\frac{M+N}{6} \delta _{\alpha \beta }\mathcal{Z}_{\eta \gamma } +\frac{2M-N}{6} ( \delta_{\alpha \eta}\mathcal{Z}_{\gamma \beta } +\delta_{\beta \eta}\mathcal{Z}_{\gamma \alpha})\right]\nonumber\\ &\quad -n{\rm \pi} \mu a^{3}(M+N) \mathcal{B}_{\alpha \beta \gamma \eta } . \end{align}

The stress depends only on the second and fourth moments of the distribution, although an approximation of the distribution in terms of these does not provide a good representation of it. Because we predict the distribution function, we are able to capture the essential role played by the fourth moment. This places our theory in the context of the work of Chacko et al. (Reference Chacko, Mari, Fielding and Cates2018), in which numerical simulations confirm the need of the fourth moment, in addition to the second moment, to describe stress reversal.

With knowledge of the distribution function, it is possible to evaluate the components of the second and fourth moments. In particular, when $F=0$, the only non-zero components are $\mathcal {B}_{1111}=\mathcal {B}_{2222}=-\mathcal {B}_{1122}$; when $F \ne 0$, then $\mathcal {Z}_{11}=-\mathcal {Z}_{22}$ are also different from zero. For example, when $c=0.64$, $k = 2.50$ and $\hat {F} = 10^{-4}$, their numerical values are $\mathcal {B}_{1111} = 0.25$ and $\mathcal {Z}_{11} = -0.57$. Then, with (6.8), the dimensionless particle pressure is

(6.10)\begin{equation} p=c \frac{b}{\bar{s}} \left(F-3\mathcal{Z}_{11} \right) = 8.51. \end{equation}

Given the force and torque balances, it is possible to characterize the role played by the torque balance in determining features of the trajectory and the distribution of particles along it. Ignoring the torque balance is equivalent to taking $S=0$, or $c_2=0$ in (2.27). This has an important influence on $\dot {\theta }$. Then, because both the distribution of near-contacting neighbours and the interval over which it is defined depend upon $\dot {\theta }$, there is a dependence of the stress upon it. For example, when $c=0.64$ and the torque balance is ignored, $p$ is $14.81,$ rather than $8.51$. That is, the value of $p$ is affected by the balance through the distribution. In contrast, $\tau$ is 18.60, rather than 18.94 and changes little, because it is independent of the shape of the distribution.

In the absence of the knowledge of the distribution function, it is possible to develop evolution equations for the approximate determination of its moments (e.g. Prantil, Jenkins & Dawson Reference Prantil, Jenkins and Dawson1993). Phan-Thien (Reference Phan-Thien1995) and Stickel et al. (Reference Stickel, Phillips and Powell2006) employ such an equation for the second moment, and break the symmetry of approach and departure by including a term in it that is proportional to $[\text {tr}(\boldsymbol {D}^{2})/2]^{1/2}$. Goddard (Reference Goddard2006) introduces a memory integral for the second moment – a representation for the solution to its evolution equation – that breaks this symmetry by incorporating a term proportional to the $\boldsymbol {D}^{2}$. Gillissen & Wilson (Reference Gillissen and Wilson2018, Reference Gillissen and Wilson2019) and Gillissen et al. (Reference Gillissen, Ness, Peterson, Wilson and Cates2019) distinguish between an anisotropic distribution of near contacts and an isotropic distribution of outer contacts, introduce a difference in association and disassociation of these contacts in the directions of compression and extension. This difference appears in the evolution equation of the second moment and provides the asymmetry necessary for normal stress differences. Theories of this type produce stress relations that are linear in the strain rate; that is, rate dependent. In contrast, we employ only a short-range repulsive force that is independent of the shear rate. Consequently, our stress relation contains contributions that are independent of rate.