2 Thermomechanics of granular suspensions

We consider a two-phase mixture consisting of a granular material and a carrier fluid. The assumptions upon which our mixture-theoretic formalism is based are the following.

1. (i) Each phase is treated as a continuum thermodynamic system.

2. (ii) The two thermodynamic continua are immiscible but occupy the same space. More specifically, they fill completely the space that they occupy (saturation condition).

3. (iii) The two thermodynamic continua are open to each other and at non-equilibrium. As such, they interact with each other in the form of momentum and energy exchanges. Herein, we assume absence of chemical reactions or phase transformations.

4. (iv) The momentum and energy exchanges between the two continua are pure, i.e. their sum vanishes.

5. (v) The postulate of phase separation holds. In other words, dissipative phenomena associated with only one phase depend exclusively on the variables of this phase.

6. (vi) Both phases are macroscopically isotropic.

Throughout this article, the indices $f$ and $s$ are reserved to denote quantities attributed to the fluid and granular phase, respectively. $\unicode[STIX]{x1D719}_{f}$ and $\unicode[STIX]{x1D719}_{s}$ stand for the volume fractions of the two phases. Then, the saturation condition of assumption (ii) reads

(2.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{f}+\unicode[STIX]{x1D719}_{s}=1. & & \displaystyle\end{eqnarray}$$

The local form of the mass, momentum and energy balance laws for each phase are derived in the standard manner, i.e. by first writing the integral balance laws over an arbitrary control volume and then by applying successively the Reynolds transport and divergence theorems. As shown in Drew & Passman (Reference Drew and Passman1999), this procedure yields the following equations,

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FC}})}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FC}}\boldsymbol{u}_{\unicode[STIX]{x1D6FC}})=0, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FC}}\boldsymbol{u}_{\unicode[STIX]{x1D6FC}})}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FC}}\boldsymbol{u}_{\unicode[STIX]{x1D6FC}}\otimes \boldsymbol{u}_{\unicode[STIX]{x1D6FC}})+\unicode[STIX]{x1D735}(p_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FC}})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D6FC}})+\boldsymbol{f}_{\unicode[STIX]{x1D6FC}}, & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FC}}e_{\unicode[STIX]{x1D6FC}}^{t})}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FC}}e_{\unicode[STIX]{x1D6FC}}^{t}\boldsymbol{u}_{\unicode[STIX]{x1D6FC}})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(p_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FC}}\boldsymbol{u}_{\unicode[STIX]{x1D6FC}})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D6FC}}\boldsymbol{\cdot }\boldsymbol{u}_{\unicode[STIX]{x1D6FC}})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FC}}\boldsymbol{q}_{\unicode[STIX]{x1D6FC}})+{\mathcal{E}}_{\unicode[STIX]{x1D6FC}}. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

In the above equations $\unicode[STIX]{x1D6FC}=f,s$ . Using standard notation, $\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}$ , $\boldsymbol{u}_{\unicode[STIX]{x1D6FC}}$ , $p_{\unicode[STIX]{x1D6FC}}$ and $e_{\unicode[STIX]{x1D6FC}}^{t}$ , denote the density, velocity vector, thermodynamic pressure and total energy, respectively, of phase $\unicode[STIX]{x1D6FC}$ . Further, $\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D6FC}}$ stands for the (pressureless) stress tensor and $\boldsymbol{q}_{\unicode[STIX]{x1D6FC}}$ for the heat-flux vector of phase $\unicode[STIX]{x1D6FC}$ . Finally, $\boldsymbol{f}_{\unicode[STIX]{x1D6FC}}$ and ${\mathcal{E}}_{\unicode[STIX]{x1D6FC}}$ represent the interphasial momentum and energy exchanges, respectively.

In view of assumption (iii) we set,

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}\triangleq \,\boldsymbol{f}_{f}=-\boldsymbol{f}_{s}, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{E}}\triangleq {\mathcal{E}}_{f}=-{\mathcal{E}}_{s}. & \displaystyle\end{eqnarray}$$

It is also convenient to decompose the interphasial energy exchange into a term related to the momentum exchange plus a remaining part. More specifically, we set

(2.7) $$\begin{eqnarray}{\mathcal{E}}=E+\boldsymbol{f}\boldsymbol{\cdot }\boldsymbol{u}_{s}.\end{eqnarray}$$

Next, we express the total energy of each phase as the sum of internal and kinetic energies, $e_{\unicode[STIX]{x1D6FC}}^{t}=e_{\unicode[STIX]{x1D6FC}}+(1/2)\boldsymbol{u}_{\unicode[STIX]{x1D6FC}}^{2}$ . Then, the above balance laws for the fluid phase are written as,

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}(\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f})}{\text{d}t_{f}}+\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f}=0, & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}\frac{\text{d}\boldsymbol{u}_{f}}{\text{d}_{f}t}+\unicode[STIX]{x1D735}(p_{f}\unicode[STIX]{x1D719}_{f})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D64B}_{f})-\boldsymbol{f}, & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}\frac{\text{d}e_{f}}{\text{d}_{f}t}+p_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f}=\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D64B}_{f}\boldsymbol{ : }\unicode[STIX]{x1D651}_{f}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}_{f}\boldsymbol{q}_{f})-E-\boldsymbol{f}\boldsymbol{\cdot }(\boldsymbol{u}_{s}-\boldsymbol{u}_{f}), & \displaystyle\end{eqnarray}$$

whereas the balance laws for the granular phase can be written as,

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}(\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s})}{\text{d}_{s}t}+\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s}=0, & \displaystyle\end{eqnarray}$$

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\frac{\text{d}\boldsymbol{u}_{s}}{\text{d}_{s}t}+\unicode[STIX]{x1D735}(p_{s}\unicode[STIX]{x1D719}_{s})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D64B}_{s})+\boldsymbol{f}, & \displaystyle\end{eqnarray}$$

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\frac{\text{d}e_{s}}{\text{d}_{s}t}+p_{s}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s}=\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D64B}_{s}\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}_{s}\boldsymbol{q}_{s})+E. & \displaystyle\end{eqnarray}$$

In the above equations, $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D6FC}}$ denotes the strain-rate tensor, $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D6FC}}=(\unicode[STIX]{x1D735}\boldsymbol{u}_{\unicode[STIX]{x1D6FC}}+\unicode[STIX]{x1D735}\boldsymbol{u}_{\unicode[STIX]{x1D6FC}}^{\text{T}})/2$ . Also, the double-dot symbol, $(\boldsymbol{ : })$ , represents the double interior product between two second-order tensors. Finally, the material time derivative with respect to the velocity of phase $\unicode[STIX]{x1D6FC}$ is defined as,

(2.14) $$\begin{eqnarray}\displaystyle \frac{\text{d}(\cdot )}{\text{d}_{\unicode[STIX]{x1D6FC}}t}={\displaystyle \frac{\unicode[STIX]{x2202}(\cdot )}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}_{\unicode[STIX]{x1D6FC}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\cdot ). & & \displaystyle\end{eqnarray}$$

One might argue that the energy equations in (2.4) are symmetric whereas equations (2.10) and (2.13) are not. However, this is merely an apparent asymmetry and is due to the particular decomposition of ${\mathcal{E}}$ that was introduced in (2.7). More specifically, it is straightforward to show that all decompositions of the form ${\mathcal{E}}=E+\boldsymbol{f}\boldsymbol{\cdot }(c_{1}\boldsymbol{u}_{s}+c_{2}\boldsymbol{u}_{f})$ , with $c_{1}$ and $c_{2}$ constants, result in the same constitutive relations. In other words, the particular choice of the constants $c_{1}$ and $c_{2}$ is inconsequential to the final results.

The thermodynamic state of each phase is specified by its state vector and thus its Gibbs relation. For the fluid phase, we assume that it consists of a simple, homogeneous fluid. As such, its state vector is $(e_{f},\unicode[STIX]{x1D70C}_{f})^{\text{T}}$ and the corresponding Gibbs relations reads,

(2.15) $$\begin{eqnarray}T_{f}\,\text{d}s_{f}=\text{d}e_{f}-{\displaystyle \frac{p_{f}}{{\unicode[STIX]{x1D70C}_{f}}^{2}}}\text{d}\unicode[STIX]{x1D70C}_{f},\end{eqnarray}$$

where $s_{f}$ stands for the entropy of the fluid.

As regards the granular phase, we observe that it can perform additional reversible work due to its microstructure and this work will depend on the volume that grains occupy and on the properties of the fluid–granular conglomerate interface. For this reason, an adequate description of the state of the granular material requires additional field variables. The exact number of such variables remains an open question but a consensus in literature is being built that the standard state vector of simple fluids, here $(e_{s},\unicode[STIX]{x1D70C}_{s})^{\text{T}}$ , should be augmented by at least two variables: the volume fraction, $\unicode[STIX]{x1D719}_{s}$ , and its gradient, $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}$ , (Goodman & Cowin Reference Goodman and Cowin1972; Savage Reference Savage1979; Passman et al. Reference Passman, Nunziato, Bailey and Reed1986; Wang & Hutter Reference Wang and Hutter1999a ,Reference Wang and Hutter b ,Reference Wang and Hutter c ; Massoudi & Mehrabadi Reference Massoudi and Mehrabadi2001; Kirchner Reference Kirchner2002; Massoudi Reference Massoudi, Antony, Hoyle and Ding2004; Ván Reference Ván2004; Papalexandris Reference Papalexandris2004b ; Fang et al. Reference Fang, Wang and Hutter2006). Then, the usual Gibbs relation is modified to accommodate the reversible work of the additional variables as follows,

(2.16) $$\begin{eqnarray}T_{s}\,\text{d}s_{s}=\text{d}e_{s}-{\displaystyle \frac{p_{s}}{{\unicode[STIX]{x1D70C}_{s}}^{2}}}\text{d}\unicode[STIX]{x1D70C}_{s}-{\displaystyle \frac{\unicode[STIX]{x1D6FD}_{s}}{\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}}}\text{d}\unicode[STIX]{x1D719}_{s}-{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}}}\boldsymbol{h}_{s}\boldsymbol{\cdot }\text{d}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}).\end{eqnarray}$$

In this equation $s_{s}$ stands for the entropy of the granular phase. The quantities $\unicode[STIX]{x1D6FD}_{s}$ and $\boldsymbol{h}_{s}$ are defined as the partial derivatives of $s_{s}$ with respect to $\unicode[STIX]{x1D719}_{s}$ and $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}$ , respectively. They can equivalently be defined with respect to the Helmholtz free energy $\unicode[STIX]{x1D713}_{s}$ of the solid phase, so that

(2.17a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{s}=\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\left({\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{s}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{s}}}\right)_{T_{s},\unicode[STIX]{x1D70C}_{s},\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}},\quad \boldsymbol{h}_{s}=\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\left({\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{s}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}}}\right)_{T_{s},\unicode[STIX]{x1D70C}_{s},\unicode[STIX]{x1D719}_{s}}. & & \displaystyle\end{eqnarray}$$

According to the Gibbs relation (2.16), $\unicode[STIX]{x1D6FD}_{s}$ is the force conjugate to the variable $\unicode[STIX]{x1D719}_{s}$ to $\unicode[STIX]{x1D719}_{s}$ . It is often referred to as intergranular stress or configuration pressure; see, for example, Goodman & Cowin (Reference Goodman and Cowin1972), Josserand, Lagrée & Lhuillier (Reference Josserand, Lagrée and Lhuillier2004) and Josserand, Lagrée & Lhuillier (Reference Josserand, Lagrée and Lhuillier2006) for dry granular materials, as well as Passman et al. (Reference Passman, Nunziato, Bailey and Reed1986), Wang & Hutter (Reference Wang and Hutter1999a ,Reference Wang and Hutter b ) and Cochran & Powers (Reference Cochran and Powers2008) for fluid-saturated ones. Physically, it describes the continuum manifestation of contact stresses resulting from strains due to intergranular interactions and, in particular, the resistance of grains to compaction. Herein, these contact stresses are modelled as purely elastic ones and thus the third term on the right-hand side of the Gibbs relation (2.16) accounts for the total work of these stresses; a more elaborate discussion can be found in Bdzil et al. (Reference Bdzil, Menikoff, Son, Kapila and Stewart1999). Quasi-static compaction experiments have shown that $\unicode[STIX]{x1D6FD}_{s}\geqslant 0$ (Elban & Chiarito Reference Elban and Chiarito1986; Mesri & Vardhanabhuti Reference Mesri and Vardhanabhuti2009). In other words, the intergranular stress is of purely repulsive nature.

As regards $\boldsymbol{h}_{s}$ , we first observe that the presumed macroscopic isotropy of the granular material implies that its free energy is an isotropic function so that $\boldsymbol{h}_{s}$ can be written as

(2.18a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{h}_{s}=\unicode[STIX]{x1D6FE}_{s}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s},\quad \unicode[STIX]{x1D6FE}_{s}=2\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\left({\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{s}}{\unicode[STIX]{x2202}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}|^{2}}}\right)_{T_{s},\unicode[STIX]{x1D70C}_{s},\unicode[STIX]{x1D719}_{s}}. & & \displaystyle\end{eqnarray}$$

In this manner, the term $\unicode[STIX]{x1D6FE}_{s}$ is the quantity conjugate to $|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}|^{2}$ . This quantity is related to stresses that are developed due to rearrangements of the interfacial area density. At this point, a couple of comments regarding the introduction of $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}$ as an additional variable are due. First we note that, as shown in Drew & Passman (Reference Drew and Passman1999), $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}$ is a measure of the interfacial area density and, as such, relates directly to the microstructure of the system. Also, as demonstrated by Dunn & Serrin (Reference Dunn and Serrin1985), the description of stresses related to variations of the interfacial area density within thermodynamically consistent theories requires that the free energy of the material is a function of the volume fraction gradient. As regards the physical meaning of these stresses, we note that macroscopically the interface between a granular material and the interstitial fluid is not a sharp discontinuity; instead it has a non-zero thickness, thus constituting a thin transition layer. Within this layer, the steep volume fraction gradients result in Korteweg stress, i.e. a weak interfacial tension that resembles surface tension, Varsakelis & Papalexandris (Reference Varsakelis and Papalexandris2016).

Finally, with respect to the Gibbs relation for the granular phase, we remark that this work is concerned with granular suspensions that are by definition non-colloidal, i.e. the Péclet numbers are high. For this reason, the Gibbs relation (2.16) does not contain any term related to Brownian effects.

Next, the mass-averaged density $\unicode[STIX]{x1D70C}$ , velocity $\boldsymbol{u}$ and entropy $s$ of the mixture are defined as follows,

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}+\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}, & \displaystyle\end{eqnarray}$$

(2.20) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}=\frac{\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\boldsymbol{u}_{s}+\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}\boldsymbol{u}_{f}}{\unicode[STIX]{x1D70C}}, & \displaystyle\end{eqnarray}$$

(2.21) $$\begin{eqnarray}\displaystyle & \displaystyle s=\frac{s_{s}\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}+s_{f}\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}}{\unicode[STIX]{x1D70C}}. & \displaystyle\end{eqnarray}$$

As is the case for any observable field, the mixture entropy obeys a balance law of the form

(2.22) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}\frac{\text{d}s}{\text{d}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{J}^{s}=\unicode[STIX]{x1D70E}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{J}^{s}$ represents the entropy flux and $\unicode[STIX]{x1D70E}$ the entropy production rate. In this relation the material derivative of the entropy is taken with respect to the mass-averaged velocity, that is,

(2.23) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\text{d}(\cdot )}{\text{d}t}}={\displaystyle \frac{\unicode[STIX]{x2202}(\cdot )}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\cdot ). & & \displaystyle\end{eqnarray}$$

Next, we note that, according to the second thermodynamic axiom, the entropy production rate is non-negative,

(2.24) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}\geqslant 0. & & \displaystyle\end{eqnarray}$$

Upon combination of equations (2.8) and (2.11) with the definitions of the mass-averaged quantities (2.19), (2.20) and (2.21), the following expression for the material derivative of the mixture entropy $s$ is obtained.

(2.25) $$\begin{eqnarray}{\displaystyle \frac{\text{d}s}{\text{d}t}}=\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}s_{f}(\boldsymbol{u}_{f}-\boldsymbol{u})+\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}s_{s}(\boldsymbol{u}_{s}-\boldsymbol{u})+\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}\frac{\text{d}s_{f}}{\text{d}t_{f}}+\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\frac{\text{d}s_{s}}{\text{d}t_{s}}.\end{eqnarray}$$

3 Derivation of the constitutive relations

In this section we derive constitutive relations for the dissipative and relaxation phenomena that occur in granular suspensions by exploiting the constraints imposed by the second thermodynamic axiom (2.24).

Our starting point is the rate equation for the mixture entropy (2.25). This equation can be expanded further upon introduction of the Gibbs relations (2.15) and (2.16) and the balance laws (2.8)–(2.13). By doing so and by collecting terms, we finally arrive at the balance equation for the mixture entropy. This equation has the form of (2.22) with the entropy flux $\boldsymbol{J}^{s}$ given by

(3.1) $$\begin{eqnarray}\displaystyle \boldsymbol{J}^{s}=\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}s_{f}(\boldsymbol{u}_{f}-\boldsymbol{u})+\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}s_{s}(\boldsymbol{u}_{s}-\boldsymbol{u})+\frac{\unicode[STIX]{x1D719}_{f}}{T_{f}}\boldsymbol{q}_{f}+\frac{\unicode[STIX]{x1D719}_{s}}{T_{s}}\boldsymbol{q}_{s}+\frac{\unicode[STIX]{x1D6FE}_{s}}{T_{s}}\frac{\text{d}\unicode[STIX]{x1D719}_{s}}{\text{d}_{s}t}\,\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}, & & \displaystyle\end{eqnarray}$$

and the entropy production rate $\unicode[STIX]{x1D70E}$ given by

(3.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E} & = & \displaystyle E\left(\frac{1}{T_{s}}-\frac{1}{T_{f}}\right)+\frac{\text{d}\unicode[STIX]{x1D719}_{s}}{\text{d}t_{s}}\left(\frac{p_{s}-\unicode[STIX]{x1D6FD}_{s}}{T_{s}}-\frac{p_{f}}{T_{f}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D6FE}_{s}}{T_{s}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right)\right)+\boldsymbol{f}\boldsymbol{\cdot }\left(\frac{\boldsymbol{u}_{f}}{T_{f}}-\frac{\boldsymbol{u}_{s}}{T_{s}}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{p_{f}}{T_{f}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\boldsymbol{\cdot }(\boldsymbol{u}_{f}-\boldsymbol{u}_{s})+\unicode[STIX]{x1D719}_{f}\boldsymbol{q}_{f}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{1}{T_{f}}\right)+\unicode[STIX]{x1D719}_{s}\boldsymbol{q}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{1}{T_{s}}\right)+\frac{\unicode[STIX]{x1D719}_{f}}{T_{f}}\unicode[STIX]{x1D64B}_{f}\boldsymbol{ : }\unicode[STIX]{x1D651}_{f}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{T_{s}}\left(\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D64B}_{s}+\unicode[STIX]{x1D63E}_{s}\right)\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}.\end{eqnarray}$$

In the last equation, the second-order tensor $\unicode[STIX]{x1D63E}_{s}$ is defined as

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63E}_{s}\triangleq \boldsymbol{h}_{s}\otimes \unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}=\unicode[STIX]{x1D6FE}_{s}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\otimes \unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}, & & \displaystyle\end{eqnarray}$$

where the symbol $\otimes$ denotes the dyadic product between two vectors.

Next, we follow the proposition of Bdzil et al. (Reference Bdzil, Menikoff, Son, Kapila and Stewart1999) who advocated to treat the differences $(T_{s}-T_{f})$ , $(p_{s}-\unicode[STIX]{x1D6FD}_{s}-p_{f})$ , and $(\boldsymbol{u}_{s}-\boldsymbol{u}_{f})$ as independent forces that drive the system toward equilibrium. This amounts to collecting the terms related to these differences in the expression for $\unicode[STIX]{x1D70E}$ above. Accordingly, we perform the following two decompositions:

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \left({\displaystyle \frac{p_{s}}{T_{s}}}-{\displaystyle \frac{p_{f}}{T_{f}}}\right)=\left({\displaystyle \frac{1}{T_{s}}}-{\displaystyle \frac{1}{T_{f}}}\right)p_{f}+\frac{1}{T_{s}}(p_{s}-p_{f}), & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \left({\displaystyle \frac{\boldsymbol{u}_{s}}{T_{s}}}-{\displaystyle \frac{\boldsymbol{u}_{f}}{T_{f}}}\right)=\left({\displaystyle \frac{1}{T_{s}}}-{\displaystyle \frac{1}{T_{f}}}\right)\boldsymbol{u}_{s}+\frac{1}{T_{f}}(\boldsymbol{u}_{s}-\boldsymbol{u}_{f}). & \displaystyle\end{eqnarray}$$

As will be shown below, these decompositions have important implications to the rate equation for $\unicode[STIX]{x1D719}_{s}$ , especially in the incompressible limit; see also Varsakelis & Papalexandris (Reference Varsakelis and Papalexandris2017). To be more explicit, this decomposition was chosen so as, in the limit, the thermomechanical model be able to predict dilatancy of the granular medium; see § 5 below. By contrast, in the model proposed by Papalexandris (Reference Papalexandris2004b ) this particular grouping of terms was not introduced. Consequently, that model cannot predict dilatancy of the granular material for constant-density flows; see Varsakelis & Papalexandris (Reference Varsakelis and Papalexandris2011).

Further, we introduce $\unicode[STIX]{x1D651}_{f}^{d}$ and $\unicode[STIX]{x1D64B}_{f}^{d}$ as the deviatoric parts of $\unicode[STIX]{x1D651}_{f}$ and $\unicode[STIX]{x1D64B}_{f}$ , respectively. They are defined by the standard decompositions

(3.6a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D651}_{f}^{d}=\unicode[STIX]{x1D651}_{f}-{\textstyle \frac{1}{3}}(\text{tr}\unicode[STIX]{x1D651}_{f})\,\unicode[STIX]{x1D644},\quad \unicode[STIX]{x1D64B}_{f}^{d}=\unicode[STIX]{x1D64B}_{f}+p_{f}^{v}\,\unicode[STIX]{x1D644}, & & \displaystyle\end{eqnarray}$$

with $p_{f}^{v}$ being the bulk viscous pressure of the fluid phase. The reason for not applying this decomposition to the granular phase will become evident below in the paper. By substituting relations (3.4), (3.5) and (3.6) to (3.2), we finally arrive at

(3.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E} & = & \displaystyle \left(E+\frac{\text{d}\unicode[STIX]{x1D719}_{s}}{\text{d}_{s}t}p_{f}\right)\left(\frac{1}{T_{s}}-\frac{1}{T_{f}}\!\right)+\frac{\text{d}\unicode[STIX]{x1D719}_{s}}{\text{d}t_{s}}\left(\frac{p_{s}-\unicode[STIX]{x1D6FD}_{s}-p_{f}}{T_{s}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D6FE}_{s}}{T_{s}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right)\right)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{\unicode[STIX]{x1D719}_{f}}{T_{f}}p_{f}^{v}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f}+\frac{1}{T_{f}}\left(\boldsymbol{f}-p_{f}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right)\boldsymbol{\cdot }(\boldsymbol{u}_{f}-\boldsymbol{u}_{s})+\unicode[STIX]{x1D719}_{f}\boldsymbol{q}_{f}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{1}{T_{f}}\right)+\unicode[STIX]{x1D719}_{s}\boldsymbol{q}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{1}{T_{s}}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x1D719}_{f}}{T_{f}}\unicode[STIX]{x1D64B}_{f}^{d}\boldsymbol{ : }\unicode[STIX]{x1D651}_{f}^{d}+\frac{1}{T_{s}}\left(\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D64B}_{s}+\unicode[STIX]{x1D63E}_{s}\right)\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}.\end{eqnarray}$$

Next we observe that, in view of (3.7), the entropy production rate is expressed as a sum of products with each product being associated with a particular dissipative process. In other words, the entropy production rate $\unicode[STIX]{x1D70E}$ can be written as

(3.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}=\mathop{\sum }_{i=1}^{8}X_{i}Y_{i}. & & \displaystyle\end{eqnarray}$$

In this expression, the quantities $X_{i}$ are identified as the thermodynamic forces, while the quantities $Y_{i}$ are identified as the thermodynamic currents. More specifically, in view of (3.7), the thermodynamic forces are

(3.9) $$\begin{eqnarray}\displaystyle & & \displaystyle X_{i}\in \left\{\left(\frac{1}{T_{s}}-\frac{1}{T_{f}}\right),\frac{1}{T_{s}}\left(p_{s}-\unicode[STIX]{x1D6FD}_{s}-p_{f}+T_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left({\displaystyle \frac{\unicode[STIX]{x1D6FE}_{s}}{T_{s}}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right)\right),\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f},\right.\nonumber\\ \displaystyle & & \displaystyle \qquad \qquad \left.(\boldsymbol{u}_{f}-\boldsymbol{u}_{s}),\unicode[STIX]{x1D735}\left(\frac{1}{T_{f}}\right),\unicode[STIX]{x1D735}\left(\frac{1}{T_{s}}\right),\unicode[STIX]{x1D651}_{f}^{d},\unicode[STIX]{x1D651}_{s}\right\},\end{eqnarray}$$

whereas the thermodynamic currents are

(3.10) $$\begin{eqnarray}\displaystyle & & \displaystyle Y_{i}\in \left\{\left(E+\frac{\text{d}\unicode[STIX]{x1D719}_{s}}{\text{d}_{s}t}p_{f}\right),\frac{\text{d}\unicode[STIX]{x1D719}_{s}}{\text{d}_{s}t},-\frac{\unicode[STIX]{x1D719}_{f}}{T_{f}}p_{f}^{v},\frac{1}{T_{f}}\left(\boldsymbol{f}-p_{f}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right),\unicode[STIX]{x1D719}_{f}\boldsymbol{q}_{f},\unicode[STIX]{x1D719}_{s}\boldsymbol{q}_{s},\frac{\unicode[STIX]{x1D719}_{f}}{T_{f}}\unicode[STIX]{x1D64B}_{f}^{d},\right.\nonumber\\ \displaystyle & & \displaystyle \qquad \qquad \left.\frac{1}{T_{s}}(\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D64B}_{s}+\unicode[STIX]{x1D63E}_{s})\right\}.\end{eqnarray}$$

It should be noted that both thermodynamic forces and currents vanish at equilibrium. Moreover, the entropy inequality (2.24) is written as

(3.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}=\mathop{\sum }_{i=1}^{8}X_{i}Y_{i}\geqslant 0. & & \displaystyle\end{eqnarray}$$

Next, we consider that the thermodynamic currents are functions of the forces. For all but the last of the currents above we assume linear relationships between currents and forces. This is equivalent to considering the Taylor expansions of these currents with respect to the forces around an equilibrium state and keeping the first-order terms of these expansions. In this manner, we arrive at

(3.12) $$\begin{eqnarray}\displaystyle Y_{i}=\mathop{\sum }_{j}L_{ij}X_{j},\quad i=1,\ldots 7. & & \displaystyle\end{eqnarray}$$

The coefficients $L_{ij}$ are the so-called phenomenological coefficients. By virtue of the presumed macroscopic isotropy of the two phases and of the representation theorem of linear isotropic tensors, the phenomenological coefficients above are scalars. Actually, the decompositions given by (3.6) were introduced so as to make more convenient the application of the representation theorem for the viscous stresses acting on the fluid.

As regards the off-diagonal elements of the matrix of phenomenological coefficients, we remark that coefficients that couple a current and a force of difference tensorial order vanish; this is a consequence of the representation theorem of linear isotropic tensors. Further, according to the principle of phase separation, coefficients that couple currents related to only one phase (such as the heat fluxes $\boldsymbol{q}_{f}$ and $\boldsymbol{q}_{s}$ ) with forces that depend on state variables of the other phase, vanish as well. Consequently, the only non-vanishing off-diagonal phenomenological coefficients are $L_{12}$ and $L_{21}$ . Also, by assuming the validity of Onsager reciprocity relations, these two coefficients are equal. Then, the only non-zero phenomenological coefficients are the following,

(3.13a,b ) $$\begin{eqnarray}\displaystyle L_{ii},\quad i=1,\ldots 7,\quad \text{and}\quad L_{12}=L_{21}. & & \displaystyle\end{eqnarray}$$

By virtue of the above equation and the linear relations (3.12), the entropy inequality (3.11) yields

(3.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}=\mathop{\sum }_{i=1}^{7}L_{ii}X_{i}^{2}+2L_{12}X_{1}X_{2}+X_{8}\,Y_{8}\geqslant 0, & & \displaystyle\end{eqnarray}$$

with

(3.15a,b ) $$\begin{eqnarray}\displaystyle X_{8}=\unicode[STIX]{x1D651}_{s},\quad Y_{8}=\frac{1}{T_{s}}(\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D64B}_{s}+\unicode[STIX]{x1D63E}_{s}). & & \displaystyle\end{eqnarray}$$

The current $Y_{8}$ has been singled out and will be treated separately because it represents the viscous stresses acting on the granular phase (divided by $T_{s}$ ). These stresses have a manifestly non-Newtonian behaviour that cannot be described by a purely linear theory. Accordingly, we provide the detailed derivation of their constitutive relations separately in the next subsection. For the time being, we remark that, by the principle of phase separation, $Y_{8}$ can only be a function of $X_{8}$ and possibly $X_{6}=\unicode[STIX]{x1D735}T_{s}^{-1}$ . If we further limit ourselves to linear dependence of $Y_{8}$ on temperature gradients then, by virtue of the representation theorem of linear isotropic tensors, the dependence of $Y_{8}$ on $X_{6}$ vanishes. Therefore, $Y_{8}$ becomes a function of $X_{8}$ only.

Furthermore, since $X_{8}$ is independent of the other thermodynamic forces, the entropy inequality (3.14) requires that the matrix of phenomenological coefficients $L_{ij}$ be positive definite and that the product $X_{8}Y_{8}$ be positive. In other words, the entropy inequality (3.14) is equivalent to the following conditions that should be independently satisfied,

(3.16) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}L_{ii}\geqslant 0,\quad i=1,\ldots 7,\\ L_{12}^{2}\leqslant L_{11}L_{22},\\ X_{8}Y_{8}\geqslant 0.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Finally, before deriving the constitutive relation for $Y_{8}$ , it is worth pointing out that the phenomenological coefficients $L_{ii}$ and $L_{12}$ can be recast in terms of transport coefficients. Those involving only one phase can be written as

(3.17a-d ) $$\begin{eqnarray}\displaystyle L_{44}=\unicode[STIX]{x1D719}_{f}T_{f}^{2}\,\unicode[STIX]{x1D705}_{f},\quad L_{55}=\unicode[STIX]{x1D719}_{s}T_{s}^{2}\,\unicode[STIX]{x1D705}_{s},\quad L_{66}={\displaystyle \frac{\unicode[STIX]{x1D719}_{f}}{T_{f}}}\,\unicode[STIX]{x1D701}_{f},\quad L_{77}={\displaystyle \frac{\unicode[STIX]{x1D719}_{f}}{T_{f}}}\,\unicode[STIX]{x1D707}_{f}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}_{f}$ and $\unicode[STIX]{x1D705}_{s}$ are the thermal conductivities of the two phases, while $\unicode[STIX]{x1D707}_{f}$ and $\unicode[STIX]{x1D701}_{f}$ are the first and second viscosity coefficients of the Newtonian stress tensor of the fluid phase, respectively.

The remaining phenomenological coefficients are related to interactions between the two phases and can be cast in the form of transport coefficients via the following definitions,

(3.18a-d ) $$\begin{eqnarray}\displaystyle L_{11}=T_{f}T_{s}\,\unicode[STIX]{x1D709},\quad L_{22}=\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D719}_{s}T_{s}{\displaystyle \frac{1}{\unicode[STIX]{x1D707}_{c}}},\quad L_{33}={\displaystyle \frac{1}{T_{f}}}\,\unicode[STIX]{x1D6FF},\quad L_{12}=L_{21}=T_{f}T_{s}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D706}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

In the above, $\unicode[STIX]{x1D709}$ represents the interphasial heat-transfer coefficient while $\unicode[STIX]{x1D6FF}$ represents the coefficient of interphasial drag. Also, $\unicode[STIX]{x1D707}_{c}$ has dimensions of viscosity and is commonly referred to as compaction viscosity. It expresses responsiveness of the granular material to compaction. Finally, $\unicode[STIX]{x1D706}$ is a coefficient for cross effect of the variation of $\unicode[STIX]{x1D719}_{s}$ due to thermal non-equilibrium between the two phases; in other words, $\unicode[STIX]{x1D706}$ expresses the capacity of the granular material to compact due to thermal loadings. According to the first of inequalities (3.16), all these transport coefficients are positive except for $\unicode[STIX]{x1D706}$ whose sign cannot be constrained on the basis of thermodynamics. Nonetheless, its amplitude is constrained by the second of inequalities (3.16) according to which,

(3.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}^{2}\leqslant {\displaystyle \frac{\unicode[STIX]{x1D709}}{\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D719}_{f}\,T_{f}\,\unicode[STIX]{x1D707}_{c}}}. & & \displaystyle\end{eqnarray}$$

3.1 Constitutive relations for the viscous stresses of the granular phase

As mentioned above, the current $Y_{8}$ is a function of $X_{8}$ only. The starting point for the derivation of a constitutive relation for $Y_{8}$ is the general representation theorem for isotropic tensors, Truesdell & Noll (Reference Truesdell and Noll2004). According to it, the expression for $Y_{8}$ reads,

(3.20) $$\begin{eqnarray}\displaystyle Y_{8}=\frac{1}{T_{s}}(\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D64B}_{s}+\unicode[STIX]{x1D63E}_{s})={\mathcal{G}}(\unicode[STIX]{x1D651}_{s})=g_{0}\unicode[STIX]{x1D644}+g_{1}\unicode[STIX]{x1D651}_{s}+g_{2}\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s}, & & \displaystyle\end{eqnarray}$$

where the coefficients $g_{0}$ , $g_{1}$ and $g_{2}$ are functions of the principal invariants of the strain-rate tensor $\unicode[STIX]{x1D651}_{s}$ and of the state variables of the granular phase.

In the derivation of the sought-after nonlinear constitutive expression for $Y_{8}$ we further require that $Y_{8}$ be non-singular and tend to zero as $\unicode[STIX]{x1D651}_{s}$ vanishes, i.e. at equilibrium. It is also noted that the component $\unicode[STIX]{x1D63E}_{s}$ of $Y_{8}$ is a function of the state variable $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}$ , cf. definition (3.3), which implies that it describes elastic stresses. As such, $\unicode[STIX]{x1D63E}_{s}$ does not vanish at equilibrium. In other words, the pressureless stress tensor $\unicode[STIX]{x1D64B}_{s}$ consists of an elastic component stemming from $\unicode[STIX]{x1D63E}_{s}$ plus a viscous component stemming from $Y_{8}$ which should tend to zero at equilibrium.

Additionally, we observe that, according to accumulated experimental and numerical evidence, the viscous stresses acting on granular suspensions admit a ‘first-order scaling’ with respect to the shear rate $\dot{\unicode[STIX]{x1D6FE}}$ . More specifically, in experiments of simple shear flow in the viscous regime, these stresses are found to be linear with respect to $|\dot{\unicode[STIX]{x1D6FE}}|$ , Boyer et al. (Reference Boyer, Guazzelli and Pouliquen2011a ), Boyer, Pouliquen & Guazzelli (Reference Boyer, Pouliquen and Guazzelli2011b ) and Dbouk, Lobry & Lemaire (Reference Dbouk, Lobry and Lemaire2013). The same has also been predicted in the simulations of Brady & Morris (Reference Brady and Morris1997), Morris & Boulay (Reference Morris and Boulay1999) and Foss & Brady (Reference Foss and Brady2000) that were based on Stokesian dynamics.

Guided by this evidence, we are seeking a constitutive relation for $Y_{8}$ that extends this scaling with respect to the shear rate to more general flows. Actually, this idea is not new and has been examined before by, among others, Fang et al. (Reference Fang, Mammoli, Brady, Ingber, Mondy and Graham2002) and Jop et al. (Reference Jop, Forterre and Pouliquen2006), who proposed a linear scaling with $\sqrt{\unicode[STIX]{x1D651}_{s}\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}}$ . We note, however, that a scaling with this quantity cannot provide by itself a tensorial representation of the stresses because this quantity is a scalar.

Herein, we propose to extend this scaling to the tensorial representation (3.20) by considering positive homogeneous functions of degree one in $\unicode[STIX]{x1D651}_{s}$ , so that

(3.21) $$\begin{eqnarray}\displaystyle {\mathcal{G}}(k\unicode[STIX]{x1D651}_{s})=k{\mathcal{G}}(\unicode[STIX]{x1D651}_{s}),\quad \forall \,k\geqslant 0. & & \displaystyle\end{eqnarray}$$

Next, we introduce the following invariants $l_{1}$ , $l_{2}$ and $l_{3}$ of $\unicode[STIX]{x1D651}_{s}$ ,

(3.22a-c ) $$\begin{eqnarray}\displaystyle l_{1}=\text{I}_{\boldsymbol{v}_{s}},\quad l_{2}=(J_{\boldsymbol{v}_{s}})^{1/2}=(-2\text{II}_{\boldsymbol{v}_{s}}+{\textstyle \frac{2}{3}}{\text{I}_{\boldsymbol{v}_{s}}}^{2})^{1/2},\quad l_{3}=(\text{III}_{\boldsymbol{v}_{s}})^{1/3}, & & \displaystyle\end{eqnarray}$$

where $\text{I}_{\boldsymbol{v}_{s}}$ , $\text{II}_{\boldsymbol{v}_{s}}$ and $\text{III}_{\boldsymbol{v}_{s}}$ are the principal invariants of $\unicode[STIX]{x1D651}_{s}$ . It can be directly verified that $l_{1}$ , $l_{2}$ and $l_{3}$ are homogeneous of degree one. Also, the following identity holds for $J_{\boldsymbol{v}_{s}}$ ,

(3.23) $$\begin{eqnarray}\displaystyle J_{\boldsymbol{v}_{s}}\triangleq \unicode[STIX]{x1D651}_{s}^{d}\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}^{d}. & & \displaystyle\end{eqnarray}$$

Since $J_{\boldsymbol{v}_{s}}$ is positive, there are bijective relations between the following invariant triplets,

(3.24) $$\begin{eqnarray}\displaystyle (\text{I}_{\boldsymbol{v}_{s}},\text{II}_{\boldsymbol{v}_{s}},\text{III}_{\boldsymbol{v}_{s}})\Longleftrightarrow (\text{I}_{\boldsymbol{v}_{s}},J_{\boldsymbol{v}_{s}},\text{III}_{\boldsymbol{v}_{s}})\Longleftrightarrow (l_{1},l_{2},l_{3}) & & \displaystyle\end{eqnarray}$$

so that the quantities $g_{0}$ , $g_{1}$ and $g_{2}$ appearing in (3.20) can be considered as functions of the invariants $l_{1}$ , $l_{2}$ and $l_{3}$ . Further, in view of the tensorial representation (3.20), equation (3.21) implies the following homogeneity constraints on these quantities,

(3.25) $$\begin{eqnarray}\displaystyle & \displaystyle g_{0}(kl_{1},kl_{2},kl_{3})=kg_{0}(l_{1},l_{2},l_{3}), & \displaystyle\end{eqnarray}$$

(3.26) $$\begin{eqnarray}\displaystyle & \displaystyle g_{1}(kl_{1},kl_{2},kl_{3})=g_{1}(l_{1},l_{2},l_{3}), & \displaystyle\end{eqnarray}$$

(3.27) $$\begin{eqnarray}\displaystyle & \displaystyle g_{2}(kl_{1},kl_{2},kl_{3})=\frac{1}{k}\,g_{2}(l_{1},l_{2},l_{3}). & \displaystyle\end{eqnarray}$$

Since $g_{2}$ is homogeneous of degree $-1$ , one has to verify that the boundedness constraint of $Y_{8}$ holds for any finite $\unicode[STIX]{x1D651}_{s}$ . To this end we remark that since

(3.28) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D651}_{s}\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}=\left({\displaystyle \frac{l_{1}^{2}}{3}}+l_{2}^{2}\right), & & \displaystyle\end{eqnarray}$$

one can define

(3.29) $$\begin{eqnarray}\displaystyle \tilde{g}_{2}=g_{2}\,\sqrt{{\displaystyle \frac{l_{1}^{2}}{3}}+l_{2}^{2}}, & & \displaystyle\end{eqnarray}$$

with $\tilde{g}_{2}$ homogeneous of degree zero. Then, the third term on the right-hand side of (3.20) is written as

(3.30) $$\begin{eqnarray}\displaystyle g_{2}\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s}=\frac{\tilde{g}_{2}}{\sqrt{{\displaystyle \frac{l_{1}^{2}}{3}}+l_{2}^{2}}}\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s}. & & \displaystyle\end{eqnarray}$$

According to this equation, if $\tilde{g}_{2}$ is bounded, then the term $g_{2}\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s}$ is always bounded and vanishes as $\unicode[STIX]{x1D651}_{s}\rightarrow 0$ .

A simple choice to satisfy the homogeneity requirements on $g_{0}$ , $g_{1}$ and $\tilde{g}_{2}$ would be to take the series expansions of these quantities up to the degree of homogeneity. This would yield

(3.31a-c ) $$\begin{eqnarray}\displaystyle g_{0}=(c_{1}l_{1}+c_{2}l_{2}+c_{3}l_{3}),\quad g_{1}=c_{4},\quad \tilde{g}_{2}=c_{5}, & & \displaystyle\end{eqnarray}$$

with the coefficients $c_{i}$ being functions of the state variables of the granular phase but independent of $\unicode[STIX]{x1D651}_{s}$ . However, as it turns out, this simple choice is unsatisfactory because it leads to the appearance of the third-order term

(3.32) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s})\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}={\textstyle \frac{3}{2}}l_{1}l_{2}^{2}+l_{3}^{3} & & \displaystyle\end{eqnarray}$$

in the product $X_{8}Y_{8}$ . The sign of this term is undetermined and so would be the sign of the product $X_{8}Y_{8}$ , since (3.32) is the highest-order term. This, in turn, implies that the entropy inequality would not be satisfied.

This conundrum can be circumvented by modifying the expressions (3.31), so as to remove the third-order term (3.32) from the entropy inequality. This is achieved by splitting the term (3.32) into two components and by subtracting one component from $g_{0}$ and the other component from $g_{1}$ , respectively. In this manner, the proposed expressions for $g_{0}$ and $g_{1}$ read,

(3.33) $$\begin{eqnarray}\displaystyle & \displaystyle g_{0}=(c_{1}l_{1}+c_{2}l_{2}+c_{3}l_{3})-\frac{3}{2}\frac{\tilde{g}_{2}}{\sqrt{{\displaystyle \frac{l_{1}^{2}}{3}}+l_{2}^{2}}}\,l_{2}^{2}, & \displaystyle\end{eqnarray}$$

(3.34) $$\begin{eqnarray}\displaystyle & \displaystyle g_{1}=c_{4}-\tilde{g}_{2}\frac{l_{3}^{3}}{\left(\displaystyle \frac{l_{1}^{2}}{3}+l_{2}^{2}\right)^{3/2}}. & \displaystyle\end{eqnarray}$$

By introducing the above expressions for $g_{0}$ and $g_{1}$ and the aforementioned expression $\tilde{g}_{2}=c_{5}$ , we arrive at the following expression for the pressureless stress tensor $\unicode[STIX]{x1D64B}_{s}$ ,

(3.35) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D64B}_{s}+\unicode[STIX]{x1D63E}_{s}) & = & \displaystyle \left(c_{1}l_{1}+c_{2}l_{2}+c_{3}l_{3}-c_{5}{\displaystyle \frac{3\,l_{2}^{2}}{2\sqrt{{\displaystyle \frac{l_{1}^{2}}{3}}+l_{2}^{2}}}}\right)\unicode[STIX]{x1D644}\nonumber\\ \displaystyle & & \displaystyle +\,\left(c_{4}-c_{5}{\displaystyle \frac{l_{3}^{3}}{\left({\displaystyle \frac{l_{1}^{2}}{3}}+l_{2}^{2}\right)^{3/2}}}\right)\unicode[STIX]{x1D651}_{s}~+c_{5}{\displaystyle \frac{1}{\sqrt{{\displaystyle \frac{l_{1}^{2}}{3}}+l_{2}^{2}}}}\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s}.\end{eqnarray}$$

It is straightforward to verify that this representation satisfies the homogeneity requirements (3.25)–(3.27). Furthermore, it is continuous with respect to its arguments, non-singular and tends to zero at equilibrium.

Also, with the above representation, the third inequality in (3.16) stemming from the second thermodynamic axiom reads,

(3.36) $$\begin{eqnarray}\displaystyle \frac{1}{T_{s}}(\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D64B}_{s}+\unicode[STIX]{x1D63E}_{s})\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}=\frac{1}{T_{s}}(c_{1}l_{1}+c_{2}l_{2}+c_{3}l_{3})l_{1}+\frac{1}{T_{s}}c_{4}\left({\displaystyle \frac{{l_{1}}^{2}}{3}}+{l_{2}}^{2}\right)\geqslant 0. & & \displaystyle\end{eqnarray}$$

This implies the following constraints on the coefficients $c_{i}$ ,

(3.37a-d ) $$\begin{eqnarray}\displaystyle c_{3}=0,\quad c_{4}\geqslant 0,\quad \left(c_{1}+{\displaystyle \frac{c_{4}}{3}}\right)\geqslant 0,\quad {c_{2}}^{2}\leqslant 4c_{4}\left(c_{1}+{\displaystyle \frac{c_{4}}{3}}\right). & & \displaystyle\end{eqnarray}$$

Next, we observe that we can decompose the second term on the right-hand side of (3.35) into a deviatoric and an isotropic component so that this equation is written in the following form,

(3.38) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D64B}_{s}+\unicode[STIX]{x1D63E}_{s}) & = & \displaystyle \left(c_{1}+{\displaystyle \frac{c_{4}}{3}}-{\displaystyle \frac{c_{5}}{3}}{\displaystyle \frac{\det \unicode[STIX]{x1D651}_{s}}{\left(\unicode[STIX]{x1D651}_{s}\mathbf{:}\unicode[STIX]{x1D651}_{s}\right)^{3/2}}}\right)(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s})\,\unicode[STIX]{x1D644}\nonumber\\ \displaystyle & & \displaystyle +\,\left(c_{4}-c_{5}{\displaystyle \frac{\det \unicode[STIX]{x1D651}_{s}}{\left(\unicode[STIX]{x1D651}_{s}\mathbf{:}\unicode[STIX]{x1D651}_{s}\right)^{3/2}}}\right)\unicode[STIX]{x1D651}_{s}^{d}\nonumber\\ \displaystyle & & \displaystyle +\,\left(c_{2}\sqrt{\unicode[STIX]{x1D651}_{s}^{d}\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}^{d}}-{\displaystyle \frac{3}{2}}c_{5}{\displaystyle \frac{\unicode[STIX]{x1D651}_{s}^{d}\mathbf{:}\unicode[STIX]{x1D651}_{s}^{d}}{\sqrt{\unicode[STIX]{x1D651}_{s}\boldsymbol{ :}\unicode[STIX]{x1D651}_{s}}}}\right)\unicode[STIX]{x1D644}\nonumber\\ \displaystyle & & \displaystyle +\,\left(c_{5}{\displaystyle \frac{1}{\sqrt{\unicode[STIX]{x1D651}_{s}\boldsymbol{ :}\unicode[STIX]{x1D651}_{s}}}}\right)\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s}.\end{eqnarray}$$

According to the above relation, two generalized Newtonian stress components can be identified by setting,

(3.39a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D701}_{s}\unicode[STIX]{x1D719}_{s}=\left(c_{1}+{\displaystyle \frac{c_{4}}{3}}-{\displaystyle \frac{c_{5}}{3}}{\displaystyle \frac{\det \unicode[STIX]{x1D651}_{s}}{\left(\unicode[STIX]{x1D651}_{s}\mathbf{:}\unicode[STIX]{x1D651}_{s}\right)^{3/2}}}\right),\quad 2\unicode[STIX]{x1D707}_{s}\unicode[STIX]{x1D719}_{s}=\left(c_{4}-c_{5}{\displaystyle \frac{\det \unicode[STIX]{x1D651}_{s}}{\left(\unicode[STIX]{x1D651}_{s}\mathbf{:}\unicode[STIX]{x1D651}_{s}\right)^{3/2}}}\right). & & \displaystyle\end{eqnarray}$$

The coefficients $\unicode[STIX]{x1D707}_{s}$ and $\unicode[STIX]{x1D701}_{s}$ are interpreted as the first and second (bulk) viscosity of the granular phase. It is interesting to observe that the representation (3.38) reduces to that of a Newtonian fluid when $c_{2}=c_{5}=0$ . In this limiting case, the constraints (3.37) translate to $\unicode[STIX]{x1D707}_{s}\geqslant 0$ and $\unicode[STIX]{x1D701}_{s}\geqslant 0$ , as expected. Nonetheless, according to (3.39), in the general case and for three-dimensional flows where $\det \unicode[STIX]{x1D651}_{s}\neq 0$ , there are additional contributions to $\unicode[STIX]{x1D707}_{s}$ and $\unicode[STIX]{x1D701}_{s}$ that depend on the strain rate.

Also, the sum of the remaining terms in (3.38),

(3.40) $$\begin{eqnarray}\displaystyle \left(c_{2}\sqrt{\unicode[STIX]{x1D651}_{s}^{d}\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}^{d}}-{\displaystyle \frac{3}{2}}c_{5}{\displaystyle \frac{\unicode[STIX]{x1D651}_{s}^{d}\mathbf{:}\unicode[STIX]{x1D651}_{s}^{d}}{\sqrt{\unicode[STIX]{x1D651}_{s}\boldsymbol{ :}\unicode[STIX]{x1D651}_{s}}}}\right)\unicode[STIX]{x1D644}+\left(c_{5}{\displaystyle \frac{1}{\sqrt{\unicode[STIX]{x1D651}_{s}\boldsymbol{ :}\unicode[STIX]{x1D651}_{s}}}}\right)\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s}, & & \displaystyle\end{eqnarray}$$

represents additional non-Newtonian viscous stresses of the granular phase. Its isotropic component is commonly referred to as the suspension pressure, see for example Morris & Boulay (Reference Morris and Boulay1999), and is responsible for the shear-induced dilatancy that is a characteristic feature of granular suspensions. Also, its deviatoric component gives rise to normal-stress differences under simple shear that are discussed in § 6.1.

We also observe that, in view of (3.39), if the bulk viscosity was assumed to be zero, then the following conditions would be true: $c_{5}=0$ and $(c_{1}+c_{4}/3)=0$ . These conditions must hold simultaneously because the coefficients $c_{i}$ do not depend on the velocity gradients. But then, the fourth constraint in (3.36), implies that $c_{2}=0$ also. In this case, the non-Newtonian viscous stresses given by (3.40) should vanish as well. However, this goes against preponderant experimental evidence that has established shear-induced dilatancy and normal-stress differences in granular suspensions. We therefore conclude that our model excludes the possibility of a Stokes-type hypothesis for the granular phase. Bulk viscosity is notoriously difficult to measure but the prediction of our model of non-zero bulk viscosity is consistent with the findings of Brady et al. (Reference Brady, Khair and Swaroop2006) who predicted non-zero $\unicode[STIX]{x1D701}_{s}$ via Stokesian dynamics simulations.

Finally, for the sake of clarity of exposition, we introduce the following short-hand notation,

(3.41a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D712}_{a}=\left(c_{2}\sqrt{\unicode[STIX]{x1D651}_{s}^{d}\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}^{d}}-{\displaystyle \frac{3}{2}}c_{5}{\displaystyle \frac{\unicode[STIX]{x1D651}_{s}^{d}\mathbf{:}\unicode[STIX]{x1D651}_{s}^{d}}{\sqrt{\unicode[STIX]{x1D651}_{s}\boldsymbol{ :}\unicode[STIX]{x1D651}_{s}}}}\right),\quad \unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D712}_{b}=c_{5}{\displaystyle \frac{1}{\sqrt{\unicode[STIX]{x1D651}_{s}\boldsymbol{ :}\unicode[STIX]{x1D651}_{s}}}}, & & \displaystyle\end{eqnarray}$$

which allows to express the stress tensor $\unicode[STIX]{x1D64B}_{s}$ in a concise manner,

(3.42) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D64B}_{s}+\unicode[STIX]{x1D63E}_{s})= & \unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D701}_{s}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s})\,\unicode[STIX]{x1D644}+2\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D707}_{s}\unicode[STIX]{x1D651}_{s}^{d}+\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D712}_{a}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D712}_{b}\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s}. & \displaystyle\end{eqnarray}$$

It should be stressed that since the coefficients $c_{i}$ are functions of the state variables, so are the coefficients $\unicode[STIX]{x1D701}_{s}$ , $\unicode[STIX]{x1D707}_{s}$ , $\unicode[STIX]{x1D712}_{a}$ and $\unicode[STIX]{x1D712}_{b}$ . In particular, they can depend on both $\unicode[STIX]{x1D719}_{s}$ and $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}$ .

4 Post-constitutive equations

The balance laws for the two phases, (2.8)–(2.10) and (2.11)–(2.13), supplemented by the constitutive equations (3.12) and (3.42), the saturation condition (2.1) and by appropriate equations of state for each phase, constitute a closed system of governing equations for granular suspensions. Further, by combining these relations, and upon introduction of the expressions for the transport coefficients (3.17) and (3.18), we finally arrive to the following post-constitutive equations.

Fluid phase

(4.1) $$\begin{eqnarray}\displaystyle \frac{\text{d}(\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f})}{\text{d}_{f}t}+\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f} & = & \displaystyle 0,\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}\frac{\text{d}\boldsymbol{u}_{f}}{\text{d}_{f}t}+\unicode[STIX]{x1D735}(p_{f}\unicode[STIX]{x1D719}_{f}) & = & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\left(\unicode[STIX]{x1D701}_{f}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f}\right)\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D644}+2\unicode[STIX]{x1D707}_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D651}_{f}^{d}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\left(\unicode[STIX]{x1D6FF}(\boldsymbol{u}_{f}-\boldsymbol{u}_{s})+p_{f}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right),\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}\frac{\text{d}e_{f}}{\text{d}_{f}t}+p_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f} & = & \displaystyle 2\unicode[STIX]{x1D707}_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D651}_{f}^{d}\boldsymbol{ : }\unicode[STIX]{x1D651}_{f}^{d}+\unicode[STIX]{x1D701}_{f}\unicode[STIX]{x1D719}_{f}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f})^{2}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D705}_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D735}T_{f})-\unicode[STIX]{x1D709}(T_{f}-T_{s})\nonumber\\ \displaystyle & & \displaystyle +\,\left(\unicode[STIX]{x1D6FF}(\boldsymbol{u}_{f}-\boldsymbol{u}_{s})+p_{f}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right)\boldsymbol{\cdot }(\boldsymbol{u}_{f}-\boldsymbol{u}_{s})+p_{f}\frac{\text{d}\unicode[STIX]{x1D719}_{s}}{\text{d}_{s}t}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D706}\,T_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D719}_{s}\left(p_{s}-\unicode[STIX]{x1D6FD}_{s}-p_{f}+T_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left({\displaystyle \frac{\unicode[STIX]{x1D6FE}_{s}}{T_{s}}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right)\right).\end{eqnarray}$$

Granular phase

(4.4) $$\begin{eqnarray}\displaystyle \frac{\text{d}(\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s})}{\text{d}_{s}t}+\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s} & = & \displaystyle 0,\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\frac{\text{d}\boldsymbol{u}_{s}}{\text{d}_{s}t}+\unicode[STIX]{x1D735}(p_{s}\unicode[STIX]{x1D719}_{s}) & = & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\left(\unicode[STIX]{x1D701}_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s}+\unicode[STIX]{x1D712}_{a}\right)\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D644}+2\unicode[STIX]{x1D707}_{s}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D651}_{s}^{d}+\unicode[STIX]{x1D712}_{b}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D63E}_{s}+\left(\unicode[STIX]{x1D6FF}(\boldsymbol{u}_{f}-\boldsymbol{u}_{s})+p_{f}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right),\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\frac{\text{d}e_{s}}{\text{d}_{s}t}+p_{s}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s} & = & \displaystyle 2\unicode[STIX]{x1D707}_{s}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D651}_{s}^{d}\!\boldsymbol{ : }\!\unicode[STIX]{x1D651}_{s}^{d}+\unicode[STIX]{x1D701}_{s}\unicode[STIX]{x1D719}_{s}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s})^{2}+\unicode[STIX]{x1D712}_{a}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D63E}_{s}\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}+\unicode[STIX]{x1D712}_{b}\unicode[STIX]{x1D719}_{s}(\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s})\boldsymbol{ : }\unicode[STIX]{x1D651}_{s}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D705}_{s}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D735}T_{s})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D709}(T_{f}-T_{s})-p_{f}\frac{\text{d}\unicode[STIX]{x1D719}_{s}}{\text{d}_{s}t}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D706}\,T_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D719}_{s}\left(p_{s}-\unicode[STIX]{x1D6FD}_{s}-p_{f}+T_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left({\displaystyle \frac{\unicode[STIX]{x1D6FE}_{s}}{T_{s}}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right)\right).\end{eqnarray}$$

Compaction equation

(4.7) $$\begin{eqnarray}\displaystyle \frac{\text{d}\unicode[STIX]{x1D719}_{s}}{\text{d}_{s}t}=\frac{\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D719}_{f}}{\unicode[STIX]{x1D707}_{c}}\left(p_{s}-\unicode[STIX]{x1D6FD}_{s}-p_{f}+T_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D6FE}_{s}}{T_{s}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right)\right)+\unicode[STIX]{x1D706}\,\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D719}_{s}(T_{f}-T_{s}). & & \displaystyle\end{eqnarray}$$

Equation (4.7) is the compaction equation and describes the evolution of the volume fraction of the granular phase. We remark that it is not derived by a conservation principle; instead it is follows directly from the constitutive relations (3.12) between thermodynamic currents and forces.

An interesting property of the system is that at the limit $\unicode[STIX]{x1D719}_{s}\rightarrow 0$ , all terms in (4.4)–(4.6) for the granular phase become zero, while (4.1)–(4.3) for the fluid reduce to the compressible Navier–Stokes equations.

It is also worth commenting on the form of the term that describes the interphasial momentum exchanges. This term reads

(4.8) $$\begin{eqnarray}\displaystyle \boldsymbol{f}=\unicode[STIX]{x1D6FF}(\boldsymbol{u}_{f}-\boldsymbol{u}_{s})+p_{f}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D6FF}$ being the interphasial drag coefficient, cf. (3.18). Jackson (Reference Jackson2000) proposed the decomposition of the interphasial momentum exchange into a drag and an Archimedean-type force that actually agrees with (4.8). This expression was subsequently employed by a number of modelling efforts, such as, for example, Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013) and Chiodi et al. (Reference Chiodi, Claudin and Andreotti2014), that further report that this decomposition is really a postulate. Our analysis, however, asserts that this decomposition is dictated by the second law of thermodynamics. More specifically, discarding this term leads to a system of governing equations that violates the positivity of the entropy production rate $\unicode[STIX]{x1D70E}$ . In other words, the particular form of (4.8) is not an ad hoc assumption but a necessary condition for thermodynamic consistency.

Actually, the inclusion of the term $p_{f}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}$ has a long history that seems to be underestimated by some contemporary studies. This term is usually referred to as nozzling term or non-conservative product; see the discussion in Powers (Reference Powers2004) and Papalexandris (Reference Papalexandris2004a ). Further, Bdzil et al. (Reference Bdzil, Menikoff, Son, Kapila and Stewart1999) remark that changes in the volume fraction can accelerate or decelerate the flow thereby acting as a nozzle. In turn, this can engender a relative motion between the two phases. The labelling of $p_{f}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}$ as a non-conservative product is also of interest. If one considers the inviscid limit of the above equations, then it is straightforward to show that, due to the presence of $p_{f}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}$ , these equations cannot be cast into a conservative form. As a consequence, even weak solutions cannot be defined and the corresponding Riemann problem is ill posed (Deledicque & Papalexandris Reference Deledicque and Papalexandris2007).

At this point, some remarks about the calibration of the various parameters that appear in the post-constitutive equations are appropriate. The parameters related to the fluid phase are the shear and bulk viscosities, $\unicode[STIX]{x1D707}_{f}$ and $\unicode[STIX]{x1D701}_{f}$ , and the thermal conductivity $\unicode[STIX]{x1D705}_{f}$ . They have been measured for a wide class of Newtonian fluids in a variety of thermodynamic conditions and their values are readily available. The situation with the parameters related to the granular phase is more complex. Due to the structure of this phase at the grain (mesoscopic) level, these parameters generally depend on $\unicode[STIX]{x1D719}_{s}$ . For constant-density flows, various expressions for the viscosities $\unicode[STIX]{x1D707}_{s}$ , $\unicode[STIX]{x1D712}_{a}$ and $\unicode[STIX]{x1D712}_{b}$ are available, mostly derived from shear-flow experiments; see, for example, Morris & Boulay (Reference Morris and Boulay1999), Jop et al. (Reference Jop, Forterre and Pouliquen2006) and Boyer et al. (Reference Boyer, Pouliquen and Guazzelli2011b ). On the other hand, data for the bulk viscosity $\unicode[STIX]{x1D701}_{s}$ are available only for very dilute mixtures, Brady et al. (Reference Brady, Khair and Swaroop2006). Also, the variation of these viscosities with temperature is largely unexplored. For the thermal conductivity $\unicode[STIX]{x1D705}_{s}$ , experimental correlations are available in the literature, cf. Carson et al. (Reference Carson, Lovatt, Tanner and Cleland2005) and Samantray, Karthikeyan & Reddy (Reference Samantray, Karthikeyan and Reddy2006). As regards $\unicode[STIX]{x1D6FD}_{s}$ , relations based on quasi-static compaction experiments have been provided by Elban & Chiarito (Reference Elban and Chiarito1986) and Mesri & Vardhanabhuti (Reference Mesri and Vardhanabhuti2009); see also the analysis in Bdzil et al. (Reference Bdzil, Menikoff, Son, Kapila and Stewart1999). By contrast, no detailed experimental data are available for $\unicode[STIX]{x1D6FE}_{s}$ . For this reason, people often assume a power-law dependence on $\unicode[STIX]{x1D719}_{s}$ based mostly on physical intuition. As regards the parameters related to phase interactions, various correlations have been proposed for the interphasial drag and heat-transfer coefficients, $\unicode[STIX]{x1D6FF}$ and $\unicode[STIX]{x1D709}$ , that are readily available in the literature. Much less is known for the compaction viscosity $\unicode[STIX]{x1D707}_{c}$ . It is often considered to be a constant whose magnitude is approximated via simplified transport models, Baer & Nunziato (Reference Baer and Nunziato1986). Similarly, no data are available for the off-diagonal coefficient $\unicode[STIX]{x1D706}$ . In many situations, however, one may assume that the evolution of the volume fraction $\unicode[STIX]{x1D719}_{s}$ depends only weakly on the thermal non-equilibrium between the two phases, especially in view of constraint (3.19). In any case, it is expected that numerical simulations based on this model, combined with available and forthcoming experimental data, will yield more accurate estimations of these parameters.

We conclude this section with a discussion on the boundary conditions that are needed to close the above system of equations. The assignment of boundary conditions for the variables related to the fluid phase, $\boldsymbol{u}_{f}$ and $p_{f}$ does not pose any problem and it can be performed in the same way as in single-phase flows. For example, at rigid walls, we can assume zero slip for the fluid velocity and Neumann condition for the fluid pressure. The same is true for the pressure of the granular phase $p_{s}$ . On the other hand, the assignment of boundary condition for $\boldsymbol{u}_{s}$ at walls requires more attention because we might have non-zero slip. Appropriate slip conditions are hard to derive and are material dependent; nonetheless, there are some works available in the literature, cf. Soltani & Yilmazer (Reference Soltani and Yilmazer1998). Similarly, the derivation of wall conditions for the volume fraction $\unicode[STIX]{x1D719}_{s}$ is quite challenging. Quite often a Neumann condition is employed but the value of the normal derivative is determined via the compaction equation or is set to zero. However, these choices are dictated by convenience and the issue is still open; additional physical evidence is required in order to derive more accurate boundary conditions for $\unicode[STIX]{x1D719}_{s}$ . Finally, as regards thermal boundary conditions at walls, quite often one may assume that the two phases are in thermal equilibrium since they are both in contact with the walls. Then, the prescription of thermal conditions can be performed in a similar manner as in single-phase flows.

5 The incompressible limit

The thermomechanical model given by (4.1)–(4.7) is valid for both compressible and incompressible flows. Nonetheless, in many practical applications, the velocities of the two phases are very small with respect to the speed of sound of the two phases and, consequently, compressibility effects are negligible. This is identified as the incompressible limit. In this case, one can employ the square of the smallest of the two Mach numbers as a perturbation parameter and perturb the governing equations with respect to this parameter. The result of this procedure is a system of equations, referred to as the ‘low-Mach-number approximation’, that is simpler than the original one in several important aspects.

In this section, we present the derivation of the low-Mach-number approximation for the model given by (4.1)–(4.7). The procedure is laborious but straightforward. For this reason, herein we present only the basic steps and the final result. For a more detailed description on the derivation of the low-Mach-number approximation for two-phase flow models, the reader is referred to Varsakelis & Papalexandris (Reference Varsakelis and Papalexandris2011).

Step 1. We consider a reference length $l_{r}$ and a reference velocity $u_{r}$ . Also, for non-dimensionalization purposes, we consider reference values of the thermodynamic variables $\unicode[STIX]{x1D70C}_{r}$ , $T_{r}$ , $p_{r}$ , specific heats ${c_{p}}_{r}$ and ${c_{v}}_{r}$ and transport coefficients $\unicode[STIX]{x1D707}_{r}$ and $\unicode[STIX]{x1D705}_{r}$ with respect to a reference state of the phase with the smallest speed of sound. In most practical cases, this is the fluid phase. Herein we make this assumption as well.

Step 2. We introduce equations of state for each phase. For example, the fluid can be considered to obey the ideal gas law, so that $p_{r}=({c_{p}}_{r}-{c_{v}}_{r})\,\unicode[STIX]{x1D70C}_{r}T_{r}$ .

Step 3. We introduce the reference Mach number $M$ as the ratio between the reference velocity and the fluid speed of sound at a reference state. Under the ideal gas law,

(5.1a,b ) $$\begin{eqnarray}\displaystyle M=\frac{u_{r}}{\sqrt{\displaystyle \frac{\hat{\unicode[STIX]{x1D6FE}}p_{r}}{\unicode[STIX]{x1D70C}_{r}}}},\quad \hat{\unicode[STIX]{x1D6FE}}=\frac{{c_{p}}_{r}}{{c_{v}}_{r}}. & & \displaystyle\end{eqnarray}$$

Then, the perturbation parameter $\unicode[STIX]{x1D716}$ is defined as

(5.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}\triangleq \hat{\unicode[STIX]{x1D6FE}}M^{2}, & & \displaystyle\end{eqnarray}$$

Step 4. The governing equations and the equations of state are non-dimensionalized with respect to the reference quantities.

Henceforth, in order to avoid the introduction of additional symbols, every variable is understood to be non-dimensionalized (except for the reference quantities).

Step 5. All variables are expanded in perturbation series of powers of $\unicode[STIX]{x1D716}$ . For example,

(5.3a,b ) $$\begin{eqnarray}\displaystyle p_{s}=p_{s}^{0}+\unicode[STIX]{x1D716}p_{s}^{1}+\cdots \,,\quad \boldsymbol{u}_{s}=\boldsymbol{u}_{s}^{0}+\unicode[STIX]{x1D716}\boldsymbol{u}_{s}^{1}+\cdots \,,\quad \text{etc}. & & \displaystyle\end{eqnarray}$$

where the superscript refers to the order of a term in the expansion.

Step 6. We insert the perturbation expansions of all variables to the governing system of equations and collect terms of the same order.

The results are presented below. Henceforth, in order to avoid unnecessary complexity in the notation, all quantities are supposed to represent the zeroth-order term of the expansions and, therefore, the superscript zero is dropped. For example, $p_{s}^{0}$ is now written as $p_{s}$ etc. The only exception to this convention is for the first-order terms of the pressures which will still be denoted as $p_{s}^{1}$ and $p_{f}^{1}$ . These are typically referred to as ‘dynamic pressures’, whereas the zeroth-order terms are referred to as ‘thermodynamic pressures’.

At order $O(\unicode[STIX]{x1D716}^{-1})$ the momentum equations of the two phases (4.2) and (4.5) yield,

(5.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}(p_{f}\unicode[STIX]{x1D719}_{f})+p_{f}\unicode[STIX]{x1D735}(\unicode[STIX]{x1D719}_{s})=0, & & \displaystyle\end{eqnarray}$$

(5.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}(p_{s}\unicode[STIX]{x1D719}_{s})-p_{f}\unicode[STIX]{x1D735}(\unicode[STIX]{x1D719}_{s})=0. & & \displaystyle\end{eqnarray}$$

Since $p_{f}$ , $p_{s}$ and $\unicode[STIX]{x1D719}_{s}$ are independent variables, equations (5.4) and (5.5) reduce to

(5.6a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}p_{s}=0=\unicode[STIX]{x1D735}p_{f},\quad \text{and}\quad p\triangleq p_{s}(t)=p_{f}(t). & & \displaystyle\end{eqnarray}$$

Also, in view of this equation, the compaction equation is identically satisfied at order $O(\unicode[STIX]{x1D716}^{-1})$ . Further, at order $O(1)$ , we arrive at the following equations.

Fluid phase

(5.7) $$\begin{eqnarray}\displaystyle \frac{\text{d}(\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f})}{\text{d}_{f}t}+\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f} & = & \displaystyle 0,\end{eqnarray}$$

(5.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}\frac{\text{d}\boldsymbol{u}_{f}}{\text{d}_{f}t}+\unicode[STIX]{x1D735}(p_{f}^{1}\unicode[STIX]{x1D719}_{f}) & = & \displaystyle \frac{1}{Re}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\left(\unicode[STIX]{x1D701}_{f}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f}\right)\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D644}+2\unicode[STIX]{x1D707}_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D651}_{f}^{d}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\left(\unicode[STIX]{x1D6FF}(\boldsymbol{u}_{f}-\boldsymbol{u}_{s})+p_{f}^{1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right),\end{eqnarray}$$

(5.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}\frac{\text{d}h_{f}}{\text{d}_{f}t} & = & \displaystyle \frac{1}{Re\,Pr}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D705}_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D735}T_{f})-\unicode[STIX]{x1D709}(T_{f}-T_{s})+\unicode[STIX]{x1D719}_{f}\frac{\hat{\unicode[STIX]{x1D6FE}}-1}{\hat{\unicode[STIX]{x1D6FE}}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}t}.\end{eqnarray}$$

Granular phase

(5.10) $$\begin{eqnarray}\displaystyle \frac{\text{d}(\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s})}{\text{d}_{s}t}+\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s} & = & \displaystyle 0,\end{eqnarray}$$

(5.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\frac{\text{d}\boldsymbol{u}_{s}}{\text{d}_{s}t}+\unicode[STIX]{x1D735}(p_{s}^{1}\unicode[STIX]{x1D719}_{s}) & = & \displaystyle \frac{1}{Re}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\left(\unicode[STIX]{x1D701}_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s}+\unicode[STIX]{x1D712}_{a}\right)\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D644}+2\unicode[STIX]{x1D707}_{s}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D651}_{s}^{d}~+\unicode[STIX]{x1D712}_{b}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D63E}_{s}+\left(\unicode[STIX]{x1D6FF}(\boldsymbol{u}_{f}-\boldsymbol{u}_{s})+p_{f}^{1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right),\end{eqnarray}$$

(5.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\frac{\text{d}h_{s}}{\text{d}_{s}t} & = & \displaystyle \frac{1}{Re\,Pr}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D705}_{s}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D735}T_{s})+\unicode[STIX]{x1D709}(T_{f}-T_{s})+\unicode[STIX]{x1D719}_{s}\frac{\hat{\unicode[STIX]{x1D6FE}}-1}{\hat{\unicode[STIX]{x1D6FE}}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}t}.\end{eqnarray}$$

Compaction equation

(5.13) $$\begin{eqnarray}\displaystyle \frac{\text{d}\unicode[STIX]{x1D719}_{s}}{\text{d}_{s}t}=Re_{c}\,\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D719}_{f}\left(p_{s}^{1}-\unicode[STIX]{x1D6FD}_{s}-p_{f}^{1}+T_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D6FE}_{s}}{T_{s}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right)\right)+\unicode[STIX]{x1D706}\,\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D719}_{s}(T_{f}-T_{s}). & & \displaystyle\end{eqnarray}$$

In the above equations, $h_{f}$ and $h_{s}$ stand for the enthalpies of the two phases. Also $Re$ and $Pr$ stand for the Reynolds and Prandtl number, respectively, based on the reference variables. Further, the non-dimensional number $Re_{c}$ is the Reynolds number based on the compaction viscosity, $\unicode[STIX]{x1D707}_{c}$ . Finally, we note that for the derivation of these results, we assumed that $\unicode[STIX]{x1D6FD}_{s}$ and $\unicode[STIX]{x1D6FE}_{s}$ scale with $\unicode[STIX]{x1D70C}_{r}u_{r}^{2}$ and $\unicode[STIX]{x1D70C}_{r}u_{r}^{2}l_{r}^{2}$ , respectively; see the relative discussion in Varsakelis & Papalexandris (Reference Varsakelis and Papalexandris2011).

The most important simplification in the incompressible limit is that the thermodynamic pressures no longer appear in the momentum equations. In this sense, the energy equation for each phase is essentially decoupled from the corresponding momentum equation. An interesting feature of the model is that it predicts thermal non-equilibrium between the two phases even at low Mach numbers, cf. (5.9) and (5.12). Moreover, when the thermodynamic pressure $p$ remains constant, as in the case of open domains, the energy equations above reduce to the standard advection–diffusion equations for $T_{f}$ and $T_{s}$ with the additional relaxation term representing interphasial heat exchange.

5.1 Special case: constant-density flows

We now consider the important special case where the densities of the two phases are constant, and the two temperatures $T_{s}$ and $T_{f}$ are also constant and equal to each other. Then, the energy equations (5.9) and (5.12) are trivially satisfied and the remaining equations simplify to the following system,

(5.14) $$\begin{eqnarray}\displaystyle \frac{\text{d}\unicode[STIX]{x1D719}_{f}}{\text{d}_{f}t}+\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f} & = & \displaystyle 0,\end{eqnarray}$$

(5.15) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D719}_{f}\frac{\text{d}\boldsymbol{u}_{f}}{\text{d}_{f}t}+\unicode[STIX]{x1D735}(p_{f}^{1}\unicode[STIX]{x1D719}_{f}) & = & \displaystyle \frac{1}{Re}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\left(\unicode[STIX]{x1D701}_{f}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f}\right)\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D644}+2\unicode[STIX]{x1D707}_{f}\unicode[STIX]{x1D719}_{f}\unicode[STIX]{x1D651}_{f}^{d}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\left(\unicode[STIX]{x1D6FF}(\boldsymbol{u}_{f}-\boldsymbol{u}_{s})+p_{f}^{1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right),\end{eqnarray}$$

(5.16) $$\begin{eqnarray}\displaystyle \frac{\text{d}\unicode[STIX]{x1D719}_{s}}{\text{d}_{s}t}+\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s} & = & \displaystyle 0,\end{eqnarray}$$

(5.17) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{s}\frac{\text{d}\boldsymbol{u}_{s}}{\text{d}_{s}t}+\unicode[STIX]{x1D735}(p_{s}^{1}\unicode[STIX]{x1D719}_{s}) & = & \displaystyle \frac{1}{Re}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\left(\unicode[STIX]{x1D701}_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s}+\unicode[STIX]{x1D712}_{a}\right)\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D644}+2\unicode[STIX]{x1D707}_{s}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D651}_{s}^{d}~+\unicode[STIX]{x1D712}_{b}\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D651}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D651}_{s}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D63E}_{s}+\left(\unicode[STIX]{x1D6FF}(\boldsymbol{u}_{f}-\boldsymbol{u}_{s})+p_{f}^{1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}\right),\end{eqnarray}$$

(5.18) $$\begin{eqnarray}\displaystyle \frac{\text{d}\unicode[STIX]{x1D719}_{s}}{\text{d}_{s}t} & = & \displaystyle Re_{c}\,\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D719}_{f}\left(p_{s}^{1}-\unicode[STIX]{x1D6FD}_{s}-p_{f}^{1}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FE}_{s}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s})\right).\end{eqnarray}$$

The above constant-density model has the important property that it does not constrain the velocity field of either phase to be divergence free. This is a direct consequence of the compaction equation (5.18) which stipulates that the volume fraction $\unicode[STIX]{x1D719}_{s}$ is not conserved along streamlines of the granular phase. We further note that, in a manner analogous to the Navier–Stokes equations which require that the dynamic pressure be determined only up to an additive constant, the above model requires that the two dynamic pressures $p_{s}^{1}$ and $p_{f}^{1}$ be determined up to the same additive constant. This is due to the structure of the compaction equation (5.18) which requires that the difference $(p_{s}^{1}-p_{f}^{1})$ be uniquely determined.

Also, the above model has the property that at the single-phase limit, $\unicode[STIX]{x1D719}_{s}\rightarrow 0$ , the equations for the fluid phase reduce to the Navier–Stokes equations, as they should. Finally, it is worth mentioning that algorithms for the numerical integration of systems of equations of this type have been recently proposed by Varsakelis & Papalexandris (Reference Varsakelis and Papalexandris2014) and Varsakelis et al. (Reference Varsakelis, Monsorno and Papalexandris2015).

6 Examples

In this section, and in order to gain insight into the structure and the predictive capacity of our model, we apply the constant-density equations (5.14)–(5.18) to two different flows of a granular suspension between two parallel plates. In particular, we consider the cases of (i) steady homogeneous shear flow with $\unicode[STIX]{x1D719}_{s}=const.$ , and (ii) pressure-driven fully developed flow.

6.1 Steady shear flow

In this example we assume that the bottom plate remains at rest whereas the upper plate is moving at a constant velocity. The (non-dimensional) shear rate, $\dot{\unicode[STIX]{x1D6FE}}$ , is then defined as the ratio between the velocity of the upper plate divided by the distance between the two plates. Assuming no-slip boundary conditions at the wall, it is easy to verify that the solution to this problem consists of linear velocity profiles given by

(6.1) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{s}=\boldsymbol{u}_{f}=\dot{\unicode[STIX]{x1D6FE}}y\hat{\boldsymbol{e}}_{x}, & & \displaystyle\end{eqnarray}$$

$\boldsymbol{e}_{x}$ being the unit vector pointing to the streamwise direction. The dynamic pressures and solid volume fraction $\unicode[STIX]{x1D719}_{s}$ are also constant and satisfy the following relation,

(6.2) $$\begin{eqnarray}p_{s}^{1}-p_{f}^{1}=\unicode[STIX]{x1D6FD}_{s}(\unicode[STIX]{x1D719}_{s})=\text{const}.\end{eqnarray}$$

By virtue of the above relation, the compaction equation (5.18) is also satisfied. Then, the strain-rate tensors are given by

(6.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D651}_{f}=\unicode[STIX]{x1D651}_{s}=\dot{\unicode[STIX]{x1D6FE}}\left(\begin{array}{@{}ccc@{}}0 & \frac{1}{2} & 0\\ \frac{1}{2} & 0 & 0\\ 0 & 0 & 0\end{array}\right) & & \displaystyle\end{eqnarray}$$

and have the following invariants, cf. (3.22).

(6.4a-c ) $$\begin{eqnarray}\displaystyle l_{1}=0,\quad l_{2}=\frac{1}{\sqrt{2}}|\dot{\unicode[STIX]{x1D6FE}}|,\quad l_{3}=0. & & \displaystyle\end{eqnarray}$$

Further, the non-dimensional stress tensors reduce to the following expressions, cf. (3.38),

(6.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64B}_{f}=\frac{\dot{\unicode[STIX]{x1D6FE}}}{Re}\left(\begin{array}{@{}ccc@{}}0 & 1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 0\end{array}\right), & \displaystyle\end{eqnarray}$$

(6.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64B}_{s}=\frac{\unicode[STIX]{x1D707}_{s}\dot{\unicode[STIX]{x1D6FE}}}{Re}\left(\begin{array}{@{}ccc@{}}0 & 1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 0\end{array}\right)-\frac{\unicode[STIX]{x1D707}_{n}|\dot{\unicode[STIX]{x1D6FE}}|}{Re}\left(\begin{array}{@{}ccc@{}}1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & \displaystyle \frac{c_{2}-\frac{3}{2}c_{5}}{c_{2}-c_{5}}\end{array}\right)\,. & \displaystyle\end{eqnarray}$$

In the last expression, the shear viscosity of the granular phase $\unicode[STIX]{x1D707}_{s}$ is given by (3.39). Moreover, $\unicode[STIX]{x1D707}_{n}$ is the so-called normal viscosity, defined by Morris & Boulay (Reference Morris and Boulay1999). In the context of our proposed thermomechanical model, $\unicode[STIX]{x1D707}_{n}$ is given by

(6.7) $$\begin{eqnarray}\unicode[STIX]{x1D707}_{n}=\frac{1}{\sqrt{2}\,\unicode[STIX]{x1D719}_{s}}(c_{5}-c_{2}).\end{eqnarray}$$

Herein, both $\unicode[STIX]{x1D707}_{s}$ and $\unicode[STIX]{x1D707}_{n}$ are understood to have been non-dimensionalized by the fluid’s shear viscosity.

The stress tensor of the granular phase $\unicode[STIX]{x1D719}_{s}\unicode[STIX]{x1D64B}_{s}$ can also be written as a function of the anisotropy coefficients, following the procedure of Morris & Boulay (Reference Morris and Boulay1999). Accordingly, the expression for $\unicode[STIX]{x1D64B}_{s}$ reads,

(6.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64B}_{s}=\frac{2\,\unicode[STIX]{x1D707}_{s}}{Re}\unicode[STIX]{x1D651}_{s}^{d}-\frac{\unicode[STIX]{x1D707}_{n}|\dot{\unicode[STIX]{x1D6FE}}|}{Re}\,\hat{\unicode[STIX]{x1D64C}},\quad \text{with}~\hat{\unicode[STIX]{x1D64C}}=\left(\begin{array}{@{}ccc@{}}1 & 0 & 0\\ 0 & \unicode[STIX]{x1D706}_{2} & 0\\ 0 & 0 & \unicode[STIX]{x1D706}_{3}\end{array}\right), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D706}_{2}$ and $\unicode[STIX]{x1D706}_{3}$ are the anisotropy coefficients. In our case, their values are the following.

(6.9a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{2}=1,\quad \unicode[STIX]{x1D706}_{3}=\frac{c_{2}-\displaystyle \frac{3}{2}c_{5}}{c_{2}-c_{5}}. & & \displaystyle\end{eqnarray}$$

Then, the first and second normal-stress differences become proportional to $(\unicode[STIX]{x1D706}_{2}-1)$ and $(\unicode[STIX]{x1D706}_{3}-\unicode[STIX]{x1D706}_{2})$ , respectively.

As regards the second normal-stress difference, which scales with $(\unicode[STIX]{x1D706}_{3}-\unicode[STIX]{x1D706}_{2})$ , we remark that in our model $\unicode[STIX]{x1D706}_{3}$ depends on the ratio $c_{5}/c_{2}$ which has to be measured experimentally. Nonetheless, if we set

(6.10) $$\begin{eqnarray}\displaystyle c_{5}={\textstyle \frac{1}{2}}\,c_{2}, & & \displaystyle\end{eqnarray}$$

then $\unicode[STIX]{x1D706}_{3}=0.5$ , which gives the same normal stresses that are predicted by the model of Brady & Morris (Reference Brady and Morris1997). There is also a general agreement that $\unicode[STIX]{x1D706}_{3}=0.5$ and that the value of $\unicode[STIX]{x1D706}_{2}$ is close to unity, which is the prediction of our model. For example, Brady & Carpen (Reference Brady and Carpen2002) and Fang et al. (Reference Fang, Mammoli, Brady, Ingber, Mondy and Graham2002) set explicitly $\unicode[STIX]{x1D706}_{2}=1$ in their computations and Stickel, Phillips & Powell (Reference Stickel, Phillips and Powell2007) reached a similar conclusion with the constitutive model that they proposed. On the other hand, the experimental measurements of Zarraga et al. (Reference Zarraga, Hill and Leighton2000), Singh & Nott (Reference Singh and Nott2003) and Boyer et al. (Reference Boyer, Pouliquen and Guazzelli2011b ) are compatible with $\unicode[STIX]{x1D706}_{2}\simeq 0.85$ and $\unicode[STIX]{x1D706}_{3}=0.5\,$ . Also, the simulations of Morris & Boulay (Reference Morris and Boulay1999) are in good agreement with these values of the anisotropy coefficients.

As regards the first normal-stress difference, since it scales with $(\unicode[STIX]{x1D706}_{2}-1)$ , our model predicts that it equals to zero. Actually, it can be shown that the representation theorem (3.20), upon which our model is based, implies a priori that is zero. We note that this has been explicitly assumed in a number of previous studies. In any case, there is a consensus that the first normal-stress difference is much smaller than the second one. However, while definite conclusions cannot be drawn yet, there is experimental evidence that the first normal-stress difference admits a small albeit non-zero value; see, for example, Zarraga et al. (Reference Zarraga, Hill and Leighton2000), Singh & Nott (Reference Singh and Nott2003) and Boyer et al. (Reference Boyer, Pouliquen and Guazzelli2011b ).

In our view, the most likely cause of it is the slow relaxation of the internal ordering of the granular phase. In fact, as is known from experimental observations and discrete-element simulations, flowing granular suspensions present a partially ordered structure because the macroscopic motion of the particles cannot be completely decoupled from their microscale configuration; a review on this issue has been provided by Morris (Reference Morris2009). In the derivation of our model, the assumption that $\unicode[STIX]{x1D64B}_{s}$ is an algebraic function of $\unicode[STIX]{x1D651}_{s}$ is equivalent to assuming that the dynamics controlling the microscale configuration is very fast with respect to the characteristic time scales of the flow. In other words, this assumption implies that any additional internal degree of freedom at the grain level relaxes very fast with respect to the flow time scales so that effectively the stresses are determined by $\unicode[STIX]{x1D719}_{s}$ , $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{s}$ and the deformation tensor $\unicode[STIX]{x1D651}_{s}$ .

This assumption, however, can be removed by the introduction of an additional tensorial order parameter that relates to the microstructure of the granular phase; say $\boldsymbol{A}_{s}$ . Possible candidates of such an order parameter include fabric tensors, Phan-Thien (Reference Phan-Thien1995), Hinrichsen & Wolf (Reference Hinrichsen and Wolf2004) and Goddard (Reference Goddard2006), a microstructure tensor based on the mean-free path length (Stickel, Phillips & Powell Reference Stickel, Phillips and Powell2006) and others. This approach allows for the description of time-dependent stresses related to relaxation phenomena. We remark that within the framework of non-equilibrium thermodynamics, the introduction of an order parameter $\boldsymbol{A}_{s}$ provides an extra term to the Gibbs equation, $-\boldsymbol{B}_{s}\boldsymbol{ : }\text{d}\boldsymbol{A}_{s}$ , with $\boldsymbol{B}_{s}$ being the partial derivative of the Helmholtz free energy of the granular phase with respect to $\boldsymbol{A}_{s}$ . Consequently, by exploiting the entropy inequality, one can derive a rate equation for $\boldsymbol{A}_{s}$ and a constitutive equation for $\unicode[STIX]{x1D64B}_{s}$ that depend both on $\unicode[STIX]{x1D651}_{s}$ and on $\boldsymbol{B}_{s}$ . This procedure has already been used by Lhuillier (Reference Lhuillier2001) and Bedeaux & Rubi (Reference Bedeaux and Rubi2002) for the case of colloidal suspensions.

In the present work, we have chosen to refrain from the introduction of a tensorial order parameter related to the granular microstructure for the following two reasons. First, currently there exists no univocal definition of such a tensorial order parameter and, to the best of our knowledge, there is not sufficient experimental guidance about how the free energy of the granular phase should depend on it. Second, from a trade-off point of view, the first normal-stress difference has been measured to be quite small, yet the introduction of a tensorial order parameter would increase considerably the complexity of the governing equations.

6.2 Pressure-driven channel flow

This test case is a convenient benchmark for suspension models because it is relatively simple to describe, yet the presence of a varying velocity gradient allows the observation of particle migration due to non-Newtonian stresses. Over the years, it has been the subject of numerous experimental studies, ranging from early works such as Karnis, Goldsmith & Mason (Reference Karnis, Goldsmith and Mason1966) to the well-known studies of Koh, Hookham & Leal (Reference Koh, Hookham and Leal1994), Hampton et al. (Reference Hampton, Mammoli, Graham, Tetlow and Altobelli1997) and Lyon & Leal (Reference Lyon and Leal1998). More recent experimental results have been obtained by the group of Gilchrist, cf. Gao, Xu & Gilchrist (Reference Gao, Xu and Gilchrist2009) and Gao et al. (Reference Gao, Kulkarni, Morris and Gilchrist2010). Experimental measurements of the group of Gilchrist and numerical results obtained via a discrete-element method have recently been reported in Yeo & Maxey (Reference Yeo and Maxey2011). Herein we compare the predictions of our model against those results.

In the experiment that we chose for our comparisons, a suspension of rigid silica microspheres in a glycerine–water solution is pumped through a rectangular channel whose width equals $40~\unicode[STIX]{x03BC}\text{m}$ . The particle diameter is $1.0~\unicode[STIX]{x03BC}\text{m}$ and the initial average concentration is set at $\langle \unicode[STIX]{x1D719}_{s}\rangle =0.411$ . Due to the small channel size and the high particle concentration, the effect of gravity is negligible and, therefore, the flow is symmetric around the centreline. The experimental data are taken when the flow has reached steady state and is fully developed. For this reason, we look explicitly for a one-dimensional and time-independent numerical solution of equations (5.14)–(5.18). In particular, we compute steady-state solutions of these equations via an iterative, under-relaxation numerical method.

In our computations, the reference values that we use for non-dimensionalisation purposes are defined as follows. The reference length is equal to the channel’s half-width, $l_{r}=20~\unicode[STIX]{x03BC}\text{m}$ . The reference velocity is the one at the centreline, $u_{r}=2\times 10^{-4}~\text{m}~\text{s}^{-1}$ . The reference density and viscosity are the ones of the liquid, $\unicode[STIX]{x1D70C}_{r}=1200~\text{kg}~\text{m}^{-3}$ and $\unicode[STIX]{x1D707}_{r}=0.04~\text{Pa}~\text{s}$ . Finally, the reference pressure is $l_{r}u_{r}^{2}$ .

As regards the material parameters, the drag coefficient $\unicode[STIX]{x1D6FF}$ is approximated by the correlation proposed in Benyahia, Syamlal & O’Brien (Reference Benyahia, Syamlal and O’Brien2005). In this particular case, the drag coefficient is so high that the fluid and solid phasial velocities practically coincide. For the viscous stresses, we remark that in this particular case the bulk stresses are zero because the flow is divergence free and the shear viscosity coefficient does not depend on $\unicode[STIX]{x1D651}_{s}$ because $\det \unicode[STIX]{x1D651}_{s}=0$ , cf. (3.39). This allows us to use the following expressions proposed in Morris & Boulay (Reference Morris and Boulay1999),

(6.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}_{s}=\frac{c_{4}}{2\unicode[STIX]{x1D719}_{s}}=\frac{2.5}{\unicode[STIX]{x1D719}_{m}-\unicode[STIX]{x1D719}_{s}}+0.1\,\frac{\unicode[STIX]{x1D719}_{s}}{(\unicode[STIX]{x1D719}_{m}-\unicode[STIX]{x1D719}_{s})^{2}}, & \displaystyle\end{eqnarray}$$

(6.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}_{n}=\frac{1}{\sqrt{2}\,\unicode[STIX]{x1D719}_{s}}(c_{5}-c_{2})=0.75\,\frac{\unicode[STIX]{x1D719}_{s}}{(\unicode[STIX]{x1D719}_{m}-\unicode[STIX]{x1D719}_{s})^{2}}. & \displaystyle\end{eqnarray}$$

In the above relations $\unicode[STIX]{x1D719}_{m}$ represents the granular volume fraction at the jamming transition, $\unicode[STIX]{x1D719}_{m}=0.6$ . The values of $\unicode[STIX]{x1D707}_{s}$ and $\unicode[STIX]{x1D707}_{n}$ , calculated at the average value $\langle \unicode[STIX]{x1D719}_{s}\rangle =0.4$ , are $14.4$ and $8.6$ , respectively. In other words, these viscosities are an order of magnitude larger than the viscosity of the liquid. The remaining coefficients $\unicode[STIX]{x1D712}_{a}$ and $\unicode[STIX]{x1D712}_{b}$ are evaluated by combining the above expression for $\unicode[STIX]{x1D707}_{n}$ with (6.10) and (3.41).

For $\unicode[STIX]{x1D6FD}_{s}$ we consider the following relation, which is a simplified version of the expression proposed in Powers (Reference Powers2004),

(6.13) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{s}=-\unicode[STIX]{x1D6FD}_{0}\ln (\unicode[STIX]{x1D719}_{m}-\unicode[STIX]{x1D719}_{s}), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}_{0}$ is a constant. Finally, as regards $\unicode[STIX]{x1D6FE}_{s}$ , we note that reliable experimental data are currently unavailable. For this reason, in our study we tested various power laws of $\unicode[STIX]{x1D719}_{s}$ . Our tests showed that high exponents gave better agreement with experiments. In particular, the following expression

(6.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}_{s}=-\unicode[STIX]{x1D6FE}_{0}\,\unicode[STIX]{x1D719}_{s}^{6}, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D6FE}_{0}$ being a constant, turned out to work well in practice. The non-dimensional values of $\unicode[STIX]{x1D6FD}_{0}$ and $\unicode[STIX]{x1D6FE}_{0}$ , are set as follows,

(6.15a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{0}=1.04\times 10^{3},\quad \unicode[STIX]{x1D6FE}_{0}=5.2\times 10^{4}. & & \displaystyle\end{eqnarray}$$

(The dimensional values of the above parameters are $5\times 10^{-2}~\text{Pa}$ and $1\times 10^{-9}~\text{N}$ , respectively. This implies, for example, that $\unicode[STIX]{x1D6FD}_{0}$ is approximately 35 times smaller than the pressure gradient along the channel multiplied by the channel half-width.) Finally, the Reynolds number of the problem is equal to $Re=1.2\times 10^{-4}$ .

Figure 1.

Pressure-driven flow: profiles of the granular volume fraction $\unicode[STIX]{x1D719}_{s}$ and velocity $u_{s}$ . The predictions of our model are plotted against the experimental results and the discrete-element simulations reported in Yeo & Maxey (Reference Yeo and Maxey2011) (denoted by Y&M in the figure legends). Average volume fraction, $\langle \unicode[STIX]{x1D719}_{s}\rangle =0.4$ .

Next, we proceed to the discussion of our results. Figure 1(a) shows numerical predictions and experimental data for the volume fraction $\unicode[STIX]{x1D719}_{s}$ . The observed particle migration towards the centre of the channel can be explained as follows. In the initial stages of the flow, the granular volume fraction is uniform. Then, the high velocity gradients in the wall regions translate into higher shear-induced normal stresses in these regions. In turn, these stresses drive the particles towards the centreline, resulting in the particle migration that is observed in the fully developed flow. As regards our model, this migration is predicted by virtue of the coupling between the momentum balance (5.17) and the compaction equation (5.18).

We can infer from figure 1(a) that the predictions of our model compare quite well with the experimental data. Nonetheless, there is a discrepancy inside a thin layer (2 to 3 grain diameters) next to the channel walls. This discrepancy may be attributed to the fact that continuum theories such as ours cannot capture flow structures that are present only at the grain scale.

Figure 1(b) shows the numerical predictions for the velocity profile. We can observe the blunting of this profile near the centreline with respect to the parabolic solution. This is due to the dependence of the viscosity coefficients on the volume fraction $\unicode[STIX]{x1D719}_{s}$ . In particular, as $\unicode[STIX]{x1D719}_{s}$ increases, so does the viscosity coefficient $\unicode[STIX]{x1D707}_{s}$ , cf. equation (6.11). For the high value of the average volume fraction $\langle \unicode[STIX]{x1D719}_{s}\rangle$ of this example, a central plug zone is predicted which exhibits practically constant velocity. As can be directly inferred from figure 1(b), the profile predicted by our model agrees quite well with the one predicted by the computations of Yeo & Maxey (Reference Yeo and Maxey2011).