1. Introduction
In a recent study, Jenkins & La Ragione (Reference Jenkins and La Ragione2015) determine the typical trajectory of a force-equilibrated pair of particles of a dense, two-dimensional suspensions of spheres subjected to a simple shearing flow. They evaluate the distribution function of such near-contacting neighbours along the trajectory and, using this distribution function and the expression for the force between the pair, they predict the particle pressure, the difference in normal stresses and the difference between the average rotation of the spheres and half the vorticity of the average velocity.
Here, we focus on extensional flows, also called pure shearing, of a dense layer of spheres and, as an extension of the previous work, also introduce the momentequilibrium. We employ a simplified Stokesian dynamics numerical simulation, perhaps more properly called lubrication dynamics (e.g. Ball & Melrose Reference Ball and Melrose1995), to guide the choice of the near-contacting pairs on a representative trajectory that contributes most to the inter-particle stress. We specify the boundary conditions on the representative trajectory, and determine the distribution of particles along it and the relationship between stress, microstructure and strain rate. We test these predictions against the results of the numerical simulations. We show that the relation between the stress and strain rate tensors involves the second and fourth moments of the particle distribution, and place this and other aspects of our approach in the context of the earlier models of Phan-Thien (Reference Phan-Thien1995), Stickel, Phillips & Powell (Reference Stickel, Phillips and Powell2006), Goddard (Reference Goddard2006), Gillissen & Wilson (Reference Gillissen and Wilson2018, Reference Gillissen and Wilson2019) and Gillissen et al. (Reference Gillissen, Ness, Peterson, Wilson and Cates2019) that focus on the second moment and that of Chacko et al. (Reference Chacko, Mari, Fielding and Cates2018), who introduce a fourth-rank tensor to describe flow reversal.
The approximate satisfaction of force and torque balances for particles in the flow plays an important role in what we do. In that regard, we operate in the spirit of Nazockdast & Morris (Reference Nazockdast and Morris2012a,Reference Nazockdast and Morrisb, Reference Nazockdast and Morris2013) or that of the statistical characterization by Thomas et al. (Reference Thomas, Ramola, Singh, Mari, Morris and Chakraborty2018) of equilibrated particles sheared in two dimensions, but in the limit of dense flows of the planar extensional flow. The analysis must be extended to three-dimensional simple shear flows before it can be placed in relation to phenomenological relations that have resulted from experiments on dense three-dimensional shearing flows (Boyer, Guazzelli & Pouliquen Reference Boyer, Guazzelli and Pouliquen2011; Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018).
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,
\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
\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
\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,
\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
\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
\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
\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
\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}$
\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
\begin{equation} F_{\alpha}^{(BA)} + \sum_{m{\neq} B}^{N^{(A)}}F_{\alpha}^{(mA)}=0; \end{equation}
and that for particle 
$B$ is
\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
\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
\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
\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
\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
\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
\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
\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
\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,
\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}
\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
\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
\begin{equation} b=-\frac{3\sqrt{3}\left(k-1\right) }{16{\rm \pi} }, \end{equation}by
\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,
\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
\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
\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
\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
\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
\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
\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)
\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,
\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
\begin{equation} \frac{\textrm{d}A}{\textrm{d}\theta }=-\frac{A}{\dot{\theta}}\frac{\textrm{d}\dot{\theta}}{\textrm{d}\theta }, \end{equation}with
\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
\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
\begin{equation} I\left( \theta \right) \equiv \int_{\theta _{0}}^{\theta }A(\theta^{\prime }) \textrm{d}\theta^{\prime}, \end{equation}as
\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),
\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
\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,
\begin{equation} \tau\equiv \tfrac{1}{2}\left({\mathsf{t}}_{11}-{\mathsf{t}}_{22}\right), \end{equation}is
\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,
\begin{equation} p\equiv -\tfrac{1}{2}\left({\mathsf{t}}_{11}+{\mathsf{t}}_{22}\right) , \end{equation}is
\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)),
\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$)
\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
\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)
\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
\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.
5. Results
We simulate monolayer systems with 1000 spheres of radius 
$a$ at area fractions, 
$c$, between 
$0.52$ and 
$0.64$. We generate initial configurations with a simple algorithm using random numbers. To reduce effects of such artificially generated initial configurations, the post-processing analyses use steady state data from 10 to 50 strain units. The repulsive force is set to be very weak 
$F_0/{\rm \pi} a^3 \mu \dot {\gamma } = 10^{-4}$ and short-ranged 
$\lambda /a = 10^{-2}$.
In the implementation of the model, we take 
$\bar {s} = 0.02$, 
$\theta _0 = 10^{-6}$, 
$s_0 = s_1 = 0.10$ and assume that 
$k$ varies linearly with 
$c$ from 2.0, at 
$c = 0.52,$ to 2.5, at 
$c = 0.64$. These values and the relation for 
$k$ are plausible and they are influenced by those measured in the simulations. The value 
$\bar {s}$ gives values of the shear stresses that are close to those measured. The value of 
$\theta _0$ is the default absolute tolerance of the solver; smaller values of 
$\theta _0$ have little influence on the shear stress, but do slightly improve the prediction of the pressure. The initial and final values of the separation were those employed in the simulations, and the variation of near-contacting neighbours with concentration, 
$k(c)$, was that measured in the simulation.
Figure 2 shows plots of the shear stress 
$\tau$ and pressure 
$p$ measured in the simulation and predicted by the model, over a range of area fraction 
$c.$ The stresses in the simulation increase in a similar manner with 
$c$, but the ratio 
$\tau / p$ decreases gradually. The predicted particle pressure is somewhat less than that measured in the simulations and the predicted shear stress is somewhat greater. The ratio of shear stress to pressure decreases with area fraction, as in the numerical simulations; but, because of under- and over-predicting, we have a greater value for the ratio.
Figure 2. Shear stress 
$\tau$ (black), pressure 
$p$ (red) and stress ratio (blue) over a range of concentrations, as measured in the simulation (solid) and predicted by the model (dashed).
Most of stress measured in the simulations is generated by closely approaching particles – defined as those with a separation less than one per cent of the particle radius. As seen in figure 3, more than 90 % of shear stress is generated from particle pairs with 
$s < 0.01$. Moreover, such near-contacting particles generates almost 100 % of the pressure 
$p$. Finally, approximately 80 % of the shear stress 
$\tau$ comes from the normal force.
Figure 3. (a) Near-contact contributions to 
$\tau$ and 
$p$. (b) The partial stress obtained using only the normal component of the pairwise force in (4.5).
We can check the concentration of stress contribution in the very narrow range of 
$s$ using distribution maps. We calculate the spatial distribution in 
$\xi \equiv \log s$. The statistics are calculated with discretized bins 
$\xi _{k} \equiv \xi _{1}+(k-1) \Delta \xi$, 
$k = 1, \dotsc , k_{{max}}$. 
$\xi _{1} = \log {10^{-7}}$ and 
$\xi _{k_{max}} = \log {10^{-1}}$. The results are plotted with 
$s$ in a logarithmic scale. Figure 4(a) displays the distribution of shear stress, indicating that the stress tends to be concentrated near the stagnation point 
$(\theta , s) \approx (0, 10^{-6})$. The region of concentration spreads until 
$\theta \sim {\rm \pi}/4$. We also separately calculate the stress of (4.5) constructed with normal forces 
$\boldsymbol {F}^{ij}_{{n }} \equiv \boldsymbol {F}^{ij}\boldsymbol {\cdot } \boldsymbol {n}^{ij} \boldsymbol {n}^{ij}$ and tangential forces 
$\boldsymbol {F}^{ij}_{{t}}\equiv \boldsymbol {F}^{ij}-\boldsymbol {F}^{ij}\boldsymbol {\cdot } \boldsymbol {n}^{ij} \boldsymbol {n}^{ij}$. As shown in figure 3(b), 80 % of the shear stress 
$\tau$ indeed comes from the normal forces. Besides systematic motions due to the shearing deformation, particle motions fluctuate due to occasional configurations of surrounding particles. Therefore, it is necessary to reconstruct averaged trajectories to compare with theoretical ones. To this end, we first calculate the averaged relative-velocity field 
$\langle \boldsymbol {U}^{(j)} - \boldsymbol {U}^{(i)}\rangle$ over all interacting pairs 
$i$ and 
$j$ in terms of the relative position coordinate 
$\Delta \boldsymbol {r}^{ij}=(2a+s, \theta )$. Owing to the symmetry of planar extension, the statistics are taken on a quadrant: 
$0 < \theta < {\rm \pi}/2$. Because we consider a situation that is very close to the singularity (Ball & Melrose Reference Ball and Melrose1995), the particles tend to approach very close to contact, i.e. a bundle of trajectories is compressed into an extremely narrow range of 
$s$. To avoid a loss of precision due to averaging, we carry out the statistical data binning with 
$\xi$ instead of 
$s$.
Figure 4. (a) Stress concentration, (b) trajectories, with colours denoting their average contribution to the stress. Both figures are for 
$c = 0.52$.
Once we evaluate the velocity field in the 
$\xi$–
$\theta$ space, i.e. 
$( \langle \dot {\xi } \rangle = \langle \dot {s}/s \rangle , \langle \dot {\theta } \rangle )$, we can obtain trajectories as streamlines of the velocity field. In figure 4(b), trajectories of the system with 
$c = 0.52$ and various initial positions are plotted. The trajectories are coloured from blue to red, according to their contribution to the stress. We identify the band of significant trajectories as those with a separation of less than 
$10^{-2}$.
In figure 5, we show the particle number density, measured in the simulations for 
$c = 0.64,$ with particle trajectories superposed. The near-contact distribution 
$A(\theta ),$ normalized by the total number of particles, is obtained by integrating across the trajectories in the region below 
$s = 10^{-2}$. The product of the integral of the particle probability distribution over the area of figure 5 below 
$s = 10^{-2}$ and the total number of particles is the coordination number.
Figure 5. Particle number density over 
$s$ and 
$\theta$, with particle trajectories superposed for 
$c = 0.64$.
In figure 6(a), we plot the trajectories from the simulation, for the smallest separation in figure 4(b), and the predicted representative trajectories of the model for concentrations of 0.52 and 0.64. The representative trajectories are located within the band and have a shape similar to the individual trajectories at smaller separations. In figure 6(b), we show distributions, 
$A(\theta )$, measured in the simulation and predicted along the representative trajectories for two values of the concentration. These share the same features and have a similar agreement as the trajectories. The asymmetry of the distributions result from the influence of the repulsive force that breaks the symmetry of approach and departure.
Figure 6. (a) Closest trajectories in the numerical simulation (solid) and the representative trajectories (dashed), for 
$c = 0.52$ (blue) and 
$c = 0.64$ (red). (b) The distribution 
$A(\theta )$ along the representative trajectory, as measured in the simulation (solid) and predicted by the model (dashed), for 
$c = 0.52$ (blue) and 
$c = 0.64$ (red).
We have employed information from the simulation on the variation in the coordination number as a function of concentration and the value of separation of near-contacting neighbours necessary to reproduce the measured particle shear stress. These, used in the model, gives it the capacity to generate particle trajectories that are representative of the stress-producing trajectories of the simulation and particle distributions along the representative trajectories with the appropriate asymmetry about 
${\rm \pi} /4$ to predict a reasonable variation of a particle pressure.
We next indicate how the structure of the model can be used as the basis for a continuum theory of dense suspensions and to provide a context for existing theories.
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
\begin{equation} {\mathsf{f}}_{\alpha \beta }=\skew4\hat{d}_{\alpha }\skew4\hat{d}_{\beta }-\tfrac{1}{2}\delta _{\alpha \beta } \end{equation}and
\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}
\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
\begin{equation} \mathcal{Z}_{\alpha \beta }=\int_{0}^{2{\rm \pi} }A(\theta){\mathsf{f}}_{\alpha \beta }\,\textrm{d}\theta \end{equation}and
\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
\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
\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
\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
\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
\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.
7. Conclusion
We have considered a planar extensional flow of a dense layer of spheres in a viscous fluid. In addition to the viscous forces associated with the flow, we assumed that there was a short-range repulsive force between the spheres. We focused on pairs of spheres, assumed that their neighbours translate with the average flow, and required that they be in force and moment equilibrium with each other and their neighbours. We then assumed that the neighbourhoods of pairs with the same orientation were equal to the average over that orientation; this permitted us to write equations for the radial and angular velocity of the relative motion of a single pair as they began and ended an interaction, and the orientation of their line of centres with the axis of greatest compression of the flow.
The possible determination of the distribution of near-contacting particles along the trajectory then leads to expressions for the particle shear stress and pressure. Stokesian dynamics simulations provided a value of particle shear stress that permitted the determination of the angle of departure for the trajectories and the separation of near-contacting neighbours. The variation of the shear stress with area fraction suggested the variation of the number of near-contacting neighbours per particle with area fraction. With this information, numerical values of the particle distribution and the particle pressure could be calculated.
The simplicity of the micro-mechanical model makes it possible to understand the influence of the normal and tangential components of the viscous force and the short-range repulsive force on the shape of the trajectory and the distribution of particles along it. It also permits the derivation of a continuum theory for dense suspensions based on the micro-mechanics of equilibrated pairs of particles and a better understanding of those based on phenomenology. However, it is important to note that correlated motions of the particles along the axis of compression in the simulations have no counterpart in the model; these are likely to force approaching particles into trajectories with significantly smaller separations. The small value of 
$\theta _0$ necessary in the model may be a way for it to incorporate the influence of such particles. Other modelling issues that remain involve the parameters 
$k$ and 
${\bar {s}}$. We believe that it is likely that measurements of the relationship between 
$k$ and 
$c$ in three dimensions will be robust and apply to different sharing flows; the determination of 
${\bar {s}}$ is less certain, although methods employed by Nazockdast & Morris (Reference Nazockdast and Morris2013) may permit its prediction.
The two-dimensional, simple shear flow that was considered earlier by Jenkins & La Ragione (Reference Jenkins and La Ragione2015) is much more complicated. In it, there are upstream and downstream trajectories, roughly, on either side of the direction of greatest rate of compression. The exact location of the limiting trajectory is influenced by the rotation of the average flow and, in the calculation, its location is difficult to determine. In the earlier study, the separation, was taken to be the average separation calculated by Torquato (Reference Torquato1995) for a dense, equilibrated system of colliding, elastic spheres in two dimensions. This varies between 0.13 and 0.06 as 
$c$ varies from 0.52 and 0.64. These values are much greater than the value 0.02 employed here. However, in the calculation for simple shear we obtain an excellent agreement between theory and simulation for the shear stress because the length of the particle interaction was taken to be much longer, and probably compensated for the greater separation. Given our experience with the simpler planar extension, we should return to the simple shear flow.
The formulation for the two-dimensional planar extensional flow can be easily extended to an axisymmetric extensional flow in three dimensions. The only essential change is in the components of the average strain rate tensor. The three-dimensional calculation can then incorporate particle elasticity and friction to permit a study of shear thickening and comparison with physical experiments.
Acknowledgements
R.S. acknowledges the support from the Wenzhou Institute, UCAS.
Funding
This research was partially supported by the National Science Foundation under Grant No. NSF PHY-1748958 and by the Gruppo Nazionale della Fisica Matematica (Italy).
Declaration of interests
The authors report no conflict of interest.
Appendix A
Force equilibrium equation for particle 
$A$ is
\begin{align} 0 &=6{\rm \pi} \mu a{\mathsf{K}}_{\alpha \beta }^{(BA) }v_{\beta }^{\left( BA\right) }-\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 ^{\left( A\right) }\hat{t}_{\alpha }^{(BA) }+{\rm \pi} \mu a^{2}\ln \left( \frac{a}{s^{(BA) }}\right) \omega ^{(B) }\hat{t}_{\alpha }^{(BA) }\nonumber\\ &\quad +\sum_{n{\neq} B}^{N^{(A) }}\left\{ \frac{3}{\bar{s}}a^{3}{\rm \pi} \mu \left({\mathsf{D}}_{\beta \xi }\skew4\hat{d}_{\xi }^{(nA) }\skew4\hat{d}_{\beta }^{(nA)}\right) \skew4\hat{d}_{\alpha }^{(nA) }+{\rm \pi} \mu a^{2}\left[ \ln \left( \frac{a}{\bar{s}}\right) -0.96\right] \omega ^{(A) } \hat{t}_{\alpha }^{(nA) }\right\} \nonumber\\ &\quad +\sum_{n{\neq} B}^{N^{(A) }}\left\{ 2a^{2}{\rm \pi} \mu \left[ \ln \left( \frac{a}{\bar{s}}\right) -0.96\right] \left({\mathsf{D}}_{\beta \xi }\hat{t} _{\xi }^{(nA) }\skew4\hat{d}_{\beta }^{(nA) }\right) \hat{t }_{\alpha }^{(nA) }-\frac{F_{0}}{\bar{s}}\skew4\hat{d}_{\alpha }^{(nA) }\right\} . \end{align} Similar expression can be written for particle 
$B.$ The difference between force equilibrium 
$A$ and 
$B$ leads to
\begin{align} 0 &=12{\rm \pi} \mu a{\mathsf{K}}_{\alpha \beta }^{(BA) }v_{\beta }^{(BA)} - 2\frac{F_{0}}{s^{(BA)}}\skew4\hat{d}_{\alpha }^{(BA) }-2\times 9.54{\rm \pi} \mu a^{2}\left(\hat{t}_{\beta }^{(BA) }{\mathsf{D}}_{\beta \xi }\skew4\hat{d}_{\xi }^{(BA) }\right) \hat{t}_{\alpha }^{(BA) } \nonumber\\ &\quad +a^{2}{\rm \pi} \mu \left[ 2\ln \left( \frac{a}{s^{(BA)}}\right) -0.96\right] S \hat{t}_{\alpha }^{(BA) }+a^{2}{\rm \pi} \mu \left[ \ln \left( \frac{a}{\bar{s}}\right) -0.96\right] S\varepsilon _{\alpha \beta }{\mathsf{Y}}_{\beta}^{(BA) } \nonumber\\ &\quad +2\frac{3}{\bar{s}}a^{3}{\rm \pi} \mu {\mathsf{D}}_{\beta \xi }{\mathsf{J}}_{\alpha \xi \beta }^{(BA)}+2a^{2}{\rm \pi} \mu \left[ 2\ln \left( \frac{a}{\bar{s}}\right) -1.92 \right] {\mathsf{D}}_{\beta \xi }{\mathsf{J}}_{\alpha \xi \beta }^{(BA)}-2\frac{F_{0}}{\bar{s}} {\mathsf{Y}}_{\alpha}^{(BA)}, \end{align}where
\begin{equation} \left.\begin{gathered} S=\omega^{(A)}+\omega^{(B) },\\ {\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{gathered}\right\} \end{equation}and
\begin{equation} {\mathsf{Y}}_{\alpha }^{(BA) }=\sum_{m{\neq} B}^{N^{(A) }}\skew4\hat{d} _{\alpha }^{(mA) }. \end{equation}
Moment equilibrium for particle 
$A$ is
\begin{align} 0 &=6{\rm \pi} \mu a{\mathsf{K}}_{\alpha \beta }^{(BA) }v_{\beta }^{(BA)} \varepsilon_{\alpha \rho }\skew4\hat{d}_{\rho }^{(BA) }~-9.54{\rm \pi} \mu a^{2}( \hat{t}_{\beta }{\mathsf{D}}_{\beta \xi }\skew4\hat{d}_{\xi }) \nonumber\\ &\quad +{\rm \pi} \mu a^{2}\left[ \ln \left( \frac{a}{s^{(BA) }} \right) -0.96\right] \omega^{(A) }+{\rm \pi} \mu a^{2}\ln \left( \frac{a}{s^{(BA) }}\right) \omega^{(B)} \nonumber\\ &\quad +\varepsilon_{\alpha \rho }\sum_{n{\neq} B}^{N^{(A)}}\left\{ \frac{3}{\bar{s}}a^{3}{\rm \pi} \mu \left({\mathsf{D}}_{\beta \xi }\skew4\hat{d}_{\xi}^{(nA)} \skew4\hat{d}_{\beta}^{(nA) }\right) \skew4\hat{d}_{\alpha}^{(nA)} \!+\!{\rm \pi} \mu a^{2}\left[ \ln \left( \frac{a}{\bar{s}}\right) \!-\!0.96\right] \omega^{(A) }\hat{t}_{\alpha }^{(nA)}\right\}\skew4\hat{d}_{\rho}^{(nA) } \nonumber\\ &\quad +\varepsilon _{\alpha \rho }\sum_{n{\neq} B}^{N^{(A) }}\left\{ 2a^{2}{\rm \pi} \mu \left[ \ln \left( \frac{a}{\bar{s}}\right) -0.96\right] \left({\mathsf{D}}_{\beta \xi }\hat{t}_{\xi }^{(nA) }\skew4\hat{d}_{\beta }^{(nA)}\right) \hat{t}_{\alpha }^{(nA)} -\frac{F_{0}}{\bar{s}}\skew4\hat{d}_{\alpha }^{(nA) }\right\} \skew4\hat{d}_{\rho }^{(nA) }. \end{align}
A similar expression holds for particle 
$B.$ The sum of moment equilibrium for particles 
$A$ and 
$B$ is
\begin{align} 0 &=12\varepsilon _{\alpha \beta }{\mathsf{K}}_{\alpha \mu }^{(BA) }v_{\mu }^{(BA) }\,\skew4\hat{d}_{\beta }^{(BA) }-2\times 9.54a\left( \hat{t}_{\mu }^{(BA) }{\mathsf{D}}_{\mu \xi }\skew4\hat{d}_{\xi }^{(BA) }\right) \nonumber\\ &\quad +a\left[ 2\ln \left( \frac{a}{s^{(BA)}}\right) -0.96\right] S+a\left[ \ln \left( \frac{a}{\bar{s}}\right) -0.96\right] S\left( k-1\right) \nonumber\\ &\quad +2a\left[ 2\ln \left( \frac{a}{\bar{s}}\right) -1.92\right] \epsilon _{\xi \nu }{\mathsf{A}}_{\nu \mu }^{(BA) }{\mathsf{D}}_{\mu \xi }, \end{align}where
\begin{equation} {\mathsf{A}}_{\nu \mu }^{(BA) }=\sum_{m{\neq} B}^{N^{(A) }}\skew4\hat{d}_{\mu}^{(mA) }\skew4\hat{d}_{\nu}^{(mA)}. \end{equation}In both the difference of force equilibrium and the sum of moment equilibrium, we have made the following approximations
\begin{equation} \left.\begin{gathered} {\mathsf{A}}_{\nu \mu }^{(BA) }={\mathsf{A}}_{\nu \mu }^{(AB) },\\ {\mathsf{J}}_{\alpha \xi \beta }^{(BA)}=-{\mathsf{J}}_{\alpha \xi \beta }^{(AB)}, \end{gathered}\right\} \end{equation}and
\begin{equation} {\mathsf{Y}}_{\alpha }^{(BA) }=-{\mathsf{Y}}_{\alpha }^{(AB) }. \end{equation} The projection of the difference in force equilibrium along the direction orthogonal to 
$t_{\alpha }^{(BA) }$ is
\begin{equation} 3{\rm \pi} \mu a\frac{a}{s^{(BA) }}\dot{s}^{(BA) }-2\frac{ F_{0}}{s^{(BA) }}+6{\rm \pi} \mu a^{2}\frac{a}{\bar{s}}\skew4\hat{d}_{\alpha }J_{\alpha \beta \gamma }{\mathsf{D}}_{\beta \gamma }-2\frac{F_{0}}{\bar{s}}{\mathsf{Y}}_{\alpha } \skew4\hat{d}_{\alpha }=0, \end{equation}
while the component along 
$t_{\alpha }^{(BA) }$is
\begin{align} 0 &=12{\rm \pi} \mu a{\mathsf{K}}_{\alpha \beta }^{(BA) }v_{\beta }^{\left( BA\right) }t_{\alpha }^{(BA) }-2\times 9.54{\rm \pi} \mu a^{2}\left( \hat{t}_{\beta }^{(BA) }{\mathsf{D}}_{\beta \xi }\skew4\hat{d}_{\xi }^{\left( BA\right) }\right) \nonumber\\ &\quad +a^{2}{\rm \pi} \mu \left[ 2\ln \left( \frac{a}{s^{(BA)}}\right) -0.96\right] S+a^{2}{\rm \pi} \mu \left[ \ln \left( \frac{a}{\bar{s}}\right) -0.96\right] S\varepsilon _{\alpha \beta }{\mathsf{Y}}_{\beta }^{(BA) }t_{\alpha }^{(BA) } \nonumber\\ &\quad +2\frac{3}{\bar{s}}a^{3}{\rm \pi} \mu {\mathsf{D}}_{\beta \xi }{\mathsf{J}}_{\alpha \xi \beta}^{(BA)}t_{\alpha}^{(BA) } +2a^{2}{\rm \pi} \mu \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}
Using (A11) and (A6), we solve for 
$S$
\begin{equation} S=-\frac{12b}{\left( 4b-k+1\right) \bar{s}[\ln ( 1/\bar{s}) -0.96] }\sin 2\theta, \end{equation}or
\begin{equation} S=-2c_{2}\sin 2\theta , \end{equation}where
\begin{equation} c_{2}=\frac{6b/\left[ \bar{s}\left( 4b-k+1\right) \right] }{\ln \left( 1/\bar{s}\right)-0.96}. \end{equation}Appendix B
The Matlab m-files for the numerical solution of the ordinary differential equations and boundary conditions follow.
