Non-Fickian macroscopic model of axial diffusion of granular materials in a long cylindrical tumbler

TL;DR

This paper introduces a non-Fickian macroscopic model for axial diffusion in granular materials within a cylindrical tumbler, accounting for particle collision-induced diffusion in a localized surface layer, with exact solutions and comparisons to classical models.

Contribution

It presents a novel non-Fickian diffusion model for granular flow in tumblers, linking parameters to measurable quantities and deriving exact solutions.

Findings

The model captures diffusion with memory effects.

Exact solutions are derived for finite-width pulse initial conditions.

The non-Fickian model differs from classical diffusion in granular flows.

Abstract

We provide new concepts for understanding transport phenomena in flows of granular materials by using a non-Fickian macroscopic model of axial diffusion of a granular material in a finite cylindrical tumbler. The model accounts for diffusion induced by particle collisions only in a thin surface flowing layer due to localization of shear within the cross-section of the drum. All model parameters are related to measurable quantities in a granular flow. It is shown that the proposed model is a member of the general class of linear constitutive relations with memory. An exact solution for the spreading of a finite-width pulse initial condition under the proposed non-Fickian model is derived and compared to the solution of a Fickian model (i.e., the ``classical'' diffusion equation).

Equations

$$\frac{\partial c_{\mathrm{d}}}{\partial t}$$

$$\frac{\partial c_{\mathrm{nd}}}{\partial t}$$

$$c_{\mathrm{d}}(z,0)=c_{\mathrm{d},0}(z),c_{\mathrm{nd}}(z,0)=c_{\mathrm{nd},0}(z),$$

$${\frac{\partial c_{\mathrm{d}}}{\partial z}|}_{z=\pm L}={\frac{\partial c_{\mathrm{nd}}}{\partial z}|}_{z=\pm L}=0.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\int_{-L}^{+L}\big{(}c_{\mathrm{d}}+c_{\mathrm{nd}}\big{)}\,\mathrm{d}z=\left.D\frac{\partial c_{\mathrm{d}}}{\partial z}\right|_{-L}^{+L}=0.$$

$$k(1+\beta )\frac{\partial c_{\mathrm{d}}}{\partial t}+\frac{{\partial }^2{c}_{\mathrm{d}}}{\partial t^2}=D\left(k\beta \frac{{\partial }^2{c}_{\mathrm{d}}}{\partial z^2}+\frac{{\partial }^3{c}_{\mathrm{d}}}{\partial t\partial z^2}\right).$$

$$\underset{t\to \infty }{\lim }{c}_{\mathrm{d}}(z,t)$$

$$\underset{t\to \infty }{\lim }{c}_{\mathrm{nd}}(z,t)$$

$$\frac{\frac{1}{2}\pi {\delta }_{0}R}{\varphi \pi R^2}=\frac{1}{2\varphi }\left(\frac{{\delta }_0}{R}\right).$$

$$\beta =\frac{{\delta }_0}{2\varphi R-{\delta }_0}\simeq \frac{1}{2\varphi }\left(\frac{{\delta }_0}{R}\right)({\delta }_0\ll R).$$

$$k\simeq \frac{\omega }{\pi }\varphi {\left(\frac{{\delta }_0}{R}\right)}^{-1}.$$

$$(1+\beta )\frac{\partial C_{\mathrm{d}}}{\partial T}+\mathcal{D}\frac{{\partial }^2{C}_{\mathrm{d}}}{\partial T^2}=\beta \frac{{\partial }^2{C}_{\mathrm{d}}}{\partial Z^2}+\mathcal{D}\frac{{\partial }^3{C}_{\mathrm{d}}}{\partial T\partial Z^2},$$

$${\tilde{C}}_{\mathrm{d}}(Z,0)$$

$$\frac{\partial {\tilde{C}}_{\mathrm{d}}}{\partial T}(Z,0)$$

$$\mathcal{D}\frac{{\mathrm{d}}^2{a}_n}{\mathrm{d}{T}^2}+\left[(1+\beta )+n^2{\pi }^2\mathcal{D}\right]\frac{\mathrm{d}{a}_n}{\mathrm{d}T}+n^2{\pi }^2\beta a_n=0.$$

$$a_n(0)$$

$${\frac{\mathrm{d}{a}_n}{\mathrm{d}T}|}_{T=0}$$

$$\overline{a_n}(s)=\frac{a_n(0)[s+(1+\beta )/\mathcal{D}+n^2{\pi }^2(1-\mathcal{D})]}{\left\{s^2+s[(1+\beta )/\mathcal{D}+n^2{\pi }^2]+n^2{\pi }^2\beta /\mathcal{D}\right\}}.$$

$$\begin{gathered} a_n(T)=\mathrm{sinc}(n\pi \ell )H(T)\exp \left(-{{\rm Y}}_{n}T\right)[\cosh \left(\frac{T\sqrt{{\Delta }_n}}{2\mathcal{D}}\right) \\ +\frac{(1+\beta )+n^2{\pi }^2\mathcal{D}(1-2\mathcal{D})}{\sqrt{{\Delta }_n}}\sinh \left(\frac{T\sqrt{{\Delta }_n}}{2\mathcal{D}}\right)], \end{gathered}$$

$${\rm Y}(\xi )$$

$$\Delta (\xi )$$

$$C_{\mathrm{d}}(Z,T)=\frac{\beta }{1+\beta }\left[\frac{1}{2}H(T)+\sum _{n=1}^{\infty }{a}_n(T)\cos (n\pi Z)\right],$$

$$\begin{gathered} C_{\mathrm{d}}(Z,T;\mathcal{D}=0)=\frac{\beta }{1+\beta }H(T) \\ \times \left[\frac{1}{2}+\sum _{n=1}^{\infty }{\mathrm{e}}^{-n^2{\pi }^2{D}_{\mathrm{F}}T}\mathrm{sinc}(n\pi \ell )\cos (n\pi Z)\right], \end{gathered}$$

$$\left[1+\frac{1}{k(1+\beta )}\frac{\partial }{\partial t}\right]q_{\mathrm{d}}=-D\frac{\beta }{1+\beta }\left[1+\frac{1}{k\beta }\frac{\partial }{\partial t}\right]\frac{\partial c_{\mathrm{d}}}{\partial z},$$

$$q_{\mathrm{d}}(z,t)=-{\int }_0^t\mathfrak{K}(t-t^{\prime })\frac{\partial c_{\mathrm{d}}}{\partial z}(z,t^{\prime })\mathrm{d}{t}^{\prime },$$

$$\mathfrak{K}(t-t^{\prime })=D\left[2\delta (t-t^{\prime })-k{\mathrm{e}}^{-k(1+\beta )(t-t^{\prime })}\right].$$

$$\mathfrak{K}(t-t^{\prime })=-D_{\alpha }\frac{1-\alpha }{\Gamma (\alpha )}(t-t^{\prime }{)}^{\alpha -2},$$

$$q_{\mathrm{d}}=-D_{\alpha }\text{prescript}0{\mathfrak{D}}_t^{1-\alpha }\frac{\partial c_{\mathrm{d}}}{\partial z}$$

$$\text{prescript}0{\mathfrak{D}}_t^{1-\alpha }\Xi (\cdots ,t)$$

$$=-\frac{1-\alpha }{\Gamma (\alpha )}{\int }_0^t\frac{\Xi (\cdots ,t^{\prime })}{(t-t^{\prime }{)}^{2-\alpha }}\mathrm{d}{t}^{\prime }$$

$$\frac{\partial c_{\mathrm{d}}}{\partial t}=D_{\alpha }\text{prescript}0{\mathfrak{D}}_t^{1-\alpha }\frac{{\partial }^2{c}_{\mathrm{d}}}{\partial z^2}.$$

Full text

1 Introduction

Granular materials at rest, in flow or under external agitation behave like various continua (liquid, solid, gases), albeit with quite “unusual” properties [32, 1]. The continuum approach has allowed successful practical modeling of granular mixing and segregation [50, 44, 63]. Yet, “[d]espite an expanding body of literature in the last ten to fifteen years, relationships between identified mechanisms are ambiguous, experimental data is scarce, and there is no accepted model[s],” as Rosato and Blackmore note in the preface to the 1999 IUTAM Symposium on Segregation in Granular Flows [52]. Arguably, this statement remains true today, even after two more decades of progress in granular mechanics.

Diffusion is a fundamental powder mixing (transport) mechanism modeled using continuum theories. A number of agitation processes lead to non-equilibrium velocity fluctuations that cause “self-diffusion” in granular flows: from tumblers [66] and vertically vibrated beds [53] (axial diffusion) to “simple” shear flows (anisotropic diffusion) [58, 49, 11] to horizontal [9] and vertical [29] Hele-Shaw cells (radial diffusion) to cylindrical Couette cells (axial and radial diffusion) [64, 26]. Then, axial mixing of granular materials in a rotating drum, as shown schematically in Fig. 1, is typically modeled using Fick’s second law, leading to the “classical” diffusion equation [38, 10, 31, 22]. For example, Cahn et al. [10] provided early experimental justification for a “diffusive” axial transport process by tracking the spread of yellow beads into identical white beads, of diameter $d_{p}=200$ $\mu$m, in a cylindrical drum of radius $R=0.1$ m (note that $d_{p}\ll R$ is required for this continuum description to hold).

Granular flow in a long drum has proven to be a simple but important system to understand [59]. Experiments in the 2000s [37, 25] reinvigorated interest in this flow because of the possibility of anomalous diffusion. Although “anomalous” scalings can also be interpreted within the theory of intermediate self-similar asymptotics of “classical” diffusion [19], the fundamentals of granular diffusion are still poorly understood. Macroscopic models based on Fick’s laws are not derived from first principles but rather fit to experimental or particle dynamics (i.e., simulation) data [61, 42]. Nevertheless, experiments have established the key features of granular flow in rotating containers [12, 8, 36, 33]. The most important conclusion drawn is that the flow in an axial cross-section of a long rotating drum (Fig. 1) is composed of a thin surface shear layer (the flowing layer) in which material flows quickly down the slope. Particles below it (the fixed bed) are in static equilibrium performing solid body rotation in unison with the container. Acknowledging these disparate flow characteristics in the cross-section, Das Gupta et al. [21] related the axial diffusivity and drift velocities in a tilted tumbler to its rotation rate, tilt angle, and the particle properties. Their model starts from the stochastic equations for a single particle and necessitates the prescription of a probability distribution of axial displacements due to interparticle collisions.

As Metcalfe and Shattuck note, “[i]t takes more than 200 revolutions to disperse the seeds … axially along the tube boundary and then mix radially. This is physically plausible given that the only motions mixing the material are the avalanches across the free surface” [46]. Thus, if the surface flow in the tumbler is the only mechanism that leads to transport and mixing (as also emphasized in Fig. 1(c,d) in [10]), then the axial spread of particles should be expected to be more complicated than the Fickian diffusion models originally proposed in [31, 21]. Specifically, the fluctuations in the axial direction are distinct from the random displacements that lead to radial mixing in the cross-section. Given that the latter dynamics are rapidly equilibrated [36, 59], the radial transport mechanisms in the tumbler are not of interest in this work. Instead, here, we would like to address a basic scientific question about the axial transport: Does taking into account the disparate flow characteristics in the tumbler cross-section lead to a non-Fickian (but not necessarily “anomalous”) axial diffusion equation? Our goal is to provide new theoretical insights into transport phenomena in flows of granular materials.

2 A macroscopic axial diffusion model accounting for localization of shear in the cross-section

Let $c(z,t)=\frac{1}{\phi A}\iint_{A}\hat{c}(r,\theta,z,t)\,\mathrm{d}A$ be the cross-sectionally averaged concentration in a partially filled cylinder with cross-sectional area $A$ and fill fraction $\phi$; $\hat{c}$ is the number of particles per unit volume. In the cross-section of the tumbler, particles continuously exchange between the surface flowing layer (in which collisions lead to random axial displacement and, thus, diffusion) and a fixed bed in which particles only rotate with the tumbler and cannot be displaced axially. The granular material is assumed to be monodisperse. Thus, to study the axial transport process, we label a certain proportion of the particles as “diffusing,” $c_{\mathrm{d}}$, and the remainder as “non-diffusing,” $c_{\mathrm{nd}}$; $c=c_{\mathrm{d}}+c_{\mathrm{nd}}$. Within the cross-sectionally averaged description, we envision that the diffusing “species” concentration, $c_{\mathrm{d}}$, can represent tagged particles being observed to spread axially. Consequently, in any axial cross-section along the tumbler, there can be both diffusing and non-diffusing particles, and their concentrations will vary in space, $z$, and time, $t$. Now, we wish to develop a transport model for the spatiotemporal evolution of $c_{\mathrm{d}}$ and $c_{\mathrm{nd}}$, accounting for the flow features thus described.

A “minimal” model of such a diffusion process (see, e.g., [45, Sec. 3.3.1]) is

[TABLE]

subject to given initial conditions (ICs) on the species:

[TABLE]

and no-flux boundary conditions (BCs) at the tumbler’s endwalls:

[TABLE]

In Eq. (1a), $D$ is an “effective” diffusivity [dimensions of (length)2/time] characterizing the mean-squared axial displacements caused by particle collisions in the flowing layer; $k$ is a rate constant (dimensions of 1/time) characterizing the flowing layer–fixed bed particle exchange using first-order kinetics; and, the dimensionless parameter $\beta$, which is $>0$ and $\neq 1$, is a partition ratio (or coefficient) related to the fact that the flowing layer and the fixed bed take up unequal areas of the cross-section. Another way to interpret $\beta$ from the mathematical structure of the model, is that there is lack of local conservation (in a given axial cross-section) between the two species, due to the fact that the diffusing species can leave the given cross-section (while the non-diffusing species cannot), which leads to unequal exchange rates (i.e., $k\neq\beta k$) between diffusing the non-diffusing species. Nevertheless, a global conservation law follows from summing Eqs. (1), integrating over space, and using the BCs from Eq. (3):

[TABLE]

More generally, Eqs. (1) represent a reaction-diffusion system in which there is local immobilization of one species [55]. Such equations also arise in the modeling of signaling and transport in certain biological systems [4, 20] wherein proteins can become “mobile” (diffusing) or “immobile” (non-diffusing) as they bind or unbind from a local substrate. A similar situation occurs in electrodeposition [40, 3], during flows in fractured reservoirs [2] and dusty gasses [56, 41], and during heat transfer between an electron gas and a metal lattice [62], among other examples (see also [14]).

Equations (1) can be written as a single partial differential equation (PDE) for $c_{\mathrm{d}}$ by taking $\partial/\partial t$ of Eq. (1a), eliminating $\partial c_{\mathrm{nd}}/\partial t$ using Eq. (1b), then eliminating $k(c_{\mathrm{d}}-\beta c_{\mathrm{nd}})$ using Eq. (1a):

[TABLE]

Observe that Eq. (5) features both second-order-in-time and mixed derivatives, thus some authors designate it as a dissipative wave equation [14]. Importantly, Eq. (5) cannot be derived from Fick’s first and second law, thus it is termed non-Fickian.

Next, we determine how the parameters $k$ and $\beta$ in Eqs. (1) relate to the physics of axial diffusion of a granular material in a tumbler (Fig. 1).

Determining $\beta$:

At steady-state ($t\to\infty$) we expect both spatial and temporal gradients to vanish for a monodisperse granular material. (Spatial gradients in concentration may exist in size-bidisperse systems even at steady state; see, e.g., [13].) Then, we may seek constant solutions to Eq. (1) and Eq. (3), and we find that

[TABLE]

where $c_{\mathrm{d},0}$ and $c_{\mathrm{nd},0}$ must now be constants given the assumption of no axial gradients. Note that $c_{\mathrm{d},0}+c_{\mathrm{nd},0}=1$ can always be enforced due to the linearity of Eqs. (1).

The cross-section of the cylinder in Fig. 1 is circular, hence the area occupied by particles is $A=\phi\pi R^{2}$. To a good approximation, the flowing layer’s shape is an ellipse with semi-major axis $R$ and semi-minor axis $\delta_{0}$ [44, §3]. The maximal depth of the flowing layer $\delta_{0}$ is typically on the order of $0.1R$ to $0.2R$ [36] (or, about 8 to 12 particles thick for beads a few mm in diameter in a $24$ cm-diameter drum [33]). Thus, the flowing layer’s fraction of the cross-sectional area is

[TABLE]

At steady state, the fractions of particles in the flowing layer, as given by Eq. (6a) and Eq. (7), must match, thus

[TABLE]

Determining $k$:

The constant $k$ quantifies the rate at which the diffusing species leaves the flowing layer. The axial diffusion process under consideration occurs on a longer time scale compared to the radial transport. This separation of time scales is supported by experiments, which have shown that radial transport equilibrium is rapidly established in a drum, and particles are exchanged between the flowing layer and the bulk in a time-periodic, or quasi-steady, manner [36, 59]. Now, the flux through the flowing layer’s interface with the bulk can be estimated as the area of the flowing layer, $\tfrac{1}{2}\pi\delta_{0}R$, times the rate-of-loss (or gain) of the diffusing species, $k$, or $Q_{\mathrm{d}}\simeq\tfrac{1}{2}\pi\delta_{0}Rk$.

At the same time, for partially filled tumblers ($\phi\lesssim 50\%$), all particles complete a transit of the cross-section in about half a revolution, specifically also entering (and leaving) the flowing layer in that time [50]. So, the flux can also be estimated as the total filled area, $A=\phi\pi R^{2}$, divided by the half period of rotation, $\pi/\omega$, yielding $Q\simeq\tfrac{1}{2}\omega\phi R^{2}$. For quasi-steady flow in the cross-section, we must have $Q=Q_{\mathrm{d}}$, hence we find that

[TABLE]

Again, typically $\delta_{0}\ll R$, so $k/\omega\gg 1$.

Nondimensionalization:

Consistent with prior literature, we use the tumbler’s axial half-length as the length scale and the axial diffusion time as the time scale: $z=LZ$, $t=(L^{2}/D)T$. Then, letting $C_{\mathrm{d}}(Z,T)=c_{\mathrm{d}}(z,t)$, Eq. (5) becomes

[TABLE]

Here, $\mathcal{D}=D/(kL^{2})=(\pi\delta_{0}D)/(\omega RL^{2})$ is a dimensionless parameter—a type of Damköhler number—that represents the ratio of the “reaction time scale” ($1/k$) to the diffusion time scale ($L^{2}/D$). We expect that $\mathcal{D}\ll 1$ since the reaction time scale is related to the exchange of particles between the flowing layer and the bulk, which occurs continuously as the tumbler rotates, while diffusion in the axial direction occurs more slowly over many revolutions. As an example, for particles with diffusivity on the order of $D\approx 1$ mm2/s [22, 57, 61] in a thin flowing layer with aspect ratio $\delta_{0}/R\approx 0.1$ in a tumbler of length $2L=600$ mm and radius $R=14.25$ mm rotating at $\omega=0.62$ rev/s $=1.24\pi$ rad/s [37], we obtain the rough estimate $\mathcal{D}\simeq 10^{-6}$.

Notice that, for $\mathcal{D}\to 0^{+}$ (e.g., $k\to\infty$, an infinitely fast flowing layer–fixed bed exchange), Eq. (10) becomes the Fickian “classical” diffusion equation postulated in previous works [38, 31, 21] but the diffusivity is multiplied by $\beta/(1+\beta)\sim\delta_{0}/R$ (for $\delta_{0}\ll R$). Therefore, when fitting solutions of the Fickian diffusion equation to data, the diffusivity may be overestimated by at least an order of magnitude, depending on the conditions under which the data is taken.

Next, we solve the axial granular diffusion initial-boundary value problem (IBVP) for Eq. (10). Specifically, we seek to highlight the effect of $\mathcal{D}>0$ on the transient axial diffusion process. As argued in [19], the time required for the concentration profile to relax away from a finite-width pulse IC can be one reason for observing “anomalous” scalings.

3 The initial-boundary value problem on $Z\in[-1,+1]$

First, we must supplement Eq. (10) with appropriate ICs, say $C_{\mathrm{d}}(z,0)=C_{\mathrm{d},0}(Z)$ and $(\partial C_{\mathrm{d}}/\partial T)_{T=0}=C_{\mathrm{d},1}(Z)$. $C_{\mathrm{d},1}(Z)$ is not arbitrary because it must satisfy a compatibility condition based on Eqs. (1). Specifically, from Eq. (1a), $C_{\mathrm{d},1}(Z)$ can be related to the non-diffusing species’ initial condition $C_{\mathrm{nd},0}(Z)$ [recall Eq. (2)] as $C_{\mathrm{d},1}(Z)=\mathcal{D}\mathrm{d}^{2}C_{\mathrm{d},0}(Z)/\mathrm{d}Z^{2}-C_{\mathrm{d},0}(Z)+\beta C_{\mathrm{nd},0}(Z)$. Now, for the problem at hand, suppose an equilibrium partition of diffusing and non-diffusing species in the cross-section at $T=0$, in which case $\beta{C}_{\mathrm{nd}}(Z,0)={C}_{\mathrm{d}}(Z,0)$. Then, a band of unit area of diffusing particles is tagged, and we seek to determine its evolution. To this end, normalize $C_{\mathrm{d}}$ as $\tilde{C}_{\mathrm{d}}(Z,T)=C_{\mathrm{d}}(Z,T)/C_{\mathrm{d},\infty}$ where $C_{\mathrm{d},\infty}=\beta/(1+\beta)$ is the constant steady-state distribution from Eq. (6a). The ICs for Eq. (10) then become

[TABLE]

where $2\ell<2$ is the band’s width ($\ell=l/L$ is dimensionless), $H(\cdot)$ is the Heaviside unit step function, $\delta(\cdot)$ is the Dirac-delta distribution and $\delta^{\prime}\equiv\mathrm{d}\delta/\mathrm{d}Z$. The BCs given in Eq. (3) carry over to $\tilde{C}_{\mathrm{d}}$.

Next by separation of variables using the eigenfunctions of $\mathrm{d}^{2}/\mathrm{d}Z^{2}$ satisfying homogeneous Neumann BCs, the solution can be written as a Fourier cosine series: $\tilde{C}_{\mathrm{d}}(Z,T)=\tfrac{1}{2}a_{0}(T)+\sum_{n=1}^{\infty}a_{n}(T)\cos(n\pi Z)$. Substituting into Eq. (10) and employing orthogonality of the eigenfunctions:

[TABLE]

The ICs from Eqs. (11) give

[TABLE]

for $n\geq 0$. Here, we have defined $\operatorname{sinc}\eta=\eta^{-1}\sin\eta$ for real $\eta\neq 0$ and $\operatorname{sinc}(0)=1$. Applying the temporal Laplace transform, $\overline{(\cdot)}=\int_{0}^{\infty}(\cdot)\,\mathrm{e}^{-sT}\,\mathrm{d}T$ ($s\in\mathbb{C}$), to Eq. (12), we obtain, for $n\geq 0$,

[TABLE]

Following [17], Eq. (14) can be inverted back to the $t$ domain by partial fractions and tables of Laplace transforms, to yield $a_{0}(t)=H(t)$, and

[TABLE]

for $n>0$, where we have defined $\Delta_{n}:=\Delta(n\pi)$ and $\Upsilon_{n}:=\Upsilon(n\pi)$, with

[TABLE]

for convenience. In general, we must consider the three cases $\Delta_{n}\gtreqqless 0$. On physical grounds, we expect $\beta<1$ (indeed, $\beta\ll 1$) as discussed in Sec. 2. Then, it can be shown $\Delta(\xi)>0$ $\forall\xi$ and $\mathcal{D}$, thus we may take $\Delta_{n}>0$.

To summarize, the exact dimensionless solution to the IBVP comprised of Eq. (10), Eq. (11) and Eq. (3) is

[TABLE]

where $a_{n}(T)$ is given in Eq. (15). To the best of our knowledge, Eq. (17) is a new exact solution [in any context in which Eq. (10) arises]. The corresponding IVP on $Z\in(-\infty,+\infty)$ is discussed in [55, Sec. V-A], including the derivation of the fundamental solution [i.e., the solution to the IVP with ICs $C_{\mathrm{d}}(Z,0)=\delta(Z)$ and $(\partial C_{\mathrm{d}}/\partial T)_{T=0}=0$]. The fundamental solution was found to consist of a singular part localized at the origin (but decaying exponentially in time) and a regular part that behaves asymptotically as the well-known Gaussian solution of the diffusion equation [i.e., Eq. (10) with $\mathcal{D}=0$]. Unlike the latter case, all terms in Eq. (17) are regular (i.e., not singular). With $C_{\mathrm{d}}(Z,T)$ now determined, $C_{\mathrm{nd}}(Z,T)$ can be easily found from the dimensionless version of Eq. (1b).

Note, however, that the term premultiplying $\sinh(\cdot)$ in Eq. (15) scales as $n$, which is expected to cause a significant Gibbs phenomenon in the final Fourier series solution in Eq. (17), when $\mathcal{D}\not\to 0$. As suggested in [35], to mitigate the Gibbs phenomenon when making plots, the Fourier series from Eq. (17) is evaluated by multiplying each term in the summation over $n$ by a cubic Lanczos $\sigma$-factor [39, pp. 221–227], i.e., $\operatorname{sinc}^{3}(n\pi/N)$, where $n=N-1$ is the final term in the truncated series. The $\sigma$-factors reduce the Gibbs phenomenon without affecting the convergence of the Fourier series.

As $T\to\infty$, the non-Fickian solution given by Eq. (17) is expected to converge to the equivalent solution of the Fickian diffusion equation [i.e., Eq. (10) with $\mathcal{D}=0$] [55], which can be computed following the same steps as above:

[TABLE]

where $D_{\mathrm{F}}:=\beta/(1+\beta)$ is dimensionless.

Figure 2 shows, respectively for $\mathcal{D}=5\times 10^{-3}$ (a,b,c) and $\mathcal{D}=5\times 10^{-2}$ (d,e,f), comparisons between the exact solutions given by Eq. (17) (solid) and Eq. (18) (dashed) at different $T$. Common values of $\beta=1/9$ (corresponding to $\delta_{0}/R=1/10$ for a half-full tumbler, $\phi=0.5$) and $\ell=1/5$ are used in all plots. The values for $\mathcal{D}$ were selected to best illustrate the features of the model. The most obvious feature of the non-Fickian model [Eq. (17)] is that, for finite $\mathcal{D}$, the initial pulse takes a nontrivial amount of time to relax, exhibiting persistent discontinuities at $Z=\pm\ell$. For $T=\mathcal{O}(1)$ (i.e., “long” times), the non-Fickian and Fickian solutions agree, but the time it takes for them to do so increases with $\mathcal{D}$. This result (also shown to hold for the IVP on $Z\in(-\infty,+\infty)$ in [55]) justifies using a Fickian diffusion model at long times. However, the non-Fickian aspects have non-trivial consequences at early times (here, in comparison to the axial diffusion time) in the transport process.

4 Discussion

Analogy to viscoelasticity.

Consider the following “constitutive law” between the axial flux $q_{\mathrm{d}}$ and concentration gradient $\partial c_{\mathrm{d}}/\partial z$ of the diffusing particle species:

[TABLE]

Then, Eq. (5) can be obtained by combining the latter with the axial conservation of mass (continuity) equation ${\partial c_{\mathrm{d}}}/{\partial t}+{\partial q_{\mathrm{d}}}/{\partial z}=0$. The significance here is that Eq. (19) is a common constitutive relation, termed the Jeffreys model [34] (see also [5, §5.2b], wherein a common sign convention for stress and flux is used), relating the stress to the rate of strain in polymeric fluids. Interestingly, in Eq. (19), the “kinematic viscosity” $D\beta/(1+\beta)$, the “relaxation time” $1/[k(1+\beta)]$, and the “retardation time” $1/(k\beta)$ are not independent (as they would, in principle, be for a polymeric fluid), since they are all functions of $\beta$.

From this analogy to viscoelasticity, an exact solution of the IBVP for Eq. (10) on $z\in[0,+\infty)$ with $c_{\mathrm{d}}(0,t)=H(t)$ can be immediately obtained by adapting any of the three representations known from non-Newtonian fluid mechanics (see, e.g., [15, 16]); similarly for the IBVP on $z\in[0,1]$ [17]. These IBVPs can be interpreted as various industrial feed problems. A time-dependent BC, $c_{\mathrm{d}}(0,t)=c_{\mathrm{d},0}(t)H(t)$, can be accommodated via the convolution theorem for the Laplace transform.

Connection to fractional (anomalous) diffusion.

A salient feature of the analogy to a Jeffreys-type viscoelastic model is that the exchange of particles between the bulk (non-diffusing) and surface (diffusing) species can clearly be interpreted as “memory” in the axial diffusion process, and quantified through $\beta$ and $k$ given explicitly in Eq. (8) and Eq. (9), respectively.

Specifically, in this analogy to linear viscoelasticity, Eq. (19) can be re-expressed as a flux–gradient relation in the form of a memory integral:

[TABLE]

with the effective kernel

[TABLE]

For the sake of argument, we have assumed that the flux $q_{\mathrm{d}}$ does not depend on the history of $\partial c_{\mathrm{d}}/\partial z$ prior to the initiation of flow ($t=0$) to set the lower limit of integration in Eq. (20). Observe that the kernel in Eq. (21) contains a Maxwell-type (exponential) “fading memory” in addition to the time-local contribution (Dirac-delta term) (see [5, §5.2b] for further details, including discussion of the factor of 2). For $\mathfrak{K}(t-t^{\prime})=2D\delta(t-t^{\prime})$ (no memory), Eq. (20) reduces to Fick’s first law [24, 6], $q_{\mathrm{d}}=-D{\partial c_{\mathrm{d}}}/{\partial z}$, and $D$ is the usual diffusivity. For a granular flow, the constitutive relation between $q_{\mathrm{d}}$ and $\partial c_{\mathrm{d}}/\partial z$ is unknown but Fick’s first law is often invoked [38, 31, 21]. A memory-integral constitutive equation for axial granular diffusion, as in Eq. (20), is reasonable because granular force networks [32] form and carry “information” across the flowing material. The latter can be interpreted as nonlocal action, leading to (possibly) highly-correlated particle distribution statistics (i.e., “memory”) in the flow [23, 43].

Now, what if $\mathfrak{K}$ decayed “slower than exponentially” [27]? Then, a generic form for $\mathfrak{K}$, which decays as $t-t^{\prime}\to\infty$, is a power law:

[TABLE]

where $\Gamma(\cdot)$ is the Gamma function, $\alpha\in(0,1)$ is a real number, and the proportionality constant $D_{\alpha}(>0)$ does not have the meaning of diffusivity. Substituting $\mathfrak{K}$ from Eq. (22) into Eq. (20), we arrive at

[TABLE]

instead of Eq. (19), where

[TABLE]

is the Riemann–Liouville fractional derivative of order $1-\alpha$ [51, 48], assuming suitable behavior of $\Xi$ as $t^{\prime}\to t$ to eliminate boundary terms upon differentiating under the integral sign. Substituting Eq. (23) into the continuity equation, we obtain a fractional diffusion equation

[TABLE]

Unlike Eq. (5), Eq. (25) has the attractive property of possessing self-similar solutions in terms of Fox’s $H$-function in the similarity variable x\big{/}t^{\alpha/2} [48]. Note that Eq. (25) is not exactly in the form of the anomalous axial granular diffusion equation posed in [37] due to technicalities regarding initial conditions for fractional PDEs. However, we interpret the intent in [37] to be the same.

Through the connection to the memory integral in Eq. (20), we see that the proposed model [Eq. (5)] and the fractional diffusion model [Eq. (25)] are two cases of a common general theory. A schematic comparison of the shapes of the kernels from Eq. (21) and Eq. (22) is shown in Fig. 3. Ultimately, only experimental measurements of the kernel $\mathfrak{K}$ would yield a definitive model for axial diffusion (accounting for various “anomalies” heretofore discussed). Such measurements remain a challenging open problem in granular mechanics.

Further refinements.

Finally, in this work, we focused on transport phenomena in a partially-filled cylindrical tumbler (fill fraction $\phi\lesssim 0.5$), which has received the most attention in the experimental studies discussed in Sec. 1. It is well known that, if $\phi>0.5$, then a “core” of quasi-static granular material forms along the cylinder axis [47, 28, 60]. Further, a non-circular tumbler cross-section [30, 18, 65] also significantly modifies the mixing and segregation dynamics. In principle, a large fill fraction or a non-circular cross-section would modify the rate constant $k$, and perhaps necessitate unequal exchange rates in each of Eqs. (1). While these details do not fundamentally change the proposed flowing layer–fixed bed exchange mechanism, it would be of interest to address these details in future work.

Acknowledgement

This paper is dedicated to memory of the late Prof. Denis L. Blackmore for his contributions to continuum models of granular flows (e.g., [7, 54]). I.C.C. would also like to acknowledge discussions (circa 2014) on wave phenomena [55] with the late Prof. A. M. Samsonov.

This work was initiated while I.C.C. was working with H.A.S. at Princeton University. At that time, I.C.C. was supported by the U.S. National Science Foundation (NSF) under grant no. DMS-1104047. I.C.C. would like to further acknowledge the hospitality of the University of Nicosia, Cyprus, where this version was completed thanks to a Fulbright U.S. Scholar award from the U.S. Department of State.

Research data availability

No data was generated as part of this study.

\printcredits