2.1. Boltzmann equation for charged particles with bidisperse size distribution

Starting from the Boltzmann equation with the probability density function, we can describe the statistical behaviour of a binary mixture of particles in the dense regime. We denote the probability density function of particles by $\,f_{pi}( \boldsymbol {x}, {\boldsymbol {c}_{pi}}, q_{pi}, t)$ at position $\boldsymbol {x}$ with velocity $\boldsymbol {c}_{pi}$ and charge $q_{pi}$ on particles for the discrete phase $i$. The number of particles in the phase $i$ with velocity between $\boldsymbol {c}_{pi}$ and $\boldsymbol {c}_{pi}+\textrm {d}\boldsymbol {c}_{ip}$ and charge between $q_{pi}$ and $q_{pi} + \textrm {d}q_{pi}$ at position $\boldsymbol {x}$ and time $t$ is then given by $\,f_{pi}\, \textrm {d} \boldsymbol {c}_{pi} \,\textrm {d}q_{pi}$. The evolution of the probability density function follows the Boltzmann equation

(2.1)\begin{equation} \frac{\partial \,f_{pi}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\boldsymbol{c}_{pi}\,f_{pi}\right) + \frac{\partial }{\partial \boldsymbol{c}_{pi}}\left(\left.\left\langle\frac{\textrm{d}\boldsymbol{c}_{pi}}{\textrm{d}t}\right\rangle\right|_{\boldsymbol{c}_{pi}, q_{pi}}\,f_{pi}\right) + \frac{\partial }{\partial q_{pi}}\left(\left.\left\langle \frac{\textrm{d}q_{pi}}{\textrm{d}t}\right\rangle\right|_{\boldsymbol{c}_{pi}, q_{pi}}\,f_{pi}\right) = \left(\frac{\partial \,f_{pi}}{\partial t}\right)_{coll}. \end{equation}

The time derivative terms in $\langle \cdot \rangle$ describe the rate of change of particle velocity and charge in the Lagrangian frame. The term on the right-hand side is the rate of change of the probability density function with particle–particle collisions.

2.2. Discrete particle equations

The rate of change of particle velocity is defined by the equation of motion as

(2.2)\begin{equation} m_{pi}\frac{\textrm{d}\boldsymbol{c}_{pi}}{\textrm{d}t}= \boldsymbol{F}_{ei}. \end{equation}

Here, $m_{pi}$ is the mass of a particle in the discrete phase i. The external forces acting on the discrete phase $i$ such as gravitational and fluid–solid interaction forces (e.g. drag and Archimedes forces) are not accounted and only the electrostatic force acting on a particle is accounted for as follows:

(2.3)\begin{equation} \boldsymbol{F}_{ei} = q_{pi}\,\boldsymbol{E}, \end{equation}

where $\boldsymbol {E}$ is the resolved electric field (higher-order terms due to polarization (Kolehmainen et al. Reference Kolehmainen, Ozel, Gu, Shinbrot and Sundaresan2018a) and magnetic forces (Genc & Derin Reference Genc and Derin2014) were neglected). The resolved electric field is computed by solving a Poisson's equation

(2.4)\begin{equation} \nabla^2 \phi ={-} \frac{\rho_q}{\epsilon}, \end{equation}

for the electrical potential, $\phi$, where $\rho _q$ is the charge density and $\epsilon$ is the electrical permittivity. Then, the resolved electric field is obtained by taking the gradient of the electrical potential

(2.5)\begin{equation} \boldsymbol{E}={-}\boldsymbol{\nabla} \phi. \end{equation}

The charge transfer occurs only by collision, therefore,

(2.6)\begin{equation} \frac{\textrm{d}q_{pi}}{\textrm{d}t} = 0. \end{equation}

2.3. Moment equations

Any macroscopic property of the discrete phase $i$ is defined using the probability density function and averaging properties over a range of velocity and charge is given as follows:

(2.7)\begin{equation} \langle \psi_{pi} \rangle = \frac{1}{n_{pi}}\int_{\mathbb{R}} \int_{\mathbb{R}^3} \psi_{pi}\,f_{pi}\,\textrm{d}\boldsymbol{c}_{pi}\,\textrm{d}q_{pi}, \end{equation}

where $n_{pi}$ is the number density of the phase $i$ particles. For each phase, mean velocity, $\boldsymbol {U}_{pi}$, granular temperature, $\varTheta _{p{i}}$, mean charge, $Q_{pi}$, and charge variance, $\mathcal {Q}_{pi}$, are then defined as

(2.8)$$\begin{gather} \boldsymbol{U}_{pi} = \frac{1}{n_{pi}}\int_{\mathbb{R}} \int_{\mathbb{R}^3} \boldsymbol{c}_{pi}\,f_{pi}\,\textrm{d}\boldsymbol{c}_{pi}\,\textrm{d}q_{pi}, \end{gather}$$

(2.9)$$\begin{gather}\varTheta_{pi} = \frac{m_{pi}}{3n_{pi}}\int_{\mathbb{R}} \int_{\mathbb{R}^3} (\boldsymbol{c}_{pi}'\boldsymbol{\cdot}\boldsymbol{c}_{pi}')\,f_{pi}\,\textrm{d}\boldsymbol{c}_{pi}\,\textrm{d}q_{pi}, \end{gather}$$

(2.10)$$\begin{gather}Q_{pi} = \frac{1}{n_{pi}}\int_{\mathbb{R}} \int_{\mathbb{R}^3} q_{pi}\,f_{pi}\,\textrm{d}\boldsymbol{c}_{pi}\,\textrm{d}q_{pi}, \end{gather}$$

(2.11)$$\begin{gather}\mathcal{Q}_{pi} = \frac{m_{pi}}{n_{pi}}\int_{\mathbb{R}} \int_{\mathbb{R}^3} q_{pi}'q_{pi}'\,f_{pi}\,\textrm{d}\boldsymbol{c}_{pi}\,\textrm{d}q_{pi}, \end{gather}$$

where $\boldsymbol {c}_{pi}'$ is the fluctuating phase velocity and $q_{pi}'$ is the fluctuating phase charge. The fluctuating phase velocity is defined as the difference between phase and mean velocities; $\boldsymbol {c}_{pi}' = \boldsymbol {c}_{pi} - \boldsymbol {U}_{pi}$. Averaging the Boltzmann equation (2.1) over a range of velocities and charges and using the relation $n_{pi}m_{pi} = \alpha _{pi}\rho _{pi}$, the Enskog equation is obtained

(2.12)\begin{equation} \frac{\partial}{\partial t}\left(\alpha_{pi}\rho_{pi}\langle\psi_{pi}\rangle\right) + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\alpha_{pi}\rho_{pi}\langle\boldsymbol{c}_{pi}\psi_{pi}\rangle\right) = \mathcal{C}(m_{pi}\psi_{pi}) + \alpha_{pi}\rho_{pi}\left\langle\frac{\textrm{d}\boldsymbol{c}_{pi}}{\textrm{d}t}\frac{\partial \psi_{pi}}{\partial \boldsymbol{c}_{pi}}\right\rangle, \end{equation}

where $\alpha _{pi}$ is the solid volume fraction, $\rho _{pi}$ is the density in the discrete phase $i$. The two terms on the left-hand side represent the transport of a quantity $\psi _{pi}$, the first term on the right-hand side represents the rate of change of the quantity averaging over collisions and the last term represents the external force (herein, it is the electrostatic force) acting on the particles. To close the system, the collisional operator, $\mathcal {C}(m_{pi}\psi _{pi})$, needs to be modelled. For a binary mixture of particles, the rate of change of a property due to collisions can be decomposed into the flux and source terms by following Jenkins & Mancini (Reference Jenkins and Mancini1987)

(2.13)\begin{align} \mathcal{C}(m_{pi}\psi_{pi}) & = \sum_{h = i,j}d_{pih}^2\int_{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w} > 0} m_{pi}(\psi_{pi}^+{-} \psi_{pi})|\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w}|\nonumber\\ & \quad \times f^*_{pih}(\boldsymbol{c}_{pi}, \boldsymbol{x}, \boldsymbol{c}_{ph}, \boldsymbol{x}+d_{pih}\boldsymbol{k})d\boldsymbol{k}\,\textrm{d}\boldsymbol{c}_{pi}\,\textrm{d}\boldsymbol{c}_{ph}\,\textrm{d}q_{pi}\,\textrm{d}q_{ph} \end{align}

(2.14)\begin{align} & = \sum_{h = i,j}\left(- \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\theta}_{ih}(m_{pi}\psi_{pi}) + \chi_{ih}(m_{pi}\psi_{pi})\right). \end{align}

The flux term, $\boldsymbol {\theta }_{ih}$, represents the redistribution of a quantity within and between phases while the source term, $\chi _{ih}$, represents the dissipation of the quantity $\psi _{pi}$ between the phases $i$ and $h$. These terms are derived by using the following integrals:

(2.15)\begin{align} \boldsymbol{\theta}_{ih}(m_{pi}\psi_{pi}) & ={-}\frac{d_{pih}^3}{2}\int_{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w} > 0} m_{pi}(\psi_{pi}^+{-} \psi_{pi})|\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w}|\boldsymbol{k}\nonumber\\ & \quad \times f^*_{pih}(\boldsymbol{c}_{pi}, \boldsymbol{x} - \tfrac{1}{2}d_{pih}\boldsymbol{k}, \boldsymbol{c}_{ph}, \boldsymbol{x}+\tfrac{1}{2}d_{pih}\boldsymbol{k})\,\textrm{d}\boldsymbol{k}\,\textrm{d}\boldsymbol{c}_{pi}\,\textrm{d}\boldsymbol{c}_{ph}\,\textrm{d}q_{pi}\,\textrm{d}q_{ph}, \end{align}

(2.16)\begin{align} \chi_{ih}(m_{pi}\psi_{pi}) & = d_{pih}^2\int_{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w} > 0} m_{pi}(\psi_{pi}^+{-} \psi_{pi})|\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w}| \nonumber\\ & \quad \times f^*_{pih}(\boldsymbol{c}_{pi}, \boldsymbol{x} - \tfrac{1}{2}d_{pih}\boldsymbol{k}, \boldsymbol{c}_{ph}, \boldsymbol{x} + \tfrac{1}{2}d_{pih}\boldsymbol{k})\,\textrm{d}\boldsymbol{k}\,\textrm{d}\boldsymbol{c}_{pi}\,\textrm{d}\boldsymbol{c}_{ph}\,\textrm{d}q_{pi}\,\textrm{d}q_{ph}. \end{align}

In these integrals, $\psi _{pi}^+$ refers to a quantity after collision, $\boldsymbol {k}$ is the unit vector from the centre of two particles at contact, $\boldsymbol {w}$ is the relative velocity between two particles and $d_{pih}$ is the mean diameter defined as $(d_{pi} + d_{ph})/2$ where $d_{pi}$ and $d_{ph}$ are diameters of phase $i$ and $h$, respectively. The symbol $f^*_{pih}$ refers to the joint pair distribution function for phases $i$ and $h$ at the contact point. With an assumption of random motion of particles, it is approximated with a Taylor's expansion at contact point as

(2.17)\begin{equation} f^{*}_{pih} = g_0\,f_{pi}\,f_{ph}\left(1 + \frac{d_{pih}}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\nabla}\ln\left( \frac{f_{ph}}{f_{pi}}\right)\right), \end{equation}

where $g_0$ is the radial distribution. To compute integrals, we presume that both charge and velocity distributions follow a Maxwellian distribution. The probability density function for the discrete phase $i$ is then defined as

(2.18)$$\begin{gather} f_{pi}(\boldsymbol{c}_{pi},q_{pi},\boldsymbol{x},t) = n_{pi} \underbrace{\left(\frac{m_{pi}}{2{\rm \pi}\varTheta_{p{i}}}\right)^{3/2} \exp\left(-\frac{m_{pi}}{2\varTheta_{p{i}}}(\boldsymbol{c}_{pi} - \boldsymbol{U}_{pi})\boldsymbol{\cdot}(\boldsymbol{c}_{pi} - \boldsymbol{U}_{pi})\right)}_{f_{pi,c}}\nonumber\\ \underbrace{\left(\frac{m_{pi}}{2{\rm \pi}\mathcal{Q}_{pi}}\right)^{1/2} \exp\left(-\frac{m_{pi}}{2\mathcal{Q}_{pi}}(q_{pi} - Q_{pi})^2\right)}_{f_{pi,q}}. \end{gather}$$

We aim to develop the constitutive closures and the charge transport equation for collisional granular flows with bidisperse charged particles where charge transfer occurs mainly via particle–particle collisions. Therefore, the correlation between charge and velocity has been neglected in the probability density function. For the dilute regime, this assumption is invalid and an additional term arising from self-diffusion of charge should be modelled. Readers are referred to a rigorous theoretical development achieved by a recent study of Montilla et al. (Reference Montilla, Ansart and Simonin2020) on the modelling of the charge-velocity correlation for monosize particles.

2.4. Revisiting hydrodynamic equations for bidisperse granular flows

Before presenting the charge transport equation, we revisit the hydrodynamic equations for the granular flows with bidisperse size distribution. The transport equations for mass ($\psi _{pi} = 1$), momentum ($\psi _{pi} = \boldsymbol {c}_{pi}$) and granular temperature ($\psi _{pi} = (\boldsymbol {c}_{pi}'\boldsymbol {\cdot }\boldsymbol {c}_{pi}')/2 = 3(\varTheta _{pi}/m_{pi})/2$) are derived from the Enskog equation (2.12). The closure relations for the collision terms are derived by following Iddir & Arastoopour (Reference Iddir and Arastoopour2005). However, there are slight differences in the derived constitutive equations discussed below. If there is no exchange of mass or breaking of particles during collisions, the mass balance for the phase $i$ is written as

(2.19)\begin{equation} \frac{\partial}{\partial t}\left(\alpha_{pi}\rho_{pi}\right) + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\alpha_{pi}\rho_{pi}\boldsymbol{U}_{pi}\right) = 0. \end{equation}

The momentum balance for the phase $i$ is written as

(2.20)$$\begin{gather} \alpha_{pi}\rho_{pi}\left[\frac{\partial}{\partial t} + \boldsymbol{U}_{pi}\boldsymbol{\cdot}\boldsymbol{\nabla}\right]\boldsymbol{U}_{pi} = \sum_{h = i,j}\left(- \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\theta}_{ih} + \boldsymbol{\chi}_{ih}\right) - \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\alpha_{pi}\rho_{pi}\langle\boldsymbol{c}_{pi}'\boldsymbol{c}_{pi}'\rangle\right) \nonumber\\ + \frac{\alpha_{pi}\rho_{pi}}{m_{pi}}Q_{pi}\boldsymbol{E}. \end{gather}$$

Here, the electric field, $\boldsymbol {E}$, is computed with (2.4) and (2.5). The first two terms on the right-hand side represent the rate of change of momentum due to collisions and redistribution due to the random velocity fluctuations, respectively. The flux term for the collisional operator is defined as

(2.21)$$\begin{gather} \boldsymbol{\theta}_{ih} = n_{pi}n_{ph}\frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})}\left(\frac{m_{pi}m_{ph}}{\varTheta_{p{i}}\varTheta_{p{h}}}\right)^{3/2}(1+e_c)g_0\frac{d_{pih}^3}{48}\left[{\rm \pi} M_1\boldsymbol{\mathsf{I}} - \frac{2d_{pih}}{5}\sqrt{\rm \pi}M_2\right.\nonumber\\ \left.\times \frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})} \sum_{l=i,h}\left(\frac{1}{\varTheta_{p{l}}}\left[(\boldsymbol{\nabla}\boldsymbol{U}_{pl})^s + \frac{5}{6} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pl}\boldsymbol{\mathsf{I}}\right] \right)\right], \end{gather}$$

with

(2.22)\begin{equation} (\boldsymbol{\nabla}\boldsymbol{U}_{pl})^s = \tfrac{1}{2}\left((\boldsymbol{\nabla}\boldsymbol{U}_{pl}) + (\boldsymbol{\nabla}\boldsymbol{U}_{pl})^{\textrm{T}}\right) -\tfrac{1}{3} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pl}\boldsymbol{\mathsf{I}}. \end{equation}

The source term is given by

(2.23) \begin{align} & \boldsymbol{\chi}_{ih} ={-} n_{pi}n_{ph}\frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})}\left(\frac{m_{pi}m_{ph}}{\varTheta_{p{i}}\varTheta_{p{h}}}\right)^{3/2}(1+e_c)g_0 \frac{d_{pih}^2}{6}\left[\vphantom{\frac{n_{ph}}{n_{pi}}}\sqrt{\frac{\rm \pi}{4}}(\boldsymbol{U}_{pi} - \boldsymbol{U}_{ph})M_3 \right.\nonumber\\ &\quad + \left. d_{pih} \frac{\rm \pi}{8}\left[M_1\left(\boldsymbol{\nabla}\ln\frac{n_{ph}}{n_{pi}} - \frac{3}{2}\boldsymbol{\nabla}\ln\frac{\varTheta_{p{h}}}{\varTheta_{p{i}}}\right) + \frac{1}{4} \left(3M_4\left(\frac{m_{ph}\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} - \frac{m_{pi}\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right)\right.\right.\right.\nonumber\\ &\quad +\left.\left.\left. 5M_5\frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})^2}\left(\frac{m_{pi}\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} - \frac{m_{ph}\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right)\right.\right.\right.\nonumber\\ &\quad +\left.\left.\left. 10B M_6\frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})}\left(\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} + \frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right)\right) \right] \right], \end{align}

where $e_c$ is the restitution coefficient with $e_c$ = 1 for fully elastic collisions. The coefficients $M_k$ ($k=1,\ldots ,6$) and $B$ are given in table 1. The derivations of these terms are not given here but the reader is referred to Iddir & Arastoopour (Reference Iddir and Arastoopour2005) for further details.

Table 1. Model coefficients of flux and source terms for the phase momentum, granular temperature and charge transport equations. The model coefficients, $M_k$ ($k=1,\ldots ,6$) and $B$, are used in (2.21) and (2.23). The model coefficients, $M_k$ ($k=1,\ldots ,14$) and $B$ are used in (2.25) and (2.26). The model coefficients, $N_k$ ($k=1,\ldots ,5$) and $B$, are used in (2.36)–(2.38), (2.40) and (2.41). The symbol $\varGamma (.)$ refers to the gamma function.

The transport equation for granular temperature of the solid phase $i$ is given by

(2.24)\begin{equation} \frac{3}{2} \frac{\partial}{\partial t}\left(\alpha_{pi}\rho_{pi}\frac{\varTheta_{p{i}}}{m_{pi}}\right) + \frac{3}{2} \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\alpha_{pi}\rho_{pi}\boldsymbol{U}_{pi}\frac{\varTheta_{p{i}}}{m_{pi}}\right) = \sum_{h = i,j}\left( -\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{q}_{ih} + \gamma_{ih}\right), \end{equation}

with the flux, $\boldsymbol {q}_{ih}$, and source, $\gamma _{ih}$, terms defined as

(2.25)\begin{align} &\boldsymbol{q}_{ih} ={-} n_{pi}n_{ph}\frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})}\left(\frac{m_{pi}m_{ph}}{\varTheta_{p{i}}\varTheta_{p{h}}}\right)^{3/2}(1+e_c)g_0d_{pih}^3\left( \frac{d_{pih}}{48}\sqrt{\rm \pi}\left[\left(\boldsymbol{\nabla}\ln\left[ \frac{n_{ph}}{n_{pi}}\right]\right.\right.\right.\nonumber\\ &\quad \left. + \frac{3}{2}\boldsymbol{\nabla}\ln\left[ \frac{\varTheta_{p{i}}}{\varTheta_{p{h}}}\right]\right) BM_7 + \frac{5}{4}\left(m_{ph}\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} - m_{pi}\frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right)BM_8 + \frac{3m_{pi}m_{ph}}{2(m_{pi} + m_{ph})^2} BM_9\nonumber\\ & \quad\left. \times \left(m_{pi}\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} - m_{ph}\frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right) + \frac{m_{pi}m_{ph}}{2(m_{pi} + m_{ph})}\left(\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} + \frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right) M_{10}\right] \nonumber\\ & \quad + \frac{(1-e_c)}{64}\frac{m_{ph}}{(m_{pi} + m_{ph})} \left[\frac{2{\rm \pi}}{3}(\boldsymbol{U}_{pi} - \boldsymbol{U}_{ph})M_1 + \sqrt{\rm \pi}d_{pih}\left[\frac{2}{3}\left(\boldsymbol{\nabla}\ln\left[ \frac{n_{ph}}{n_{pi}}\right] + \frac{3}{2}\boldsymbol{\nabla}\ln\left[ \frac{\varTheta_{p{i}}}{\varTheta_{p{h}}}\right]\right)M_2\right.\right.\nonumber\\ & \quad + \frac{1}{2}\left(m_{ph}\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} - m_{pi}\frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right)M_{11} + \frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})^2}\left(m_{pi}\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} - m_{ph}\frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right)M_{12}\nonumber\\ &\quad \left.\left.\left. + 2B \frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})}\left(\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} + \frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right)M_{13}\right]\right]\right), \end{align}

and

(2.26)\begin{align} & \gamma_{ih} = n_{pi}n_{ph}\frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})}\left(\frac{m_{pi}m_{ph}}{\varTheta_{p{i}}\varTheta_{p{h}}}\right)^{3/2}(1+e_c)g_0d_{pih}^2\left[\frac{\sqrt{\rm \pi}}{4}BM_7 - \frac{{\rm \pi} d_{pih}}{160} \right.\nonumber\\ & \quad \times\left[\left(m_{ph}\frac{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{ph}}{\varTheta_{p{h}}} - m_{pi}\frac{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pi}}{\varTheta_{p{i}}}\right)M_{14} + 5B\frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})}\left(\frac{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{ph}}{\varTheta_{p{h}}} + \frac{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pi}}{\varTheta_{p{i}}}\right)M_6\right]\nonumber\\ & \quad- \frac{1}{8}\frac{m_{ph}}{(m_{pi} + m_{ph})}(1-e_c)\left[\sqrt{\rm \pi}M_2 - \frac{{\rm \pi} d_{pih}}{8}\left(m_{ph}\frac{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{ph}}{\varTheta_{p{h}}} - m_{pi}\frac{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pi}}{\varTheta_{p{i}}}\right)BM_6\right.\nonumber\\ &\left.\left. \quad - \frac{{\rm \pi} d_{pih}}{8} \frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})}\left(\frac{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{ph}}{\varTheta_{p{h}}} + \frac{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pi}}{\varTheta_{p{i}}}\right)M_5\right]\right]. \end{align}

The derived models are very similar to the ones proposed by Iddir & Arastoopour (Reference Iddir and Arastoopour2005) but there are differences in the high-order terms of the model coefficients (see table 1). The differences might result from the Taylor expansion of the relative velocity and the centre of mass velocity multiplication (see (A13)) for integral of the collision operator. Our approximation is detailed in Appendix A. Unfortunately, Iddir & Arastoopour (Reference Iddir and Arastoopour2005) did not give an explicit explanation about how they treated this term. The assessment benchmark of these revisited constitutive equations through hard-sphere simulation results is given in § 3.

2.5. Transport equation for phase mean charge

In this section, we present the transport equation for mean charge for each phase. Assuming that charge and velocity distributions are uncorrelated and using (2.12), the transport equation for the mean charge is given by

(2.27)\begin{equation} \frac{\partial}{\partial t}\left(\frac{\alpha_{pi}\rho_{pi}}{m_{pi}}Q_{pi}\right) + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\alpha_{pi}\rho_{pi}}{m_{pi}}Q_{pi}\boldsymbol{U}_{pi}\right) = \mathcal{C}(q_{pi}). \end{equation}

The charge transfer during collision between particles is based on the model proposed by Laurentie et al. (Reference Laurentie, Traoré and Dascalescu2013) with the phase effective work function, $\varphi _{ph=i,j}$, and the electric field, $\boldsymbol {E}$. The transfer charge between particles $l$ and $m$ within the same phase (i.e. the phase $i$) is given by

(2.28)$$\begin{gather} q_{pi}^{(l)+} = q_{pi}^{(l)} + dq = q_{pi}^{(l)} + \mathcal{A}_{max}(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w}) \epsilon_0\left(- \boldsymbol{E}\boldsymbol{\cdot}\boldsymbol{k} + \frac{q_{pi}^{(m)} - q_{pi}^{(l)}}{{\rm \pi}\epsilon_0d_{pi}^2}\right), \end{gather}$$

(2.29)$$\begin{gather}q_{pi}^{(m)+} = q_{pi}^{(m)} - dq = q_{pi}^{(m)} - \mathcal{A}_{max}(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w})\epsilon_0\left(- \boldsymbol{E}\boldsymbol{\cdot}\boldsymbol{k} + \frac{q_{pi}^{(m)} - q_{pi}^{(l)}}{{\rm \pi}\epsilon_0d_{pi}^2}\right). \end{gather}$$

For between different phases;

(2.30)$$\begin{gather} q_{pi}^{(l)+} = q_{pi}^{(l)} + dq = q_{pi}^{(l)} + \mathcal{A}_{max}(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w}) \epsilon_0\left(\frac{\varphi_{pi} - \varphi_{pj}}{\delta_c e} - \boldsymbol{E}\boldsymbol{\cdot}\boldsymbol{k} + \frac{1}{{\rm \pi}\epsilon_0}\left(\frac{q_{pj}^{(m)}}{d_{pj}^2} - \frac{q_{pi}^{(l)}}{d_{pi}^2}\right)\right), \end{gather}$$

(2.31)$$\begin{gather}q_{pj}^{(m)+} = q_{pj}^{(m)} - dq = q_{pj}^{(m)} - \mathcal{A}_{max}(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w})\epsilon_0\left(\frac{\varphi_{pi} - \varphi_{pj}}{\delta_c e} - \boldsymbol{E}\boldsymbol{\cdot}\boldsymbol{k} + \frac{1}{{\rm \pi}\epsilon_0}\left(\frac{q_{pj}^{(m)}}{d_{pj}^2} - \frac{q_{pi}^{(l)}}{d_{pi}^2}\right)\right). \end{gather}$$

In (2.30) and (2.31), $\delta _c$ is the cutoff distance of electron transfer, $e$ is the elementary charge and $\epsilon _0$ is the electrical permittivity in a vacuum; $\mathcal {A}_{max}$ is the maximum overlapping area computed with help of the contact Hertz theory as (Kolehmainen et al. Reference Kolehmainen, Ozel, Boyce and Sundaresan2017a)

(2.32)\begin{equation} \mathcal{A}_{max} = 2{\rm \pi} r_p^*\left(\frac{15m_p^*}{16Y_p^*\sqrt{r_p^*}}\right)^{2/5}|\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w}|^{4/5} = \mathcal{A}^*|\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{w}|^{4/5}. \end{equation}

The maximum overlapping area is approximated by a collision of two elastic particles that follows a conversion of the particle kinetic energy into the potential energy of a Hertzian spring (it is assumed that the electric potential energy is negligible as compared with the potential energy in the spring). In (2.32), $\mathcal {A}^*$ is the effective area which is only a function of particle physical properties such as the effective Young modulus, $Y_p^*$, the effective mass, $m_p^*$, and the effective radius, $r_p^*$, that are defined as

(2.33a–c)\begin{equation} \frac{1}{Y_p^*} = \frac{1 - \nu_{pi}^2}{Y_{pi}} + \frac{1 - \nu_{pj}^2}{Y_{pj}},\quad \frac{1}{r_p^*} = \frac{1}{r_{pi}} + \frac{1}{r_{pj}},\quad \frac{1}{m_p^*} = \frac{1}{m_{pi}} + \frac{1}{m_{pj}}. \end{equation}

The maximum contact area could be more accurately computed for granular material with viscoelastic particles by following Schwager & Pöschel (Reference Schwager and Pöschel2008) or Brilliantov & Pöschel (Reference Brilliantov and Pöschel2010). However, due to the complex nature of contact electrification, an agreement on the charge transfer model in granular material has not yet been reached, as discussed by Matsusaka et al. (Reference Matsusaka, Maruyama, Matsuyama and Ghadiri2010b). The contact electrification also depends on many other particle properties such as shape, roughness, surface functionalization etc. Therefore, we preferred to use the elastic spheres in the phenomenological charge transfer model.

The closure of the collisional operator in (2.27) is defined as follows:

(2.34)\begin{equation} \mathcal{C}(q_{pi}) = \sum_{h = i,j}\left(- \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\theta}_{ih}^q(q_{pi}) + \chi_{ih}^q(q_{pi})\right), \end{equation}

where the flux term, $\boldsymbol {\theta }_{ih}^q(q_{pi})$, represents the spatial redistribution of charge and the source term, $\chi _{ih}^q(q_{pi})$, represents the charge transfer between phases. The derivations of flux and source terms for the phase charge transport equation are discussed in Appendix A and their final forms are given in (A31) and (A32), respectively. Here, we present these equations in a compact form by following Kolehmainen et al. (Reference Kolehmainen, Ozel and Sundaresan2018b). The flux term, $\boldsymbol {\theta }_{ih}^q$, is then written as

(2.35)\begin{align} \boldsymbol{\theta}_{ih}^q ={-} \boldsymbol{\sigma}^{\theta}_{ih}\boldsymbol{\cdot}\boldsymbol{E} - \kappa^{\theta}_{ih}\left( \frac{\boldsymbol{\nabla} Q_{ph}}{d_{ph}^2} + \frac{\boldsymbol{\nabla} Q_{pi}}{d_{pi}^2} \right) - \boldsymbol{D}^{\theta}_{ih}\left(\frac{\varphi_{pi} - \varphi_{ph}}{\delta_c e} + \frac{1}{{\rm \pi}\epsilon_0}\left(\frac{Q_{ph}}{d_{ph}^2} - \frac{Q_{pi}}{d_{pi}^2}\right)\right), \end{align}

with the triboelectric conductivity tensor, $\boldsymbol {\sigma }^{\theta }_{ih}$, the triboelectric diffusivity, $\kappa ^{\theta }_{ih}$, and the triboelectric phase coupling coefficient, $\boldsymbol {D}^{\theta }_{ih}$. These coefficients are defined below. The first term on the-right-hand side in (2.35) represents a current density due the electric field resulting from the charge on the particles. The second term results from the dispersion of charge while the third term arises due to charge difference between particles during a collision. The triboelectric phase coupling coefficient appears due to the non-equipartition of granular temperature (see (2.38)). These terms are defined in explicit forms as

(2.36)\begin{align} \boldsymbol{\sigma}^{\theta}_{ih} & = n_{pi}n_{ph}\left(\frac{m_{pi}m_{ph}}{\varTheta_{p{i}}\varTheta_{p{h}}}\right)^{3/2}\mathcal{A}^*\epsilon_0g_0\frac{d_{pih}^3}{8}\sqrt{\rm \pi}\left[ -\frac{5}{21}N_1\boldsymbol{\mathsf{I}} + \frac{3}{1102}d_{pih}\frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})}N_5 \right.\nonumber\\ &\left. \quad \times\left[\sum_{l = i,j}\frac{1}{\varTheta_{p{l}}}\left( (\boldsymbol{\nabla}\boldsymbol{U}_{pl}) + (\boldsymbol{\nabla}\boldsymbol{U}_{pl})^{\textrm{T}} + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pl}\boldsymbol{\mathsf{I}}\right)\right]\right], \end{align}

(2.37)\begin{align} \kappa^{\theta}_{ih} & = n_{pi}n_{ph}\left(\frac{m_{pi}m_{ph}}{\varTheta_{p{i}}\varTheta_{p{h}}}\right)^{3/2}\mathcal{A}^*g_0\frac{d_{pih}^4}{\sqrt{\rm \pi}}\frac{5}{336} N_1, \end{align}

(2.38)\begin{align} \boldsymbol{D}^{\theta}_{ih} & = n_{pi}n_{ph}\left(\frac{m_{pi}m_{ph}}{\varTheta_{p{i}}\varTheta_{p{h}}}\right)^{3/2}\mathcal{A}^*\epsilon_0g_0\frac{5}{112}d_{pih}^4\sqrt{\rm \pi}\left[\frac{1}{3}\left(\boldsymbol{\nabla}\ln\left( \frac{n_{ph}}{n_{pi}}\right) + \frac{3}{2}\boldsymbol{\nabla}\ln\left( \frac{\varTheta_{p{i}}}{\varTheta_{p{h}}}\right)\right) N_1\right.\nonumber\\ & \quad + \frac{1}{8}\left(m_{ph}\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} - m_{pi}\frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right)N_2 + \frac{1}{6}\frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})^2} \left(m_{pi}\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} - m_{ph}\frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right)N_3 \nonumber\\ &\left. \quad + \frac{1}{3}B\frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})}\left(\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} + \frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right)N_4\right]. \end{align}

The source term, $\chi ^q_{ih}$, is written as

(2.39)\begin{equation} \chi^q_{ih} ={-} \boldsymbol{\sigma}^{\chi}_{ih}\boldsymbol{\cdot}\boldsymbol{E} + D^{\chi}_{ih} \left(\frac{\varphi_{pi} - \varphi_{ph}}{\delta_c e} + \frac{1}{{\rm \pi}\epsilon_0}\left(\frac{Q_{ph}}{d_{ph}^2} - \frac{Q_{pi}}{d_{pi}^2}\right)\right) , \end{equation}

with the triboelectric source conductivity vector, $\boldsymbol {\sigma }^{\chi }_{ih}$, and the triboelectric source phase coupling coefficient, $D^{\chi }_{ih}$, that are defined as

(2.40) \begin{align} D^{\chi}_{ih} & = n_{pi}n_{ph}\left(\frac{m_{pi}m_{ph}}{\varTheta_{p{i}}\varTheta_{p{h}}}\right)^{3/2}\notag\\ &\quad \times \mathcal{A}^*\epsilon_0g_0d_{pih}^2\frac{5\sqrt{\rm \pi}}{28}\left[ N_1 - \frac{7d_{pih}}{57}\frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})}\left(\frac{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{ph}}{\varTheta_{p{h}}} + \frac{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pi}}{\varTheta_{p{i}}}\right)N_5\right], \end{align}

(2.41)\begin{align} \boldsymbol{\sigma}^{\chi}_{ih} & = n_{pi}n_{ph}\left(\frac{m_{pi}m_{ph}}{\varTheta_{p{i}}\varTheta_{p{h}}}\right)^{3/2}\mathcal{A}^*\epsilon_0g_0d_{pih}^3\frac{5\sqrt{\rm \pi}}{168}\left[ \left[\boldsymbol{\nabla}\ln\left( \frac{n_{ph}}{n_{pi}}\right) + \frac{3}{2}\boldsymbol{\nabla}\ln\left( \frac{\varTheta_{p{i}}}{\varTheta_{p{h}}}\right)\right]N_1 \right.\nonumber\\ & \quad + \frac{3}{4}\left(m_{ph}\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} - m_{pi}\frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right) N_2 + \frac{m_{pi}m_{ph}}{2(m_{pi} + m_{ph})^2}\left(m_{pi}\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} - m_{ph}\frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right) N_3 \nonumber\\ &\left.\left. \quad + \frac{m_{pi}m_{ph}}{(m_{pi} + m_{ph})}\left(\frac{\boldsymbol{\nabla}\varTheta_{p{h}}}{\varTheta_{p{h}}^2} + \frac{\boldsymbol{\nabla}\varTheta_{p{i}}}{\varTheta_{p{i}}^2}\right)B N_4\right]\right] . \end{align}

The coefficients $N_k$ ($k=1,\ldots ,5$) and $B$ in the flux and source terms are listed in table 1.

2.6. Constitutive relations with linear departures from mixture properties

The model derived above has been given under a general form assuming the probability function to be a Maxwellian distribution with a distinct granular temperature and mean velocity for each solid phase. These assumptions about the probability function are valid as long as the flow is nearly elastic and close to equilibrium. To extend the range of application for the model, it should include a non-Maxwellian distribution (i.e. Galvin et al. Reference Galvin, Dahl and Hrenya2005; Garzó et al. Reference Garzó, Dufty and Hrenya2007b) to take into account the inelastic collisions and their consequences on the collisional terms as well as the deviation from equilibrium. It is important to underline that our charge model assumed there is no correlation between the charge and velocity probability function, therefore analytical derivations of charge collision integrals given in Appendix A are independent of the hydrodynamic model for either nearly elastic or inelastic collisions.

In this study, the assessment for the model including the charge transport equation for a binary mixture is done assuming that the particles are fully elastic and the system deviated slightly from equilibrium. To be consistent with this underlying assumption, the model could be simplified with the granular temperature and mean velocity of each solid phase undergoing a linear variation from the mixture value by following Jenkins & Mancini (Reference Jenkins and Mancini1987). The granular temperature and mean velocity for phase $h$ ($h = i,j$) can be defined as

(2.42a,b)\begin{equation} \varTheta_{p{h}} = \varTheta_{p{m}} + \delta_{\varTheta_{p{h}}}, \quad \boldsymbol{U}_{ph} = \boldsymbol{U}_{pm} + \delta_{\boldsymbol{U}_{ph}}, \end{equation}

where $\delta _{\varTheta _{p{h}}}$ and $\delta _{\boldsymbol {U}_{ph}}$ are the linear deviations for the granular temperature and mean velocity, respectively, and the subscript $m$ refers to the mixture property. As these deviations being small, the coefficients $M_k$ ($k = [1,14]$) and $N_l$ ($l = [1,5]$) listed in table 1 can be simplified by accounting for the first term only; the higher-order terms being proportional to the difference between the granular temperature via the coefficient $B$ can be neglected ($B \propto (\delta _{\varTheta _{p{j}}} - \delta _{\varTheta _{p{i}}}) \approx 0$). Assuming that the gradient term multiplied by the linear deviations turns to zero ($\delta _{\varTheta _{p{h}}}\boldsymbol {\nabla } \approx 0$), all collisional terms given above can be reduced in new expressions as a function of the mixture properties. For the momentum transport equation, the flux and source expressions become, respectively,

(2.43) \begin{align} \boldsymbol{\theta}_{ih}^{{m}} & = {\rm \pi} n_{pi}n_{ph}(1+e_c)g_0\frac{d_{pih}^3}{3}\left[\left(\varTheta_{p{m}} + \frac{m_{pi}\delta_{\varTheta_{p{j}}} + m_{ph}\delta_{\varTheta_{p{i}}}}{(m_{pi} + m_{ph})}\right)\boldsymbol{\mathsf{I}}\right.\nonumber\\ &\quad - \left. \frac{2d_{pih}}{5}\left( \frac{(m_{pi}m_{ph})}{{\rm \pi}(m_{pi} + m_{ph})}\varTheta_{p{m}}\right)^{1/2} \times\left[\left((\boldsymbol{\nabla}\boldsymbol{U}_{pm}) + (\boldsymbol{\nabla}\boldsymbol{U}_{pm})^{\textrm{T}}\right) + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pm}\boldsymbol{\mathsf{I}} \right]\vphantom{\frac{m_{pi}\delta_{\varTheta_{p{j}}} + m_{ph}\delta_{\varTheta_{p{i}}}}{(m_{pi} + m_{ph})}}\right], \end{align}

(2.44)\begin{align} \boldsymbol{\chi}_{ih}^{{m}} & = {\rm \pi}n_{pi}n_{ph}(1+e_c)g_0\frac{d_{pih}^3}{3}\varTheta_{p{m}}\left[\frac{4}{d_{pih}}(\delta_{\boldsymbol{U}_{pj}} - \delta_{\boldsymbol{U}_{pi}})\left(\frac{m_{pi}m_{ph}}{2{\rm \pi}(m_{pi} + m_{ph})\varTheta_{p{m}}}\right)^{1/2} \right.\nonumber\\ &\left. \quad + \boldsymbol{\nabla}\ln\frac{n_{pi}}{n_{ph}} + \frac{(m_{ph} - m_{pi})}{(m_{pi} + m_{ph})}\boldsymbol{\nabla}\ln\varTheta_{p{m}} \right]. \end{align}

For the granular temperature transport equation, the flux and source terms are given by

(2.45) \begin{align} \boldsymbol{q}_{ih}^{{m}} & = \frac{n_{pi}n_{ph}}{(m_{pi} + m_{ph})}(1+e_c)g_0 d_{pih}^4\varTheta_{p{m}}\left[ -\frac{2}{3}\left(\frac{2{\rm \pi}\varTheta_{p{m}}m_{pi}m_{ph}}{(m_{pi}+m_{ph})}\right)^{1/2}\boldsymbol{\nabla}\ln\varTheta_{p{m}} + m_{ph}(1-e_c)\right. \nonumber\\ & \quad \times\left[\frac{\rm \pi}{6d_{pih}}(\delta_{u,h} - \delta_{u,i}) + \frac{1}{4}\left(\frac{2{\rm \pi}(m_{pi} + m_{ph})\varTheta_{p{m}}}{m_{pi}m_{ph}}\right)^{1/2}\right.\nonumber\\ &\left. \left. \quad \times\left(\frac{2}{3}\boldsymbol{\nabla}\ln\left(\frac{n_{pi}}{n_{ph}}\right) + \frac{(m_{ph} - m_{pi})}{(m_{pi} + m_{ph})}\boldsymbol{\nabla}\ln\varTheta_{p{m}}\right)\right]\right], \end{align}

(2.46)\begin{align} \gamma_{ih}^{{m}} & = 2n_{pi}n_{ph}g_0d_{pih}^2(1+e_c)\varTheta_{p{m}}\left[2\frac{(\delta_{\varTheta,h} - \delta_{\varTheta,i})}{\varTheta_{p{m}}(m_{pi}+m_{ph})}\left(\frac{2{\rm \pi} m_{pi}m_{ph}\varTheta_{p{m}}}{(m_{pi}+m_{ph})}\right)^{1/2}\right. \nonumber\\ & \quad - \frac{{\rm \pi} d_{pih}(m_{ph} - m_{pi})}{10(m_{pi}+m_{ph})}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pm} - \frac{m_{ph}}{(m_{pi}+m_{ph})}(1-e_c)\left[\left(\frac{2{\rm \pi}(m_{pi} + m_{ph})\varTheta_{p{m}}}{m_{pi}m_{ph}}\right)^{1/2}\right.\nonumber\\ &\left. \quad - \frac{{\rm \pi} d_{pih}}{2} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pm}\right] + \frac{{\rm \pi} d_{pih}}{2(m_{pi}+m_{ph})}\left[ \frac{1}{5}\left(m_{pi}\boldsymbol{\nabla}\boldsymbol{\cdot}\delta_{u,i} - m_{ph}\boldsymbol{\nabla}\boldsymbol{\cdot}\delta_{u,h}\right)\right.\nonumber\\ &\quad + \left.\left. \frac{m_{ph}}{2}(1-e_c) \left(\boldsymbol{\nabla}\boldsymbol{\cdot}\delta_{u,h} + \boldsymbol{\nabla}\boldsymbol{\cdot}\delta_{u,i}\right)\vphantom{\frac{1}{5}}\right]\vphantom{\frac{1}{5}}\right]. \end{align}

In addition in our generic model, the collisional terms in the mean charge transport equation can also be expressed in term of small deviations from the mixture parameter for the granular temperature and the mean velocity. Following (2.35), the triboelectric conductivity tensor, diffusivity and phase coupling coefficient in the flux term for unlike particles collision expression can be shortened as

(2.47)\begin{align} \boldsymbol{\sigma}^{\theta {m}}_{ih} & = \sqrt{\rm \pi}g_0d_{pih}^3n_{pi}n_{ph}\mathcal{A}^*\epsilon_0 \left(2\varTheta_{p{m}}\frac{(m_{pi} + m_{ph})}{m_{pi}m_{ph}}\right)^{9/10}\varGamma\left(\frac{12}{5}\right)\left[ -\frac{5}{21}\boldsymbol{\mathsf{I}}\right. \nonumber\\ &\quad + \left. d_{pih}\frac{6}{551}\left(\frac{6m_{pi}m_{ph}}{5(m_{pi}+m_{ph})\varTheta_{p{m}}}\right)^{1/2}\left[ (\boldsymbol{\nabla}\boldsymbol{U}_{pm}) + (\boldsymbol{\nabla}\boldsymbol{U}_{pm})^{\textrm{T}} + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pm}\boldsymbol{\mathsf{I}}\right] \right], \end{align}

(2.48)\begin{align} \kappa^{\theta {m} }_{ih} & = \frac{5}{42}\sqrt{\rm \pi}g_0d_{pih}^4n_{pi}n_{ph}\mathcal{A}^*\epsilon_0 \left(2\varTheta_{p{m}}\frac{(m_{pi} + m_{ph})}{m_{pi}m_{ph}}\right)^{9/10}\varGamma\left(\frac{12}{5}\right), \end{align}

(2.49)\begin{align} \boldsymbol{D}^{\theta {m} }_{ih} & = \frac{1}{14}\sqrt{\rm \pi}g_0d_{pih}^4n_{pi}n_{ph}\mathcal{A}^*\epsilon_0 \left(2\varTheta_{p{m}}\frac{(m_{pi} + m_{ph})}{m_{pi}m_{ph}}\right)^{9/10}\varGamma\left(\frac{12}{5}\right)\left[\frac{5}{3}\boldsymbol{\nabla}\ln\left( \frac{n_{ph}}{n_{pi}}\right)\right. \nonumber\\ &\left. \quad -\frac{11}{4}\frac{(m_{ph} - m_{pi})}{(m_{pi} + m_{ph})}\boldsymbol{\nabla}\ln\varTheta_{p{m}} \right]. \end{align}

For similar particle collision ($h = i$), the triboelectric phase coupling coefficient turns to zero and the other coefficients become

(2.50)\begin{align} \boldsymbol{\sigma}^{\theta {m} }_{ii} & = 2^{9/5}\sqrt{\rm \pi}g_0d_{pi}^3n_{pi}^2\mathcal{A}^*\epsilon_0 \left(\frac{\varTheta_{p{m}}}{m_{pi}}\right)^{9/10}\varGamma\left(\frac{12}{5}\right)\left[ -\frac{5}{21}\boldsymbol{\mathsf{I}} \right.\nonumber\\ &\quad + \left. d_{pi}\frac{6}{551}\left(\frac{3m_{pi}}{5\varTheta_{p{m}}}\right)^{1/2}\left[ (\boldsymbol{\nabla}\boldsymbol{U}_{pm}) + (\boldsymbol{\nabla}\boldsymbol{U}_{pm})^{\textrm{T}} + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pm}\boldsymbol{\mathsf{I}}\right] \right], \end{align}

(2.51)\begin{align} \kappa^{\theta {m} }_{ii} & = 2^{4/5}\frac{5}{21}\sqrt{\rm \pi}g_0d_{pi}^4n_{pi}^2\mathcal{A}^*\epsilon_0 \left(\frac{\varTheta_{p{m}}}{m_{pi}}\right)^{9/10}\varGamma\left(\frac{12}{5}\right). \end{align}

For the source term, the coefficients in (2.39) for unlike particle collision expression become

(2.52)\begin{align} \boldsymbol{\sigma}^{\chi {m} }_{ih} & = g_0\mathcal{A}^*\epsilon_0n_{pi}n_{ph}d_{pih}^3\frac{5\sqrt{\rm \pi}}{21}\left(2\varTheta_{p{m}}\frac{(m_{pi} + m_{ph})}{m_{pi}m_{ph}}\right)^{9/10}\varGamma\left(\frac{12}{5}\right)\left[\boldsymbol{\nabla}\ln\left( \frac{n_{ph}}{n_{pi}}\right) \right.\nonumber\\ &\quad - \left. \frac{9}{10}\frac{(m_{ph} - m_{pi})}{(m_{pi} + m_{ph})}\boldsymbol{\nabla}\ln\varTheta_{p{m}} \right], \end{align}

(2.53)\begin{align} D^{\chi {m} }_{ih} & = g_0\mathcal{A}^*\epsilon_0n_{pi}n_{ph}d_{pih}^2\frac{10\sqrt{\rm \pi}}{7}\left(2\varTheta_{p{m}}\frac{(m_{pi} + m_{ph})}{m_{pi}m_{ph}}\right)^{9/10}\varGamma\left(\frac{12}{5}\right)\nonumber\\ & \quad \times\left[ 1 - \frac{28d_{pih}}{57}\left(\frac{6m_{pi}m_{ph}}{5(m_{pi}+m_{ph})\varTheta_{p{m}}}\right)^{1/2} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}_{pm} \right]. \end{align}

For similar particle collisions, the source term for the mean charge transport equation turns to zero.

3.1. Spatially homogeneous bidisperse granular gas simulations

We performed hard-sphere simulations of elastic granular gases in a fully periodic cubic box with a dimension of $32d_{pj}\times 32d_{pj}\times 32d_{pj}$. Here, $d_{pj}$ refers to the larger particle with a diameter of $300\ \mathrm {\mu }\textrm {m}$. The particle density and diameter, the domain-averaged solid volume fraction, the granular temperature, the initial charge, the number of particles ($N_p$) and the number of collisions per particles ($n_c/N_p$) for each phase are listed for three simulation cases in table 2. For all simulations, the particles were randomly distributed in the domain and velocities were initialized with a Maxwellian distribution with a zero mean velocity for each solid phase. The initial charges followed a bimodal distribution imposing a mixture charge, $Q_{pm}$, equal to zero. The difference of the work function was set to zero as well. As the mixture charge was zero all at times, the macroscopic electric field turns to zero, therefore, there was no electrostatic force acting on particles. For the phase $i$, we set the initial mean charge equal to the reference charge, $Q_p^0$, while the initial mean charge for the phase $j$ was imposed by the ratio of particle number density as

(3.9)\begin{equation} \left. \begin{gathered} Q_{pi}(x, 0) = Q_p^0 \\ Q_{pj}(x, 0) ={-} \frac{n_{pi}}{n_{pj}}Q_p^0 \end{gathered} \right\}. \end{equation}

Table 2. Particle properties and flow parameters for three spatially homogeneous flow configurations. Case-A refers to a monodisperse case (two solid phase classes with the same particle properties but with opposite initial charges); Case-B refers to a bidisperse case with a particle diameter ratio of $R_d$ = 5 and the same particle density; Case-C refers to a bidisperse case with a particle diameter ratio of $R_d$ = 3 and a particle density ratio of $R_{\rho } = 10$. For all cases, a bimodal distribution of charge is imposed by following (3.9) with a total charge equal to zero, $Q_{pm} = 0$ (electrically neutral condition) and the restitution coefficient is set to unity, $e_c = 1$.

For spatially homogeneous flow with elastic collisions ($e_c=1$) and without mean convection, the set of equations given in the previous sections will be simplified to ordinary differential equations for granular temperature and charge evolution for each phase as

(3.10)$$\begin{gather} \frac{\textrm{d}\varTheta_{p{i}}}{\textrm{d}t} = n_{pj}\frac{m_{pi}m_{pj}}{m_{pi} + m_{pj}}\left(\frac{m_{pi}m_{pj}}{\varTheta_{p{i}}\varTheta_{p{j}}}\right)^{3/2}\frac{\sqrt{\rm \pi}}{3}g_0d_{pij}^2BM_7, \end{gather}$$

(3.11)$$\begin{gather}\frac{\textrm{d}\varTheta_{p{j}}}{\textrm{d}t} ={-} n_{pi}\frac{m_{pi}m_{pj}}{m_{pi} + m_{pj}}\left(\frac{m_{pi}m_{pj}}{\varTheta_{p{i}}\varTheta_{p{j}}}\right)^{3/2}\frac{\sqrt{\rm \pi}}{3}g_0d_{pij}^2BM_7, \end{gather}$$

(3.12)$$\begin{gather}\frac{\textrm{d}Q_{pi}}{\textrm{d}t} = n_{pj}\left(\frac{m_{pi}m_{pj}}{\varTheta_{p{i}}\varTheta_{p{j}}}\right)^{3/2}\frac{\mathcal{A}^*}{\sqrt{\rm \pi}}g_0\frac{5d_{pij}^2}{14}\left(\frac{Q_{pj}}{d_{pj}^2} - \frac{Q_{pi}}{d_{pi}^2}\right)N_1, \end{gather}$$

(3.13)$$\begin{gather}\frac{\textrm{d}Q_{pj}}{\textrm{d}t} ={-} n_{pi}\left(\frac{m_{pi}m_{pj}}{\varTheta_{p{i}}\varTheta_{p{j}}}\right)^{3/2}\frac{\mathcal{A}^*}{\sqrt{\rm \pi}}g_0\frac{5d_{pij}^2}{14}\left(\frac{Q_{pj}}{d_{pj}^2} - \frac{Q_{pi}}{d_{pi}^2}\right)N_1. \end{gather}$$

We first started with the monodisperse flow configuration represented by Case-A in table 2 where we had two solid classes with the same particle properties and domain-averaged solid volume fraction but opposite initial charges. The granular temperature for each phase was also identical, therefore the charge only evolved due to the mean charge difference between phases. The total solid volume fraction was set to $\langle \alpha _p\rangle = 0.2$ and half of the particles were assigned initial charge of $Q_p$ = $-1\ \textrm {fC}$ while the other half were assigned opposite charge of $Q_p = 1\ \textrm {fC}$. The scaled charge evolutions by simulation and model prediction for each identical phase are shown in figure 1. The orange line shows the solutions of (3.12) and (3.13) while the symbols show hard-sphere simulation results for phases $i$ (${\odot }$ in blue) and $j$ (${\varDelta }$ in blue). The mean charge of each phase follows an exponential trend in time until they reach the total charge which is equal to zero. In Case-B, we performed a simulation of bidisperse solid mixtures with a particle diameter ratio of $R_d = 5$ and compared results with solutions of an equation set given in (3.10), (3.11), (3.12) and (3.13). The total solid volume fraction was set to $\langle \alpha _p\rangle = 0.2$ with a large-to-small particle solid volume fraction ratio of $R_{\alpha }$ = 3. We also imposed different initial granular temperatures for each phase. The initial mean charges for the phases $i$ and $j$ were equal to $Q_{pi} = -1\ \textrm {fC}$ and $Q_{pj} = 42.3\ \textrm {fC}$, respectively. Granular temperature and charge equations were coupled with sequential solutions. Figure 2 shows the time evolution of the scaled granular temperature and the scaled charge for each phase. The granular temperature for each phase rapidly reaches the equilibrium state (dimensionless granular temperature which is equal to one) at $t^{*}=10$ (figure 2a). In contrast, the mean charge goes to zero with a slower trend (figure 2b). One can argue that the mean charge for the phase $i$ follows an exponential decay after the granular temperature reaches equilibrium, which is similar to the monodisperse case (Case-A).

Figure 1. Evolution of scaled charge of the phase $h\,(h=i,j)$ for Case-A. Particle properties and flow parameters are given in table 2. The orange lines show the solutions of (3.12) and (3.13). Here, ${\odot }$ in blue: hard-sphere simulation results for the phase $i$ and ${\varDelta }$ in blue: hard-sphere simulation results for the phase $j$. Charge was scaled by using (3.6).

Figure 2. Evolution of (a) scaled granular temperature and (b) scaled charge of the phase $h\,(h=i,j)$ for Case-B. Particle properties and flow parameters are given in table 2. The orange lines show the solutions of (3.10), (3.11) in (a) and show the solutions of (3.12) and (3.13) in (b). Here, ${\odot }$ in blue: hard-sphere simulation results for the phase $i$ and ${\varDelta }$ in blue: hard-sphere simulation results for the phase $j$. Granular temperature and charge were scaled by using (3.1a–c) and (3.6), respectively.

In Case-C, we imposed the particle diameter ratio of $R_d = 3$ and the particle density ratio of $R_{\rho } = 10$ at the same time. The total solid fraction and the large-to-small particle solid fraction ratio were identical to those of Case-B. The initial mean charges were $Q_{pi} = -1\ \textrm {fC}$ and $Q_{pj} = 9\ \textrm {fC}$. The granular temperature and mean charge evolution are shown in figure 3. A similar pattern as Case-B was observed with a quick evolution to the equilibrium temperature and slower evolution for the mean charge. It can be observed that, before $t^* = 10$, the mean charge evolution does not follow the exponential decrease due to the variation of the granular temperature. The Supplementary Movies 1 and 2 available at https://doi.org/10.1017/jfm.2021.739 show particle motions and charge evolution by the hard-sphere simulation for Case-C. For all spatially homogeneous granular gas cases, the model predictions are in excellent agreement with hard-sphere simulation results.

Figure 3. Evolution of (a) scaled granular temperature and (b) scaled charge of the phase $h\,(h=i,j)$ for Case-C. Particle properties and flow parameters are given in table 2. The orange lines show the solutions of (3.10), (3.11) in (a) and the solutions of (3.12) and (3.13) in (b). Here, ${\odot }$ in blue: hard-sphere simulation results for the phase $i$ and ${\varDelta }$ in blue: hard-sphere simulation results for the phase $j$. Granular temperature and charge were scaled by using (3.1a–c) and (3.6), respectively.

These simulation cases also allowed us to probe the truncation terms in $M_7$ and $N_1$ (see table 1). We computed the mean errors between hard-sphere simulation results and model predictions for phase granular temperature and phase mean charge for different truncation orders and truncated the expansion if the mean error is lower than 5 %.

3.2. Quasi-one-dimensional bidisperse granular gas simulations with spatial gradients

As a second assessment benchmark, we performed hard-sphere simulations of bidisperse elastic granular gases in a fully periodic rectangular box with a dimension of $384\,d_{pj}\times 12\,d_{pj}\times 12\,d_{pj}$ ($d_{pj}$ is the diameter of the larger particle). For these simulations, we imposed the solid volume fraction for each phase with a step function through the domain and granular temperature difference between phases to validate gradient terms in the derived models. Similar to homogeneous granular gas cases, the velocities were initiated with a Maxwellian distribution with a mean velocity equal to zero for each phase. The total charge was equal to zero with an initial mean charge of each phase imposed as follows:

(3.14)\begin{equation} \left(n_{pi}^L + n_{pi}^R\right)Q_{pi}(x, 0) + \left(n_{pj}^L + n_{pj}^R\right)Q_{pj}(x, 0) = 0. \end{equation}

Here, $n_{ph}^L$ and $n_{ph}^R$ refer to particle number densities of the phase $h$ for left and right sides of the domain. If the solid phase is initially charged, the charge follows a Maxwellian distribution with a pre-defined non-zero mean value and varies with a step function in the domain. The phase charge variance is then equal to a small value $\mathcal {Q}_p/m_{ph} = 1 \times 10^{-32}\ \textrm {C}^2$. The simulation campaign for quasi-one-dimensional bidisperse granular gases where we varied particle properties, domain-averaged solid volume fractions, granular temperatures and initial mean charges are listed in table 3. For all these cases, the collisions were elastic ($e_c=1$) and the electrostatic force was not taken into account.

Table 3. Particle properties and flow parameters for quasi-one-dimensional granular gas simulations. Case-D refers to a monodisperse case; Case-E refers to a bidisperse case with a particle diameter ratio of $R_d = 3$; Case-F refers to a bidisperse case with a particle diameter ratio of $R_d$ = 6; Case-G refers to a bidisperse case with a particle diameter ratio of $R_d$ = 3 and a particle density ratio of $R_{\rho }$ = 10. For all cases, a bimodal distribution of charge is imposed by following (3.9) with a total charge equal to zero, $Q_{pm}$ = 0 (electrically neutral condition) and the restitution coefficient is set to unity, $e_c$ = 1.

By following Fox (Reference Fox2014), instead of solving the granular temperature balance equation for each phase, we solve the total kinetic energy, $E_{ph}$, for a solid phase $h$, which is a conserved quantity. The total kinetic energy tensor for a solid phase $h$ is defined as

(3.15)\begin{equation} \boldsymbol{\mathsf{E}}_{ph} = \tfrac{1}{2}\left(\boldsymbol{\sigma}_{ph} + \boldsymbol{U}_{ph}\otimes\boldsymbol{U}_{ph}\right), \end{equation}

with the fluctuating kinetic energy tensor, $\boldsymbol {\sigma }_{ph}$. The granular temperature is given by the trace of $\boldsymbol {\sigma }_{ph}$ as

(3.16)\begin{equation} \frac{\varTheta_{p{h}}}{m_{ph}} = \frac{1}{3}\mathrm{tr}{(\boldsymbol{\sigma}_{ph})}. \end{equation}

Hence, the total kinetic energy is given by

(3.17)\begin{equation} E_{ph} = \frac{3}{2}\frac{\varTheta_{p{h}}}{m_{ph}} + \frac{1}{2}\mathrm{tr}(\boldsymbol{U}_{ph}\otimes\boldsymbol{U}_{ph}). \end{equation}

The complete set of mass, momentum, total kinetic energy and charge transport equation for a phase $h$ ($h = i,j$) is written in a conservative form as

(3.18)\begin{align} \left. \begin{aligned} & \frac{\partial}{\partial t} [ \alpha_{ph}] + \boldsymbol{\nabla}\boldsymbol{\cdot} [\alpha_{ph} \boldsymbol{U}_{ph}] = 0\\ & \frac{\partial}{\partial t} [ \alpha_{ph} \boldsymbol{U}_{ph}] + \boldsymbol{\nabla}\boldsymbol{\cdot} \left[\alpha_{ph} \left(\boldsymbol{U}_{ph}\otimes\boldsymbol{U}_{ph}) + \frac{1}{\rho_{ph}}(P_{ph}^{kin}\boldsymbol{\mathsf{I}} + \sum_{l = i,j}\boldsymbol{\theta}_{hl} \right) \right] = \frac{1}{\rho_{ph}}\sum_{l = i,j}\boldsymbol{\chi}_{hl}\\ & \frac{\partial}{\partial t} [ \alpha_{ph} E_{ph}] + \boldsymbol{\nabla}\boldsymbol{\cdot} \left[\alpha_{ph}E_{ph}\boldsymbol{U}_{ph} + \frac{\boldsymbol{U}_{ph}}{\rho_{ph}}\boldsymbol{\cdot}\left(P_{ph}^{kin}\boldsymbol{\mathsf{I}} + \sum_{l = i,j}\boldsymbol{\theta}_{hl}\right) + \frac{1}{\rho_{ph}}\sum_{l = i,j}\boldsymbol{q}_{hl}\right] \\ & \quad = \frac{1}{\rho_{ph}}\sum_{l = i,j}(\boldsymbol{\chi}_{hl}\boldsymbol{\cdot}\boldsymbol{U}_{ph} + \gamma_{hl})\\ & \frac{\partial}{\partial t} [ \alpha_{ph}Q_{ph}] + \boldsymbol{\nabla}\boldsymbol{\cdot} \left[\alpha_{ph}Q_{ph}\boldsymbol{U}_{ph} + \frac{m_{ph}}{\rho_{ph}} \sum_{l = i,j}\boldsymbol{\theta}_{hl}^{q}\right] = \frac{m_{ph}}{\rho_{ph}} \sum_{l = i,j}\chi_{hl}^q \end{aligned} \right\}. \end{align}

In the first block of the equation set, we have time derivative terms for the conserved quantities, $\alpha _{ph}, \alpha _{ph}\boldsymbol {U}_{ph}, \alpha _{ph} E_{ph}$ and $\alpha _{ph}Q_{ph}$. In the second block, the spatial fluxes with the collisional flux terms representing variable exchange within and between phases are given, as well as the kinetic granular pressure, $P_{ph}^{kin}$, which is defined as (e.g. Gidaspow Reference Gidaspow1994)

(3.19)\begin{equation} P_{ph}^{kin}=\alpha_{ph}\,\rho_{ph}\,\frac{\varTheta_{p{h}}}{m_{ph}}. \end{equation}

On the right-hand side, we have non-conservative source terms representing the dissipation of a quantity between different phases. The phase solid volume fraction, $\alpha _{ph}$, the phase velocity, $\boldsymbol {U}_{ph}$, the phase total kinetic energy, $E_{ph}$, and the phase mean charge, $Q_{ph}$, are then found from the conserved quantities. The phase granular temperature, $\varTheta _{p{h}}$, is computed by (3.17). The given conservative forms can be solved using any finite-volume method. Here, the time derivatives were discretized with the Euler method whereas the spatial fluxes were computed with a Lax–Friedrichs scheme with van Albada slope limiter.

We first compared hard-sphere simulation results with the developed model predictions for a quasi-one-dimensional granular gas simulation of monodisperse particles named Case-D in table 3. The charge and solid volume fraction were imposed with a step function and the total charge inside the system was equal to zero. The initial conditions for hard-sphere simulation with blue dots are shown in figure 4. This figure also shows the evolution of solid volume fraction, velocity, granular temperature and charge by the simulation and the model predictions (solid lines) at $t^* = 33.33$ (top) and $t^* = 66.6$ (bottom). The phase velocity rapidly evolves due to the solid volume fraction gradient and the wave-like behaviour is seen through the periodic domain in time (figure 4b). The non-uniform granular temperature profile also develops as a wave-like shape (figure 4c). Initially, the granular pressure induces a solid flux that leads to spatial redistribution of the granular temperature due to the stress production term in the fluctuating energy equation and ‘thermal’ diffusion. Wave-like solid fluxes from left and right boundaries also lead to local increases and decreases in granular temperature. As expected, the charge distribution dissipates and goes to zero (figure 4d). All these behaviours are very well captured by our model predictions. However, the Eulerian model overestimated the solid volume fraction and the granular temperature at $x/L=0.2$ and $t^* = 33.33$ (figure 4a,c). Additionally, the location of velocity sharp gradient $t^* = 66.66$ was slightly mispredicted by the model.

Figure 4. Evolution of (a) phase solid volume fraction, (b) phase scaled velocity, (c) phase scaled granular temperature and (d) phase scaled mean charge for Case-D at two time instants. Top and bottom figures refer to the simulation results and the model predictions at $t^*=33.33$ and $t^*=66.66$, respectively. Here, ${\boldsymbol {\odot }}$ in blue: initial conditions for the hard-sphere simulation. black solid line: Eulerian model predictions and ${\boldsymbol {\odot }}$ in orange: hard-sphere simulation results. Variables were scaled by using (3.1a–c), (3.2) and (3.6).

Case-E presents a quasi-one-dimensional simulation case of a bidisperse solid mixture with a particle diameter ratio of $R_d=3$. The total solid volume fraction $\langle \alpha _p\rangle$ was equal to 0.237 and, initially, only the phase $i$ particles were charged. The simulation results and the model predictions for the phases $i$ and $j$ at $t^* = 25$ and $t^* = 50$ are shown in figures 5 and 6, respectively. One can see that both solid phases rapidly reach the mixture mean velocity and granular temperature. After phases reach the equilibrium state, the flow is mainly driven by the gradient of solid volume fraction, which is a slow process for this specific case (figures 5a and 6a). The charge evolution shows the charge repartition between the solid phases and the phase $j$ picks up a large amount of charge before decreasing towards zero charge in time (figures 5d and 6d). This second stage is slower as the granular temperature of each phase has reached the mixture value; the contribution from the electric field in (2.35) tends to zero and the exchange of charge is driven mainly by the difference of charge between phases. To conclude, the hard-sphere simulation results are very well predicted by the Eulerian model.

Figure 5. Evolution of (a) phase solid volume fraction, (b) phase scaled velocity, (c) phase scaled granular temperature and (d) phase scaled mean charge of the phase $h\,(h=i,j)$ for Case-E at $t^*=25$. Here, ${\boldsymbol {\odot }}$ in blue: initial conditions for the hard-sphere simulation. black solid line: Eulerian model predictions and ${\boldsymbol {\odot }}$ in orange: hard-sphere simulation results. Variables were scaled by using (3.1a–c), (3.2) and (3.6).

Figure 6. Evolution of (a) phase solid volume fraction, (b) phase scaled velocity, (c) phase scaled granular temperature and (d) phase scaled mean charge of the phase $h\,(h=i,j)$ for Case-E at $t^*=50$. Here, black solid line: Eulerian model predictions and ${\boldsymbol {\odot }}$ in orange: hard-sphere simulation results. Variables were scaled by using (3.1a–c), (3.2) and (3.6).

We compared the simulation results and the model predictions for a larger particle diameter ratio and the charged particles for both phases in Case-F. The particle diameter ratio, $R_d$, was set to 6 and the average solid fraction was slightly less than that of Case-E. The simulation results and the model predictions for the phases $i$ and $j$ at $t^{*}=28.5$ and $t^*=85.5$ are shown in figures 7 and 8, respectively. Similar to Case-E, the phase granular temperatures and the phase velocities quickly saturate to the mixture values. The predictions of hydrodynamic variables evolution are in very good agreement with the simulation results. However, there is a discrepancy between the results for charge evolution and particularly, the phase charges are overestimated at $t^*=85.5$ by the model predictions. This difference might be explained by the self-diffusion of charge with velocity-charge correlation (Montilla et al. Reference Montilla, Ansart and Simonin2020) in the dilute regime which is not modelled for this case.

Figure 7. Evolution of (a) phase solid volume fraction, (b) phase scaled velocity, (c) phase scaled granular temperature and (d) phase scaled mean charge of the phase $h\,(h=i,j)$ for Case-F at $t^*=28.5$. Here, ${\boldsymbol {\odot }}$ in blue: initial conditions for the hard-sphere simulation. black solid line: Eulerian model predictions and ${\boldsymbol {\odot }}$ in orange: hard-sphere simulation results. Variables were scaled by using (3.1a–c), (3.2) and (3.6).

Figure 8. Evolution of (a) phase solid volume fraction, (b) phase scaled velocity, (c) phase scaled granular temperature and (d) phase scaled mean charge of the phase $h\,(h=i,j)$ for Case-F at $t^*=85.5$. Here, black solid line: Eulerian model predictions and ${\boldsymbol {\odot }}$ in orange: hard-sphere simulation results. Variables were scaled by using (3.1a–c), (3.2) and (3.6).

In Case-G, we set the particle diameter ratio, $R_d$, to 3 and set the density ratio, $R_{\rho }$, to 10. The average solid volume fraction and the number of particles are identical to Case-E. Each phase is initially charged following (3.14). The simulation results and the model predictions for Case-G at $t^{*}=25$ and $t^*=50$ are shown in figures 9 and 10, respectively. The charge for both phases shows a ‘wavy’ pattern in the domain and these trends are very well captured with our model. The animation of charge evolution by the hard-sphere simulation for Case-G is given in Supplementary Movie 3.

Figure 9. Evolution of (a) phase solid volume fraction, (b) phase scaled velocity, (c) phase scaled granular temperature and (d) phase scaled mean charge of the phase $h\,(h=i,j)$ for Case-G at $t^*=25$. Here, ${\boldsymbol {\odot }}$ in blue: initial conditions for the hard-sphere simulation. black solid line: Eulerian model predictions and ${\boldsymbol {\odot }}$ in orange: hard-sphere simulation results. Variables were scaled by using (3.1a–c), (3.2) and (3.6).

Figure 10. Evolution of (a) phase solid volume fraction, (b) phase scaled velocity, (c) phase scaled granular temperature and (d) phase scaled mean charge of the phase $h\,(h=i,j)$ for Case-G at $t^*=50$. Here, black solid line: Eulerian model predictions and ${\boldsymbol {\odot }}$ in orange: hard-sphere simulation results. Variables were scaled by using (3.1a–c), (3.2) and (3.6).

3.2.1. Quasi-one-dimensional bidisperse granular gas simulation with a work function difference

In this short section, we show how the work function difference between phases generates charge with the Eulerian model. Starting from Case-E given in table 3, the work functions of 3.9 and 4.2 eV were imposed for the phases $i$ and $j$, respectively. The Eulerian simulation was performed with a negative work function difference $\Delta \varphi _p =\varphi _{pi}-\varphi _{pj}$ for a duration of $t^*$ = 13 410 (equal to 60 s in physical time). It is worth noting that the hard-sphere simulation is not computationally affordable for this duration but it is well established that the charge build-up is a slow process (i.e. it takes more than a few minutes to reach the saturated charge in the vibrated beds shown by Kolehmainen et al. Reference Kolehmainen, Sippola, Raitanen, Ozel, Boyce, Saarenrinne and Sundaresan2017b). Therefore, the simulations with longer duration or steady-state solutions are necessary to have a better understanding of effects of the tribocharging on the hydrodynamics.

The evolution of hydrodynamic variables and charges for Case-E with a work function difference by the Eulerian predictions is shown in figure 11. The solid volume fraction for each phase slowly reaches a flat profile (figure 11a) and the phase velocities dissipate quickly and become zero at $t^*>13\,410$ (figure 11b). The granular temperatures reach an equilibrium value which is slightly lower than that of Case-E (figure 6c). Due to the work function difference, a bipolar charge distribution occurs and each phase reaches an equilibrium value (figure 11d). For this case, we also accounted for the electrostatic force in the phase momentum equations. After a short duration ($t^* > 223.5$), the electric field became very small in the domain, therefore, the electrostatic force has a limited effect on the momentum and energy evolution. We also performed the same case with a positive function difference (not shown here) and obtained a very similar evolution for hydrodynamic variables with an inverse bipolar charge distribution (the phase $i$ has a negative charge whereas the phase $j$ has a positive charge).

Figure 11. Evolution of (a) phase solid volume fraction, (b) phase scaled velocity, (c) phase scaled granular temperature and (d) phase scaled mean charge of the phase $h\,(h=i,j)$ for Case-E with a negative work function difference between phases by the Eulerian predictions at various time instants. Work functions for phases $i$ and $j$ are 3.9 and 4.2 eV, respectively. Black, blue, red and green lines refer the Eulerian predictions at $t^* = 0, 223.5, 2\,235, 13\,410$, respectively. Variables were scaled by using (3.1a–c), (3.2) and (3.6).

3.2.2. Knudsen number analysis

We solved the hydrodynamic equations along the $x$-direction for quasi-one-dimensional simulations and, here, the validity of these computations in the continuum regime was tested with a Knudsen number analysis. The Knudsen number, $Kn$, is defined as

(3.20)\begin{equation} Kn = \frac{\lambda}{L}, \end{equation}

where $\lambda$ is the mean free path and $L$ is a macroscopic length scale. In a dense granular fluid, the mean free path is given by Garzó (Reference Garzó2005) as

(3.21)\begin{equation} \lambda= \frac{d_p}{6\sqrt{2}\alpha_p\,g_0(\alpha_p)}. \end{equation}

Fullmer et al. (Reference Fullmer, Liu, Yin and Hrenya2017) and Wang et al. (Reference Wang, Chen, Bian, Zhao and Wang2019) chose $\alpha _p/|\boldsymbol {\nabla }\alpha _p|$, $U_p/|\boldsymbol {\nabla } U_p|$ and $\varTheta _{p{}}/|\boldsymbol {\nabla } \varTheta _{p{}}|$ for the characteristic length scale, $L$. By following these studies, we define the three Knudsen numbers based on the gradients of solid volume fraction, velocity and granular temperature for the phase $i$ as follows:

(3.22)$$\begin{gather} {Kn}_{\alpha_{pi}} = \frac{5}{6\sqrt{2}}\frac{d_{pi}\,|\boldsymbol{\nabla} \alpha_{pi}|}{\alpha_{pi}^2\,g_0(\alpha_{pi})}, \end{gather}$$

(3.23)$$\begin{gather}{Kn}_{U_{pi}} = \frac{5}{12}\frac{d_{pi}\,|\boldsymbol{\nabla} U_{pi}|}{\alpha_{pi}\,g_0(\alpha_{pi})\sqrt{\varTheta_{p{i}}/m_{pi}}}, \end{gather}$$

(3.24)$$\begin{gather}{Kn}_{\varTheta_{p{i}}} = \frac{5}{6\sqrt{2}}\frac{d_{pi}\,|\boldsymbol{\nabla} \varTheta_{p{i}}/m_{pi}|}{\alpha_{pi}\,g_0(\alpha_{pi})\varTheta_{p{i}}/m_{pi}}. \end{gather}$$

For the quasi-one-dimensional simulations, instead of computing a $Kn$ value averaged over the computation domain that would be misleading, the Knudsen profiles were computed along the $x$-direction at various time instants. An example of the profiles of the three Knudsen numbers defined by (3.22), (3.23) and (3.24) for each phase ($i$ and $j$) and the mixture is shown in figure 12 for Case-E at $t^*=25$ (top) and $t^*=50$ (bottom). One can see that all three Knudsen numbers are in the range of several orders of magnitude from approximately $10^{-6}$ to $10^{-1}$. So that, the low Knudsen number assumption is valid for our simulation cases.

Figure 12. Profiles of Knudsen number based on (a) phase solid volume fraction (3.22), (b) phase velocity (3.23) and (c) phase granular temperature (3.24), for Case-E at $t^*=25$ (top) and $t^*=50$ (bottom). Blue, red and green lines refer to hydrodynamic variables of phases $i$, $j$ and mean, respectively.

3.3. Wall-bounded segregating bidisperse granular flow

In the previous validation cases, the granular temperature for each phase evolves to the mixture granular temperature. In this section, we further validate the proposed model with a three-dimensional steady segregating granular flow simulation where the non-equipartition of granular temperature persists. Figure 13 shows the computational domain which is a cubic box with a size of $L_x=L_y=L_z=24\,d_{pj}$ ($d_{pj}$ is the larger particle diameter, see table 4). A granular temperature gradient along the $x$-direction is imposed by left (cold) and right (hot) stationary conducting walls with different constant granular temperatures.

Figure 13. Computational domain of segregating granular flow simulation. Different granular temperatures were imposed to left (cold – blue) and right (hot – red) surfaces. Both surfaces are conducting wall with an effective work function difference between particles and walls of $0.001\ \text {eV}$. Electric potential was equal to zero at walls. Green spheres refer to larger particles while brown particles refer to smaller particles.

Table 4. Particle properties and flow parameters for a three-dimensional steady segregating flow configuration. Case-H refers to a bidisperse case with a particle diameter ratio of $R_d$ = 2 and the same particle density (the mass ratio is equal 8). Initial charge for each phase is imposed as $1\ \textrm {fC}$. The restitution coefficient is set to $e_c = 0.9$.

By following Galvin et al. (Reference Galvin, Dahl and Hrenya2005), for the hard-sphere simulation, we imposed post-collision velocities of a particle colliding with a wall as

(3.25)$$\begin{gather} c_{post,x} = \sqrt{-\frac{4}{3}\frac{\varTheta_{wall}}{m_{pi}}\ln{z_1}}\cos{(2{\rm \pi} z_2)}, \end{gather}$$

(3.26)$$\begin{gather}c_{post,y} = \sqrt{-\frac{4}{3}\frac{\varTheta_{wall}}{m_{pi}}\ln{z_3}}\cos{(2{\rm \pi} z_4)}, \end{gather}$$

(3.27)$$\begin{gather}c_{post,z} = \sqrt{-\frac{4}{3}\frac{\varTheta_{wall}}{m_{pi}}\ln{z_5}}\cos{(2{\rm \pi} z_6)}, \end{gather}$$

where $\varTheta _{wall}$ is the imposed granular temperature at the wall and $z_1-z_6$ are random numbers generated from a uniform distribution within the interval of $[0, 1]$. The post-collision velocity along the $x$-direction, $c_{post,x}$, is reversed according to the inward normal vector of the wall. We also imposed the effective work function difference between particles and conducting walls, $\varphi _{w-p}$, to maintain charge in the domain. Periodic boundary conditions are imposed in the $y$- and $z$-directions. We do not account for any external force, therefore the time-averaged mean velocity for each phase is equal to zero.

The particle densities and diameters and the domain-averaged solid volume fraction for each phase are listed for a simulation case (Case-H) in table 4. The granular temperature ratio between cold and hot walls was set to $\varTheta _{hot\text {-}wall}/\varTheta _{cold\text {-}wall}=6$. The particles were randomly distributed in the domain and velocities were initialized with a Maxwellian distribution with a zero mean velocity for each solid phase. The total number of particles was approximately 12 500. The particle velocities were updated when they contacted walls based on (3.25), (3.26) and (3.27) and the kinetic energy dissipated through inelastic collisions between particles (the restitution coefficient, $e_c$, was 0.9). The effective work function difference between particles and walls was set to $\Delta \varphi _{w-p} = \varphi _w - \varphi _{p} = 0.001\ \textrm {eV}$. The initial charge on particles was equal to $Q_{pi} = Q_{pj} = 1\ \textrm {fC}$.

For the hard-sphere simulation, we monitored the time evolution of domain-averaged kinetic energy and charge to ensure that the flow reached a statistically steady state. We computed the Eulerian variables such as granular temperature, mean charge for each phase and mixture granular temperature for each time step and the time averaging of these variables were carried out during a duration while an additional 4000 collisions per particle occurred. To compute the Eulerian variables, we divided the computational domain into cubic cells with a length of $2\,d_{pj}$ (the mesh configuration is $12\times 12\times 12$). As we determined the length of the unit cell, we ensured that a further mesh refinement had no effect on the computed variables. The particle Lagrangian quantities such as velocity and charge were mapped by using a simple injection method where any Lagrangian properties were averaged into a cell without using any smoothing function. A further averaging over periodic directions ($y$ and $z$) was undertaken to generate the profiles along the $x$-direction. The same mesh configuration was also used to solve Poisson's equation to compute the electric field at the contact point for charge transfer by a finite difference approach during the simulation (see Appendix B for the Lagrangian hard-sphere method).

For the charge transfer model, a boundary condition for a conducting wall is necessary. Following our earlier work (Kolehmainen et al. Reference Kolehmainen, Ozel and Sundaresan2018b), the charge balance for each phase ($h = i,j$) at a conducting wall is computed by

(3.28)\begin{equation} \sigma_{qh,w} \left( \frac{\Delta \varphi_{w-p}}{\delta_c e} - \frac{2Q_{ph}}{{\rm \pi} \varepsilon_0 d_{ph}^2} \right) + \kappa_{qh} \boldsymbol{n}_w \boldsymbol{\cdot} \boldsymbol{\boldsymbol{\nabla}} Q_{ph} - \left( \sigma_{qh} - \sigma_{qh,w} \right) \boldsymbol{E} \boldsymbol{\cdot} \boldsymbol{n}_w = 0. \end{equation}

Here, $\boldsymbol {n}_w$ is the inward normal unit vector of wall and triboelectric conductivity at the wall, $\sigma _{qh,w}$, is given by

(3.29)\begin{align} \sigma_{qh,w} = \varepsilon_0 \alpha_{ph} \left( 1 + 2 (1+e_w) \alpha_{ph} g_0 \right)\left( \frac{6}{ \sqrt[10]{2} (1+e_w) d_{ph}^2 } { \left(\frac{15\,m_{ph}(1-\nu_{ph}^2)}{8Y_{ph}\sqrt{d_{ph}}}\right)^{{2}/{5}} \varGamma\left(\frac{9}{10}\right) } \right) \varTheta_{ph}^{9/10}, \end{align}

where $e_w$ is the restitution coefficient of wall–particle collisions which is set to the particle–particle restitution coefficient, $e_c$. The symbols $\sigma _{qh}$ and $\kappa _{ph}$ are the triboelectric conductivity and diffusivity defined in (A34) and (A35), respectively.

Equation (3.30) shows the complete set of the equations simplified for a one-dimensional steady wall-bounded flow configuration

(3.30) \begin{equation} \left. \begin{aligned} & \frac{\textrm{d}}{\textrm{d}\kern0.7pt x}\left(\frac{1}{\rho_{pi}}\left[P_{pi}^{kin} + \sum_{l = i,j}\theta_{il}\right] + \frac{1}{\rho_{pj}}\left[P_{pj}^{kin} + \sum_{l = i,j}\theta_{jl}\right] \right) = 0\\ & \frac{\textrm{d}}{\textrm{d}\kern0.7pt x}\sum_{l = i,j}q_{hl} = \sum_{l = i,j}\gamma_{hl} \\ & \frac{\textrm{d}}{\textrm{d}\kern0.7pt x}\sum_{l = i,j}\theta_{hl}^{q} = \sum_{l = i,j}\chi_{hl}^q \end{aligned} \right\}. \end{equation}

The first line is the mixture momentum balance for which the flux term and the kinetic pressure are defined in (2.21) and (3.19). The energy balance for phase $h$ is given in the second line with the flux and the source terms defined in (2.25) and (2.26). The last line represents the charge balance for phase $h$ with the flux and source terms defined in (A31) and (A32). All terms have been simplified for a one-dimensional flow configuration (see figure 13) where the mean phase velocities are equal to zero.

Figure 14 shows profiles of the phase solid volume fraction, the scaled phase granular temperature and the scaled phase mean charge for Case-H along the $x$-direction. The solid volume fraction distributions (figure 14a,d) show a clear segregation of particles. The larger particles (the phase $j$) are mainly located around $x/L\approx 0.3$ and very few of them stay close to the hot wall. The smaller particles (the phase $i$) tend to group around the larger particles ($x/L\approx 0.1$ and $x/L\approx 0.4$) and, similarly to the larger particles, the number of smaller particles decreases close to the hot wall. The model predictions (black solid line) of both solid volume fraction profiles are good in agreement with hard-sphere simulation results (${\boldsymbol {\odot }}$ in orange). The scaled granular temperature (the granular temperature is divided by cold wall temperature) for the small particles (figure 14b) has a v-shape profile with a minimum value at the highest total solid concentration location ($x/L\approx 0.3$). The scaled temperature for the larger particles (figure 14e) has also a v-shape profile up to the half of the domain but it decreases to zero at $x/L>0.5$ due to the absence of the larger particles. At $x/L>0.5$, the model overestimated the granular temperature, which points out the limitation of the proposed model. It is worth noting that the larger particles are located away from walls and mainly interacted with small particles. Therefore, the imposed wall temperatures have a very limited effect on the fluctuating energy of the larger particles and in the very dilute region ($x/L>0.5$). A small number of particles stay nearly still.

Figure 14. Profiles of phase solid volume fractions (a,d), phase scaled granular temperatures (b,e) and phase scaled mean charges (c,f) for Case-H. Black solid line: Eulerian model predictions and ${\boldsymbol {\odot }}$ in orange: hard sphere simulation results. Variables were scaled by using $\varTheta _{cold\text {-}wall}$ and $Q_{eq,h}$ from (3.31).

Figure 14(c,f) presents the mean charge profiles scaled by the equilibrium charge for both phases. The equilibrium charge, $Q_{eq,h}$, is defined as

(3.31)\begin{equation} Q_{eq,h} = \frac{1}{2}\frac{{\rm \pi}\,\varepsilon_{0}}{\delta_c\,e}\Delta \varphi_{w-p}d_{ph}^2. \end{equation}

We refer to the equilibrium charge as the maximum charge that a particle can acquire by interactions of an isolated particle with a wall (without the electric field effect on the charge transfer) (Kolehmainen et al. Reference Kolehmainen, Ozel, Boyce and Sundaresan2017a). One can see that the maximum value of the scaled charge for both phases is smaller than one and this shows how the electric field hinders the total charge in the domain. As we discussed in previous sections, the proposed model has been assessed for granular flows with a mixture solid volume fraction between $0.2$ and $0.4$. For dilute flows ($\alpha _{pi}<0.1$), the correlation between charge and velocity, or namely the self-diffusion of charge, should be accounted for, which is crucial to accurately predicting the charge distribution at $x/L>0.5$ ($\alpha _{pi}<0.05$ and $\alpha _{pj}<0.005$). By following Kolehmainen et al. (Reference Kolehmainen, Ozel and Sundaresan2018b), an ad hoc model for the self-diffusion of charge has been included as

(3.32)\begin{equation} \kappa_{qh}^{+} = \frac{d_{ph} \sqrt{\varTheta_{ph}/m_{ph}}}{9 \sqrt{\rm \pi} g_0 V_{ph}}, \end{equation}

where $V_{ph}$ is the volume of a particle within the solid phase $h$. The proposed model prediction (black solid line) is in very good agreement with the hard-sphere simulation results (${\boldsymbol {\odot }}$ in orange) for the larger particles (figure 14f), even for dilute region ($x/L>0.5$). However, the predicted charge profile for the smaller particles is less accurate and the total amount of charge is underestimated (figure 14c).

To discuss why the non-equipartition of granular temperature is crucial for the charge distribution, we computed the mean charge with the proposed model but, instead of using different granular temperatures for each phase, we used the mixture granular temperature (3.2) only for the constitutive equations (3.30) by imposing $\varTheta _{p{i}} = \varTheta _{p{j}}$ = $\varTheta _{p{m}}$ (that leads $B = 0$). The flux and the source terms in the charge transfer equation are simplified for the equipartitioning of granular temperature as follows:

(3.33)\begin{align} & \theta_{ih}^q ={-} n_{pi}n_{ph}\left(\frac{m_{pi}m_{ph}}{\varTheta_{p{i}}\varTheta_{p{h}}}\right)^{3/2}\mathcal{A}^*\epsilon_0g_0\frac{d_{pih}^3}{8}\frac{5\sqrt{\rm \pi}}{21}N_1 \nonumber\\ &\quad \times \left[{-}E + \frac{d_{pih}}{2} \frac{1}{{\rm \pi}\epsilon_0}\left[ \left(\frac{Q_{ph}}{d_{ph}^2} - \frac{Q_{pi}}{d_{pi}^2}\right)\times\frac{\textrm{d}}{\textrm{d}\kern0.7pt x}\left(\ln\left( \frac{n_{ph}}{n_{pi}}\right)\right) + \frac{\textrm{d}}{\textrm{d}\kern0.7pt x}\left(\frac{Q_{ph}}{d_{ph}^2}\right) + \frac{\textrm{d}}{\textrm{d}\kern0.7pt x}\left(\frac{Q_{pi}}{d_{pi}^2}\right)\right]\right], \end{align}

(3.34)\begin{align} & \chi^{q}_{ih} = n_{pi}n_{ph}\left(\frac{m_{pi}m_{ph}}{\varTheta_{p{i}}\varTheta_{p{h}}}\right)^{3/2}\mathcal{A}^*\epsilon_0g_0d_{pih}^2\frac{5\sqrt{\rm \pi}}{28}N_1\left[ \frac{1}{{\rm \pi}\epsilon_0}\left(\frac{Q_{ph}}{d_{ph}^2} - \frac{Q_{pi}}{d_{pi}^2}\right) - \frac{d_{pih}}{6}E\frac{\textrm{d}}{\textrm{d}\kern0.7pt x}\left(\ln\left( \frac{n_{ph}}{n_{pi}}\right)\right)\right]. \end{align}

To emphasize how the non-equipartition of granular temperature takes place, the phase granular temperature with hard-sphere simulation is scaled by the mixture granular temperature and shown in figure 15(a). The larger particles contribute a big portion of the mixture granular temperature at $0.1< x/L<0.3$ whereas the small particles generate the mixture granular temperature at $x/L>0.6$. Figure 15(b,c) shows profiles of the scaled mean charge for each phase. One can see that using the mixture granular temperature leads to a nearly uniform charge distribution. Particularly, the charge conductivity, $\boldsymbol {\sigma }_{ih}$, is overestimated for the larger particles at $x/L>0.6$ due to overestimation of the granular temperature of the phase $j$ (orange dots in figure 15) by using the mixture granular temperature (dash red line in figure 15), therefore it leads to a smoother distribution with an underestimation of mean charge. As shown in (3.33)–(3.34), the charge distribution along the $x$-direction is driven only by the charge difference between the solid phases and the gradient of natural logarithm of solid volume fractions for the case with the equipartition of granular temperature. The underestimation of mean charge further reduces charge transfer due to the charge differences between phases. The slight variation on the left side ($0.2< x/L<0.4$) results from the segregation of particles and the formation of a band of smaller particles (phase $i$) around the larger particles (phase $j$). These results show that the non-equipartition of fluctuating kinetic energy should be accounted for in bidisperse granular flows with charged particles such as wall-bounded flows commonly associated with charged particle transport in powder technology.

Figure 15. (a) Phase granular temperatures scaled by the mixture granular temperature, (3.2), with hard-sphere simulations (b,c) scaled phase mean charges for Case-H. Black solid line: Eulerian model predictions with the non-equipartition of granular temperature, blue solid line: Eulerian model predictions with the mixture granular temperature and ${\boldsymbol {\odot }}$: hard-sphere simulation results. Phase mean charges are scaled by (3.31).