Charge Transport Equation for Bidisperse Collisional Granular Flows with Nonequipartitioned Fluctuating Kinetic Energy

Starting from the Boltzmann-Enskog kinetic equations, the charge transport equation for bidisperse granular flows with contact electrification is derived with separate mean velocities, total kinetic energies, charges and charge variances for each solid phase. To close locally-averaged transport equations, a Maxwellian distribution is presumed for both particle velocity and charge. The hydrodynamic equations for bidisperse solid mixtures are first revisited and the resulting model consisting of the transport equations of mass, momentum, total kinetic energy, which is the sum of the granular temperature and the trace of fluctuating kinetic tensor, and charge is then presented. The charge transfer between phases and the charge build-up within a phase are modelled with local charge and effective work function differences between phases and the local electric field. The revisited hydrodynamic equations and the derived charge transport equation with constitutive relations are assessed through hard-sphere simulations of three-dimensional spatially homogeneous, quasi-onedimensional spatially inhomogeneous bidisperse granular gases and a three-dimensional segregating bidisperse granular flow with conducting walls.

The governing physics of charge transfer are still under debate (Williams 2011;Lacks & Sankaran 2011). There are three main mechanisms suggested in the literature: (i) electron transfer (Harper 1967), (ii) ion transfer (McCarty & Whitesides 2008), and (iii) bulk material transfer (Williams 2012). All three have experimental evidence supporting them (Matsusaka et al. 2010a). In the electron transfer model, the driving force for electron transfer between contacting materials is the difference between the work functions of the materials. The electron transfer model probes the charge transfer between the conducting materials very well, but it is not applicable for insulators which have low charge mobility (Duke & Fabish 1978;Bailey 2001). The ion transfer mechanism proposes that insulators mainly exchange ions located on their surfaces during contact (McCarty & Whitesides 2008). The ions are not necessarily part of the material, but can be tied to the environment properties (e.g. humidity) and during a mechanical contact between two surfaces, some of the ions may transfer from surface-to-surface that leads to different overall charges on the surfaces (Wiles et al. 2004;McCarty & Whitesides 2008;Waitukaitis et al. 2014;Schella et al. 2017). When particles come into contact, they may also exchange material with one another. The material exchanged can have a non-zero charge difference that leads a charge transfer through a mechanism referred to as the bulk material transfer. While this possible mechanism has been known for some time (Salaneck et al. 1976), its predictability and reproducibility is questionable (Lowell & Rose-Innes 1980).
One can ask whether tribocharging via electron or ion transfer mechanisms can be captured through a simple modelling framework, which is suitable for easy integration to large-scale granular flows. Recently, we developed a Computational Fluid Dynamics-Discrete Element Method (CFD-DEM) approach for gas-solid flows that accurately predicts the effects of tribocharging on flow hydrodynamics (Kolehmainen et al. 2016(Kolehmainen et al. , 2017a. In this approach, charge transfer between particles and charge build-up in the overall system are accounted for short-range electrostatic forces using the Coulomb force with neighbouring particles and long-range electrostatic forces via Poisson's equation (Kolehmainen et al. 2016). The charge accumulation on particles is modelled by an effective-work-function based model (Laurentie et al. 2013). The effective work function is a lumped parameter that can be used to quantify charging rates and extents observed in specific experimental studies (Laurentie et al. 2013;Naik et al. 2015Naik et al. , 2016Kolehmainen et al. 2017b;Sippola et al. 2018) or quantum calculations (Naik et al. 2015). Similar CFD-DEM approaches were also developed by Pei et al. (2016) and Grosshans & Papalexandris (2017). We validated the computational framework against experimental measurements of charge on monodisperse particles in vibrated and fluidised beds (Kolehmainen et al. 2017a;Sippola et al. 2018). The studies show that the total charge in the system is well-predicted with the developed models. CFD-DEM simulations, however, are limited to relatively small systems in cm-range and not affordable for large industrial-scale systems due to the highly demanding computational effort. To achieve simulations of gas-solid flows with charged particles in larger systems, the kinetic-theory based Eulerian-Eulerian models (also called two-fluid) with tribocharging have been recently developed for monodisperse particles (Kolehmainen et al. 2018b;Ray et al. 2019)(the readers are referred to a rich literature on the two-fluid model without charge transfer, e.g. Savage & Jeffrey (1981); Jenkins & Savage (1983); Lun et al. (1984); Garzó & Dufty (1999a)). Singh & Mazza (2019) has also developed hydrodynamic equations to study homogeneous and quasi-monodisperse aggregation of charged granular gases. In two-fluid models with tribocharging, the mean charge transport equation was derived from the Boltzmann equation with an assumption of Maxwellian distributions for particle velocities and charges, and coupled with the two-fluid hydrodynamic equations. These models allow for conduction of mean charge through collisions in the presence of electric field and the boundary condition capturing charge generation at the solid boundary. Ray et al. (2019) and Montilla et al. (2020) further proposed a model for the velocity-charge covariance that accounts for the self-diffusion of charge. Kolehmainen et al. (2018b) validated the constitutive equations for mean charge transfer through hard-sphere simulation results whereas Ray et al. (2019) validated the developed models through gas-solid fluidised bed experimental data (Sowinski et al. 2012).
The recent charge transport models are only applicable for particles with a uniform size distribution. The gas-solid systems and granular flows containing particles with a variety of sizes and masses (polydisperse particles) experience specific clustering, deposition dynamics due to tribocharging that are not well understood. Furthermore, there is no consensus on the charge distribution based on particle size. As an example, Salama et al. (2013) and Schella et al. (2017) studied tribocharging of particles with bidisperse size distribution and concluded that larger particles tended to obtain a more negative charge than smaller particles. In contrast, Forward et al. (2009);Zhao et al. (2003); Lee et al. (2018) and Liu et al. (2020) observed the opposite behaviour. Very recently, Ray et al. (2020) extended their monodisperse charge model for bidisperse particles to study steady-state solution of a bipolar charging of the particles with different sizes but the same material. The charge transport closures were derived by following the kinetic theory of Jenkins & Mancini (1987) for bidisperse granular flows with assuming that deviations of phase granular temperature from the equipartitioned granular temperature are small. However, several studies showed the importance of non-equipartition of granular temperature for bidisperse segregated granular flows (Alam & Luding 2003Galvin et al. 2005;Liu et al. 2007;Serero et al. 2008). The non-equipartition of granular temperature was also shown by Wildman & Parker (2002) and Feitosa & Menon (2002) experiments where binary mixtures of solid particles were agitated in vibrating fluidised beds. It was discussed that the non-equipartition of granular temperature further increased the driving forces associated to size segregation with the gradient terms of phase granular temperatures. The extensions of kinetic theory of granular flows with bidisperse particles and non-equipartitioned granular temperatures were proposed by Garzó & Dufty (1999b); Huilin et al. (2001); Iddir & Arastoopour (2005); ;  for dilute and dense granular flows. In this study, we develop the charge transport equation for collisional bidisperse granular flows with separate mean velocities, charges, charge variances and fluctuating kinetic energies for each phase without accounting for the interstitial fluid effect. The developed model predictions are assessed through a set of hard-sphere simulations of bidisperse granular flows with various particle sizes, particle mass ratios and mixture solid volume fractions in a range from 0.2 to 0.4. The structure of the paper is as follows; we revisit mass, momentum and granular temperature transport equations for bidisperse granular flows in §2. In the latter part of this section, we present the charge transport equation with constitutive relations for binary solid mixtures. Hard-sphere simulations of spatially homogeneous and inhomogeneous elastic granular flows are introduced in §3 and their results are compared with the developed model predictions. In §3.2.1, we discuss how the work function difference within a binary mixture generates charge in inhomogeneous flow. In §3.3, we present the hard-sphere simulation results with the proposed model predictions for a segregating bidisperse granular flow bounded with conducting walls. In §4, we summarise our findings and discuss further developments to the proposed model.

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 (x, c pi , q pi , t) at position x with velocity c pi and charge q pi on particles for the discrete phase i. The number of particles in the phase i with velocity between c pi and c pi + dc ip ; and charge between q pi and q pi + dq pi at position x and time t is then given by f pi dc pi dq pi . The evolution of the probability density function follows the Boltzmann equation: (2.1) The time derivative terms in . 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.

Discrete Particle Equations
The rate of change of particle velocity is defined by the equation of motion as: Here, the external forces acting on the discrete phase i such as gravitational and fluidsolid interactions forces (e.g. drag and Archimedes forces) are not accounted and only the electrostatic force acting on a particle is accounted for as follows: where E is the resolved electric field (higher order terms due to polarization (Kolehmainen et al. 2018a) and magnetic forces (Genc & Derin 2014) were neglected). The resolved electric field is computed by solving a Poisson's equation for the electrical potential, φ, where ρ q is the charge density; and ǫ is the electrical permittivity. Then, the resolved electric field is obtained by taking the gradient of the electrical potential: (2.5) The charge transfer occurs only by collision, therefore,

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: where n pi is the number density of the phase i particles. For each phase, mean velocity, U pi , granular temperature, Θ pi , mean charge, Q pi , and charge variance, Q pi , are then defined as: where 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; Averaging the Boltzmann Equation (2.1) over a range of velocities and charges and using the relation n pi m pi = α pi ρ pi , the Enskog equation is obtained: where α pi is the solid volume fraction, ρ pi is the density and m pi is the mass of a particle in the discrete phase i. The two terms on the left-hand-side represent the transport of a quantity ψ 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, C(m pi ψ 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 (1987): (2.14) The flux term, θ ih , represents the redistribution of a quantity within and between phases while the source term, χ ih , represents the dissipation of the quantity ψ pi between the phases i and h. These terms are derived by using the following integrals: In these integrals, ψ + pi refers a quantity after collision, k is the unit vector from the centre of two particles at contact, 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 contact point. With an assumption of random motion of particles, it is approximated with a Taylor's expansion at contact point as: 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) 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 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. (2020) on the modelling of the charge-velocity correlation for monosize particles.

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 (ψ pi = 1), momentum (ψ pi = c pi ), granular temperature (ψ pi = 1 2 c ′ pi · c ′ pi = 3 2 Θ pi /m pi ) are derived from the Enskog equation (2.12). The closure relations for the collision terms are derived by following Iddir & Arastoopour (2005). 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: The momentum balance for the phase i is written as: (2.20) Here, the electric field, E , is computed with (2.4) and (2.3). 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: The source term is given by: where e c is the restitution coefficient with e c = 1 for fully elastic collisions. The coefficients M k (k = 1..6) and B are given in table 1. The derivation of these terms are not given here but the reader is referred to Iddir & Arastoopour (2005) for further details. The transport equation for granular temperature of the solid phase i is given by with the flux, q ih , and source, γ ih , terms defined as and The derived models are very similar to the ones proposed by Iddir & Arastoopour (2005) 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 (A 13)) for integration of the collision operator. Our approximation is detailed in Appendix A. Unfortunately, Iddir & Arastoopour (2005) 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.

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: The charge transfer during collision between particles is based on the model proposed by Laurentie et al. (2013) with the phase effective work function, ϕ ph=i,j , and the electric field, E . The transfer charge between particles l and m within the same phase (i.e. the phase i) is given by (2.29) For between different phases; (2.31) In (2.30) and (2.31), δ c is the cutoff distance of electron transfer, e is the elementary charge and ǫ 0 is the electrical permittivity in a vacuum. A max is the maximum overlapping area computed with help of the contact Hertz theory as (Kolehmainen et al. 2017a): 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 to the potential energy in the spring). In (2.32), 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.33) The maximum contact area could be more accurately computed for granular material with viscoelastic particles by following Schwager & Pöschel (2008) or Brilliantov & Pöschel (2010). However, due to the complex nature of contact electrification, an agreement on the charge transfer model in granular material has not been yet reached as discussed by Matsusaka et al. (2010b). The contact electrification also depends on many other particle properties as shape, roughness, surface functionalisation 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: where the flux term, θ q ih (q pi ), represents the spatial redistribution of charge and the source term, χ q ih (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 (A 31) and (A 32), respectively. Here, we present these equations in a compact form by following Kolehmainen et al. (2018b). The flux term, θ q ih , is then written as: with the triboelectric conductivity tensor, σ θ ih , the triboelectric diffusivity, κ θ ih , and the triboelectric phase coupling coefficient, D θ 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 charge on 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: The source term, χ q ih , is written as: with the triboelectric source conductivity vector, σ χ ih , and the triboelectric source phase coupling coefficient, D χ ih that are defined as: The coefficients N k (k = 1..5) and B in the flux and source terms are listed in table 1.

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 Maxwellian distribution with 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. (2005); ) 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 (1987). The granular temperature and mean velocity for phase h (h = i, j) can be defined as: where δ Θ ph and δ 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 Assuming that the gradient term multiplied by the linear deviations turns to zero (δ Θ ph ∇ ≈ 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: For the granular temperature transport equation, the flux and source terms are given by: In addition in our generic model, the collisional terms in the mean charge transport equation can be also 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: For similar particles collision (h = i), the triboelectric phase coupling coefficient turns to zero and the others coefficients become: ). (2.51) For the source term, the coefficients in (2.39) for unlike particle collision expression become: (2.53) For similar particles collisions, the source term for the mean charge transport equation turns to zero.

Model assessment
We assess the generic models through Lagrangian hard-sphere simulations for three case studies: (i) spatially homogeneous elastic granular gases with a random initial distribution of monodisperse and bidisperse solids with a bimodal charge distribution, (ii) quasi-1D bidisperse elastic granular gases with spatial gradients, and (iii) a threedimensional segregating inelastic bidisperse granular flow with conducting walls. For these simulations, we varied the particle diameter ratio, R d = d pj /d pi , the solid volume fraction ratio, R α = α pj /α pi , the particle density ratio, R ρ = ρ pj /ρ pi of phases i and j, and phase initial charges. Recalling that the proposed model is applicable only for moderately dense and dense flows where charge transfer is mainly driven by collisions, the model was assessed for granular flows with the mixture solid volume fraction in a range from 0.2 to 0.4. The readers are referred to Appendix B for details of the Lagrangian hard-sphere method. Briefly, the Lagrangian hard-sphere solver is based on the timestepped algorithm where the predicted locations are computed for each particle to ensure no overlapping occurs. In case of overlapping between two particles, the positions of these particles are reversed to the previous time step and the velocities are corrected by following the collision rule. The new locations of particles are then updated using the corrected velocity. If the overlapping distance is more than 10% percent of particle radius, the time-step is decreased. The hard-sphere simulations with a mixture volume fraction larger than 0.4 produced several collisions where the ratio of overlapping was greater than 5% of the particle diameter even for very small time-step. This was deemed unacceptable, therefore these simulations were excluded from this study. To compute the electric field at contact point, we first map particle charges into the Eulerian cells to compute charge densities and then solve the Poisson's equation with a spectral method for fully periodic simulations (Sections 3.1 and 3.2). For bounded simulations with conducting walls, we use a finite difference method to discretize the Poisson equation with the electric potential set to zero at walls.
As we present hard-sphere simulation results and Eulerian model predictions, we use the following dimensionless quantities for the phase h (h = i, j): The subscript, m, refers to the mixture quantities defined as: For spatially homogeneous and quasi-1D granular gas simulations, the mean charge of the system is set to zero, Q pm = 0. We scale the mean phase charge as with a reference charge, Q 0 p = 1 fC. For all simulations below, the Poisson's ratio and the Young's modulus were kept constant (ν ph = 0.42 and Y ph = 0.5 MPa). The radial distribution function, g 0 , proposed by Jenkins & Mancini (1987) with the coefficients µ and ξ  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ρ = 10. For all cases, a bimodal distribution of charge is imposed by following (3.9) with a total charge equal to zero, Qpm = 0 (electrically neutral condition) and the restitution coefficient is set to unity, ec = 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 32 d pj × 32 d pj × 32 d pj . Here, d pj refers to the larger particle with a diameter of 300 µ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 initialised with a Maxwellian distribution with a zero mean velocity for each solid phase. The initial charges followed a bimodal distribution with imposing the mixture charge, Q pm , equal to zero. The difference of 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 0 p , while the initial mean charge for the phase j was imposed by the ratio of particle number density as: (3.9) For spatially homogeneous flow with elastic collision (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: 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 the domainaveraged 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 α p = 0.2 and half of the particles were assigned initial charge of Q p = −1 fC while the other half were assigned opposite charge of Q p = 1 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 (⊙) and j (∆). 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 α p = 0.2 with a large-to-small particle solid volume fraction ratio of R α = 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 fC and Q pj = 42.3 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 2(a)). In contrast, the mean charge goes to zero with a slower trend ( figure 2(b)). 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). PSfrag replacements PSfrag replacements ⊙: hard-sphere simulation results for the phase i and ∆: hard-sphere simulation results for the phase j. Granular temperature and charge were scaled by using (3.1) and (3.6), respectively.
In Case-C, we imposed the particle diameter-ratio of R d = 3 and the particle densityratio of R ρ = 10 at the same time. The total solid fraction and the large-to-small particle solid fraction ratio were identical to of Case B. The initial mean charges were Q pi = −1 fC and Q pj = 9 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 Material 1 and 2 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.
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%.  Table 3. Particle properties and flow parameters for quasi-1D 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ρ = 10. For all cases, a bimodal distribution of charge is imposed by following (3.9) with a total charge equal to zero, Qpm = 0 (electrically neutral condition) and the restitution coefficient is set to unity, ec = 1.

Quasi-1D 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 × 12 d pj × 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: n L pi + n R pi Q pi (x, 0) + n L pj + n R pj Q pj (x, 0) = 0. (3.14) Here, n L ph and n R ph 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; Q p /m ph = 1 × 10 −32 C 2 . The simulation campaign for quasi-1D 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.
By following Fox (2014), 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 Hence, the total kinetic energy is given by 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: In the first block of the equation set, we have time derivative terms for the conserved quantities, α ph , α ph U ph , α ph E ph and α 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 kin ph , which is defined as (e.g. Gidaspow (1994)) On the right-hand-side, we have non-conservative source terms representing the dissipation of quantity between different phases. The phase solid volume fraction, α ph , the phase velocity, 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, Θ ph , 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-1D 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 dot-points (⊙) 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 4(b)). The non-uniform granular temperature profile also develops as a wave-like shape (figure 4(c)). Initially, the granular pressure induces solid flux that leads to spatial redistribution of the granular temperature due to 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 4(d)). 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 (figures 4(a) and c). Additionally, the location of velocity sharp gradient t * = 66.66 was slightly mispredicted by the model. Case-E presents a quasi-1D simulation case of a bidisperse solid mixture with a particle diameter ratio of R d = 3. The total solid volume fraction α p 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 5(a) and 6(a)). 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 5(d) and 6(d)). 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 phase. To conclude, the hard-sphere simulation results are very well predicted by the Eulerian model.
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 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. 2020) in the dilute regime which is not modelled for this case.
In Case-G, we set the particle diameter-ratio, R d , to 3 and set the density ratio, R ρ , 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 Material 3.

Quasi-1D 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 ∆ϕ p = ϕ pi − ϕ pj for duration of t * = 13 410 (equal to 60 s in physical time). It is worth to note 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 (e.g. it takes more than few minutes to reach the saturated charge in the vibrated beds shown by Kolehmainen et al. (2017b)). Therefore, the simulations with longer duration or steady-state solutions are necessary to have 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 11-(a)) and the phase velocities dissipate quickly and become zero at t * > 13 410 ( figure 11-(b)). The granular temperatures reach a equilibrium value which is slightly lower than of Case-E ( figure 6-(c)). Due to the work function difference, a bipolar charge distribution occurs and each phase reaches an equilibrium value ( figure 11-(d)). 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).

Knudsen number analysis
We solved the hydrodynamic equations along x-direction for quasi-one dimensional simulations and here, the validity of these computations in the continuum regime was  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.1), (3.2) and (3.6).
tested with the Knudsen number analysis. The Knudsen number, Kn, is defined as where λ 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ó (2005) as: (3.21) Fullmer et al. (2017) and Wang et al. (2019) chose α p /|∇α p |, U p /|∇U p | and Θ p /|∇Θ 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: For the quasi 1-D simulations, instead of computing a Kn value averaged over the computation domain that would be misleading, the Knudsen profiles were computed along 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 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 assumption is valid for our simulation cases.

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 3D 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 x-direction is imposed by left (cold) and right (hot) stationary conducting walls with different constant granular temperatures.
By following Galvin et al. (2005), for the hard-sphere simulation, we imposed postcollision velocities of a particle colliding with a wall as where Θ wall is the imposed granular temperature at wall, z 1 − z 6 are random numbers generated from a uniform distribution within the interval of [0, 1]. The post-collision velocity along 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, ϕ w−p , to maintain charge in the domain. Periodic boundary conditions are imposed in 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 y x z 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 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 3D 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 was imposed to 1 fC. The restitution coefficient was set to ec = 0.9.
ratio between cold and hot walls was set to Θ hot-wall /Θ cold-wall = 6. The particles were randomly distributed in the domain and velocities were initialised 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 with 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 ∆ϕ w−p = ϕ w − ϕ p = 0.001 eV. The initial charge on particles was equal to Q pi = Q pj = 1 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 number of 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×12×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 x-direction. The same mesh configuration was also used to solve the Poisson's equation to compute the electric field at 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. 2018b), charge balance for each phase (h = i, j) at a conducting wall is computed by: Here, n w is the inward normal unit vector of wall and triboelectric conductivity at the wall, σ qh,w , is given by: (3.29) where e w is the restitution coefficient of wall-particle collisions which is set to the particle-particle restitution coefficient, e c . The symbols σ qh and κ ph are the triboelectric conductivity and diffusivity defined in (A 34) and (A 35), respectively.
(3.30) shows the complete set of the equations simplified for a 1-D steady wall-bounded flow configuration: 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 (A 31) and (A 32). All terms have been simplified for a 1-D 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 x-direction. The solid volume fraction distributions (figures 14(a) and (d)) show a clear segregation of particles. The larger particles (the phase j) are mainly located around x/L ≈ 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 ≈ 0.1 and x/L ≈ 0.4) and, similarly to the larger particles, the number of the smaller particles decreases close to the hot wall. The model predictions ( ) of both solid volume fraction profiles are good in agreement with hard-sphere simulation results (⊙). The scaled granular temperature (the granular temperature is divided by cold wall temperature) for the small particles ( figure 14(b)) has a v-shape profile with a minimum value at the highest total solid concentration location (x/L ≈ 0.3). The scaled temperature for the larger particles (figure 14(e)) 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 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 to note 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. Figures 14(c) and (f) present the mean charge profiles scaled by the equilibrium charge for both phases. The equilibrium charge, Q eq,h , is defined as: We refer 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. 2017a). One can see that the maximum value of the scaled charge for both phases is smaller than one and that 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 (α 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 predict charge distribution at x/L > 0.5 (α pi < 0.05 and α pj < 0.005). By following Kolehmainen et al. (2018b), an ad-hoc model for the self-diffusion of charge has been included as: where V ph is the volume of a particle within the solid phase h. The proposed model prediction ( ) is in very good agreement with the hard-sphere simulation results (⊙) for the larger particle (figure 14(f)), even for dilute region (x/L > 0.5). However, the predicted charge profile for the smaller particle is less accurate and the total amount of charge is underestimated ( figure 14(c)).
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 Θ pi = Θ pj = Θ pm (that leads B = 0). The flux and the source terms in the charge transfer equation are simplified for the equal-partitioning of granular temperature as follows: (3.34) To emphasise 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. Figures 15(b) and (c) show 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, σ ih , is overestimated for the larger particles at x/L > 0.6 due to overestimation of the granular temperature of the phase j (orange points 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 x-direction is driven only by charge difference between 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 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 bidisperse granular flows with charge particles such as wall bounded flows commonly associated with charged particle transport in powder technology.

Conclusion
In this study, we have revisited kinetic-theory based hydrodynamic equations and derived charge transport equation for bidisperse granular flows with tribocharging. Each solid phase has separate mean velocity, total fluctuating kinetic energy, which is the sum of the granular temperature and the trace of fluctuating kinetic tensor, charge variance and mean charge. To close mass, momentum, total kinetic energy and mean charge balance equations, a Maxwellian distribution for particle velocity and charge without a cross-correlation (an assumption of both velocity and charge are independent variables) has been used for local-averaging of the Boltzmann equation. The constitutive relations of collisional flux and source terms for momentum, granular temperature and charge balance equations, which account for the rate of change of the quantities between phases and within a phase, are then presented. We introduced a finite-volume scheme to discretize and solve the transport equations and assessed the proposed models through various hard-sphere simulations of three-dimensional spatially homogeneous and quasione-dimensional spatially inhomogeneous bidisperse granular gases. For these cases, the mixture solid volume fraction were set into a range from 0.2 to 0.4. However, the hardsphere algorithm cannot handle granular flows with mixture solid volume fraction higher than 0.4 due to time-step limitation. We will improve our hard-sphere time-stepping algorithm to study these flows in a future study. In these simulations, we varied particle diameter and density ratios, initial phase charges and phase granular temperatures. The proposed model predictions were in very good agreement with the hard-sphere simulation results. Finally, a segregating bidisperse granular flow in a wall-bounded domain where the non-equipartition of granular temperature persisted was studied with a steady state solution of the proposed model and hard-sphere simulation. This steady-state solution had an acceptable level of accuracy with the hard-sphere simulation results. This case also showed the importance of accounting for the charge-velocity correlation in the dilute region and the non-equipartition of granular temperature to have an accurate charge distribution.
In this study, we limited hard-sphere simulations and the proposed model predictions to elastic granular flows (except the three-dimensional segregating bidisperse flow). For a further study, we will extend the proposed models by the Chapman-Enskog expansion by following Iddir & Arastoopour (2005) or couple the revisited Enskog theory (e.g. ) with the developed charge transfer models to consider a wider range of inelastic bidisperse granular flows.
We will also extend the proposed model with accounting for the charge-velocity correlation which is significant in dilute granular flows (Montilla et al. 2020). In this study, we only focus on the granular flows without interstitial fluid effect but the interstitial fluid has a huge impact on granular material hydrodynamics in particle technology applications such as fluidised bed and pneumatic conveying. Introducing the fluid phase will be also a topic of a future study. Additionally, the charge variance is assumed to be a constant thorough this study but as discussed by Singh & Mazza (2019), the charge variance plays a role in agglomeration of particles in homogeneous granular gases. It is necessary to develop the transport equation for charge variance and study its effects on the charge transport properties.

Appendix B. Hard sphere modelling
The hard sphere model in this study is based on the time-stepped algorithm (Kolehmainen et al. 2018b). The particles are moved along their trajectories with small steps according to their velocity: where ∆t is the time step; x (n) pi and v (n) pi are the particle position and velocity at time n∆t; and x (n+1) * pi is the predicted particle location. With the predicted locations, we ensure that particles are not overlapping according to their respective radius: x (n+1) * pi − x (n+1) * pj > d pi +d pj . When overlapping occurs, the time is reversed for these two particles to first time of contact. The velocities of particles are updated according to hard sphere model: where k ij is unit vector pointing from particle i to particle j, e c is the coefficient of restitution and m p is the mass of particle. The predicted locations in (B 1) are computed using the update velocities due to collision. During a particle-particle contact, the charge transfer occurs by following (2.30) and (2.31). The electric field at particle positions is computed by mapping all the particle charges into the Eulerian cells to compute charge densities ρ q . For fully periodic simulations, the electric field at cell each is resolved with a spectral method as: where λ is the wave-number vector; and F and F −1 refer to Fourier transform and inverse Fourier transform respectively. When a collision between a particle and the wall occurs, the particle velocity is updated following (3.25)-(3.27). During a particle-wall contact, the transfer of charge is modelled by Kolehmainen et al. (2018b) as: where α q is a geometrical factor depending on the type of collision; α q = 2 πd 2 p . A max is the maximum area of contact defined in (2.32) with the following effective parameters; Y * p = Yp 1−ν 2 p , r * p = dp 2 and m * p = m p . β q is defined as: with the first term representing the work function difference between particle and wall, ∆ϕ w−p , the second term is due to the charge carried by particle and the last term is due to the electric field resulting from charge on surrounding particles. For the wall bounded flow configuration, the electric field is solved by finite difference based on (2.4) and (2.5). First, the linear system for the electrical potential is solved using a second-order central difference with the discrete three-dimensional Laplacian matrix L: where ρ s q,cell is the interpolated charge density at the surface of the cell. At the walls, we impose the electric potential equal be zero. Finally, the electric field is then computed following (2.5) with a first-order finite difference method in each direction.