2.1 Population balance

In the population balance, each agglomerate size is treated as its own species, with a separate balance equation on the number density of that species. In its most general form, the population balance is given by (Ramkrishna Reference Ramkrishna2000)

for an agglomerate of volume $V$ with a number density of $n$ and velocity $\boldsymbol{u}$ . The left-hand side represents the transient and convective contributions, while the right-hand side represents the source and/or sinks of agglomerates of volume $V$ . More specifically, the net generation, $S$ , is defined as (Ramkrishna Reference Ramkrishna2000)

where $B_{a}$ is the volumetric birth rate due to aggregation, $D_{a}$ is the death rate due to aggregation, $B_{b}$ is the birth rate due to breakage and $D_{b}$ is the death rate due to breakage. It should be noted that birth and death due to nucleation (coalescence due to a phase change or chemical reaction) are not considered here, as the focus of this work is on the aggregation and breakage kernels.

The birth rate due to aggregation used here is similar to that generally given in the literature (Ramkrishna Reference Ramkrishna2000), and is defined as

where $\unicode[STIX]{x1D6FD}(V^{\prime },V-V^{\prime },t)$ is the ‘aggregation kernel’, which represents the frequency of aggregation events (i.e. when two agglomerates collide and form a larger agglomerate) between agglomerates of sizes $V^{\prime }$ and $V-V^{\prime }$ , per unit volume and per number density of each agglomerate size undergoing aggregation. Here, $\unicode[STIX]{x1D6FF}_{V-V^{\prime },V^{\prime }}$ is the Kronecker delta function,

It is worthwhile to note that $1+\unicode[STIX]{x1D6FF}_{V-V^{\prime },V}$ ensures the proper counting of birth due to aggregation events during symmetric aggregation (i.e. the case of $V^{\prime }=V-V^{\prime }$ is only summed over by the integrand in (2.3) once, whereas all other combinations $V^{\prime }\neq V-V^{\prime }$ are summed over twice), and does not appear in prior works. It should also be noted that the birth rate due to aggregation given in (2.3) only accounts for aggregation events in which all of the particles in the agglomerates of sizes $V^{\prime }$ and $V-V^{\prime }$ form an agglomerate of size $V$ ; collisions involving simultaneous aggregation and breakage in which an agglomerate of size $V$ is formed are accounted for in the birth rate due to breakage ( $B_{b}$ ), as detailed below.

Similarly, the rate of death rate due to aggregation is taken as

where $\unicode[STIX]{x1D6FF}_{V,V^{\prime }}$ is the corresponding Kronecker delta function. Again, the $1+\unicode[STIX]{x1D6FF}_{V,V^{\prime }}$ term is not included in prior works. This term is included here to represent the events in which two agglomerates of size $V$ collide and aggregate, thereby leading to the death of two agglomerates of size $V$ .

Before presenting the equations for the birth and death rates due to breakage ( $B_{b}$ and $D_{b}$ ), it is worth noting that an agglomerate can break due to (i) collision with another agglomerate, (ii) collision with a wall, (iii) shear forces from a surrounding fluid or (iv) centrifugal forces arising from its own rotation (Ramachandran et al. Reference Ramachandran, Immanuel, Stepanek, Litster and Doyle2009; Donahue et al. Reference Donahue, Davis, Kantak and Hrenya2012). Mechanism (i) is the focus of this work, although the approach taken here could easily be extended to mechanism (ii). Consequently, the birth and death rates due to breakage will take a form that is proportional to the frequency of collisions between two agglomerates, and hence the concentration of each (similar to the birth and death rates due to aggregation, $B_{a}$ and $D_{a}$ ).

The birth rate due to breakage is defined similarly to previous works (Fan et al. Reference Fan, Marchisio and Fox2004) as

where $\unicode[STIX]{x1D6FC}(V^{\prime },V)$ is the breakage kernel representing the frequency of breakage events between agglomerates of sizes $V^{\prime }$ and $V^{\prime \prime }$ . The quantity $b(V\mid T_{s},V^{\prime },V^{\prime \prime })$ is the daughter distribution function, which gives the number of agglomerates of size $V$ produced by a collision resulting in any breakage of agglomerates of size $V^{\prime }$ and $V^{\prime \prime }$ for a system of mixture granular temperature $T_{s}$ (i.e. a measure of the fluctuating kinetic energy of particles and agglomerates). The term $(1+\unicode[STIX]{x1D6FF}_{V^{\prime },V^{\prime \prime }})/2$ is included to prevent double counting of breakage events due to the second integral.

Finally, the death rate due to breakage (Fan et al. Reference Fan, Marchisio and Fox2004) is given by

In this work, the aggregation kernel $\unicode[STIX]{x1D6FD}$ used in (2.3) and (2.5) and the breakage kernel $\unicode[STIX]{x1D6FC}$ used in (2.6) and (2.7) are defined in a similar manner to previous literature (Goldschmidt Reference Goldschmidt2001). A notable exception is that the breakage kernel $\unicode[STIX]{x1D6FC}$ is now also dependent on the frequency of collisions (Fan et al. Reference Fan, Marchisio and Fox2004),

where $\unicode[STIX]{x1D713}_{a}(V,V^{\prime })$ and $\unicode[STIX]{x1D713}_{b}(V,V^{\prime })$ are the ‘success factors’ of aggregation and breakage respectively, which represent the fraction of collisions that result in agglomeration or breakage. Here, $\unicode[STIX]{x1D714}(V,V^{\prime })$ is the collision frequency per unit volume of agglomerates of size $V$ and $V^{\prime }$ , given by

$T_{s}$ is the granular temperature (which represents the variance of the velocities of individual particles relative to the local motion), $m_{V}$ is the mass of an agglomerate of size $V$ , $g_{0}$ is the radial distribution function (e.g. that given by Carnahan Starling (Mansoori Reference Mansoori1971)) and $d_{V,V^{\prime }}$ is the collision cylinder diameter of agglomerates of sizes $V$ and $V^{\prime }$ . A collision cylinder diameter interpretation is adopted here as agglomerates are aspherical, and thus whether or not they are in contact (undergoing a collision) is no longer a function of just position, but also orientation. The term $\text{dof}_{V}$ represents the sum of translational and rotational degrees of freedom of an agglomerate of size $V$ (e.g. for frictionless particles, a singlet has three translational degrees of freedom and $\text{dof}_{1}=3$ , a doublet has three translational and two rotational degrees of freedom and $\text{dof}_{2}=5$ ). It should be noted that this expression for the collision frequency does not take into account interphase slip, and does not assume the system to be homogeneous. In the validation systems explored here, there is no interphase slip and the system is homogeneous, so this expression is applicable.

To recap, the population-balance equations for each agglomerate size (2.1) include sink and source terms for the rates of birth and death due to aggregation and breakage (2.3), (2.5)–(2.7). The fundamental closures for the sink and source terms depend on the aggregation and breakage kernels, $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FC}$ respectively, i.e. the frequency of aggregation and breakage events (2.8)–(2.9). The kernels now require a closure for the success factors of aggregation and breakage ( $\unicode[STIX]{x1D713}_{a}$ and $\unicode[STIX]{x1D713}_{b}$ ).

2.2 Population-balance closures: success factors of aggregation $(\unicode[STIX]{x1D713}_{a})$ and breakage $(\unicode[STIX]{x1D713}_{b})$

As mentioned above, the closures for the aggregation $(\unicode[STIX]{x1D6FD}(V,V^{\prime }))$ and breakage $(\unicode[STIX]{x1D6FC}(V,V^{\prime }))$ kernels are proportional to the fraction of collisions that result in aggregation or breakage, or the success factors $\unicode[STIX]{x1D713}_{a}$ and $\unicode[STIX]{x1D713}_{b}$ . Numerous previous works (e.g. Lian, Thornton & Adams Reference Lian, Thornton and Adams1998; Jackson Reference Jackson2000; Calvert et al. Reference Calvert, Ghadiri and Tweedie2009; Donahue et al. Reference Donahue, Hrenya, Davis, Nakagawa, Zelinskaya and Joseph2010, Reference Donahue, Davis, Kantak and Hrenya2012; Ringl et al. Reference Ringl, Bringa, Bertoldi and Urbassek2012) have demonstrated that the result of a collision between two cohesive bodies (single particles and/or agglomerates) depends on their impact (relative) velocity. Moreover, a series of numerical simulations of two-particle collisions of cohesive particles indicates that there is a normal component of the impact velocity, $v_{n}$ , which demarcates the agglomeration and breakage outcomes, as illustrated in figure 1 (black solid line), which is consistent with past experimental observations for cohesion due to liquid (Donahue et al. Reference Donahue, Hrenya, Davis, Nakagawa, Zelinskaya and Joseph2010, Reference Donahue, Davis, Kantak and Hrenya2012). Accordingly, $\unicode[STIX]{x1D713}_{a}$ and $\unicode[STIX]{x1D713}_{b}$ should depend on the distribution of $v_{n}$ in a given system. In terms of continuum variables, the distribution of individual collision impact velocities is related to the (continuum) granular temperature, since the latter is a measure of the particle velocity distribution relative to the local mean velocity (i.e. a measure of particle velocity fluctuations). Therefore, the closures for the success factors $\unicode[STIX]{x1D713}_{a}$ and $\unicode[STIX]{x1D713}_{b}$ should depend on $T_{s}$ .

Figure 1.

(a) Collision schematic; $\boldsymbol{v}_{rel}=\boldsymbol{v}_{1}-v_{2}$ . (b) Regime map of binary collisions between two cohesive primary particles (singlets). Here, $v^{\ast }$ is a non-dimensionalized impact velocity ( $=v/(A/(\unicode[STIX]{x1D70C}d_{1}^{3}))^{1/2}$ , where $A$ is the Hamaker constant). The black line and points show the magnitude of the normal component of the critical impact velocity $(v_{a,1,1,crit}^{\ast })$ , which demarcates aggregation and rebound from DEM simulations of cohesive particles described in § 3.2.1. The red dashed line and points show the magnitude of the critical total relative impact velocity $(v_{a,tot,1,1,crit}^{\ast })$ . It should be noted that the normal component $v_{a,1,1,crit}^{\ast }$ is approximately independent of the impact angle, whereas the total component $v_{a,tot,1,1,crit}^{\ast }$ is not.

Expressions similar to previous works (Kantak et al. Reference Kantak, Hrenya and Davis2009; Liu Reference Liu2011) are used to calculate the success factors of both agglomeration and breakage. Namely, the success factors are defined here as

where $v_{a,V,V^{\prime },crit}$ and $v_{b,V,V^{\prime },crit}$ are the critical normal impact velocities of agglomeration and breakage for collisions between agglomerates of sizes $V$ and $V^{\prime }$ respectively. For the purposes of generality, the critical velocities of agglomeration and breakage are allowed to take on different values for a given pair of colliding objects (i.e. the critical velocities of agglomeration and breakage are different from each other), as well as for different sizes of agglomerates undergoing collisions. Here, $f(v_{n})$ is the normalized probability density function that a collision between agglomerates will have a normal impact velocity of $v_{n}$ . This method for calculating the success factors is depicted graphically in figure 2 for $\unicode[STIX]{x1D713}_{a}$ . For a given system, the critical velocities of agglomeration (vertical dashed line in figure 2) and breakage depend on the type of cohesion (van der Waals, liquid bridging, etc.) and materials involved, and can be determined mathematically from a particle–particle cohesion model (i.e. analytical solutions to two-body collisions or the DEM (e.g. Bryan & Rober Pfeffer Reference Bryan and Rober Pfeffer1991; Weber et al. Reference Weber, Hoffman and Hrenya2004)) or from experiments (e.g. Donahue et al. Reference Donahue, Davis, Kantak and Hrenya2012). Accordingly, general expressions for the success factors as functions of these critical velocities can be determined if $f(v_{n})$ is known. Kantak et al. (Reference Kantak, Hrenya and Davis2009) determined $f(v_{n})$ from a Gaussian fit of DEM data, which confirmed the accuracy of the expressions above but is not useful in a general continuum framework since the impact velocities are quantities specific to an individual collision (i.e. not continuum quantities) and thus remain unknown. Instead, we derive an analytical prediction of $f(v_{n})$ as a function of the granular temperature, a flow variable that is solved for within the continuum framework. Liu (Reference Liu2011) analytically determined the distribution of total relative velocity $(v_{tot})$ as a function of the granular temperature $T_{s}$ . However, figure 1 demonstrates that it is more appropriate to consider only the normal component of the impact velocity, for which the critical value is constant across various impact angles ( $\unicode[STIX]{x1D703}$ ). Put another way, considering $\unicode[STIX]{x1D713}$ as a function of $v_{n}$ only, instead of $\unicode[STIX]{x1D713}$ as a function of $v_{tot}$ and $\unicode[STIX]{x1D703}$ , is much cleaner since the $\unicode[STIX]{x1D703}$ distribution (also comprising quantities specific to individual collisions and not solved for in the continuum framework) does not have to be assumed or estimated. In the following, a novel analytical closure for $f(v_{n})$ is derived as a function of $T_{s}$ .

Figure 2.

Illustrative schematic of prediction of the success factor of aggregation for particles. The solid black line shows the normal impact velocity distribution $f(v_{n})$ . The dashed black line shows the critical velocity of aggregation $v_{a,V,V^{\prime },crit}$ . The grey area, which is shown expanded in the inset, is the success factor of aggregation $\unicode[STIX]{x1D713}_{a}$ .

Assuming a Maxwellian distribution of particle velocities and independence of each velocity component for agglomerates of size $V$ , the magnitude of the velocity distribution of the agglomerate (or single particle) is given by

where $\boldsymbol{v}$ is the agglomerate translational velocity, $\boldsymbol{u}$ is the local mean velocity of the agglomerates of the same size and $v_{x},v_{y},v_{z}$ are the $x$ , $y$ and $z$ components of the agglomerate translational velocity respectively. The term $\unicode[STIX]{x1D70E}_{\boldsymbol{u}}^{2}$ is the variance of the magnitude of velocities of agglomerates of the same size. It is related to the granular energy $E$ and granular temperature by

The assumption of isotropy gives

It should be noted that the incorporation of anisotropy is highly non-trivial (Sela & Goldhirsch Reference Sela and Goldhirsch1998), and it may not be necessary to incorporate it in the theory even when such higher-order effects are present (Mitrano et al. Reference Mitrano, Zenk, Benyahia, Galvin, Dahl and Hrenya2013; Fullmer & Hrenya Reference Fullmer and Hrenya2017). Further, the SSF and HCS explored in the current work are isotropic. Application of (2.15) to (2.13) yields

where $u_{x}$ is the mean velocity of the agglomerate size in the $x$ direction. By solving (2.14) for $\unicode[STIX]{x1D70E}_{\boldsymbol{u}}^{2}$ , and substituting $\unicode[STIX]{x1D70E}_{\boldsymbol{u}}^{2}$ and (2.16) into (2.13), the velocity distribution for a single agglomerate (or particle) in one direction becomes

where $p$ is the normalized probability distribution function, $v_{x}$ is the $x$ component of the velocity of an individual agglomerate and $u_{x}$ is the $x$ component of the local mean species velocity. Assuming that particles and/or agglomerates of different masses display an approximate equipartition of energy ( $T_{s,V}=T_{s,V^{\prime }}\equiv T_{s}$ ), (2.17) can be rewritten for two agglomerates of different masses $V$ and $V^{\prime }$ as

By manipulating (2.19) (redefining the coordinate system to replace $v_{x,V^{\prime }}$ with $-v_{x,V^{\prime }}$ and $u_{x,V^{\prime }}$ with $-u_{x,V^{\prime }}$ throughout), the following is obtained:

Combining the two Gaussians (2.18) and (2.20), and assuming molecular chaos (velocities of colliding particles are not correlated) gives

It should be noted that this expression includes interphase slip, and assumes that the system is homogeneous. These assumptions are appropriate in the validation systems explored here as there is no interphase slip and these systems are homogeneous. It should also be noted that the direction $x$ is arbitrary, so it is defined as the normal direction between the centres of mass of agglomerates of sizes $V$ and $V^{\prime }$ ,

In (2.21), the difference in mean velocities in the $x$ (normal) direction is the magnitude of the relative velocity between the agglomerate sizes,

where $\boldsymbol{u}_{\boldsymbol{V}}$ and $\boldsymbol{u}_{\mathbf{V}^{\prime }}$ are the mean velocities of agglomerates of sizes $V$ and $V^{\prime }$ . The distribution of normal velocities between agglomerates of sizes $V$ and $V^{\prime }$ is therefore given by

In order to get the impact velocity distribution $f(v_{n})$ , three additional terms ( $C,H(v_{n})$ and $\unicode[STIX]{x1D714}(v_{n})$ ) must be multiplied by the distribution of normal velocities,

where

Inclusion of the Heaviside function, $H(v_{n})$ , removes the relative velocities of particles moving away from one another (i.e. negative impact velocities) from the distribution as there is no chance of a collision occurring. Here, $C$ is a normalization constant which ensures that the distribution is normalized,

Therefore, $C$ is found by substituting (2.25) into (2.27) and rearranging,

The term $\unicode[STIX]{x1D6FA}(v_{n})$ accounts for the fact that impacts between particles with higher normal velocities will occur more frequently,

where $P(v_{n})$ is the time it takes for two particles destined to collide to approach one another, given by

Here, $l_{arb}$ is the arbitrary distance between the particles when they begin to approach one another. By substituting (2.26), (2.28), (2.29) and (2.30) into (2.25), the distribution of impact velocities becomes

where

and (2.28) is used to obtain the normalization constant,

Further, it should be noted that in the special case where the mean velocities of the two particle types in question are the same (i.e.  $|\boldsymbol{u}_{\boldsymbol{V}}-\boldsymbol{u}_{\boldsymbol{V}^{\prime }}|=0$ , and also the case when $V^{\prime }=V$ ), then $C_{norm}=1$ and $f(v_{n})=(v_{n}/T_{s}\unicode[STIX]{x1D70D})\exp [-v_{n}^{2}/2T_{s}\unicode[STIX]{x1D70D}]H(v_{n})$ .

The derivation of $f(v_{n})$ presented above is validated in figure 3 via DEM simulations of an SSF system (see § 3.2.2 for details on the method), in which the granular temperature is constant and the local mean velocities of each agglomerate size are equal. The impact velocities were measured just prior to collision in the DEM. Moreover, to isolate any error associated with the new derivation, the $T_{s}$ measured in the DEM was used as an input to the predictions for the impact velocity distribution (instead of solving for $T_{s}$ using continuum theory) for both the current work (2.34) and the previous prediction for the total impact velocity given by Liu (Reference Liu2011). The results indicate that the prediction of $f(v_{n})$ from the current work adequately captures the normal component of the impact velocity distribution observed in the DEM simulations. It is also worth noting that the prediction for the total component of the impact velocity distribution derived by Liu (Reference Liu2011) also matches well with the corresponding DEM data. However, as demonstrated below in figure 9 and alluded to above, the use of the total component of the impact velocity results in an inaccurate prediction of the success factor due to the dependence of its critical velocity on the impact angle, whereas the normal component of the critical velocity is independent of the impact angle (figure 1).

Figure 3.

Comparison of normal impact velocity distributions $f(v_{n})$ predicted by continuum theory (2.13) and the DEM for singlet–singlet collisions (see § 3.2.2 for more details) in SSF (solid volume fraction $\unicode[STIX]{x1D700}_{s}=0.1$ , shear rate $\unicode[STIX]{x1D6FE}^{\ast }=\unicode[STIX]{x1D6FE}d_{1}/v_{a,1,1,crit}=1.08$ ). It should be noted that Liu (Reference Liu2011) predicts the total impact velocity $(v_{tot}^{\ast })$ rather than the normal impact velocity $(v_{n}^{\ast })$ .

The distribution of impact velocities, $f(v_{n})$ , in (2.31) can now be integrated with respect to $v_{n}$ according to (2.11) and (2.12) in order to give the success factors. This integration will eliminate $v_{n}$ , a variable not present in the continuum theory, and yield a success factor that depends on $T_{s}$ , $\boldsymbol{u}_{\boldsymbol{V}}$ and $\boldsymbol{u}_{\boldsymbol{V}^{\prime }}$ , with each of these determined from kinetic-theory-based balances present in the continuum theory.

In addition to the constitutive quantities and their related variables described above $(B_{a},D_{a},B_{b},D_{b},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D713}_{a},\unicode[STIX]{x1D713}_{b},b,f(v_{n}))$ , several parameters are still needed to fully close the population balance and kinetic-theory-based balances; namely, the cohesion-specific parameters include the critical velocities of agglomeration $(v_{a,V,V^{\prime },crit})$ and breakage $(v_{b,V,V^{\prime },crit})$ appearing in (2.11)–(2.12), the collision cylinder diameter $d_{V,V^{\prime }}$ (2.8)–(2.9) and the daughter distributions $b(V\mid T_{s},V^{\prime },V^{\prime \prime })$ appearing in (2.6). The closures for the population balance presented in this section are dependent on $T_{s}$ , $\boldsymbol{u}_{\boldsymbol{V}}$ and $\boldsymbol{u}_{\boldsymbol{V}^{\prime }}$ , and thus appropriate balances are needed to determine these flow-field variables, as detailed below.

2.3 Momentum and granular energy balances

The solids-phase momentum balance for monodisperse inelastic spheres takes the general form (Gidaspow Reference Gidaspow1994)

where $P$ is the granular pressure, $\unicode[STIX]{x1D749}$ is the solids stress tensor and $\boldsymbol{F}_{body}$ is the body force vector.

Equation (2.34) describes the momentum balance of a monodisperse system. Since a distribution of agglomerate sizes is expected in a given system, a polydisperse kinetic theory (Garzó, Dufty & Hrenya Reference Garzó, Dufty and Hrenya2007a ; Garzó, Hrenya & Dufty Reference Garzó, Hrenya and Dufty2007b ; Hrenya Reference Hrenya2011) is expected to be most appropriate. Polydispersity is not relevant for the demonstration case presented here, however, as we target a lightly cohesive system in order to isolate the new aspects of our theory and thus more easily identify the source of any discrepancies. Accordingly, a monodisperse approximation is appropriate for our current purposes.

Similarly, the granular energy balance takes the following form:

where $-(PI+\unicode[STIX]{x1D749})\boldsymbol{ : }(\unicode[STIX]{x1D735}\boldsymbol{u})$ gives the generation of granular energy due to shear, $\boldsymbol{q}$ is the flux of granular energy and $\unicode[STIX]{x1D701}$ is the dissipation rate (all of which are dependent on the granular temperature and coefficient of restitution). It should be recalled that no fluid effects are given in (2.34) or (2.35), as the current focus is on a granular system, with the intent of isolating the new aspects of our theory.

To close the momentum and granular energy balances, closures are required for $P$ , $\unicode[STIX]{x1D749}$ , $\boldsymbol{q}$ and $\unicode[STIX]{x1D701}$ . These closures are generally obtained using an analogy with dense-gas kinetic theory – i.e. based on the Chapman–Enskog expansion of the Enskog equation, which is modified to account for the inelasticity of the particles (Hrenya Reference Hrenya2011). These closures depend on particle properties such as $e_{int}$ , as well as flow variables such as $T_{s}$ , $\unicode[STIX]{x1D700}_{s}$ and $\boldsymbol{u}$ . Explicit closures can be found in detail elsewhere (e.g. Garzó & Dufty (Reference Garzó and Dufty1999) for a monodisperse system). In a cohesive system, agglomerate shape is expected to have an effect on these closures, which are based on an assumption of spherical particles. Rigorous inclusion of these effects is outside the scope of this work, and is not relevant for the lightly cohesive approximately monodisperse system examined here. Accordingly, the focus of kinetic-theory modifications considered here will remain on the effects of cohesion on the coefficient of restitution.

2.4 Closure of kinetic-theory-based balances: effective restitution coefficient

As discussed earlier, collisions between cohesive particles are more dissipative than collisions between non-cohesive particles due to the added acceleration experienced by approaching particles with attractive forces. For non-cohesive particles, the (intrinsic) coefficient of restitution is defined as the absolute value of the ratio of the post-collision normal relative velocity $v_{n,f}$ and the pre-collision normal relative velocity $v_{n,0}$ of two colliding particles $(e_{int}=|v_{n,f}/v_{n,0}|)$ , and is an input to the kinetic-theory-based closures for the transport coefficients. For elastic particles, $e_{int}=1$ , and for inelastic particles, $0<e_{int}<1$ . For cohesive particles, the additional dissipation decreases the value of $|v_{n,f}/v_{n,0}|$ to some effective value, $e_{eff}$ , which depends on the velocity at which the particles collide (Davis et al. Reference Davis, Rager and Good2002; Donahue et al. Reference Donahue, Hrenya, Davis, Nakagawa, Zelinskaya and Joseph2010; Gollwitzer et al. Reference Gollwitzer, Rehberg, Kruelle and Huang2012; Abbasfard, Evans & Moreno-Atanasio Reference Abbasfard, Evans and Moreno-Atanasio2016); at infinitely high impact velocities, $e_{eff}\rightarrow e_{int}$ .

To account for cohesion, we seek to use an impact-velocity-dependent coefficient of restitution, ultimately cast in terms of its continuum counterpart $T_{s}$ . First, the dependence of $e_{eff}$ on $v_{n}$ can be determined from experiments or DEM simulations; an illustrative example is shown in the dotted line of figure 4 for the case of van der Waals forces. Aggregation occurs when $e_{eff}=0$ and rebound occurs when $e_{eff}>0$ . This dependence can be captured via a numerical fit of the data from experiments or the DEM. Alternatively, an analytical approximation based on the square-well model is used here for purposes of simplicity. In particular, an analytical prediction for $e_{eff}$ is available for cohesive hard spheres described by a square well (Weber et al. Reference Weber, Hoffman and Hrenya2004),

where $e_{int}$ is the intrinsic (non-cohesive) coefficient of restitution, assumed to be constant for this work, and $D_{well}$ is the well depth, given by Weber et al. (Reference Weber, Hoffman and Hrenya2004),

Figure 4 shows a comparison of $e_{eff}$ predicted by the square-well model with that observed in two-particle simulations (DEM) of particles subject to van der Waals cohesion. The results indicate that the analytical prediction of $e_{eff}$ based on the square-well model agrees quantitatively well with the DEM data. Further, the expression for $e_{eff}(v_{n})$ given by the square-well model (2.36) is expected to be extendible to other types of cohesion where the $e_{eff}$ versus $v_{n}$ plot takes on the same general shape – e.g. liquid bridging experiments (Davis et al. Reference Davis, Rager and Good2002; Weber et al. Reference Weber, Hoffman and Hrenya2004; Liu et al. Reference Liu, van Wachem, Mudde, Chen and van Ommen2016c ).

Figure 4.

Effective coefficient of restitution measured in DEM simulations of collisions between two particles with van der Waals forces described in § 3.2.1 (red dashed line) compared with the prediction from the square-well model ((2.36), blue solid line), as well as the intrinsic property value (black dashed line).

Rigorous incorporation of an impact-velocity-dependent coefficient of restitution, as illustrated in figure 4, into the derivation of the transport coefficients is non-trivial at best even for simple systems (Murphy & Subramaniam Reference Murphy and Subramaniam2015), and is outside the current scope. Another approach to account for the additional dissipation due to cohesion is by using a mean value of the effective coefficient of restitution $e_{eff}$ , or $\langle e_{eff}\rangle$ , in the transport coefficients and dissipation rate. This value of $\langle e_{eff}\rangle$ , like $e_{eff}$ itself, lies in the range $0<\langle e_{eff}\rangle <e_{int}$ . Whereas Murphy & Subramaniam (Reference Murphy and Subramaniam2015) used $e_{eff}$ in the dissipation term in the granular energy balance, here, $\langle e_{eff}\rangle$ is instead calculated and utilized in all of the kinetic-theory-based closures. The mean effective coefficient of restitution is calculated via

where $f(v_{n})$ is defined above in terms of $T_{s}$ (see (2.31)) and $e_{eff}(v_{n})$ is dictated by the specific type of cohesion under consideration. It should be noted that this integrand has no analytical solution and must be integrated numerically for each numerical cell in the domain of a given system.

4.1 Extraction of cohesion-specific quantities needed for continuum-theory closures: collision cylinder diameter and critical velocities

Three quantities were extracted from DEM simulations for use in the continuum theory: the collision cylinder diameter for singlet–doublet collisions, $d_{V,V^{\prime }}$ (via SSF), the critical velocity of aggregation for singlet–singlet collisions, $v_{a,1,1,crit}$ (via controlled two-particle collisions), and the critical velocity of breakage for singlet–doublet collisions, $v_{b,1,2,crit}$ (via SSF). The results for these cohesion-specific quantities are summarized in table 2. For the lightly cohesive system explored here, doublet–doublet collisions were neglected because the low concentration of doublets led to a very low frequency of doublet–doublet collisions relative to singlet–doublet collisions, thus making the breakage of doublets due to these types of collisions negligible.

The quantity $d_{V,V^{\prime }}$ was determined by averaging the distance between the centre of the doublet and the centre of the singlet for collisions of singlets and doublets detected in SSF simulations. The quantity $v_{a,1,1,crit}$ was measured using controlled two-particle DEM simulations (see figure 1).

Table 2.

Cohesion-specific quantities extracted from DEM SSF simulations and used as inputs to continuum SSF and HCS simulations.

As described in § 3.2.2, $v_{b,1,2,crit}$ was determined from the measured success factors of breakage $\unicode[STIX]{x1D713}_{b}$ in DEM SSF and using the closure for the success factor of breakage given in (2.12). The value of $v_{b,1,2,crit}$ shows a slight downward trend as the shear rate $\unicode[STIX]{x1D6FE}^{\ast }$ increases; see figure 8. A plausible explanation for this observation is the following: as the shear rate increases, agglomerate rotation increases, thus making it easier for agglomerates to break (the rotational motion adds to the normal impact velocity due to translational motion alone). Additional simulations for different values of $v_{b,1,2,crit}$ across the range measured in figure 8 were performed, and $\unicode[STIX]{x1D713}_{b,1,2}$ was found to be insensitive to $v_{b,1,2,crit}$ in this range, so use of a mean value in the continuum theory was deemed to be acceptable.

4.2 Verification and validation of continuum theory: SSF

Simple shear flow was first used to validate the closure for the success factor of aggregation $\unicode[STIX]{x1D713}_{a}$ (2.11) and then to validate the closure for $\langle e_{eff}\rangle$ (2.38). As detailed below, the corresponding comparisons with quantities extracted from DEM simulations constitute validation (theory, or closure, accuracy) rather than verification (code accuracy) (Grace & Taghipour Reference Grace and Taghipour2004) since $\unicode[STIX]{x1D713}_{a}$ and $\langle e_{eff}\rangle$ depend only on the critical velocity of aggregation $v_{a,1,1,crit}$ determined from DEM simulations of controlled two-particle collisions (rather than SSF). Next, SSF was used for verification (code accuracy) of the aggregation and breakage kernels (2.8)–(2.9). This last comparison is verification (code accuracy) rather than validation (theory accuracy), as the cohesion-specific inputs of the critical velocity of breakage $v_{b,1,2,crit}$ and the collision cylinder diameter $d_{V,V^{\prime }}$ , which were determined from the SSF system, were subsequently used in the continuum-theory predictions for SSF. Further, these latter verifications utilized $n_{1}$ , $n_{2}$ and $T_{s}$ extracted from the DEM simulations, rather than the predicted values, to ensure that any resulting discrepancies were traceable directly to the aggregation and breakage kernels themselves.

Figure 8.

Critical velocity of breakage calculated from DEM SSF simulations with $\unicode[STIX]{x1D700}_{s}=0.1$ (black points). The dashed line shows the mean value used in continuum theory.

It should be recalled that the success factor $\unicode[STIX]{x1D713}_{a,1,1}$ is used to determine the rate of aggregation in the system (2.3), and is therefore a key quantity for the accurate prediction of the agglomerate size distributions $n_{1}$ and $n_{2}$ . For simple shear flow, $\unicode[STIX]{x1D713}_{a,1,1}$ was extracted directly from the DEM, and compared with the current theory (2.11), as well as a previous theoretical expression (Liu Reference Liu2011). Both of these theoretical expressions depend on $v_{a,1,1,crit}$ and $T_{s}$ . It should be recalled that the value of $v_{a,1,1,crit}$ for these predictions was determined from controlled two-particle DEM simulations, and the value of $T_{s}$ used in the prediction was extracted directly from the DEM SSF simulations in order to isolate any discrepancies as stemming from the closure for $\unicode[STIX]{x1D713}_{a,1,1}$ , rather than uncertainty in $T_{s}$ . Figure 9 demonstrates that the new theory outlined here correctly predicts the success factor measured via the DEM, whereas the previous theory of Liu (Reference Liu2011) significantly underpredicts $\unicode[STIX]{x1D713}_{a,1,1}$ . In the current theory, $\unicode[STIX]{x1D713}_{a,1,1}$ is determined from the distribution of normal impact velocities. The disagreement between $\unicode[STIX]{x1D713}_{a,1,1}$ predicted by the Liu theory (Reference Liu2011) and the DEM simulation results is not surprising since the value of $\unicode[STIX]{x1D713}_{a,1,1}$ in the Liu theory (Reference Liu2011) is determined from the distribution of total impact velocities (including tangential components). These results are consistent with figure 3. That is, at lower impact velocities, the probability of an impact having a given total relative velocity as predicted by the Liu theory (Reference Liu2011) is much smaller than the probability of a collision having the same given normal relative impact velocity as predicted by the current theory. Integration of these impact velocity distributions from 0 to a constant $v_{a,1,1,crit}$ to determine $\unicode[STIX]{x1D713}_{a,1,1}$ (2.11) causes the Liu theory (Reference Liu2011) to underpredict $\unicode[STIX]{x1D713}_{a,1,1}$ , as evidenced in figure 9. It should also be noted that for the limiting case of aggregation between singlets, the prediction of $\unicode[STIX]{x1D713}_{a,1,1}$ given by the current model is consistent with the ratio of the frequency of collisions that successfully result in agglomeration to the total collision frequency predicted by Murphy & Subramaniam (Reference Murphy and Subramaniam2015).

Figure 9.

Success factor of aggregation for $\unicode[STIX]{x1D700}_{s}=0.1$ . The solid line shows the prediction by Liu (Reference Liu2011).

Another closure that can be validated directly via SSF simulations is the kinetic-theory treatment involving $\langle e_{eff}\rangle$ . For SSF, the momentum balances are solved identically such that only the energy balance and population balance require solution. It should be recalled that the granular energy balance has one-way coupling with the population balance, such that the granular energy balance affects the population balance, but not vice versa. Therefore, the granular energy balance can be solved in isolation for $T_{s}$ , and any discrepancies between the predicted value and that extracted from the DEM will arise from the kinetic-theory treatment itself (i.e.  $\langle e_{eff}\rangle$ ), and not due to the new closures in the population balance. A comparison between the $T_{s}$ measured in the DEM and that predicted by the granular energy balance is given in figure 10. The kinetic theory for a non-cohesive system (using the intrinsic property for the coefficient of restitution $e_{int}$ in the dissipation and shear viscosity closures in the granular energy balance (3.14)) is included for comparison with the current (cohesive) theory and DEM simulations. The kinetic theory using $e_{int}$ gives an overestimation of the granular temperature compared with when $\langle e_{eff}\rangle$ is used. The excellent agreement between the DEM results and the kinetic-theory treatment using $\langle e_{eff}\rangle$ provides initial validation for this approach, and also provides support for the monodisperse assumption. Further, this approach gives an excellent prediction of the $T_{s}$ for a system outside of that explored by Murphy & Subramaniam (Reference Murphy and Subramaniam2015), which included the effect of $e_{eff}$ (rather than $\langle e_{eff}\rangle$ used here) in the dissipation term only, and was thus restricted to the HCS. Further, the approach used by Murphy & Subramaniam (Reference Murphy and Subramaniam2015) does not include the shear term (i.e. granular shear viscosity $\unicode[STIX]{x1D707}$ , which also depends on $e_{eff}$ ) included here.

Figure 10.

Comparison of the non-dimensional granular temperature $(T_{s}^{\ast }=T_{s}/(\unicode[STIX]{x1D70C}d_{1}^{3}v_{a,1,1,crit}^{2}))$ predicted in SSF by the continuum theory and the DEM for various solids volume fractions and shear rates. The continuum predictions are given by solid lines. The dashed line corresponds to a simulation using the intrinsic, rather than effective, coefficient of restitution. Points represent DEM data.

Figure 11.

Comparison of the fraction of particles in doublets predicted in SSF by the continuum theory (solid lines) and the DEM (points) for various solids volume fractions and shear rates. The error bars show the standard deviation of the fraction of particles at incremental times after steady state was reached for one simulation.

Figure 12.

Terms associated with aggregation and breakage rates for $\unicode[STIX]{x1D700}_{s}=0.1$ for (a–c) singlet–singlet interactions and (d–f) singlet–doublet interactions. Collision frequencies and aggregation and breakage rates are non-dimensionalized by multiplying by $d_{1}^{4}/v_{a,1,1,crit}$ .

Above, the closures of the new continuum theory that can be validated using SSF ( $\unicode[STIX]{x1D713}_{a}$ and $\langle e_{eff}\rangle$ ) are compared successfully with the DEM simulations (theory accuracy); we will now turn our attention to the remaining closures for which verification (code accuracy) can be carried out in SSF, and then move on to the HCS for validation of the same quantities. In particular, for verification purposes, the fraction of particles in doublets $(2n_{2}/n_{particles})$ predicted by the continuum theory (solution of both the kinetic theory and population balances) for SSF is compared with the DEM results. This comparison is verification since the critical velocity of breakage $v_{b,1,2,crit}$ was determined from the SSF DEM simulations, and this value is then used in the calculation of the fraction of particles in doublets (2.1). The resulting comparisons are shown in figure 11 for a range of concentrations and shear rates. The DEM simulations display a trend of a decreasing fraction of particles in doublets with an increasing shear rate. The same behaviour is captured by the continuum theory. This behaviour is expected, as an increased shear rate leads to an increased (steady state) $T_{s}$ . An increased $T_{s}$ (i.e. impact velocity) decreases the likelihood of agglomeration ( $\unicode[STIX]{x1D713}_{a,1,1}$ ) and increases the likelihood of breakage ( $\unicode[STIX]{x1D713}_{b,1,2}$ ), leading to less aggregation and more breakage, and thus fewer doublets.

Quantitatively, a non-systematic deviation exists between the fraction of particles in doublets observed in the DEM and the continuum-theory prediction (figure 11). More specifically, the continuum theory overpredicts the amount of aggregation in some cases (for all shear rates for $\unicode[STIX]{x1D700}_{s}=0.05$ , for most shear rates at $\unicode[STIX]{x1D700}_{s}=0.1{-}0.2$ ) and underpredicts for other conditions (the lowest shear rates at $\unicode[STIX]{x1D700}_{s}=0.1{-}0.25$ and for all shear rates at $\unicode[STIX]{x1D700}_{s}=0.3$ ). To identify the source of this deviation, the individual terms in the population balance are considered. In particular, the success factors of aggregation and breakage, the collision frequencies of singlet–singlet and singlet–doublet collisions, and the aggregation and breakage rates are considered individually in figure 12 for the case of $\unicode[STIX]{x1D700}_{s}=0.1$ . Namely, the product of the success factors (a,d) with the collision frequency (b,e) gives the rates of aggregation and breakage (c,f). Panels (a–c) represent singlet–singlet collisions and (d–f) represent singlet–doublet collisions.

The results for the different terms in the population balance shown in figure 12 explain the overprediction of the fraction of particles in doublets observed in figure 11 for $\unicode[STIX]{x1D700}_{s}=0.1$ . More specifically, figure 12(a) demonstrates that the continuum theory overpredicts the success factor $\unicode[STIX]{x1D713}_{a,1,1}$ (except at shear rates greater than $\unicode[STIX]{x1D6FE}^{\ast }\sim 1.6$ ) and consequently the singlet–singlet aggregation rate in figure 12(c). Further, the continuum theory slightly underpredicts the collision frequency of singlet–doublet collisions $N_{1,2}$ (figure 12 e) for all shear rates, which causes a small underprediction in the breakage rate of doublets (figure 12 f). Collectively, these effects lead to an overprediction of particles in doublets for higher shear rates (figure 11 b), which explains the deviation between the continuum theory and the DEM observed for shear rates greater than $\unicode[STIX]{x1D6FE}^{\ast }\sim 0.98$ in figure 11(b). In figure 11(b), however, the continuum theory appears to accurately predict the fraction of particles in doublets at the lowest shear rates of $\unicode[STIX]{x1D6FE}^{\ast }\sim 0.87{-}0.98$ . One possible explanation for the absence of the expected mismatch is that the continuum theory is missing a sink of doublets, which would counteract the overprediction expected from figure 12(c,f).

To explain the discrepancy noted above, it is hypothesized that the birth of larger agglomerates (larger than doublets) due to aggregation may provide the missing net sink of doublets. Indeed, figure 13 indicates that larger agglomerates are present, especially at lower shear rates. More specifically, the doublet fraction in the DEM was observed to remain at a statistical steady state in time, requiring that the net source of doublets be zero. However, calculations confirm that an estimate of the net source of doublets observed in the DEM $(\unicode[STIX]{x1D700}_{s}=0.1)$ , based only on the aggregation rate of singlets less the breakage rate of doublets, is indeed non-zero for $\unicode[STIX]{x1D6FE}^{\ast }<0.2$ (not shown for the sake of brevity). Therefore, agglomerates involving three particles or more, which were ignored in the application of the continuum theory to these lightly cohesive systems, are needed for the net source of doublets at steady state to equal zero at such low shear rates. Moreover, the results shown in figures 9–11 have better agreement at higher shear rates (higher $T_{s}$ ) due to the lack of larger agglomerates (figure 13). Accordingly, for purposes of verification, attention should be restricted to the higher shear rates of $\unicode[STIX]{x1D6FE}^{\ast }>0.98$ , since the continuum theory was applied assuming that only singlets and doublets are present (i.e. population balances for $n_{3}$ , etc., which are beyond the scope of this work, were not included). Within this range, the results in figures 11 and 12 provide sound verification of the aggregation and breakage kernels contained in the new continuum theory.

Figure 13.

Distribution of larger agglomerates in DEM SSF for $\unicode[STIX]{x1D700}_{s}=0.1$ . The error bars show the standard deviation of the fraction of particles in agglomerates for different times after steady state was reached in one simulation.

Figure 14.

Granular temperature in the HCS predicted by the current work and Murphy & Subramaniam (Reference Murphy and Subramaniam2015) compared with the DEM.

4.3 Validation of continuum theory: HCS

As detailed above, SSF provides a simple steady-state validation of the use of $\langle e_{eff}\rangle$ in kinetic-theory closures as well as the success factor $\unicode[STIX]{x1D713}_{a}$ appearing in the population balance. The HCS, by contrast, provides a simple transient validation case for the terms mentioned above, as well as the aggregation and breakage kernels ( $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FC}$ respectively) in the population balance. First, the granular temperature $T_{s}$ predicted by the continuum theory for the HCS is compared with $T_{s}$ measured in the DEM, as well as the prediction given by Murphy & Subramaniam (Reference Murphy and Subramaniam2015), in figure 14. The results are all in very good agreement, providing initial validation of the kinetic-theory treatment using $\langle e_{eff}\rangle$ . The slight deviation of the granular temperature predicted by Murphy & Subramaniam’s (Reference Murphy and Subramaniam2015) theory from the DEM results and the new theory presented here is associated with the fourth moment of $f^{(0)}$ , or equivalently the $c^{\ast }$ constant in the theory of Garzó & Dufty (Reference Garzó and Dufty1999). This constant takes into account the deviation of elastic particles from an inelastic gas. Such effects are inherent in the DEM simulations, and are taken into account in the new theory through the use of Garzó and Dufty’s constitutive relations (1999), but were not taken into account in the work of Murphy & Subramaniam (Reference Murphy and Subramaniam2015).

It should be noted that the model given by Murphy & Subramaniam (Reference Murphy and Subramaniam2015) does not incorporate the effects of $e_{eff}$ into the terms other than the dissipation rate of granular energy, and therefore can only be applied to the HCS. Again, it should be noted that there is one-way coupling between the granular energy balance and the population balance, such that the granular energy balance affects the population balance, but not vice versa. Accordingly, the excellent predictions for $T_{s}$ resulting from the $\langle e_{eff}\rangle$ treatment indicate that any potential discrepancy in predictions for $n_{1}$ and $n_{2}$ can be traced to the population balance, and do not stem from the kinetic-theory treatment.

Finally, the fraction of particles in agglomerates in the HCS predicted by the continuum theory is compared with the DEM data in figure 15. The fraction of particles in agglomerates is observed to increase in time. More specifically, as time increases, $T_{s}$ decreases due to dissipation from inelastic collisions (figure 14). As $T_{s}$ decreases, the success factor of agglomeration $\unicode[STIX]{x1D713}_{a}$ increases (due to lower impact velocities) while the collision frequency $N_{1,1}$ decreases. The former leads to an increase in doublet formation, while the latter has the opposite effect. Thus, at earlier times when $T_{s}$ is high, the higher rate of doublet formation can be traced to the high collision frequency despite the lower success factor. As $T_{s}$ decreases with time, the reduced collision frequency leads to a lower rate of doublet formation (smaller slope), although doublet formation still increases with time since the likelihood of agglomeration is greater for smaller impact velocities. The continuum theory captures the inflection at $t^{\ast }\sim 0.2{-}0.5$ caused by these competing mechanisms of aggregation.

Figure 15.

Fraction of particles in doublets for the HCS with non-dimensional time, $t^{\ast }=tv_{a,crit,1,1}/d_{1}$ . The lines show the continuum prediction, the points show the DEM results and the error bars show the standard deviation of five repeats.