2.1 The 10-moment system

Going beyond the Navier–Stokes (NS)-order hydrodynamics of five field variables, here we will work with an ‘extended’ set of 10 hydrodynamic fields: (i) the mass density

(2.2) $$\begin{eqnarray}{\it\rho}(\boldsymbol{x},t)\equiv mn(\boldsymbol{x},t)=m\int f(\boldsymbol{c},\boldsymbol{x},t)\,\text{d}\boldsymbol{c},\end{eqnarray}$$

(ii) the coarse-grained velocity

(2.3) $$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)\equiv \langle \boldsymbol{c}\rangle =\frac{1}{n(\boldsymbol{x},t)}\int \boldsymbol{c}f(\boldsymbol{c},\boldsymbol{x},t)\,\text{d}\boldsymbol{c},\end{eqnarray}$$

and (iii) the second-moment tensor

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D648}(\boldsymbol{x},t)\equiv \langle \boldsymbol{C}\boldsymbol{C}\rangle =\frac{1}{n(\boldsymbol{x},t)}\int \boldsymbol{C}\boldsymbol{C}f(\boldsymbol{c},\boldsymbol{x},t)\,\text{d}\boldsymbol{c},\end{eqnarray}$$

where $n(\boldsymbol{x},t)$ is the number density and $\boldsymbol{C}\equiv \boldsymbol{c}-\boldsymbol{u}$ is the peculiar/fluctuation velocity of particles. The last hydrodynamic field (2.4) is required to account for normal stress differences (Jenkins & Richman Reference Jenkins and Richman1988; Chou & Richman Reference Chou and Richman1998; Saha & Alam Reference Saha and Alam2014) which are the primary focus of the present work. The balance equations for the mass, momentum and second moment, respectively, can be obtained by taking the appropriate moment of (2.1):

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\partial }{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\right){\it\rho}=-{\it\rho}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\left(\frac{\partial }{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\right)\boldsymbol{u}=-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}, & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\left(\frac{\partial }{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\right)\unicode[STIX]{x1D648}=-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}-\unicode[STIX]{x1D64B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}-(\unicode[STIX]{x1D64B}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u})^{\text{T}}+\boldsymbol{\aleph }, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D64B}$ is the total stress, a second-rank tensor, given by

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D64B}\equiv {\it\rho}\unicode[STIX]{x1D648}+{\it\bf\Theta}(m\boldsymbol{C}),\end{eqnarray}$$

$\unicode[STIX]{x1D64C}$ is the flux of the second moment, a third-rank tensor, given by

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D64C}\equiv {\it\rho}\langle \boldsymbol{C}\boldsymbol{C}\boldsymbol{C}\rangle +{\it\bf\Theta}(m\boldsymbol{C}\boldsymbol{C}),\end{eqnarray}$$

and $\boldsymbol{\aleph }$ is the collisional source of the second moment, a second-rank tensor, given by

(2.10) $$\begin{eqnarray}\boldsymbol{\aleph }\equiv \boldsymbol{\aleph }(m\boldsymbol{C}\boldsymbol{C}).\end{eqnarray}$$

In (2.8)–(2.9), the first and second terms refer to the corresponding kinetic and collisional contributions, respectively. The integral expression for the collisional flux term ${\it\bf\Theta}({\it\psi})$ of any particle-level property ${\it\psi}$ is given by

(2.11) $$\begin{eqnarray}\displaystyle {\it\bf\Theta}({\it\psi}) & = & \displaystyle -\frac{{\it\sigma}^{3}}{2}\iiint _{\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{k}>0}({\it\psi}_{1}^{\prime }-{\it\psi}_{1})\boldsymbol{k}\int _{0}^{1}f^{(2)}(\boldsymbol{c}_{1},\boldsymbol{x}-{\it\omega}{\it\sigma}\boldsymbol{k},\boldsymbol{c}_{2},\boldsymbol{x}+{\it\sigma}\boldsymbol{k}-{\it\omega}{\it\sigma}\boldsymbol{k})\nonumber\\ \displaystyle & & \displaystyle \times \,(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{g})\,\text{d}{\it\omega}\,\text{d}\boldsymbol{k}\,\text{d}\boldsymbol{c}_{1}\,\text{d}\boldsymbol{c}_{2},\end{eqnarray}$$

where ${\it\omega}$ is an integration variable (Jenkins & Richman Reference Jenkins and Richman1985), and that for the collisional source of ${\it\psi}$ , $\boldsymbol{\aleph }({\it\psi})$ , is

(2.12) $$\begin{eqnarray}\boldsymbol{\aleph }({\it\psi})=\frac{{\it\sigma}^{2}}{2}\iiint _{\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{k}>0}({\it\psi}_{1}^{\prime }+{\it\psi}_{2}^{\prime }-{\it\psi}_{1}-{\it\psi}_{2})f^{(2)}(\boldsymbol{c}_{1},\boldsymbol{x}-{\it\sigma}\boldsymbol{k},\boldsymbol{c}_{2},\boldsymbol{x})(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{g})\,\text{d}\boldsymbol{k}\,\text{d}\boldsymbol{c}_{1}\,\text{d}\boldsymbol{c}_{2}.\end{eqnarray}$$

Note that the origin of the collisional flux term (2.11) is tied to the ‘macroscopic’ nature of particles (and hence to the ‘denseness’ of the matter) and this term is zero in a dilute gas of point particles.

Defining the granular temperature as $T\equiv \unicode[STIX]{x1D614}_{{\it\alpha}{\it\alpha}}/3$ and taking the trace of (2.7), we obtain the balance equation for the granular energy

(2.13) $$\begin{eqnarray}\frac{3}{2}{\it\rho}\left(\frac{\partial }{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\right)T=-\frac{\partial q_{{\it\alpha}}}{\partial x_{{\it\alpha}}}-\unicode[STIX]{x1D617}_{{\it\alpha}{\it\beta}}\frac{\partial u_{{\it\beta}}}{\partial x_{{\it\alpha}}}-\mathscr{D},\end{eqnarray}$$

and that of the deviator of the second moment

(2.14) $$\begin{eqnarray}\displaystyle \frac{1}{2}{\it\rho}\left(\frac{\partial }{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\right)\widehat{\unicode[STIX]{x1D614}}_{{\it\alpha}{\it\beta}} & = & \displaystyle -\frac{\partial }{\partial x_{{\it\gamma}}}\left(\unicode[STIX]{x1D618}_{{\it\gamma}{\it\alpha}{\it\beta}}-\frac{2}{3}q_{{\it\gamma}}{\it\delta}_{{\it\alpha}{\it\beta}}\right)\nonumber\\ \displaystyle & & \displaystyle -\left\{\frac{1}{2}\left(\unicode[STIX]{x1D617}_{{\it\gamma}{\it\alpha}}\frac{\partial u_{{\it\beta}}}{\partial x_{{\it\gamma}}}+\unicode[STIX]{x1D617}_{{\it\gamma}{\it\beta}}\frac{\partial u_{{\it\alpha}}}{\partial x_{{\it\gamma}}}\right)-\frac{1}{3}\unicode[STIX]{x1D617}_{{\it\gamma}{\it\xi}}\frac{\partial u_{{\it\xi}}}{\partial x_{{\it\gamma}}}{\it\delta}_{{\it\alpha}{\it\beta}}\right\}+\frac{1}{2}\widehat{\aleph }_{{\it\alpha}{\it\beta}}.\qquad \quad\end{eqnarray}$$

In the above equations,

(2.15) $$\begin{eqnarray}q_{{\it\alpha}}\equiv {\textstyle \frac{1}{2}}\unicode[STIX]{x1D618}_{{\it\alpha}{\it\beta}{\it\beta}}={\textstyle \frac{1}{2}}{\it\rho}\unicode[STIX]{x1D614}_{{\it\alpha}{\it\beta}{\it\beta}}+{\textstyle \frac{1}{2}}{\it\Theta}_{{\it\alpha}{\it\beta}{\it\beta}}\end{eqnarray}$$

is the heat flux vector, and

(2.16) $$\begin{eqnarray}\mathscr{D}\equiv -{\textstyle \frac{1}{2}}\aleph _{{\it\beta}{\it\beta}}=-{\textstyle \frac{1}{2}}\boldsymbol{\aleph }(mC^{2})\end{eqnarray}$$

is the rate of dissipation of kinetic energy per unit volume. The balance equations (2.5)–(2.6) and (2.13), along with constitutive relations for (2.8), (2.15) and (2.16), constitute the NS-order hydrodynamics for which the equation for the deviatoric part of the second-moment tensor (2.14) is satisfied identically.

For an extended hydrodynamic description of granular fluids, incorporating normal stress differences, the balance equation (2.7) for the full second-moment tensor along with the mass and momentum balances (2.5)–(2.6) are needed. For a closure of (2.7), the deviatoric part of the third order $\unicode[STIX]{x1D618}_{{\it\gamma}{\it\alpha}{\it\beta}}$ ,

(2.17) $$\begin{eqnarray}\widehat{\unicode[STIX]{x1D618}}_{{\it\gamma}{\it\alpha}{\it\beta}}=\unicode[STIX]{x1D618}_{{\it\gamma}{\it\alpha}{\it\beta}}-{\textstyle \frac{1}{5}}(\unicode[STIX]{x1D618}_{{\it\gamma}{\it\xi}{\it\xi}}{\it\delta}_{{\it\alpha}{\it\beta}}+\unicode[STIX]{x1D618}_{{\it\alpha}{\it\xi}{\it\xi}}{\it\delta}_{{\it\gamma}{\it\beta}}+\unicode[STIX]{x1D618}_{{\it\beta}{\it\xi}{\it\xi}}{\it\delta}_{{\it\alpha}{\it\gamma}}),\end{eqnarray}$$

is assumed to be zero, leaving only its isotropic part, the heat flux vector (2.15), to be evaluated as a constitutive relation. In this paper, we specifically focus on the uniform shear flow (see § 3) for which (2.15) is zero; therefore we are left to determine constitutive relations for the stress tensor (2.8) and the source of the second moment (2.10) as discussed in §§ 3–5.

2.2 Anisotropic Gaussian distribution function and the maximum entropy principle

To evaluate (2.11) and (2.12), we adopt the molecular chaos ansatz (Chapman & Cowling Reference Chapman and Cowling1970)

(2.18) $$\begin{eqnarray}f^{(2)}(\boldsymbol{c}_{1},\boldsymbol{x}-{\it\sigma}\boldsymbol{k},\boldsymbol{c}_{2},\boldsymbol{x})=g_{0}({\it\nu})f(\boldsymbol{c}_{1},\boldsymbol{x}-{\it\sigma}\boldsymbol{k})f(\boldsymbol{c}_{2},\boldsymbol{x}),\end{eqnarray}$$

that relates the two-particle distribution function as a product of two single-particle velocity distribution functions. The positional correlation is taken care of via the well-known contact radial distribution function of Carnahan & Starling (Reference Carnahan and Starling1969),

(2.19) $$\begin{eqnarray}g_{0}({\it\nu})=\frac{(2-{\it\nu})}{2(1-{\it\nu})^{3}},\end{eqnarray}$$

with

(2.20) $$\begin{eqnarray}{\it\nu}=n{\rm\pi}{\it\sigma}^{3}/6\end{eqnarray}$$

being the local volume fraction of particles. Finally we assume that the single-particle velocity distribution is an anisotropic Maxwellian/Gaussian

(2.21) $$\begin{eqnarray}f(\boldsymbol{c},\boldsymbol{x},t)=\frac{n}{(8{\rm\pi}^{3}|\unicode[STIX]{x1D648}|)^{1/2}}\exp \left(-\frac{1}{2}\boldsymbol{C}\boldsymbol{\cdot }\unicode[STIX]{x1D648}^{-1}\boldsymbol{\cdot }\boldsymbol{C}\right),\end{eqnarray}$$

with $|\unicode[STIX]{x1D648}|=\det (\unicode[STIX]{x1D648})$ , which contains complete information about the second-moment tensor $\unicode[STIX]{x1D648}$ . This form was originally used by Holway (Reference Holway1966) to improve certain problems in the BGK (Bhatnagar–Gross–Krook) model of gas dynamics, resulting in what is popularly known as the ‘Ellipsoidal’ BGK model. In fact it is straightforward to show that (2.21) follows from the maximum entropy principle (Jaynes Reference Jaynes1957; Holway Reference Holway1966) when an extended set of hydrodynamic fields ( ${\it\rho}$ , $\boldsymbol{u}$ and $\unicode[STIX]{x1D648}$ ) is assumed to hold.

While (2.21) is an appropriate leading-order distribution function for a non-equilibrium steady state, such as steady uniform shear flow (Goldreich & Tremaine Reference Goldreich and Tremaine1978; Araki & Tremaine Reference Araki and Tremaine1986; Araki Reference Araki1988; Jenkins & Richman Reference Jenkins and Richman1988; Richman Reference Richman1989; Chou & Richman Reference Chou and Richman1998; Lutsko Reference Lutsko2004; Saha & Alam Reference Saha and Alam2014), the isotropic Maxwellian/Gaussian

(2.22) $$\begin{eqnarray}f(\boldsymbol{c},\boldsymbol{x},t)=\frac{n}{(2{\rm\pi}T)^{3/2}}\exp \left(-\frac{C^{2}}{2T}\right),\end{eqnarray}$$

(i.e. (2.21) with $\unicode[STIX]{x1D614}_{{\it\alpha}{\it\beta}}=T{\it\delta}_{{\it\alpha}{\it\beta}}$ ) holds for the rest state of a gas at equilibrium.

3.1 Construction of the second-moment tensor from its eigenvectors

Recalling that the granular temperature $T=\unicode[STIX]{x1D614}_{{\it\alpha}{\it\alpha}}/3$ is the isotropic measure of the second-moment tensor $\unicode[STIX]{x1D648}$ , we can decompose $\unicode[STIX]{x1D648}$ into an isotropic tensor and a traceless deviatoric tensor:

(3.5) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D648}}{T}=\unicode[STIX]{x1D644}+\frac{\widehat{\unicode[STIX]{x1D648}}}{T},\end{eqnarray}$$

where $\widehat{\unicode[STIX]{x1D648}}/T$ is the dimensionless counterpart of its deviatoric/traceless tensor whose eigenvalues ${\it\xi}$ , ${\it\varsigma}$ and ${\it\zeta}$ satisfy ${\it\xi}+{\it\varsigma}+{\it\zeta}=0$ . From (3.5) it follows that the eigenvalues of $\unicode[STIX]{x1D648}$ are $T(1+{\it\xi})$ , $T(1+{\it\varsigma})$ and $T(1+{\it\zeta})$ , and let us assume that the corresponding orthonormal set of eigendirections are $|\unicode[STIX]{x1D614}_{1}\rangle$ , $|\unicode[STIX]{x1D614}_{2}\rangle$ and $|\unicode[STIX]{x1D614}_{3}\rangle$ (see figure 2), respectively. Therefore, we can express the second-moment tensor $\unicode[STIX]{x1D648}$ as

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D648}=T(1+{\it\xi})|\unicode[STIX]{x1D614}_{1}\rangle \langle \unicode[STIX]{x1D614}_{1}|+T(1+{\it\varsigma})|\unicode[STIX]{x1D614}_{2}\rangle \langle \unicode[STIX]{x1D614}_{2}|+T(1+{\it\zeta})|\unicode[STIX]{x1D614}_{3}\rangle \langle \unicode[STIX]{x1D614}_{3}|.\end{eqnarray}$$

Referring to figure 2, we assume that the shear-plane eigenvectors $|\unicode[STIX]{x1D614}_{1}\rangle$ and $|\unicode[STIX]{x1D614}_{2}\rangle$  can be obtained by rotating the system of axes at an angle $({\rm\pi}/4+{\it\phi})$ , with ${\it\phi}$ being unknown, in the anticlockwise sense about the $z$ -axis which coincides with $|\unicode[STIX]{x1D614}_{3}\rangle$ :

(3.7a-c ) $$\begin{eqnarray}|\unicode[STIX]{x1D614}_{1}\rangle =\left[\begin{array}{@{}c@{}}\cos \left({\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right)\\ \sin \left({\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right)\\ 0\end{array}\right],\quad |\unicode[STIX]{x1D614}_{2}\rangle =\left[\begin{array}{@{}c@{}}-\text{sin}\left({\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right)\\ \cos \left({\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right)\\ 0\end{array}\right]\quad \text{and}\quad |\unicode[STIX]{x1D614}_{3}\rangle =\left[\begin{array}{@{}c@{}}0\\ 0\\ 1\end{array}\right].\end{eqnarray}$$

We further assume that the contact vector $\boldsymbol{k}$ makes an angle ${\it\varphi}$ with $|\unicode[STIX]{x1D614}_{3}\rangle \,$ , and ${\it\theta}$ is the angle between $|\unicode[STIX]{x1D614}_{1}\rangle$ and $\boldsymbol{k}-(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{z})\boldsymbol{z}$ , the projection of $\boldsymbol{k}$ on the shear plane, as shown in figure 2. Inserting (3.7) into (3.6), we obtain the following expression for the second-moment tensor

(3.8) $$\begin{eqnarray}\unicode[STIX]{x1D648}=T[{\it\delta}_{{\it\alpha}{\it\beta}}]+\widehat{\unicode[STIX]{x1D648}},\end{eqnarray}$$

and its deviatoric part is

(3.9) $$\begin{eqnarray}\widehat{\unicode[STIX]{x1D648}}=T\left[\begin{array}{@{}ccc@{}}{\it\lambda}^{2}+{\it\eta}\sin 2{\it\phi} & -{\it\eta}\cos 2{\it\phi} & 0\\ -{\it\eta}\cos 2{\it\phi} & {\it\lambda}^{2}-{\it\eta}\sin 2{\it\phi} & 0\\ 0 & 0 & -2{\it\lambda}^{2}\end{array}\right],\end{eqnarray}$$

where we have introduced the following notation

(3.10a,b ) $$\begin{eqnarray}{\it\eta}\equiv \frac{1}{2}({\it\varsigma}-{\it\xi})\geqslant 0\quad \text{and}\quad {\it\lambda}^{2}\equiv \frac{1}{2}({\it\varsigma}+{\it\xi})=-\frac{{\it\zeta}}{2}\geqslant 0.\end{eqnarray}$$

Note that while ${\it\eta}\sim (T_{x}-T_{y})$ is a measure of the anisotropy of the second-moment tensor in the shear plane, ${\it\lambda}^{2}$ is a measure of the excess temperature ( ${\sim}(T-T_{z})$ , where $T_{i}$ is the temperature along the $i$ th direction, see § 6.1) along the mean vorticity direction. It is straightforward to verify that the eigenvalues in the shear plane can be expressed in terms of ${\it\eta}$ and ${\it\lambda}$ via

(3.11a,b ) $$\begin{eqnarray}{\it\xi}={\it\lambda}^{2}-{\it\eta}\quad \text{and}\quad {\it\varsigma}={\it\lambda}^{2}+{\it\eta}>{\it\xi},\end{eqnarray}$$

with the eigenvalue, ${\it\zeta}$ , along the vorticity direction ( $z$ ), being given by (3.10).

Let us define the dimensionless shear rate (Savage & Jeffrey Reference Savage and Jeffrey1981)

(3.12) $$\begin{eqnarray}R\equiv \frac{\dot{{\it\gamma}}}{4\sqrt{T/{\it\sigma}^{2}}}=\left(\frac{\sqrt{T}}{{\it\sigma}\dot{{\it\gamma}}/4}\right)^{-1}\equiv \frac{v_{sh}}{v_{th}},\end{eqnarray}$$

which can be interpreted as the inverse of the square root of dimensionless temperature. Equation (3.12) is called the Savage–Jeffrey parameter (Savage & Jeffrey Reference Savage and Jeffrey1981) which is a measure of the mean shear velocity ( $v_{sh}={\it\sigma}\dot{{\it\gamma}}/2$ over a particle diameter) relative to the thermal velocity ( $v_{th}\propto \sqrt{T}$ ) associated with the random motion of particles. The second-moment tensor (3.8) in USF, constructed from its eigenbasis, is therefore completely determined when $R$ , ${\it\eta}$ , ${\it\phi}$ and ${\it\lambda}^{2}$ are specified.

3.2 The balance of second moment in USF

In steady uniform shear flow, the number density $n$ , the velocity gradient $\boldsymbol{{\rm\nabla}}\boldsymbol{u}$ and the components of the second-moment tensor $\unicode[STIX]{x1D648}$ are constants and the contracted third moment vanishes. Therefore, the mass and momentum balance equations, (2.5) and (2.6), are identically satisfied. The remaining balance equation (2.7), $\unicode[STIX]{x1D617}_{{\it\delta}{\it\beta}}u_{{\it\alpha},{\it\delta}}+\unicode[STIX]{x1D617}_{{\it\delta}{\it\alpha}}u_{{\it\beta},{\it\delta}}=\aleph _{{\it\alpha}{\it\beta}}$ , for the second-moment tensor simplifies to

(3.13) $$\begin{eqnarray}{\it\rho}\unicode[STIX]{x1D614}_{{\it\delta}{\it\beta}}(\unicode[STIX]{x1D60B}_{{\it\alpha}{\it\delta}}+\unicode[STIX]{x1D61E}_{{\it\alpha}{\it\delta}})+{\it\rho}\unicode[STIX]{x1D614}_{{\it\delta}{\it\alpha}}(\unicode[STIX]{x1D60B}_{{\it\beta}{\it\delta}}+\unicode[STIX]{x1D61E}_{{\it\beta}{\it\delta}})+{\it\Theta}_{{\it\delta}{\it\beta}}\unicode[STIX]{x1D60B}_{{\it\alpha}{\it\delta}}+{\it\Theta}_{{\it\delta}{\it\alpha}}\unicode[STIX]{x1D60B}_{{\it\beta}{\it\delta}}=\unicode[STIX]{x1D608}_{{\it\alpha}{\it\beta}}+\widehat{\unicode[STIX]{x1D60C}}_{{\it\alpha}{\it\beta}}+\widehat{\unicode[STIX]{x1D60E}}_{{\it\alpha}{\it\beta}},\end{eqnarray}$$

where we made use of $u_{{\it\alpha},{\it\delta}}=(\unicode[STIX]{x1D60B}_{{\it\alpha}{\it\delta}}+\unicode[STIX]{x1D61E}_{{\it\alpha}{\it\delta}})$ , $\unicode[STIX]{x1D617}_{{\it\alpha}{\it\beta}}={\it\rho}\unicode[STIX]{x1D614}_{{\it\alpha}{\it\beta}}+{\it\Theta}_{{\it\alpha}{\it\beta}}$ and the following decomposition for the collisional source of second moment (Jenkins & Richman Reference Jenkins and Richman1988; Chou & Richman Reference Chou and Richman1998; Saha & Alam Reference Saha and Alam2014):

(3.14) $$\begin{eqnarray}\aleph _{{\it\alpha}{\it\beta}}=\unicode[STIX]{x1D608}_{{\it\alpha}{\it\beta}}+\widehat{\unicode[STIX]{x1D60C}}_{{\it\alpha}{\it\beta}}+\widehat{\unicode[STIX]{x1D60E}}_{{\it\alpha}{\it\beta}}+{\it\Theta}_{{\it\alpha}{\it\delta}}\unicode[STIX]{x1D61E}_{{\it\beta}{\it\delta}}+{\it\Theta}_{{\it\beta}{\it\delta}}\unicode[STIX]{x1D61E}_{{\it\alpha}{\it\delta}}.\end{eqnarray}$$

The integral expression for the collisional stress tensor is

(3.15) $$\begin{eqnarray}{\it\Theta}_{{\it\alpha}{\it\beta}}=\frac{3(1+e){\it\rho}{\it\nu}g_{0}({\it\nu})}{{\rm\pi}^{3/2}}\int k_{{\it\alpha}}k_{{\it\beta}}(\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D648}\boldsymbol{\cdot }\boldsymbol{k})\mathfrak{G}({\it\chi})\,\text{d}\boldsymbol{k},\end{eqnarray}$$

the tensor $\unicode[STIX]{x1D608}_{{\it\alpha}{\it\beta}}$ is

(3.16) $$\begin{eqnarray}\unicode[STIX]{x1D608}_{{\it\alpha}{\it\beta}}=-\frac{6(1-e^{2}){\it\rho}{\it\nu}g_{0}({\it\nu})}{{\it\sigma}{\rm\pi}^{3/2}}\int k_{{\it\alpha}}k_{{\it\beta}}(\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D648}\boldsymbol{\cdot }\boldsymbol{k})^{3/2}\mathfrak{F}({\it\chi})\,\text{d}\boldsymbol{k},\end{eqnarray}$$

and the traceless tensors (denoted by a hat), $\widehat{\unicode[STIX]{x1D60C}}_{{\it\alpha}{\it\beta}}$ and $\widehat{\unicode[STIX]{x1D60E}}_{{\it\alpha}{\it\beta}}$ , are

(3.17) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{\unicode[STIX]{x1D60C}}_{{\it\alpha}{\it\beta}}=-\frac{12(1+e){\it\rho}{\it\nu}g_{0}}{{\it\sigma}{\rm\pi}^{3/2}}\int (k_{{\it\alpha}}j_{{\it\beta}}+j_{{\it\alpha}}k_{{\it\beta}})(\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D648}\boldsymbol{\cdot }\boldsymbol{j})(\unicode[STIX]{x1D660}\boldsymbol{\cdot }\unicode[STIX]{x1D648}\boldsymbol{\cdot }\boldsymbol{k})^{1/2}\mathfrak{F}({\it\chi})\,\text{d}\boldsymbol{k}, & \displaystyle\end{eqnarray}$$

(3.18) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{\unicode[STIX]{x1D60E}}_{{\it\alpha}{\it\beta}}=\frac{6(1+e){\it\rho}{\it\nu}g_{0}}{{\rm\pi}^{3/2}}\!\int (k_{{\it\alpha}}j_{{\it\beta}}+j_{{\it\alpha}}k_{{\it\beta}})[(\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D648}\boldsymbol{\cdot }\boldsymbol{k})(\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\boldsymbol{j})-(\boldsymbol{k}\boldsymbol{\cdot }\widehat{\unicode[STIX]{x1D648}}\boldsymbol{\cdot }\boldsymbol{j})(\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\boldsymbol{k})]\mathfrak{G}({\it\chi})\,\text{d}\boldsymbol{k}.\hspace{-3.60004pt} & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The integrals (3.15)–(3.18) are to be evaluated over $\text{d}\boldsymbol{k}$ such that $\text{d}\boldsymbol{k}=\sin {\it\varphi}\,\text{d}{\it\varphi}\,\text{d}{\it\theta}$ , with the limits of the integrations being ${\it\theta}\in (0,2{\rm\pi})$ and ${\it\varphi}\in (0,{\rm\pi})$ .

In (3.17)–(3.18), $\boldsymbol{j}$ represents an unit vector perpendicular to the contact vector $\boldsymbol{k}$ that lies in the plane formed by $\boldsymbol{g}$ and $\boldsymbol{k}$ such that

(3.19a ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{k}=\left[\begin{array}{@{}c@{}}\cos \left({\it\theta}+{\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right)\sin {\it\varphi}\\ \sin \left({\it\theta}+{\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right)\sin {\it\varphi}\\ \cos {\it\varphi}\end{array}\!\!\!\right] & \displaystyle\end{eqnarray}$$

and

(3.19b ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}=\frac{1}{\sqrt{2}}\left[\begin{array}{@{}c@{}}\cos {\it\varphi}\cos \left({\it\theta}+{\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right)-\sin \left({\it\theta}+{\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right)\\ \cos {\it\varphi}\sin \left({\it\theta}+{\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right)+\cos \left({\it\theta}+{\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right)\\ -\text{sin}\,{\it\varphi}\end{array}\!\!\!\right]. & \displaystyle\end{eqnarray}$$

The following two analytic functions (Araki & Tremaine Reference Araki and Tremaine1986) appear in the integrand of (3.15)–(3.18):

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle \mathfrak{F}({\it\chi})\equiv -\sqrt{{\rm\pi}}({\textstyle \frac{3}{2}}{\it\chi}+{\it\chi}^{3})\text{erfc}({\it\chi})+(1+{\it\chi}^{2})\exp (-{\it\chi}^{2}), & \displaystyle\end{eqnarray}$$

(3.21) $$\begin{eqnarray}\displaystyle & \displaystyle \mathfrak{G}({\it\chi})\equiv \sqrt{{\rm\pi}}({\textstyle \frac{1}{2}}+{\it\chi}^{2})\text{erfc}({\it\chi})-{\it\chi}\exp (-{\it\chi}^{2}), & \displaystyle\end{eqnarray}$$

where

(3.22) $$\begin{eqnarray}{\it\chi}(R,{\it\eta},{\it\phi},{\it\lambda};{\it\theta},{\it\varphi})\equiv \frac{{\it\sigma}(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{k})}{2\sqrt{(\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D648}\boldsymbol{\cdot }\boldsymbol{ k})}}=\frac{2R\sin ^{2}{\it\varphi}\cos (2{\it\phi}+2{\it\theta})}{\sqrt{(1-{\it\eta}\sin ^{2}{\it\varphi}\cos 2{\it\theta}+{\it\lambda}^{2}(3\sin ^{2}{\it\varphi}-2))}}.\quad\end{eqnarray}$$

The origin of ${\it\chi}$ can be traced to the excluded volume effects of macroscopic particles (Jenkins & Richman Reference Jenkins and Richman1988; Saha & Alam Reference Saha and Alam2014), and hence ${\it\chi}=0$ in the dilute limit and, consequently,

(3.23a,b ) $$\begin{eqnarray}\mathfrak{F}({\it\chi})=1\quad \text{and}\quad \mathfrak{G}({\it\chi})=\frac{\sqrt{{\rm\pi}}}{2},\quad \text{as }{\it\nu}\rightarrow 0.\end{eqnarray}$$

3.3 Second-moment balance in rotated coordinate frame

Let us now rewrite (3.13) in a new coordinate system $x^{\prime }y^{\prime }z^{\prime }$ , formed by the orthonormal triad of eigenvectors of $\unicode[STIX]{x1D648}$ , i.e., with respect to the coordinate system whose axes coincide with $|\unicode[STIX]{x1D614}_{1}\rangle \,$ , $|\unicode[STIX]{x1D614}_{2}\rangle$ and $|\unicode[STIX]{x1D614}_{3}\rangle \,$ , respectively. This amounts to a transformation, see figure 2, via the following rotation matrix,

(3.24) $$\begin{eqnarray}\mathscr{R}=\left[\begin{array}{@{}ccc@{}}\cos \left({\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right) & -\!\sin \left({\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right) & 0\\ \sin \left({\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right) & \cos \left({\it\phi}+{\displaystyle \frac{{\rm\pi}}{4}}\right) & 0\\ 0 & 0 & 1\end{array}\right],\end{eqnarray}$$

that transforms the second-moment tensor,

(3.25) $$\begin{eqnarray}\unicode[STIX]{x1D648}^{\prime }=T\left[\begin{array}{@{}ccc@{}}1+{\it\lambda}^{2}-{\it\eta} & 0 & 0\\ 0 & 1+{\it\lambda}^{2}+{\it\eta} & 0\\ 0 & 0 & 1-2{\it\lambda}^{2}\end{array}\right],\end{eqnarray}$$

into a diagonal matrix. With a prime over a quantity denoting its value in the new coordinate frame, the following relations hold:

(3.26) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{k}^{\prime }=\cos {\it\theta}\sin {\it\varphi}|\unicode[STIX]{x1D614}_{1}\rangle +\sin {\it\theta}\sin {\it\varphi}|\unicode[STIX]{x1D614}_{2}\rangle +\cos {\it\varphi}|\unicode[STIX]{x1D614}_{3}\rangle , & \displaystyle\end{eqnarray}$$

(3.27) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}^{\prime }=\frac{1}{\sqrt{2}}[(\cos {\it\varphi}\cos {\it\theta}-\sin {\it\theta})|\unicode[STIX]{x1D614}_{1}\rangle +(\cos {\it\varphi}\sin {\it\theta}+\cos {\it\theta})|\unicode[STIX]{x1D614}_{2}\rangle -\sin {\it\varphi}|\unicode[STIX]{x1D614}_{3}\rangle ], & \displaystyle\end{eqnarray}$$

(3.28) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}^{\prime }=2\dot{{\it\gamma}}\left[x^{\prime }\sin \left({\it\phi}+\frac{{\rm\pi}}{4}\right)+y^{\prime }\cos \left({\it\phi}+\frac{{\rm\pi}}{4}\right)\right]\!\left[\cos \left({\it\phi}+\frac{{\rm\pi}}{4}\right)\!|\unicode[STIX]{x1D614}_{1}\rangle -\sin \left({\it\phi}+\frac{{\rm\pi}}{4}\right)\!|\unicode[STIX]{x1D614}_{2}\rangle \right],\hspace{-3.60004pt} & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.29a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}^{\prime }=\dot{{\it\gamma}}\left[\begin{array}{@{}ccc@{}}\cos 2{\it\phi} & -\text{sin}\,2{\it\phi} & 0\\ -\text{sin}\,2{\it\phi} & -\text{cos}\,2{\it\phi} & 0\\ 0 & 0 & 0\end{array}\right]\quad \text{and}\quad \unicode[STIX]{x1D652}^{\prime }=\unicode[STIX]{x1D652}. & \displaystyle\end{eqnarray}$$

The last equation confirms that the spin tensor $\unicode[STIX]{x1D652}$ is invariant under the planar rotation (3.24).

With the aid of (3.25)–(3.29), the second-moment balance equation (3.13) transforms into four independent equations in the rotated coordinate frame: (i) the trace of (3.13),

(3.30) $$\begin{eqnarray}-4{\it\eta}{\it\rho}T\dot{{\it\gamma}}\cos 2{\it\phi}+2\dot{{\it\gamma}}[({\it\Theta}_{x^{\prime }x^{\prime }}-{\it\Theta}_{y^{\prime }y^{\prime }})\cos 2{\it\phi}-2{\it\Theta}_{x^{\prime }y^{\prime }}\sin 2{\it\phi}]=\unicode[STIX]{x1D608}_{x^{\prime }x^{\prime }}+\unicode[STIX]{x1D608}_{y^{\prime }y^{\prime }}+\unicode[STIX]{x1D608}_{z^{\prime }z^{\prime }},\end{eqnarray}$$

(ii) the $z^{\prime }$ – $z^{\prime }$ component of its deviatoric part

(3.31) $$\begin{eqnarray}-4{\it\eta}{\it\rho}T\dot{{\it\gamma}}\cos 2{\it\phi}+2\dot{{\it\gamma}}[({\it\Theta}_{x^{\prime }x^{\prime }}-{\it\Theta}_{y^{\prime }y^{\prime }})\cos 2{\it\phi}-2{\it\Theta}_{x^{\prime }y^{\prime }}\sin 2{\it\phi}]=-3\widehat{{\it\Gamma}}_{z^{\prime }z^{\prime }},\end{eqnarray}$$

(iii) the difference between the $x^{\prime }$ – $x^{\prime }$ and $y^{\prime }$ – $y^{\prime }$ components

(3.32) $$\begin{eqnarray}4(1+{\it\lambda}^{2}){\it\rho}T\dot{{\it\gamma}}\cos 2{\it\phi}+2\dot{{\it\gamma}}({\it\Theta}_{x^{\prime }x^{\prime }}+{\it\Theta}_{y^{\prime }y^{\prime }})\cos 2{\it\phi}={\it\Gamma}_{x^{\prime }x^{\prime }}-{\it\Gamma}_{y^{\prime }y^{\prime }},\end{eqnarray}$$

and, finally, (iv) the off-diagonal $x^{\prime }$ – $y^{\prime }$ component

(3.33) $$\begin{eqnarray}2{\it\rho}T\dot{{\it\gamma}}[{\it\eta}-(1+{\it\lambda}^{2})\sin 2{\it\phi}]-({\it\Theta}_{x^{\prime }x^{\prime }}+{\it\Theta}_{y^{\prime }y^{\prime }})\dot{{\it\gamma}}\sin 2{\it\phi}={\it\Gamma}_{x^{\prime }y^{\prime }},\end{eqnarray}$$

where

(3.34) $$\begin{eqnarray}{\it\Gamma}_{{\it\alpha}{\it\beta}}=\unicode[STIX]{x1D608}_{{\it\alpha}{\it\beta}}+\widehat{\unicode[STIX]{x1D60C}}_{{\it\alpha}{\it\beta}}+\widehat{\unicode[STIX]{x1D60E}}_{{\it\alpha}{\it\beta}},\end{eqnarray}$$

see (3.14)–(3.18). Various collision integrals ( ${\it\Theta}_{{\it\alpha}^{\prime }{\it\beta}^{\prime }}$ , $\unicode[STIX]{x1D608}_{{\it\alpha}^{\prime }{\it\beta}^{\prime }}$ and ${\it\Gamma}_{{\it\alpha}^{\prime }{\it\beta}^{\prime }}$ ) appearing in (3.30)–(3.33) can be expressed in terms of the following three integrals (see appendix A in supplementary materials available at http://dx.doi.org/10.1017/jfm.2016.237):

(3.35) $$\begin{eqnarray}\displaystyle \mathscr{H}_{{\it\alpha}{\it\beta}{\it\gamma}}^{{\it\delta}p}({\it\eta},R,{\it\phi},{\it\lambda}) & \equiv & \displaystyle \int _{{\it\theta}=0}^{2{\rm\pi}}\int _{{\it\varphi}=0}^{{\rm\pi}}\sin ^{{\it\alpha}}2{\it\theta}\cos ^{{\it\beta}}2{\it\theta}\sin ^{{\it\delta}}{\it\varphi}\cos ^{p}{\it\varphi}\nonumber\\ \displaystyle & & \displaystyle \times \,(1-{\it\eta}\sin ^{2}{\it\varphi}\cos 2{\it\theta}+{\it\lambda}^{2}(3\sin ^{2}{\it\varphi}-2))^{{\it\gamma}/2}\nonumber\\ \displaystyle & & \displaystyle \times \,\mathfrak{F}({\it\chi}[{\it\eta},R,{\it\phi},{\it\lambda};{\it\theta},{\it\varphi}])\,\text{d}{\it\varphi}\,\text{d}{\it\theta},\end{eqnarray}$$

(3.36) $$\begin{eqnarray}\displaystyle \mathscr{J}_{{\it\alpha}{\it\beta}{\it\gamma}}^{{\it\delta}p}({\it\eta},R,{\it\phi},{\it\lambda}) & \equiv & \displaystyle \int _{{\it\theta}=0}^{2{\rm\pi}}\int _{{\it\varphi}=0}^{{\rm\pi}}\sin ^{{\it\alpha}}2{\it\theta}\cos ^{{\it\beta}}2{\it\theta}\sin ^{{\it\delta}}{\it\varphi}\cos ^{p}{\it\varphi}\nonumber\\ \displaystyle & & \displaystyle \times \,(1-{\it\eta}\sin ^{2}{\it\varphi}\cos 2{\it\theta}+{\it\lambda}^{2}(3\sin ^{2}{\it\varphi}-2))^{{\it\gamma}/2}\nonumber\\ \displaystyle & & \displaystyle \times \,\mathfrak{G}({\it\chi}[{\it\eta},R,{\it\phi},{\it\lambda};{\it\theta},{\it\varphi}])\,\text{d}{\it\varphi}\,\text{d}{\it\theta},\end{eqnarray}$$

(3.37) $$\begin{eqnarray}\displaystyle \mathscr{K}_{{\it\alpha}{\it\beta}}^{{\it\delta}p}({\it\eta},R,{\it\phi},{\it\lambda}) & \equiv & \displaystyle \int _{{\it\theta}=0}^{2{\rm\pi}}\int _{{\it\varphi}=0}^{{\rm\pi}}\sin ^{{\it\alpha}}2{\it\theta}\cos ^{{\it\beta}}2{\it\theta}\sin ^{{\it\delta}}{\it\varphi}\cos ^{p}{\it\varphi} [\,\!(1-2{\it\lambda}^{2})(\sin (2{\it\phi}+2{\it\theta})\nonumber\\ \displaystyle & & \displaystyle -\cos {\it\varphi}\times \cos (2{\it\phi}+2{\it\theta}))+\sin ^{2}{\it\varphi}(3{\it\lambda}^{2}\sin (2{\it\phi}+2{\it\theta})-{\it\eta}\sin 2{\it\phi})]\nonumber\\ \displaystyle & & \displaystyle \times \,\mathfrak{G}({\it\chi}[{\it\eta},R,{\it\phi},{\it\lambda};{\it\theta},{\it\varphi}])\,\text{d}{\it\varphi}\,\text{d}{\it\theta},\end{eqnarray}$$

with $\mathfrak{F}({\it\chi})$ and $\mathfrak{G}({\it\chi})$ being given by (3.20) and (3.21) respectively. Of course, the integrals (3.35)–(3.37) can be evaluated numerically via any quadrature method.

In § 4, we outline an approximate method to evaluate the integrals (3.35)–(3.37) analytically via a power series expansion which helps to reduce the integro-algebraic equations (3.30)–(3.33) into a set of algebraic equations for four unknowns ${\it\eta}$ , $R$ , ${\it\phi}$ and ${\it\lambda}$ . More importantly, we derive ‘closed-form’ analytical solutions of the second-moment balance at the Burnett order (see § 4.1) and beyond (§ 4.2).

4.1 Exact solution at Burnett order

Retaining terms up to second order $O({\it\eta}^{m}{\it\lambda}^{n}R^{p}\sin ^{q}(2{\it\phi}),m+n+p+q\leqslant 2)$ in the resulting infinite series for each integral (3.35)–(3.37) and substituting them into (3.30)–(3.33), we obtain the following set of coupled nonlinear algebraic equations

(4.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\begin{array}{@{}l@{}}20\sqrt{{\rm\pi}}\{1+{\textstyle \frac{4}{5}}(1+e){\it\nu}g_{0}\}{\it\eta}R\cos 2{\it\phi}+128(1+e){\it\nu}g_{0}R^{2}\\ \quad -\,3(1-e^{2}){\it\nu}g_{0}(10+{\it\eta}^{2}+32R^{2}+8\sqrt{{\rm\pi}}{\it\eta}R\cos 2{\it\phi})=0\end{array}\\ \begin{array}{@{}l@{}}35\sqrt{{\rm\pi}}{\it\eta}R\cos 2{\it\phi}+(1+e){\it\nu}g_{0}\{32(1+3e)R^{2}-3(3-e)({\it\eta}^{2}+21{\it\lambda}^{2})\\ \quad -\,8\sqrt{{\rm\pi}}(4-3e){\it\eta}R\cos 2{\it\phi}\}=0\end{array}\\ 5\sqrt{{\rm\pi}}R\cos 2{\it\phi}-(1+e){\it\nu}g_{0}\{3(3-e){\it\eta}+2(1-3e)\sqrt{{\rm\pi}}R\cos 2{\it\phi}\}=0\\ 5({\it\eta}-\sin 2{\it\phi})+2(1+e)(1-3e){\it\nu}g_{0}\sin (2{\it\phi})=0.\end{array}\right\}\end{eqnarray}$$

These equations represent the second-moment balance equations at ‘Burnett order’ (Burnett Reference Burnett1935) since all terms up to the second order in the shear rate have been retained.

Equations (4.3) admit an exact solution

(4.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\it\eta}={\displaystyle \frac{\{5-2(1+e)(1-3e){\it\nu}g_{0}\}}{5}}\sin 2{\it\phi}\\ \begin{array}{@{}rcl@{}}{\it\lambda}^{2}\, & =\, & {\displaystyle \frac{10(1-e)}{21(3-e)}}+{\displaystyle \frac{[(7-3e)\{5-2(1+e)(1-3e){\it\nu}g_{0}\}-18(1+e)^{2}(3-e){\it\nu}g_{0}]}{525(3-e)}}\\ \, & \, & \times \,\{5-2(1+e)(1-3e){\it\nu}g_{0}\}\sin ^{2}2{\it\phi}\end{array}\\ R={\displaystyle \frac{3(1+e)(3-e)}{5\sqrt{{\rm\pi}}}}{\it\nu}g_{0}\tan 2{\it\phi}\\ {\displaystyle \frac{{\it\eta}}{R}}\cos (2{\it\phi})={\displaystyle \frac{\sqrt{{\rm\pi}}}{3(1+e)(3-e)}}\cos ^{2}(2{\it\phi})\left({\displaystyle \frac{5}{{\it\nu}g_{0}}}+2(1+e)(3e-1)\right),\end{array}\right\}\end{eqnarray}$$

where $\sin ^{2}(2{\it\phi})=\mathscr{Y}$ is the real positive root of the quadratic equation

(4.5) $$\begin{eqnarray}\displaystyle & \hspace{-12.0pt} & \displaystyle (11-3e)\{5-2(1+e)(1-3e){\it\nu}g_{0}\}^{2}{\rm\pi}\mathscr{Y}^{2}-\mathscr{Y}[(11-3e)\{5-2(1+e)(1-3e){\it\nu}g_{0}\}^{2}{\rm\pi}\nonumber\\ \displaystyle & \hspace{-12.0pt} & \displaystyle \quad +\,96(1+3e)(1+e)^{2}(3-e)^{2}{\it\nu}^{2}g_{0}^{2}+250{\rm\pi}(1-e)]+250{\rm\pi}(1-e)=0.\end{eqnarray}$$

For specified values of ${\it\nu}$ and $e$ , the non-coaxiality angle ${\it\phi}$ is determined from (4.5) and the remaining quantities are from (4.4): this is dubbed the ‘Burnett-order’ solution for ${\it\phi}$ , ${\it\eta}$ , ${\it\lambda}^{2}$ and $R$ as functions of ${\it\nu}$ and  $e$ .

4.2 Beyond Burnett order: perturbation solution

To obtain solutions beyond the Burnett order, i.e. at $O({\it\eta}^{m}{\it\lambda}^{n}R^{p}\sin ^{q}(2{\it\phi}),m+n+p+q>2)$ , we must solve related nonlinear algebraic equations as given by (B 1) in appendix B. We could not find an ‘exact’ solution of (B 1) either at super-Burnett or super–super Burnett order (except in the dilute limit, see appendix B.2). Therefore we look for perturbation solution of (B 1) by taking our Burnett-order solution (4.4)–(4.5) as the leading solution:

(4.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\it\eta}={\it\eta}^{(2)}+{\it\varepsilon}{\it\eta}^{(3)}+{\it\varepsilon}^{2}{\it\eta}^{(4)}+\cdots \\ {\it\lambda}={\it\lambda}^{(2)}+{\it\varepsilon}{\it\lambda}^{(3)}+{\it\varepsilon}^{2}{\it\lambda}^{(4)}+\cdots \\ R=R^{(2)}+{\it\varepsilon}R^{(3)}+{\it\varepsilon}^{2}R^{(4)}+\cdots \\ \sin 2{\it\phi}=\sin 2{\it\phi}^{(2)}+{\it\varepsilon}\sin 2{\it\phi}^{(3)}+{\it\varepsilon}^{2}\sin 2{\it\phi}^{(4)}+\cdots \,.\end{array}\right\}\end{eqnarray}$$

In the above expressions ${\it\varepsilon}\sim \dot{{\it\gamma}}$ and the superscript ‘2’ corresponds to the ‘Burnett-order’ solution and the superscripts ‘3’ and ‘4’ correspond to the corrections at the third and fourth order, respectively, in the shear rate. Plugging (4.6) into the corresponding third- and fourth-order equations and after performing some cumbersome algebra, we obtain the following solution for the correction terms at third order:

(4.7) $$\begin{eqnarray}{\it\eta}^{(3)}=0={\it\lambda}^{(3)}=R^{(3)}=\sin 2{\it\phi}^{(3)}.\end{eqnarray}$$

The fourth-order correction terms are found to be non-zero and are given in appendix B.1 (supplementary materials).

Figure 3. Variations of (a)  ${\it\eta}$ , (b)  ${\it\lambda}^{2}$ , (c)  $R$ and (d)  ${\it\phi}$ (degrees) with density ( ${\it\nu}$ ). While the solid black lines denote the exact numerical solution of second-moment equations, the blue dashed and the red dot-dashed lines denote the exact Burnett-order solution and the perturbation solution at fourth order, respectively.

4.3 Comparison between analytical and numerical solutions for ${\it\eta}$ , ${\it\phi}$ , ${\it\lambda}$ and $R$

Figure 3(a–d) shows a comparison of the ‘exact’ numerical solution of second-moment equations with its (i) exact Burnett-order solution and (ii) perturbation solution at the fourth order for the variations of the shear-plane anisotropy ${\it\eta}$ , the excess temperature ${\it\lambda}^{2}$ , the dimensionless shear rate $R$ and the non-coaxiality angle ${\it\phi}$ (degrees) with density ( ${\it\nu}$ ) for three values of the restitution coefficient $e=0.9$ , 0.7 and 0.5. In each panel, while the solid black lines denote the exact numerical solution, the blue dashed lines and red dot-dashed lines denote the analytical solutions at the second and fourth order, respectively. It is seen that while the Burnett-order solution provides a good agreement for ${\it\eta}$ , ${\it\lambda}^{2}$ , ${\it\phi}$ and $R$ up to a restitution coefficient of $e\geqslant 0.9$ , the fourth-order solution is required for more dissipative particles ( $e=0.5$ ) for a reasonable agreement over the whole range of density (for ${\it\lambda}^{2}$ , $R$ and ${\it\phi}$ , see panels $b$ , $c$ and $d$ respectively). On the other hand, panel $a$ indicates that even the fourth-order perturbation solution is not enough to quantitatively predict the shear-plane anisotropy ${\it\eta}$ at $e=0.5$ in the dense limit (for  ${\it\nu}>0.3$ ).

All transport coefficients in USF can be expressed as functions of ${\it\eta}$ , ${\it\lambda}^{2}$ , ${\it\phi}$ and $R$ as discussed in §§ 5.1–5.4. We will further check the ability of the Burnett- and fourth-order solutions of the second-moment equation to predict $p$ , ${\it\mu}$ , $\unicode[STIX]{x1D615}_{1}$ and $\unicode[STIX]{x1D615}_{2}$ in § 5.5 as functions of ${\it\nu}$ and  $e$ .

5.1 Shear stress and viscosity

Retaining all terms up to the fourth order $O({\it\eta}^{m}{\it\lambda}^{n}R^{p}\sin ^{q}(2{\it\phi}),m+n+p+q\leqslant 4)$ , the dimensionless shear stress can be written as (see appendix C)

(5.4) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D617}_{xy}^{\ast }}{{\it\nu}T^{\ast }} & = & \displaystyle -{\it\eta}\cos 2{\it\phi}-\frac{4(1+e){\it\nu}g_{0}}{105\sqrt{{\rm\pi}}}\left[21R\left\{8+\sqrt{{\rm\pi}}\frac{{\it\eta}\cos 2{\it\phi}}{R}\right\}+48{\it\lambda}^{2}R\right.\nonumber\\ \displaystyle & & \displaystyle +\left.4R^{3}\left\{32-\frac{{\it\eta}^{2}}{R^{2}}(2+(1+2\cos ^{2}2{\it\phi}))\right\}\right],\end{eqnarray}$$

with the dimensionless temperature $T^{\ast }$ being given by (5.3). The expression for the dimensionless shear viscosity, ${\it\mu}^{\ast }=-\unicode[STIX]{x1D617}_{xy}/{\it\rho}_{p}U_{R}^{2}=-\unicode[STIX]{x1D617}_{xy}^{\ast }$ , follows from (5.4):

(5.5) $$\begin{eqnarray}\displaystyle {\it\mu}^{\ast } & = & \displaystyle \frac{{\it\nu}\sqrt{T^{\ast }}}{8}\left[\frac{{\it\eta}\cos 2{\it\phi}}{R}+\frac{4(1+e){\it\nu}g_{0}}{105\sqrt{{\rm\pi}}}\left(21\left\{8+\sqrt{{\rm\pi}}\frac{{\it\eta}\cos 2{\it\phi}}{R}\right\}\right.\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\left.\underbrace{48{\it\lambda}^{2}+128R^{2}-4{\it\eta}^{2}\{2+(1+2\cos ^{2}2{\it\phi})\}}_{}\right)\right],\end{eqnarray}$$

where the under-braced terms represent nonlinear contributions beyond the NS order.

Neglecting quadratic- and higher-order terms in (5.5), we obtain the NS-order expression for the shear viscosity:

(5.6) $$\begin{eqnarray}{\it\mu}_{NS}^{\ast }=\frac{{\it\nu}\sqrt{T^{\ast }}}{8}\left[\frac{{\it\eta}\cos 2{\it\phi}}{R}+\frac{4(1+e){\it\nu}g_{0}}{5}\left(\frac{8}{\sqrt{{\rm\pi}}}+\frac{{\it\eta}\cos 2{\it\phi}}{R}\right)\right]+O(R^{2}).\end{eqnarray}$$

The elastic limit of the Burnett-order solution (4.4) for ${\it\eta}\cos 2{\it\phi}/R$ , with ${\it\phi}\rightarrow 0$ (which holds at NS order),

(5.7) $$\begin{eqnarray}\frac{{\it\eta}\cos 2{\it\phi}}{R}\stackrel{e=1}{\stackrel{{\it\phi}=0}{\longrightarrow }}\frac{5\sqrt{{\rm\pi}}}{12}\left(\frac{1}{{\it\nu}g_{0}}+\frac{8}{5}\right),\end{eqnarray}$$

can be substituted into (5.6) to arrive at

(5.8) $$\begin{eqnarray}{\it\mu}_{NS}^{\mathit{elastic}}=\sqrt{T^{\ast }}\left[\frac{5\sqrt{{\rm\pi}}}{96g_{0}}\left(1+\frac{8}{5}{\it\nu}g_{0}\right)^{2}+\frac{8}{5\sqrt{{\rm\pi}}}{\it\nu}^{2}g_{0}\right].\end{eqnarray}$$

This expression (5.8) matches exactly with the shear viscosity for an elastic hard-sphere system (Chapman & Cowling Reference Chapman and Cowling1970).

5.2 Normal stress components and the pressure

The diagonal components of the stress tensor, correct up to $O({\it\eta}^{m}{\it\lambda}^{n}R^{p}\sin ^{q}(2{\it\phi}),m+n+p+q\leqslant 4)$ , have the following expressions:

(5.9) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D617}_{xx}^{\ast }}{{\it\nu}T^{\ast }} & = & \displaystyle (1+{\it\lambda}^{2}+{\it\eta}\sin 2{\it\phi})+\frac{2(1+e){\it\nu}g_{0}}{1155}\left[\vphantom{\displaystyle \frac{8}{\sqrt{{\rm\pi}}}}33(35+96R^{2}+14{\it\eta}\sin 2{\it\phi}+14{\it\lambda}^{2})\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{8}{\sqrt{{\rm\pi}}}{\it\eta}R\cos 2{\it\phi}\{3(66+5{\it\eta}^{2}-22{\it\lambda}^{2})-160R^{2}-22{\it\eta}\sin 2{\it\phi}\}\right],\end{eqnarray}$$

(5.10) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D617}_{yy}^{\ast }}{{\it\nu}T^{\ast }} & = & \displaystyle (1+{\it\lambda}^{2}-{\it\eta}\sin 2{\it\phi})+\frac{2(1+e){\it\nu}g_{0}}{1155}\left[\vphantom{\displaystyle \frac{8}{\sqrt{{\rm\pi}}}}33(35+96R^{2}-14{\it\eta}\sin 2{\it\phi}+14{\it\lambda}^{2})\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{8}{\sqrt{{\rm\pi}}}{\it\eta}R\cos 2{\it\phi}\{3(66+5{\it\eta}^{2}-22{\it\lambda}^{2})-160R^{2}+22{\it\eta}\sin 2{\it\phi}\}\right],\end{eqnarray}$$

(5.11) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D617}_{zz}^{\ast }}{{\it\nu}T^{\ast }} & = & \displaystyle (1-2{\it\lambda}^{2})+\frac{2(1+e){\it\nu}g_{0}}{1155}\nonumber\\ \displaystyle & & \displaystyle \times \left[33(35+32R^{2}-28{\it\lambda}^{2})+\frac{8}{\sqrt{{\rm\pi}}}{\it\eta}R\cos 2{\it\phi}\{(66+3{\it\eta}^{2}-32R^{2})\}\right].\end{eqnarray}$$

The dimensionless mean pressure, correct up to $O({\it\eta}^{m}{\it\lambda}^{n}R^{p}\sin ^{q}(2{\it\phi}),m+n+p+q\leqslant 4)$ , is given by

(5.12) $$\begin{eqnarray}p^{\ast }={\it\nu}T^{\ast }\left[1+\frac{2(1+e){\it\nu}g_{0}}{315}\left\{315+\underbrace{672R^{2}+\frac{8}{\sqrt{{\rm\pi}}}{\it\eta}R\cos 2{\it\phi}(42+3{\it\eta}^{2}-32R^{2}-12{\it\lambda}^{2})}_{}\right\}\right]\!.\end{eqnarray}$$

Neglecting the ‘under-braced’ nonlinear terms in (5.12), we obtain the well-known expression for pressure,

(5.13) $$\begin{eqnarray}p_{NS}^{\ast }={\it\nu}T^{\ast }(1+2(1+e){\it\nu}g_{0}),\end{eqnarray}$$

at the NS order.

5.3 First and second normal stress differences

Subtracting (5.10) from (5.9) the expression for the first normal stress difference (third equation in (5.2)) is found to be

(5.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D615}_{1}^{\ast }=2{\it\eta}\sin (2{\it\phi})\left(1+\frac{4(1+e){\it\nu}g_{0}}{105}\left[21-\frac{8}{\sqrt{{\rm\pi}}}{\it\eta}R\cos (2{\it\phi})\right]\right){\it\nu}T^{\ast }. & & \displaystyle\end{eqnarray}$$

Similarly, the expression for the second normal stress difference (last equation in (5.2)) is

(5.15) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D615}_{2}^{\ast } & = & \displaystyle [3{\it\lambda}^{2}-{\it\eta}\sin (2{\it\phi})]{\it\nu}T^{\ast }+\frac{32(1+e){\it\nu}^{2}T^{\ast }g_{0}}{1155}\left[264\left(\frac{1}{2}+\frac{7}{{\it\nu}}\unicode[STIX]{x1D615}_{2}^{k\ast }\right)R^{2}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{1}{\sqrt{{\rm\pi}}}{\it\eta}R\cos 2{\it\phi}\{66+6{\it\eta}^{2}-64R^{2}-33{\it\lambda}^{2}+11{\it\eta}\sin 2{\it\phi}\}\right],\end{eqnarray}$$

where $T^{\ast }$ is the dimensionless temperature (5.3). Both (5.14) and (5.15) are correct at fourth order $O({\it\eta}^{m}{\it\lambda}^{n}R^{p}\sin ^{q}(2{\it\phi}),m+n+p+q\leqslant 4)$ .

5.4 Source of the second moment and the collisional dissipation

In USF the collisional source of the second moment (3.14) takes the following form:

(5.16) $$\begin{eqnarray}\boldsymbol{\aleph }=\left(\begin{array}{@{}ccc@{}}\aleph _{xx} & \aleph _{xy} & 0\\ \aleph _{yx} & \aleph _{yy} & 0\\ 0 & 0 & \aleph _{zz}\end{array}\right),\end{eqnarray}$$

with its non-zero components being given by

(5.17) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\aleph _{xx}=\unicode[STIX]{x1D608}_{xx}+\widehat{\unicode[STIX]{x1D60C}}_{xx}+\widehat{\unicode[STIX]{x1D60E}}_{xx}+2\dot{{\it\gamma}}{\it\Theta}_{xy}\\ \aleph _{yy}=\unicode[STIX]{x1D608}_{yy}+\widehat{\unicode[STIX]{x1D60C}}_{yy}+\widehat{\unicode[STIX]{x1D60E}}_{yy}-2\dot{{\it\gamma}}{\it\Theta}_{xy}\\ \aleph _{zz}=\unicode[STIX]{x1D608}_{zz}+\widehat{\unicode[STIX]{x1D60C}}_{zz}+\widehat{\unicode[STIX]{x1D60E}}_{zz}\\ \aleph _{xy}=\unicode[STIX]{x1D608}_{xy}+\widehat{\unicode[STIX]{x1D60C}}_{xy}+\widehat{\unicode[STIX]{x1D60E}}_{xy}+\dot{{\it\gamma}}({\it\Theta}_{yy}-{\it\Theta}_{xx}).\end{array}\right\}\end{eqnarray}$$

The integral expressions for ${\it\Theta}_{{\it\alpha}{\it\beta}}$ (3.15), $\unicode[STIX]{x1D608}_{{\it\alpha}{\it\beta}}$ (3.16), $E_{{\it\alpha}{\it\beta}}$ (3.17) and $G_{{\it\alpha}{\it\beta}}$ (3.18) have been evaluated (appendix A) in terms of ${\it\eta}$ , ${\it\lambda}$ and $R$ , correct up to $O({\it\eta}^{m}{\it\lambda}^{n}R^{p}\sin ^{q}(2{\it\phi}),m+n+p+q\leqslant 4)$ , and the resulting truncated series for each term in (5.17) is written down in appendix D (supplementary materials).

5.4.1 Collisional dissipation

The constitutive expression for the collisional dissipation rate (2.16) follows directly from the trace of (5.16):

(5.18) $$\begin{eqnarray}\displaystyle \mathscr{D} & = & \displaystyle -\frac{1}{2}\aleph _{{\it\beta}{\it\beta}}=-\frac{1}{2}(\unicode[STIX]{x1D608}_{xx}+\unicode[STIX]{x1D608}_{yy}+\unicode[STIX]{x1D608}_{zz})=\frac{3(1-e^{2}){\it\rho}{\it\nu}g_{0}T^{3/2}}{{\it\sigma}{\rm\pi}^{3/2}}\mathscr{H}_{003}^{10}({\it\eta},{\it\lambda}^{2},R,{\it\phi})\nonumber\\ \displaystyle & = & \displaystyle \frac{{\it\rho}{\it\nu}g_{0}(1-e^{2})T^{3/2}}{70{\it\sigma}\sqrt{{\rm\pi}}}\left[840+32\left\{\left(84+21\sqrt{{\rm\pi}}\frac{{\it\eta}}{R}\cos 2{\it\phi}\right)\right.\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\!\left.32R^{2}-2{\it\eta}^{2}(2+\cos 4{\it\phi})+24{\it\lambda}^{2}\vphantom{\left\{\!\left(\displaystyle \frac{{\it\eta}}{R}\right)\right\}}\right\}R^{2}+3(28{\it\eta}^{2}+{\it\eta}^{4}-8{\it\eta}^{2}{\it\lambda}^{2}+84{\it\lambda}^{4})\right],\qquad\end{eqnarray}$$

where we have retained terms up to $O({\it\eta}^{m}{\it\lambda}^{n}R^{p}\sin ^{q}(2{\it\phi}),m+n+p+q\leqslant 4)$ . In the isotropic limit ( ${\it\eta}\rightarrow 0$ and ${\it\lambda}^{2}\rightarrow 0$ ), we obtain

(5.19) $$\begin{eqnarray}\mathscr{D}=\mathscr{D}^{0}(1+{\textstyle \frac{16}{5}}R^{2}+{\textstyle \frac{128}{105}}R^{4}),\end{eqnarray}$$

where

(5.20) $$\begin{eqnarray}\mathscr{D}^{0}=\frac{12{\it\rho}_{p}{\it\nu}^{2}g_{0}(1-e^{2})T^{3/2}}{{\it\sigma}\sqrt{{\rm\pi}}}\end{eqnarray}$$

is the corresponding bare part valid at NS order (Jenkins & Richman Reference Jenkins and Richman1985). The quadratic-order shear-rate dependence in (5.18) agrees qualitatively with that calculated by Sela & Goldhirsch (Reference Sela and Goldhirsch1998) via a Burnett-order Chapman–Enskog expansion of the inelastic Boltzmann equation. The dependence of (5.18) on ${\it\eta}$ and ${\it\lambda}$ indicates that the Grad-level collisional dissipation depends on both normal stress differences,

(5.21) $$\begin{eqnarray}\mathscr{D}\equiv \mathscr{D}(\cdots \,;\unicode[STIX]{x1D615}_{1},\unicode[STIX]{x1D615}_{2}),\end{eqnarray}$$

since $({\it\eta},{\it\lambda})\sim (\unicode[STIX]{x1D615}_{1},\unicode[STIX]{x1D615}_{2})$ as we demonstrate below.

5.4.2 Dilute limit: dependence on NSDs

To clarify the functional dependence of (5.18) on normal stress differences, we focus on the dilute limit. In appendix B.2 (supplementary materials) we showed that the second-moment equations admit an exact solution in the dilute limit ( ${\it\nu}\rightarrow 0$ ):

(5.22) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\it\eta}^{2}={\textstyle \frac{1}{12}}(6+\unicode[STIX]{x1D615}_{1}^{k}+2\unicode[STIX]{x1D615}_{2}^{k})\unicode[STIX]{x1D615}_{1}^{k}\\ {\it\lambda}^{2}={\textstyle \frac{1}{6}}(\unicode[STIX]{x1D615}_{1}^{k}+2\unicode[STIX]{x1D615}_{2}^{k}),\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D615}_{1}^{k}=(\unicode[STIX]{x1D617}_{xx}^{k}-\unicode[STIX]{x1D617}_{yy}^{k})/p$ and $\unicode[STIX]{x1D615}_{2}^{k}=(\unicode[STIX]{x1D617}_{yy}^{k}-\unicode[STIX]{x1D617}_{zz}^{k})/p$ are the kinetic parts of the ‘scaled’ first and second normal stress differences, respectively. Substituting (5.22) into (5.18) and retaining terms up to quadratic order in $\unicode[STIX]{x1D615}_{1}^{k}$ and $\unicode[STIX]{x1D615}_{2}^{k}$ , we obtain the following expression for the collisional dissipation rate,

(5.23) $$\begin{eqnarray}\displaystyle \mathscr{D} & = & \displaystyle \frac{3{\it\rho}_{p}{\it\nu}^{2}(1-e^{2})T^{3/2}}{70{\it\sigma}\sqrt{{\rm\pi}}}[280+{\it\eta}^{2}(28+{\it\eta}^{2}-8{\it\lambda}^{2})+84{\it\lambda}^{4}]+\text{h.o.t.},\nonumber\\ \displaystyle & = & \displaystyle \frac{3{\it\rho}_{p}{\it\nu}^{2}(1-e^{2})T^{3/2}}{70{\it\sigma}\sqrt{{\rm\pi}}}\left[280+\frac{\unicode[STIX]{x1D615}_{1}^{k}}{2}\left(28+\frac{1}{2}\unicode[STIX]{x1D615}_{1}^{k}+\frac{10}{3}(\unicode[STIX]{x1D615}_{1}^{k}+2\unicode[STIX]{x1D615}_{2}^{k})\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{7}{3}(\unicode[STIX]{x1D615}_{1}^{k}+2\unicode[STIX]{x1D615}_{2}^{k})^{2}\vphantom{\displaystyle \frac{\unicode[STIX]{x1D615}_{1}^{k}}{2}}\right]+O((N_{i}^{k})^{3}),\end{eqnarray}$$

which holds in the dilute limit. That the Grad-level dissipation rate depends on the normal stress difference was pointed out previously (Saha & Alam Reference Saha and Alam2014) for the case of granular shear flow in two dimensions.

5.5 Validation of non-Newtonian constitutive relations

For the ‘exact’ numerical solution of (3.30)–(3.33), first we evaluate the integrals $\mathscr{H}_{{\it\alpha}{\it\beta}{\it\gamma}}^{{\it\delta}p}$ (3.35), $\mathscr{J}_{{\it\alpha}{\it\beta}{\it\gamma}}^{{\it\delta}p}$ (3.36), and $\mathscr{K}_{{\it\alpha}{\it\beta}}^{{\it\delta}p}$ (3.37), that appear in (A 1)–(A 2), numerically using the standard quadrature rule, for specified values of ${\it\nu}$ and  $e$ . Substituting the numerically evaluated integrals into (3.30)–(3.33) results in a system of nonlinear algebraic equations which again is solved by the same Newton’s method. The values of ${\it\eta}$ , $R$ , ${\it\phi}$ and ${\it\lambda}$ thus obtained are inserted into the expressions for pressure ( $p$ , (C 6)), viscosity ( ${\it\mu}$ , (C 8)) and the normal stress differences ( $\unicode[STIX]{x1D615}_{1}$ and $\unicode[STIX]{x1D615}_{2}$ , (C 9)) as given in appendix C. Such numerically obtained transport coefficients are dubbed ‘exact’ numerical solutions since a very high accurate solution can be obtained, limited only by the truncation error of the quadrature rule and the machine precision. In the following, such exact numerical solutions for $p$ , ${\it\mu}$ , $T$ , $\unicode[STIX]{x1D615}_{1}$ and $\unicode[STIX]{x1D615}_{2}$ are compared with those obtained from (i) the (exact) Burnett-order solution (§ 4.1) and (ii) the perturbation solution at fourth order (§ 4.2) for ${\it\eta}$ , $R$ , ${\it\phi}$ and  ${\it\lambda}$ .

Figure 4. Comparison between the ‘Burnett-order’ analytical solution (blue dashed lines), fourth-order perturbation solution (red dot-dashed lines) and the ‘exact’ numerical solution (black solid lines) for the variations of the (a) pressure, (b) shear viscosity and (c) granular temperature with volume fraction  $({\it\nu})$ . The filled circles denote the simulation data of Alam & Luding (Reference Alam, Luding, Garcia-Rojo, Herrmann and McNamara2005) for $e=0.9$ .

Figure 4(a–c) displays the density variations of (a) the pressure  $p$ , (b) the shear viscosity ${\it\mu}$ and (c) granular temperature $T$ for three values of the restitution coefficient 0.9, 0.7 and 0.5. In each panel, the ‘exact’ numerical solution (denoted by the black solid line) is compared with (i) Burnett order (blue dashed line) and (ii) the perturbation solution at fourth order (red dot-dashed line). It is seen that the Burnett-order solutions for $p$ , ${\it\mu}$ and $T$ are almost indistinguishable from their exact numerical value at small dissipation ( $e=0.9$ ); moreover, this agreement seems to hold uniformly for the whole range of density. On the other hand, retaining the fourth-order terms yields a better agreement for $p$ , ${\it\mu}$ and $T$ at large dissipation ( $e=0.5$ ).

The ability of the fourth-order series solution to quantitatively predict $p$ and ${\it\mu}$ at any density also holds for both the first and second normal stress differences, see figure 5(a,b). Note that the plotted quantities in figure 5 are the ‘scaled’ first and second normal stress differences defined via

(5.24) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D615}_{1}=\frac{\unicode[STIX]{x1D617}_{xx}-\unicode[STIX]{x1D617}_{yy}}{p} & \displaystyle\end{eqnarray}$$

(5.25) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D615}_{2}=\frac{\unicode[STIX]{x1D617}_{yy}-\unicode[STIX]{x1D617}_{zz}}{p} & \displaystyle\end{eqnarray}$$

respectively, with the expressions for $\unicode[STIX]{x1D617}_{xx}$ , $\unicode[STIX]{x1D617}_{yy}$ , $\unicode[STIX]{x1D617}_{zz}$ and $p$ given in (5.9)–(5.12); equations (5.24) and (5.25) are measures of two normal stress differences with respect to the mean pressure. In figure 5(b) we find that $\unicode[STIX]{x1D615}_{2}$ undergoes a sign reversal at some finite density. The location ( ${\it\nu}={\it\nu}_{cr}$ ) of the sign reversal of $\unicode[STIX]{x1D615}_{2}$ appears to be independent of the restitution coefficient, see figure 5(b).

The event-driven simulation data for the uniform shear flow of inelastic hard spheres, previously carried out by Alam & Luding (Reference Alam, Luding, Garcia-Rojo, Herrmann and McNamara2005) for a restitution coefficient of $e=0.9$ , are superimposed in figures 4 and 1. Note that Alam & Luding (Reference Alam, Luding, Garcia-Rojo, Herrmann and McNamara2005) used event-driven techniques to conduct these simulations of smooth inelastic hard spheres in a cubic box by implementing the Lees–Edwards boundary condition (Lees & Edwards Reference Lees and Edwards1972) along the gradient ( $y$ ) direction with periodic boundary conditions along the streamwise ( $x$ ) and spanwise ( $z$ ) directions: the other details of simulation can be obtained from the original paper. Figure 4(a–c) indicates that our theoretical predictions are in good agreement with the simulation data for $p$ , ${\it\mu}$ and $T$ as well as for two normal stress differences $\unicode[STIX]{x1D615}_{1}$ and $\unicode[STIX]{x1D615}_{2}$ (figure 1) over a large range of density up to the freezing density ( ${\it\nu}\sim 0.5$ ).

We may conclude from the comparative study in figures 3–5 that the terms retained up to the fourth order ( $\text{super}^{2}$ -Burnett solutions) in the series expansion (4.1)–(4.2) (that appear in various collision integrals (3.35)–(3.37) in the second-moment balance equations (3.30)–(3.33)) provide an adequate accuracy to predict all transport coefficients ( $p$ , ${\it\mu}$ , $\unicode[STIX]{x1D615}_{1}$ and $\unicode[STIX]{x1D615}_{2}$ ) in a sheared granular fluid over a large range of ${\it\nu}$ and  $e$ .

Figure 5. Same as figure 4, but for the variations of (a) the scaled first $(\unicode[STIX]{x1D615}_{1})$ and (b) second $(\unicode[STIX]{x1D615}_{2})$ normal stress differences with volume fraction  $({\it\nu})$ .

6.1 The Boltzmann/dilute limit

To see whether the anisotropy of $\unicode[STIX]{x1D648}$ can be tied with the non-Newtonian nature of the stress tensor, we recall that the exact expression for the stress tensor (2.8) is $\unicode[STIX]{x1D64B}\equiv {\it\rho}\unicode[STIX]{x1D648}$ in the dilute limit. The first normal stress difference is given by

(6.2) $$\begin{eqnarray}(\unicode[STIX]{x1D617}_{xx}-\unicode[STIX]{x1D617}_{yy})={\it\rho}(\unicode[STIX]{x1D614}_{xx}-\unicode[STIX]{x1D614}_{yy})\sim {\it\eta}\sin 2{\it\phi}\geqslant 0,\end{eqnarray}$$

while the second normal stress difference is

(6.3) $$\begin{eqnarray}(\unicode[STIX]{x1D617}_{yy}-\unicode[STIX]{x1D617}_{zz})={\it\rho}(\unicode[STIX]{x1D614}_{yy}-\unicode[STIX]{x1D614}_{zz})\sim 3{\it\lambda}^{2}-{\it\eta}\sin 2{\it\phi}\leqslant 0\end{eqnarray}$$

and both are non-zero, and hence the stress tensor is non-Newtonian if ${\it\phi}\neq 0$ and/or ${\it\eta}\neq 0$ and ${\it\lambda}\neq 0$ . Therefore, the anisotropy of $\unicode[STIX]{x1D648}$ gives rise to finite normal stress differences, resulting in the stress tensor being non-Newtonian.

It is clear from (6.2)–(6.3) that the normal stress differences and the anisotropy of the second-moment tensor are intertwined with each other. If the eigendirections of $\unicode[STIX]{x1D648}$ are coaxial with those of the shear tensor $\unicode[STIX]{x1D63F}$ , i.e.  ${\it\phi}=0$ , the first normal stress difference (6.2) vanishes; however, the second normal stress difference (6.3) can still be non-zero ( $=3{\it\lambda}^{2}$ ) even if ${\it\phi}=0$ since the ‘excess’ temperature (6.1) along the vorticity direction could differ from zero in USF (see figure 6). Therefore, the origin of $\unicode[STIX]{x1D615}_{1}$ is tied to the ‘non-coaxiality’ between the eigendirections of $\unicode[STIX]{x1D648}$ and $\unicode[STIX]{x1D63F}$ , but that of $\unicode[STIX]{x1D615}_{2}$ is the non-zero ‘excess’ temperature along the vorticity direction.

6.2 Effect of density

In this section we discuss the effects of finite density on the origins of two NSDs.

6.2.1 First normal stress difference

The kinetic and collisional contributions ( $\unicode[STIX]{x1D615}_{1}^{\ast }=\unicode[STIX]{x1D615}_{1}^{k\ast }+\unicode[STIX]{x1D615}_{1}^{c\ast }$ ) to the first normal stress difference (5.14) are given by

(6.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D615}_{1}^{k\ast }=2{\it\eta}\sin (2{\it\phi}){\it\nu}T^{\ast }, & \displaystyle\end{eqnarray}$$

(6.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D615}_{1}^{c\ast }=\frac{8(1+e){\it\nu}g_{0}}{105}\left[21-\frac{8}{\sqrt{{\rm\pi}}}{\it\eta}R\cos (2{\it\phi})\right]{\it\eta}\sin (2{\it\phi}){\it\nu}T^{\ast }. & \displaystyle\end{eqnarray}$$

Note that both (6.4) and (6.5) vanish in the limits of ${\it\eta}\rightarrow 0$ and/or ${\it\phi}\rightarrow 0$ : while the former represents the limit of vanishing temperature anisotropy in the shear plane (3.10), the latter corresponds to the eigendirections of the second-moment tensor $\unicode[STIX]{x1D648}$ and the shear tensor $\unicode[STIX]{x1D63F}$ being coaxial (viz. figure 2). Therefore we conclude that the origin of first normal stress difference is tied to (i) the ‘finite’ temperature anisotropy and/or (ii) the ‘non-coaxiality’ between the eigendirections of $\unicode[STIX]{x1D648}$ and $\unicode[STIX]{x1D63F}$ at any density.

The leading terms in both (6.4) and (6.5) are

(6.6) $$\begin{eqnarray}{\it\eta}\sin 2{\it\phi}=O(\dot{{\it\gamma}}^{2}),\end{eqnarray}$$

of Burnett order in the shear rate. The leading-order corrections in (6.5) are

(6.7) $$\begin{eqnarray}R^{2}{\it\eta}\sin (2{\it\phi})\left(\frac{{\it\eta}}{R}\cos (2{\it\phi})\right)=O(\dot{{\it\gamma}}^{4}).\end{eqnarray}$$

It is noteworthy that the excess temperature (6.1), $T_{z}^{ex}\propto {\it\lambda}^{2}$ , does not affect the kinetic part of the first NSD, but it affects the collisional part of the first NSD at sixth order and beyond in the shear rate.

6.2.2 Second normal stress difference and its sign reversal

The kinetic and collisional components of the second normal stress difference (5.15) at $O({\it\eta}^{m}{\it\lambda}^{n}R^{p}\sin ^{q}(2{\it\phi}),m+n+p+q\leqslant 4)$ are given by

(6.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D615}_{2}^{k\ast }=[3{\it\lambda}^{2}-{\it\eta}\sin (2{\it\phi})]{\it\nu}T^{\ast }, & \displaystyle\end{eqnarray}$$

(6.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D615}_{2}^{c\ast } & = & \displaystyle \frac{32(1+e){\it\nu}^{2}T^{\ast }g_{0}}{1155}\left[264\left(\frac{1}{2}+\frac{7}{{\it\nu}}\unicode[STIX]{x1D615}_{2}^{k\ast }\right)R^{2}+\frac{1}{\sqrt{{\rm\pi}}}{\it\eta}R\cos 2{\it\phi}\right.\nonumber\\ \displaystyle & & \displaystyle \times \left.(66+6{\it\eta}^{2}-64R^{2}-33{\it\lambda}^{2}+11{\it\eta}\sin 2{\it\phi})\vphantom{\left(\displaystyle \frac{1}{\sqrt{{\rm\pi}}}\right)}\right],\end{eqnarray}$$

where $T^{\ast }$ is the dimensionless temperature (5.3). In the limit of vanishing of the ‘shear-plane’ temperature anisotropy ( ${\it\eta}\rightarrow 0$ ) and/or the coaxiality ( ${\it\phi}\rightarrow 0$ ) between the eigendirections of $\unicode[STIX]{x1D648}$ and $\unicode[STIX]{x1D63F}$ , we can simplify (6.8) and (6.9) into

(6.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D615}_{2}^{k\ast }=3{\it\lambda}^{2}{\it\nu}T^{\ast }\propto T_{z}^{ex}\geqslant 0, & \displaystyle\end{eqnarray}$$

(6.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D615}_{2}^{c\ast }=\frac{4(1+e){\it\nu}^{2}T^{\ast }g_{0}}{35}[32R^{2}+21{\it\lambda}^{2}]\geqslant 0. & \displaystyle\end{eqnarray}$$

Both (6.10) and (6.11) hold strictly in the dense limit since ${\it\eta}$ and ${\it\phi}$ approach zero as ${\it\nu}\rightarrow {\it\nu}_{\mathit{max}}$ . Note that even the kinetic part of the second normal stress difference (6.10) remains positive in the dense limit as long as ${\it\lambda}^{2}\propto T_{z}^{ex}>0$ (see figure 6).

It is evident from (6.10) and (6.11) that the second normal stress difference in the dense limit (6.11) remains finite and positive, which is in contrast to the zero first normal stress difference in the same limit. Moreover, even if ${\it\lambda}^{2}=0$ and ${\it\eta}=0$ , $\unicode[STIX]{x1D615}_{2}\equiv \unicode[STIX]{x1D615}_{2}^{c\ast }$ remains finite as ${\it\nu}\rightarrow {\it\nu}_{\mathit{max}}$ as long as the shear rate is finite ( $R^{2}>0$ ) and hence this is ‘shear induced’. Therefore, the origin of non-zero second normal stress difference in the dense limit is tied to the imposed shear field.

Equation (6.8) indicates that the $\unicode[STIX]{x1D615}_{2}^{k\ast }$ can be positive or negative depending on the relative magnitudes of $3{\it\lambda}^{2}$ and ${\it\eta}\sin (2{\it\phi})$ and can undergo a sign change at a finite density if

(6.12) $$\begin{eqnarray}3{\it\lambda}^{2}-{\it\eta}\sin (2{\it\phi})=0.\end{eqnarray}$$

(Recall from figure 5 b that $\unicode[STIX]{x1D615}_{2}$ undergoes a sign reversal at a finite density.) The density variation of (6.12), obtained by solving (3.30)–(3.33) numerically, as depicted in figure 7(a), clarifies the above point: $3{\it\lambda}^{2}<{\it\eta}\sin (2{\it\phi})$ in the dilute limit and $3{\it\lambda}^{2}>{\it\eta}\sin (2{\it\phi})$ in the dense limit for any value of the restitution coefficient. The variation of the critical density ${\it\nu}_{cr}^{k}$ , at which (6.12) holds, with the restitution coefficient is shown in figure 7(b) – clearly, ${\it\nu}_{cr}^{k}(e)$ is a decreasing function of $e$ and can be fitted via the following linear function:

(6.13) $$\begin{eqnarray}{\it\nu}_{cr}^{k}(e)=0.27-0.086e,\end{eqnarray}$$

denoted by the solid line in figure 7(b). Since $\unicode[STIX]{x1D615}_{2}^{c}=0$ at ${\it\nu}=0$ and is a monotonically increasing function of ${\it\nu}$ , the critical density, ${\it\nu}_{cr}$ , at which total second normal stress difference ( $\unicode[STIX]{x1D615}_{2}=\unicode[STIX]{x1D615}_{2}^{k}+\unicode[STIX]{x1D615}_{2}^{c}=0$ ) changes sign from negative to positive would be slightly lower than (6.13).

In conclusion, the second normal stress difference is negative and positive in the dilute and dense limits, respectively, and the ‘sign reversal’ of $\unicode[STIX]{x1D615}_{2}$ at some finite density is directly tied to the sign reversal of its kinetic component $\unicode[STIX]{x1D615}_{2}^{k}$ . The above analysis further confirms that the origin of $\unicode[STIX]{x1D615}_{2}$ is tied to the ‘excess’ temperature ( $T_{z}^{ex}\propto {\it\lambda}^{2}$ , viz. (6.1)) along the mean vorticity direction in the dilute limit, but its origin in the dense limit is tied to the imposed shear field.

Figure 7. (a) Variation of $(3{\it\lambda}^{2}-{\it\eta}\sin 2{\it\phi})$ with density ${\it\nu}$ for two values of restitution coefficient $e=0.9$ (thick solid line) and $e=0.1$ (dashed line); the arrows locate the respective critical density at which the above quantity is zero. (b) Variation of the critical density ${\it\nu}_{cr}^{k}$ , at which ${\it\eta}\sin 2{\it\phi}=3{\it\lambda}^{2}$ , with  $e$ ; the circles denote numerical solution and the solid line is a linear fit (6.13).