Minerals Engineering 119 (2018) 244–262

Contents lists available at ScienceDirect

Minerals Engineering
journal homepage: www.elsevier.com/locate/mineng

A hierarchical simulation methodology for rotary kilns including granular T

flow and heat transfer
⁎
Peter J. Witta, Matthew D. Sinnottb, Paul W. Clearyb, M. Philip Schwarza,
a
CSIRO Mineral Resources, Private Bag 10, Clayton South, VIC 3169, Australia
b
CSIRO Data61, Private Bag 10, Clayton South, VIC 3169, Australia

A R T I C L E I N F O A B S T R A C T

Keywords: Rotary kilns are used in several minerals processing operations, as well as related industries such as cement
Rotary kiln manufacture. Since they involve complex multiphase heat and mass transfer processes, optimization would
Simulation ideally be guided by numerical modelling. A new generic hierarchical approach is proposed in this paper, se-
Computational fluid dynamics (CFD) quentially combining 2D Discrete Element Method (DEM) simulations of slices of the bed and 3D two-phase
Discrete element method (DEM)
computational fluid dynamics (CFD) models of the entire kiln, encompassing both bed and gaseous freeboard.
Minerals processing
DEM simulations have been used to derive mean solids velocities in a bed for various rotation speeds and particle
Granular flow
Coulombic friction sizes, and these were then used to calibrate and validate the solids rheology model used in the subsequent CFD
Drucker-Prager model model. A modified Coulombic friction for the bed in the CFD model was shown to give satisfactory agreement
Hierarchical multiscale modelling with the DEM bed results over a range of rotation speeds. Similarly, thermal mixing simulations carried out using
the 2D DEM model for the granular bed were used to validate a CFD model of heat transfer after calibration of
the small-scale diffusive term. The capability of the resultant CFD model was demonstrated by 3D simulation of
the bed of a large rotary kiln.

1. Introduction or more flames will typically be located. In the Becher kiln, air injection
into the kiln is staged along the length of the reactor so as to control
Many minerals processing operations utilize rotary kilns for pro- temperature and oxidation degree. Such complexity of heat and mass
cesses such as roasting, calcining and nodulization. Examples include transfer processes cannot be accounted for in simple models such as
alumina calcining, the ilmenite reduction stage of the Becher process, analytical expressions, empirical correlations, or even one-dimensional
iron ore reduction, and pet coke calcination. models. Detailed measurements are extremely difficult to make within a
Industrial applications of rotary kilns can be classified into four kiln, and even measurements of the solids and gas as they leave the kiln
main types: are very time consuming and expensive. Temperature measurements at
the shell are similarly expensive, and all plant measurement attempts
• Drying and pre-heating (including pellet induration) are problematic because of the inherent variability of kiln operations,
• Calcination (eg. alumina, limestone, pet coke, TiO , etc) 3 and the long time-scales involved. Consequently, a numerical modelling
• Mineral oxide reduction and other similar reactions approach is needed if optimization is to be carried out in a predictive
• Cement manufacture way.
One dimensional models have been developed and applied for many
Though there has been a trend in some applications towards the use years with some success in predicting global behaviour and quantities
of fluidized beds and suspension systems, kilns are still widely used. such as overall conversion (see, for example, Nicholson, 1995). How-
Though rotary kilns are robust and the processes generally stable, ever detailed information, such as would be required to avoid ring
issues such as ring accretion formation and minimization of energy formation, requires a multi-dimensional modelling approach. Compu-
consumption have been difficult to address decisively. Temperature tational fluid dynamics (CFD) models are naturally suited to dealing
distribution can be difficult to predict because reaction kinetics and with the gaseous freeboard including combustion processes, and can
heat transfer depend both on mixing of the solids in the bed, and on also be readily extended to incorporate heat transfer through the shell.
transfer between the bed and the gaseous freeboard space in which one Discrete Element Method (DEM) models are most naturally applied to

⁎
Corresponding author.
E-mail address: [email protected] (M.P. Schwarz).

https://doi.org/10.1016/j.mineng.2018.01.035
Received 29 July 2017; Received in revised form 22 January 2018; Accepted 25 January 2018
Available online 22 February 2018
0892-6875/ © 2018 Elsevier Ltd. All rights reserved.
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

modelling the granular bed, but the great length of many kilns, com- advantage is short run times and that they can be tuned for a specific
bined with the small particle size (ilmenite feed is typically less than plant, thus allowing numerous changes in operating conditions to be
1 mm) and long timescales mean that practical design calculations quickly assessed. They can also be used in real time by operators by
using DEM to track every particle in an industry scale kiln is still being incorporated into process control systems. A disadvantage of 1D
computationally intractable for design purposes. models is that they are based heavily on plant operating data, making
A new hybrid approach is described in this paper, hierarchically them only accurate over a limited range of operating parameters or for
combining 2D DEM simulations and 3D two-phase CFD models of the a specific application. Furthermore, they lack the ability to account for
entire kiln, encompassing coupled dynamics and thermal behaviour of 3D details, and these are often important to kiln operation.
both bed and gaseous freeboard. The approach has been demonstrated Thornton and Batterham (1982) describe a 1D phenomenological
in this work: DEM simulations were used to derive mean solids velo- model of a rotary kiln used for processing iron-ore pellets. This model
cities in a bed for various rotation speeds and particle sizes, which were accounts for heat transfer by radiation and convection in the gas phase
then used to confirm the validity of the solids rheology model used in and solid conduction through the kiln wall. Motion of the solids is
the CFD model, and calibrate the model parameters. This approach considered with heat transfer within the bed and from the bed to the
allows account to be taken of processes occurring at the particle scale as gas and walls. This model was further developed by Thurlby (1988) and
well as their contribution to processes occurring at the overall kiln scale associated researchers. Boateng and Barr (1996) extended this concept
over very long timescales. This is an example of a multi-scale approach, of a 1D model by coupling it to a 2D transverse model that accounts for
as described by Schwarz and Liu (2005), Liu and Schwarz (2009) and (transverse) bed mixing and segregation. Kritzinger and Kingsley
Schwarz et al. (2015). It should be emphasized that this hierarchical (2015) showed how this type of 1D model can also be extended to
approach is quite different from the common coupled CFD-DEM tech- account for multiple air jets (or flames) along the length of a rotary kiln,
nique in which the CFD component is used simulate the interstitial gas as is often used in pyrometallurgical minerals processing, e.g. in iron
and DEM to simulate the granular solids. In the present hierarchical ore reduction. While such 1D models are clearly useful for design in-
approach, the CFD and DEM modelling are carried out sequentially, vestigations, CFD simulations will undoubtedly be much more accurate
rather than in the close-coupled concurrent way of the method com- for such complex multi-dimensional combustion processes where tem-
monly known as “CFD-DEM”. Furthermore, both CFD and DEM stages perature and species concentrations are not uniform on transverse
simulate the granular phase (at different scales), though the CFD cross-sections, and more importantly, have transverse profiles that are
component can also be used to simulate the interstitial gas when used in dependent on their longitudinal position along the kiln. Similarly, DEM
a two-fluid mode. The present approach is intended to be applicable to models will be able to better capture more complex granular bed
a wide range of different kiln processes, and takes advantage of the fact phenomena such as size and density segregation. In some rotary kiln
that CFD models are readily customized to incorporate specific flame processes, gases are generated by reaction of solids in the bed, causing
types, specific gas-solids reactions, kiln internals, and so on. fluidisation of the particulate bed in some zones. Such phenomena
cannot be captured in any way by a 1D model.
2. Previous work on modelling rotary kilns
2.2. DEM models
Only a limited amount of information appears in the literature de-
scribing detailed measurements of granular bed dynamics in plant or DEM models typically focus on the dynamics of particle motion,
pilot plant operations, Nicholson’s (1995) work being an exception. On simulating individual particle collisions, and are therefore valuable for
the other hand, there is a rich literature on solids motion in idealized elucidating the small-scale physics affecting bed dynamics. As a result,
cylinders, both horizontal and inclined (see, for example, Metcalfe DEM simulations are expensive in computer time, particularly for in-
et al., 1999). The axial motion closely follows plug flow, with sig- dustrial scale kilns with fine particles. Explicit integration of particle
nificant mixing occurring in the transverse plane. Axial transport is well motion can also be a problem for systems with very long time-scales:
characterised both experimentally and theoretically, with formulas for thermal time-scales are extremely long for industrial rotary kilns. The
residence time available, even for the situation with dam rings in- Discrete Element Method is described by Cleary et al. (1998), Cleary
stalled. While the motion in the transverse plane is well understood and Sawley (2002) and Cleary (2004, 2009) in relation to various
qualitatively, no simple methods to predict mixing and segregation particle processing operations. It has particularly been used to explore
exist. size, density and shape segregation in slowly rotating 3D drums which
Behaviour of a collection of particles in the transverse plane of a have very similar bed dynamics to rotary kilns (Pereira et al. (2011,
rotating cylinder can be categorized into one of six regimes: slipping, 2014) and Pereira and Cleary (2013, 2017)).
slumping/avalanching, rolling, cascading, cataracting or centrifuging. Finnie et al. (2005) studied longitudinal and transverse mixing in
Of the six regimes only the first four are usually found in industrial horizontal isothermal rotary kilns, using two and three-dimensional
rotary kilns. The mixing rate and potential for segregation is strongly DEM. The focus was on the effect of the main operating conditions on
influenced by the operating regime of the kiln bed. Parameters such as quantitative measures of longitudinal and transverse mixing in short
rotational speed, kiln diameter, percentage fill, friction between the horizontal (i.e. non-inclined) rotating drums with a ratio of kiln radius
wall and bed, and particle diameter and density determine the oper- to average particle radius of 40.0. While this is a valuable study for
ating regime of the bed. DEM models have had considerable success in understanding trends, the particle size is very much larger than that
predicting these complex phenomena. used in industrial full-scale kilns, and the kiln is short: while it is long
Mathematical models of kilns can conveniently be divided into two enough to capture some 3D effects, the aspect ratio is much smaller
categories: one-dimensional process models which divide the kiln than that used in industrial kilns. The size of the simulation is thus very
lengthwise into several slices or zones and carry out heat and mass far from what would be needed for a realistic simulation of a full-scale
balances on each zone; and multi-dimensional CFD and DEM models industrial kiln. Kwapinska et al. (2006) analysed the details of trans-
that attempt to solve for the gas and solids motion and chemical com- verse mixing between two populations of particles in a rotary drum
positions at each point in the kiln. using two-dimensional DEM simulations, comparing the results with
experimental data, and pointing out that accurate prediction of inter-
2.1. One-dimensional models particle mixing is needed for the models to be able to predict thermal
behaviour in kilns. Chaudhuri et al. (2010) has demonstrated how such
In general, the 1D models can provide axial profiles of temperature, particle dynamics simulations can be extended to simulating heat
residence time, solids fill and gas composition along the kiln. Their transfer in rotary calciners.

245
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

2.3. CFD models questionable, as discussed further below.

Amongst other published techniques for rotary kiln simulation is the
While DEM models are ideal for simulating the details of particle work by Küssel et al. (2009) on kilns for limestone calcination. A CFD
dynamics in a granular bed, CFD models have been used extensively for model of the gas-space flame was used to provide heat transfer in-
heat and mass transfer in complex fluid reaction systems, with exten- formation to a classical one-dimensional process model of the kiln bed,
sions to multiple phases. The first published rotary kiln application of a in a manner analogous to the technique used by Mastorakos et al.
three-dimensional CFD model was that developed by Bui et al. (1995) (1999). While this is an advance on a simulation solving 1D equations
for a coke calcining rotary kiln. Four separate sub-models were used, for both gaseous freeboard and bed, it suffers from all the drawbacks
one each for the kiln bed, gas space (with reactions), kiln walls and gas mentioned above for 1D bed models.
space radiation. Coupling protocols link the four models together with a
number of iterations over the four models required to converge the 3. Model formulation
entire solution. By treating the bed as a viscous liquid, the CFD plat-
form, PHOENICS, could be used to predict the granular flow behaviour 3.1. A hierarchical hybrid model
of the bed. A body-fitted grid was applied, with the bed shape specified
by assuming a typical roll regime bed profile as seen in laboratory ex- The hybrid CFD-DEM model described here aims to exploit the
periments, rather than predicting the bed profile. Different viscosities strengths of DEM modelling in terms of simulating the detailed inter-
were used for the active layer, stagnant layer and in the axial direction. actional dynamics of particles and of CFD modelling in terms of well-
Values for the viscosities were selected to match the bed behaviour seen established techniques for simulating reacting turbulent flow and heat
in the laboratory experiments. The authors explain that real kiln ob- transfer at the large scale. The basic principle is to use DEM simulations
servations and small-scale model studies yield values of active-layer to determine appropriate constitutive relations to allow incorporation
surface velocities of 1.42 ± 0.4 m/s, and to obtain these velocity va- of a continuum model for the granular flow and heat transfer within a
lues in the model, the kinematic viscosities required were found to be two-phase (gas-solids) 3D CFD formulation.
1.0 × 10−5 m2/s for the active layer and 1.5 m2/s for the underlayer, The reader is referred to work by Cleary et al. (1998), Cleary and
which constitutes a difference of 5 orders of magnitude. They assume Sawley (2002) and Cleary (2004, 2009) which give details of the DEM
zero slip at the kiln wall. Conduction in the kiln refractory walls was method and the CSIRO DEM software that is used in these simulations
modelled by solution of the energy equation with boundary conditions described here. Similarly, the techniques used for solving the CFD
obtained from the other sub-models. Heat transfer within the bed is equations (conservation equations for mass, momentum and, if needed,
dominated by advection, which was accounted for by the three-di- energy) are not detailed here, but can be found in the CFX User manuals
mensional fluid flow model of the bed. In addition, the thermal con- (AEA Technology).
ductivity of the bed is taken to be 0.1 W/m K, which, the authors claim, For gas or liquid single-phase flow, the conservation equations are
takes into account the porosity of the coke. This yields a thermal dif- the Navier-Stokes equations. This set of continuum equations, derived
fusivity of the bed of the order of 10−7 m2/s, which is clearly negligible from the conservation of mass, momentum and energy, are widely ac-
compared with advection. Assuming valid boundary conditions, the cepted as the governing equations for the flow of a homogeneous
accuracy of the simulated bed heat transfer will depend entirely on the Newtonian fluid. In two-phase gas-particle flows, the well-accepted
accuracy of the predicted bed velocities. The overall simulation scheme two-fluid technique solves conservation equations similar to the single-
was used to investigate the operation of a 60 m long and 2.4 m diameter phase Navier-Stokes equations, obtained by space and time averaging of
coke calcining kiln. The calculation scheme was inventive and prag- mass, momentum and energy for each phase. Additional constitutive
matic given the computational capabilities of the early 1990 s, but the equations for the transfer of momentum and energy between each
approximation that the bed profile shape can be assumed the same as phase are required. For the gas phase, the stress terms are the laminar
observed in laboratory experiments clearly has its limitations. In rea- (viscous) stresses, or in turbulent flows the so-called Reynolds stresses.
lity, the bed surface shape is expected to depend on rotation rate, For the particle phase, the stress terms are more complex, with accepted
particle size and particle shape, so that a single bed surface profile, as models only just starting to emerge. One means of determining the
assumed in the model approach, is unlikely to be adequate. Assumption particle phase stress terms is from the kinetic theory of granular flows
of effective solids viscosity for the bed is even more problematic, given (Lun et al., 1984). In this approach, stresses and viscosity for the solid
that bed movement will also depend on these operating parameters. phase are derived from a granular temperature in a manner analogous
Clinker formation in a coal fired rotary kiln was modelled by to that used in gases to obtain the relationship between viscosity and
Mastorakos et al. (1999) using a two-dimensional axi-symmetric geo- gas temperature from a consideration of molecular collisions. For par-
metry CFD model of the free-board, in which coal combustion, radiative ticle flows, the granular temperature is a measure of fluctuations in the
transfer, and heat conduction were accounted for. This gas phase model velocity of particles about their local mean velocity. Further details on
was coupled to a one-dimensional model of the bed in which the clinker the derivation and implementation of the terms for gas-particle flows
charge moved axially at constant speed and was assumed to be well are given by Witt et al. (1998). These techniques have been successfully
mixed (in temperature and composition) in each cross-section by action used for simulation of fluidized bed processes (e.g. Witt et al., 1999,
of the kiln rotation. Bed reactions are included subject to the 1D as- Rahaman et al., 2003, Feng et al., 2012, Hassan et al., 2016).
sumption that concentrations are uniform at each axial cross-section. Solids material in rotary kilns exhibits complex flow behaviour, and
Heat transfer is also simulated through the kiln walls, including an physical models needed for the CFD equations have not yet been de-
assumed uniform-thickness layer of solid clinker. veloped and validated for such flows. In principle, the same two-fluid
Eghlimi et al. (1999) performed CFD modelling of the three-di- formulation used for fluidized beds can be applied to the granular beds
mensional behaviour and mixing of raw oil shale and recycled com- of kilns, but what is missing from the current models are the con-
busted oil shale particles in a rotary kiln for a retorting process for the stitutive relations governing the flow behaviour of the solids phase.
extraction of oil. An Eulerian-Eulerian three-phase model was used with Specifically, the major challenge is to include the complex rheological
one phase being gas and separate solid phases used for the two different behaviour of bed solids in the model. Experimental investigations of
oil shale particles represented as pseudo-fluids. A granular flow model solids flow in rotary kilns clearly show two distinct solid flow regimes
(Syamlal and O’Brien, 1989) based on kinetic theory was used to de- (Boateng and Barr, 1997): fluidized flow in the active region and
termine particle viscosities. Only limited details of the model and re- sheared granular flow in the bed below. More fundamental studies of
sults were presented, with no validation described. Applicability of the shearing particle flow using DEM simulations (Zhang and Campbell,
kinetic theory of granular flow to a bed of closely packed particles is 1992) have also shown this transition between rapid granular flow and

246
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

packed dense beds. τ = μ fr σ + c (7)

Kinetic theory models are thus potentially suitable for predicting the
where μ fr is the coefficient of friction, c the cohesion and σ the normal
behaviour in the active surface layer. However, as discussed by Zhang
stress. For convenience, an internal angle of friction, ϕ, is defined as:
and Campbell (1992), kinetic theory models are unlikely to be suitable
for the more stagnant zone below, or the transition between fluid-like tan φ = μ fr (8)
and solid-like behaviour. For these regions, a stress model using a yield The critical angle of friction, ϕr, which applies at the boundary
criterion and shear stress model based on the normal stress would need between flowing and immobile (static) regions is somewhat different
to be developed, incorporated into the CFD solution scheme, and then from the dynamic angle of friction, ϕD, which applies to a closely
validated. For fluidized beds this has been undertaken by Vun et al. packed flowing granular bed. In the shear cell apparatus used by Hanes
(2010) where there are typically only small areas of defluidization and and Inman (1985), tan ϕD increases with height in the apparatus, which
thus flows involving yielding granular flow. Application of this ap- corresponds to increase in shear rate. Hanes and Inman explain this
proach to the granular flow in kilns would require further extension and change as resulting “either from the angular distribution of collisions in
validation. the flow or a relative change in the importance of multiple collisions
and sliding friction between the grains”.
3.2. Fluid-like behaviour For soil and mechanical and granular flows in silos the viscous-
plastic Drucker-Prager model has been used (Karlsson, 1995, Karlsson
When a force is applied to a fluid, it deforms or flows with re- et al., 1998). In the Drucker-Prager model the yield criterion is given
arrangement of the fluid molecules. While the fluid is in motion a ve- by:
locity gradient exists, which gives rise to a shear stress. To maintain a
0 = q + mp (9)
shear stress the fluid requires a continued motion or a velocity gradient.
For fluids the simplest constitutive stress model is that of a Newtonian with m being a material property (a friction coefficient) defined as:
fluid. For the Newtonian case the stress terms are given by:
18sinφ
m=
∂P 4 ∂u 4 ∂v ∂w ⎞ 9−sin2 φ (10)
σx = − + ⎛κ + μ⎞ + ⎛κ− μ⎞ ⎛ + ⎜ ⎟
∂x ⎝ 3 ⎠ ∂x ⎝ 3 ⎠ ⎝ ∂y ∂z ⎠ (1)
The variables p and q are the stress invariants:
∂P 4 ∂v 4 ∂u ∂w ⎞ 1
σy = − + ⎛κ + μ⎞ + ⎛κ− μ⎞ ⎛ + p= σii
∂y ⎝ 3 ⎠ ∂y ⎝ 3 ⎠ ⎝ ∂x ∂z ⎠ (2) 2 (11)

∂P 4 ∂w 4 ∂u ∂v ⎞ 2
σz = − + ⎛κ + μ⎞ + ⎛κ− μ⎞ ⎛ + q= sij sij
∂z 3 ⎠ ∂z
⎜ ⎟
3 (12)
⎝ ⎝ 3 ⎠ ⎝ ∂x ∂y ⎠ (3)
where sij is the deviatoric stress tensor:
∂u ∂v ⎞
τxy = τyx = μ ⎛ +
⎜ ⎟
sij = σij−pδij (13)
⎝ ∂y ∂x ⎠ (4)
At walls, fluids are assumed to be attached to the wall and some
where σ is the normal stress, P the pressure, τ the shear stress, κ the bulk
form of velocity profile between the wall and moving fluid is assumed
viscosity and μ the shear viscosity.
to obtain the shear stress. Granular materials can either move at the
wall velocity or slip along the wall surface with the shear stress given by
3.3. Solid-like behaviour
a function of the normal stress:

With elastic solid materials, when a force is applied the solid will τwall = μ wall σn (14)
deform elastically creating a shear stress in the material. If the load is where μ wall is the wall friction coefficient and σn the normal stress at the
removed the material will return to its original shape. As long as the wall. This gives the upper limit for stress at the wall in the absence of
force is applied the deformation will remain (assuming that the max- slip and once exceeded, slip between the material and the wall occurs.
imum stress the material can withstand is not exceeded). For a solid the
shear stress is given by:
4. Determination of constitutive parameters from DEM
τ = Gγ (5)
The purpose of the DEM component of the kiln model is to collect
where G is the shear modulus and γ the strain. Note that the strain rate
the numerical data required for the rheological model of the bed to be
is given by:
used in the full scale CFD model of the kiln, namely average velocity
∂u ∂v ⎞ fields, collisional stresses, wall slip velocity and void fractions. These
γ̇ = ⎛⎜ + ⎟
continuum quantities are calculated using suitable spatial and temporal
⎝ ∂y ∂x ⎠ (6)
averaging from the particle information predicted using the “pure”
DEM method via the CSIRO DEM software mentioned in Section 3.1. A
3.4. Constitutive model for granular material used in the CFD model similar approach in a more general context has been followed by Luding
(2009a, 2009b) and others.
Granular materials exhibit both fluid-like and solid-like behaviour To this end, DEM simulations of avalanching granular material were
depending on the applied load. Under light or moderate loads a gran- performed using the CSIRO DEM software for a 2-metre diameter cy-
ular material will deform elastically as the particles are locked together. linder rotating at low speeds, and a parametric study of variation in
Once the applied load exceeds a critical value the particles are able to stress distributions and solids fraction for various parameters was car-
slide relative to one another, this results in a slip or failure plane oc- ried out to provide “data” for the CFD model development and cali-
curring: many such slippages may manifest as continuous plastic de- bration and to assess their influence on the bed behaviour. Effective
formation. The thickness of the slip plane is commonly of the order of shear viscosities and bed-shell friction coefficients were also derived
ten particle diameters (Nedderman, 1992). The shear stress behaviour from the simulations for various operating parameter combinations.
of many granular materials can be described by the Coulomb yield Most simulations were carried out for circular disc particles (in 2D), but
criterion: the effect of particle “blockiness” on constitutive properties has also

247
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

been investigated for base case operating parameters. For all the DEM values vary more continuously across the depth of the bed.
simulations, the values of the normal and tangential spring stiffness The bed-shell friction coefficient is required for the continuum
were taken to be 100,000 N/m and 50,000 N/m respectively, and the model in order to predict bed-shell slip velocities and needs to be cal-
friction coefficient for both particle-particle and particle-wall contacts culated from the DEM simulations. This has been determined from the
was taken to be 0.5. These values are based on previous work and are spatial-average along the shell boundary of the ratios of time-averaged
not intended to represent any particular material, since the present shear and normal collisional stresses. This friction coefficient is in-
work is a demonstration of the approach rather than a simulation of a tended for use in the full-scale CFD model of the kiln.
particular kiln application.
4.2. Influence of various parameters on the bed behaviour
4.1. Base DEM simulation
Several other ‘pure’ DEM simulations were performed, varying dif-
The base DEM simulation was performed with a 2 m diameter cy- ferent parameters to assess their influence on the bed behaviour and the
linder filled to 20% of the volume and rotating at 3 rpm (which is 10% average rheological and flow quantities. All the simulations were run
of the critical speed at which particles centrifuge or a Froude number of for 60 s (corresponding to 3 full revolutions at 3 rpm in the base case,
0.1). The particles had a density of 2700 kg/m3 and were chosen to be and 1.5 and 6 revolutions for the other rotation rates), and the figures
20 mm diameter discs (i.e. 100 particles across the diameter of the show the cumulative results at the end of each simulation. In each case
cylinder) with a narrow uniform PSD (particle size distribution) range it was verified that the simulations were performed over a sufficient
of ± 5%. There are around 1700 particles in such a configuration. The duration to ensure convergence of the average bed quantities. All si-
simulation was carried out for 60 s (3 full revolutions of the cylinder) mulations have been filtered by volume fraction when deriving the
and the various bed quantities were averaged cumulatively every 10 s shear viscosity to eliminate noise from spatial discontinuities in the
to ensure that their values had converged by the end of the simulation. strain rate distribution at the bed free surface.
The fine grained cylindrical grid used for this spatiotemporal averaging
used 30 cells across the kiln diameter with a cell size of radial length 4.2.1. Influence of the fill level
0.07 m. The averaged quantities consist of the velocity distributions, the Alternative fill levels of 10% and 30% were chosen (to compare
collisional and streaming stresses, the granular temperature and the with a level of 20% in the base case) with particles being 20 mm dia-
volume fractions (see Cleary (2008) for definitions of these quantities meter discs ( ± 5%). Fig. 2 shows the bed viscosity variation for the
and their calculation from DEM data). three fill levels. At 10% fill, this region covers much of the bed-shell
Fig. 1 shows the shear viscosity distribution for the base simulation interface. With increasing fill level, the high viscosity region moves to
case derived from the computed shear stress and strain rate (and shown the right and becomes more concentrated below the bottom of the
on a cylindrical grid). We see that the shape of the viscosity pattern is avalanching surface. The magnitude of the shear viscosity varies line-
dominated by the shear stress distribution with the strain rate being arly with bed depth (and hence bed weight) over this range.
relatively uniform within the bed. The regions of high viscosity appear The wall friction coefficient was found to be fairly insensitive to the
to lie at the bottom of the bed close to the kiln wall. The parts of the bed depth. Table 1 indicates, at most, a subtle increase in bed-shell
space where there is no data (because they are above the surface of the friction as a function of fill level (bed depth). The packing of particles in
bed) are shown as black. The granular free surface passes through the the dense regions of the bed (near the shell) is also fairly insensitive to
upper coloured cells and cannot be identified precisely. Any dis- fill level, and a similar number of particle/shell contacts are maintained
continuities in the strain rate appear close to the bed surface where the irrespective of the volume of material in the bed.
shear stress is zero and therefore plays little part in influencing the
viscosities. The presence of a transient free surface leads to some noise 4.2.2. Influence of the rotation rate
in the estimates of the average strain rate at these surfaces, so all The rotation rate in the base case was set as 3 rpm. Simulations with
viscosities are filtered according to volume fraction to ignore data with rates half as fast and twice as fast were carried out for the same fill level
insufficient material contributing to the calculation. of 20%. Fig. 3 shows the shear viscosities for rates of 1.5 rpm (5% of the
It is interesting to compare the viscosities with those assumed in the critical speed), 3 rpm (10% critical) and 6 rpm (20% critical).
Bui et al. (1995) model described above. To obtain velocities similar to As the angular speed of the kiln increases, the regions of high
those seen in laboratory experiments, Bui et al. assumed a two zone viscosity become more pronounced. They spread more evenly over
viscosity distribution, with the active layer value of order 0.1 N s/m2, progressively larger areas and there is a matching decline in the max-
and the value for the underlayer taken to be order 1500 N s/m2. This is imum viscosity. The shear stress is largely invariant to rotation rate and
broadly consistent with the values plotted in Fig. 1, though the DEM the system remains in the same regime (discrete avalanching) for all

Fig. 1. Base DEM simulation - spatial distributions of: shear stress (left), strain rate (centre) and shear viscosity (right).

248
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

Fig. 2. Dependence of shear viscosity determined from the DEM simulations on fill level: 10% (left), 20% (middle) and 30% (right).

Table 1 fairly insensitive to rotation rate. There is a small variation but no clear
Dependence of bed-shell friction coefficient on parameters investigated in the DEM si- trend is evident.
mulations.

Parameter description Parameter Value Bed-shell friction 4.2.3. Influence of the particle size distribution (PSD) range
coefficient The PSD range was increased to ± 10%, ± 20% and ± 100% (to
Fill level (% of cylinder volume) 10 0.64 compare against ± 5% in the base case). Fig. 4 shows the variation of
20 0.66 viscosity pattern with PSD range. As the PSD becomes wider, the high
30 0.72 viscosity region becomes more compact and concentrated at the base of
Rotation rate (% of critical 5 0.70 the bed with commensurate increases in its maximum value.
velocity) 10 0.66 Table 1 indicates that the bed-shell friction is sensitive to PSD range.
30 0.67 For the largest value of PSD range, small particles are able to fill the
PSD range (in% around the mean 5 0.66 gaps between larger ones and the number of particle contacts along the
value) 10 0.67 shell is strongly increased. This strengthens the wall-bed contact and
20 0.68
leads to a higher bed-shell friction coefficient.
100 0.84

# of particles (across the cylinder 25 0.58

diam.) 50 0.63 4.2.4. Influence of average particle size
100 0.66 The average particle diameter was also varied. The values chosen
200 0.70 were 80 mm, 40 mm, 10 mm, 6.6 mm, 5 mm and 2.5 mm (i.e. 25, 50,
300 0.90
200, 300, 400 and 600 particles, respectively across the diameter of the
400 1.10
600 1.42 cylinder) in contrast to the base value of 20 mm (corresponding to 100
particles across). The numbers of particles in these simulations were
Blockiness index 2 0.66
3 0.67 respectively about 100, 400, 7000, 15,000, 27,000 and 60,000.
4 0.72 Fig. 5 shows that for particle diameters larger than 40 mm, the
6 0.91 granular material becomes too coarse to adequately assign a continuous
viscosity distribution. This is primarily due to the high resolution of the
data collection grid in our DEM simulations. One could use a coarser
these three different speeds. However, the rate of strain in the bed does
data grid, but this is not useful for calculating the continuous quantities
increase naturally with rotation rate. This results in a decrease in the
needed for the CFD model since this has minimum resolution require-
magnitude of the peak viscosity as the strain becomes progressively
ments for the viscosity distribution. With decreasing particle size, the
more evenly distributed. Table 1 shows that the bed-shell friction is
flow of material becomes more fluid-like and there is a drop in

Fig. 3. Dependence of shear viscosity determined from the DEM simulations on the rotation rate: 1.5 rpm (left), 3 rpm (middle) and 6 rpm (right).

249
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

Fig. 4. Dependence of shear viscosity determined from the DEM simulations on PSD range: ± 5% (top left), ± 10% (top right), ± 20% (bottom left) and ± 100% (bottom right).

magnitude of the peak viscosity. This can be seen in Fig. 6, which plots restricts their ability to rotate against the cylinder wall, resulting in an
the maximum viscosity determined from the DEM simulations against increase in the bed-shell friction.
the number of particles across the diameter of the cylinder (inversely
proportional to the average particle size). Fig. 5 shows that the high
viscosity region moves to the right and becomes concentrated in the 4.3. Summary of DEM study on viscosity and bed-shell friction
lowest part of the bed as the particle size decreases.
From Table 1, we can see a clear trend between friction coefficient The DEM simulations have identified important qualitative and
and particle size. For smaller particles, the number of shell/particle quantitative features of the effective bed viscosity and bed-shell friction
contacts is increased, resulting in a larger bed-shell friction as both the that need to be incorporated into a continuum model. The shear visc-
bed and bed-wall connections becomes stronger and are more able to osity distribution is qualitatively similar for all cases studied in that its
resist shear. maximum values lie at the bottom of the bed close to the cylinder wall.
The shear viscosity distribution closely follows the pattern of shear
4.2.5. Influence of the particle shape (blockiness/angularity) stress with modification by the strain rate distribution. It varies subtly
The effect of varying the particle shape was also considered for with some of the parameters varied.
varying blockiness index (with the particle shape described mathema- The maximum viscosities were found to be highly sensitive to all the
tically as a super-quadric) of the particles from 2 (circular) to 3 (fairly parameters studied. The highest viscosity case would occur for a bed of
round), 4 and 6 (relatively square). Fig. 7 shows that the blockiness of non-circular particles with a large fill level and low rotation rate. Bed-
the particles also influences the shear viscosity. ‘Blocky’ particles pack shell friction is strongly dependent on PSD range, particle size and
more densely in the bed than do the circular particles and lock together shape, but is largely insensitive to macroscopic quantities of the system
better (Cleary, 2008). Thus, the volume fraction increases with the such as fill level and rotation rate. This suggests the details of the bed
blockiness index, as does the ability of the particle bed to resist shear. microstructure are important for understanding the friction force be-
This results in larger peak values of shear stress leading to an increase in tween the bed and kiln wall. Shear viscosity and bed-shell friction are
the maximum viscosity. Note that the free surface becomes steeper with both highly sensitive to particle shape and further work in under-
increasing blockiness, reflecting the increasing angle of repose of the standing the effect of non-circular shape is desirable.
material. This correlates with the increase in the shear strength of the
material (Cleary, 2009). Increasing blockiness of the particles also

250
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

Fig. 5. Dependence of shear viscosity determined from the DEM simulations on the number of particles across the diameter of the cylinder (inversely proportional to the average particle
size): (a) 25, (b) 50, (c) 100 (base case), (d) 200, and (e) 300.

6000 with the black line showing ideal Coulomb behaviour given by Eq. (7).
Thus, the solid material is behaving in a manner that is generally
5000
consistent with a Coulomb material, although with a reasonable
Viscosity (N.s/m2)

4000 amount of spread about the solid line. However, implementation of this
model does not always predict bed behaviour precisely in accord with
3000 that found in the DEM simulations. More detailed investigation of the
rheology predicted by the DEM simulations, via a plot of the viscosity
2000
obtained from the DEM simulations against pressure (Fig. 10), shows
1000 that there is an approximately linear relationship between viscosity and
pressure for each case, with the slope being dependent on the specific
0 operating parameters for each case.
0 50 100 150 200 250 300 350 Plotting the coefficient of friction, m, (defined in Eq. (9)) against
Number of particles strain rate (Fig. 11) shows that m is not constant but is dependent on the
shear rate. By fitting a line through the grouped data points, a re-
Fig. 6. Dependence of maximum shear viscosity determined from the DEM simulations on
the number of particles across the diameter of the cylinder (inversely proportional to the
lationship between m and the strain rate has been obtained and is used
average particle size). as a modification to the ideal model of Eq. (7). The dependence of m on
strain rate agrees with the finding of Hanes and Inman (1985) men-
tioned above, that the dynamic internal angle of friction increases with
5. Implementation of CFD solids rheology model
strain rate in their shear cell apparatus. It is worth noting that an al-
ternative approach could be to fit the data to a non-linear form of Eq.
5.1. Analysis of DEM stress and strain data
(9), similar to that used in the Modified Cam Clay model.
At this point, a comment should be made on the procedure that
A plot of the shear stress against strain rate data for all points from
would be followed to determine the specific constitutive parameters for
the DEM model is shown below in Fig. 8. The solid black line indicates
a particular granular solid for an actual industrial situation. The DEM
the behaviour for a Newtonian fluid given in Eq. (4). This plot indicates
simulations described here have been used to determine the general
that the solid granular material is highly non-Newtonian, as one may
rheological behaviour in the kiln granular flow regime (and will be used
have expected.
below to similarly characterize heat transfer at a particle scale).
Fig. 9 shows a plot of shear stress against pressure or normal stress,
Parameters in the model will depend on details of particle size

251
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

Fig. 7. Dependence of shear viscosity determined from the DEM simulations on the particle shape (blockiness index): 2 (top left), 3 (top right), 4 (bottom left) and 6 (bottom right).

10% 3 RPM 3500

3500
20% 3 RPM 3000
3000
30% 3 RPM
Shear Stress [Pa]
Shear Stress [Pa]

2500
2500 20% 1.5 RPM
20% 6 RPM 2000 10% 3 RPM
2000
μ = 500 20% 3 RPM
1500
1500
30% 3 RPM
1000
1000 20% 1.5 RPM
500 20% 6 RPM
500
m = 0.25
0
0 0 2000 4000 6000 8000 10000 12000
0 0.5 1 1.5 2 2.5 3 3.5
Strain Rate [1/s] Pressure [Pa]

Fig. 8. Plot of shear stress against strain rate for the DEM data at all fill levels and ro- Fig. 9. Plot of shear stress against normal stress (pressure) for all the DEM data from all
tational speeds. Percentage in the legend refers to fill level. the calibration cases. Percentage in the legend refers to fill level.

distribution, particle shape, particle roughness, and potentially also (2017) and used by them to validate the use of a fluid simulation code
cohesiveness. In principle, individual particles could be characterized in (Smoothed Particle Hydrodynamics) to model a granular flow (as de-
detail, with a DEM model then used to determine the rheological scribed by Jonsén et al., 2012).
parameters, but it is more likely that an empirical approach would be
used. A triaxial shear test is the standard method for characterizing 5.2. Comparison of two-dimensional CFD with DEM kiln results
strength of granular materials, but a well-defined characterization ex-
periment in which the granular materials flow in a regime similar to To establish a level of confidence in results predicted from the ca-
that occurring in the kiln may also be needed for this purpose. An ex- librated CFD bed rheology model, a two-dimensional CFD model based
ample is the column collapse experiment described by Larsson et al. on the DEM simulations presented in Section 4 has been run. The solid

252
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

3500 20000
18000
3000
16000 Pressure
2500 14000
Pressure DEM
Viscosity [Pa-s]

12000
2000 shear
10% 3 RPM 10000
Shear DEM
1500 20% 3 RPM 8000
6000 visc
1000 30% 3 RPM
4000 Visc DEM
20% 1.5 RPM
500 2000
20% 6 RPM
0
0
0 2000 4000 6000 8000 10000 12000

PM

PM

PM
10

20

30
Pressure

5R

0R

0R
1.

3.

6.
Fig. 10. Plot of viscosity against normal stress (pressure) for all the DEM data from all the
2
calibration cases. Percentage in the legend refers to fill level.
1.8
1.6
0.6 strain
1.4
1.2 Strain DEM
0.5
1 Xvel
0.4 0.8 Xvel DEM
0.6 Slip
0.3 0.4
m

Slip DEM
0.2
0.2 10% 3 RPM 0
20% 3 RPM
0.1

PM

PM

PM
%
30% 3 RPM

10

20

30

5R

0R

0R
1.

3.

6.
0 20% 1.5 RPM
0 0.5 1 1.5 2 2.5 3 20% 6 RPM
3.5
Fig. 12. Comparison of the CFD and DEM data near the bottom of the kiln for 10, 20 and
Strain Rate [1/s] Model 30% fill level cases and 1.5, 3.0 and 6.0 rpm cases. Top plot compares pressure, shear and
viscosity, and bottom plot compares strain, x-component of velocity and wall slip.
Fig. 11. Plot of friction coefficient, m, (defined in Eq. (9)) against strain rate for all the
DEM data from all the calibration cases. Percentage in the legend refers to fill level.
Table 2
Comparison of average wall slip between the DEM and CFD models.
rheology model employed is given by Eqs. (9)–(14), together with the
modification from Fig. 11 described in Section 5.1. The cylinder dia- 10% fill 20% fill 30% fill 20% fill 10% fill
meter was 2 m, and was rotated at 1.5, 3 and 6 rpm with a fill of 20% 3 RPM 3 RPM 3 RPM 1.5 RPM 6 RPM
(matching the conditions for the original DEM simulation configura-
DEM 36.3% 16.2% 11.1% 17.2% 18.0%
tions). Fill levels of 10 and 30% were also run at 3 rpm to assess the CFD 20% 18.9% 12.6% 11.3% 37.8%
sensitivity to rotation rate.
Values predicted by the CFD model for the key variables char-
acterizing the kiln near its base are compared with DEM data at the difference is that the CFD model predicts a convex bed surface, while
same location in Fig. 12. These results show that the spot values are in the DEM bed shape is flat or slightly concave. Another difference is that
generally good agreement with the DEM data. Both the trend and ab- the volume fraction obtained from the DEM model near the bed surface
solute values are reasonably well predicted. The most significant dif- varies quite substantially near where the bed surface meets the kiln
ference is the viscosity value in the 1.5 rpm case. Wall slip velocities are walls, whereas the CFD model gives a sharp change in volume fraction
given in Table 2 along with the DEM values from Section 4. The CFD across the surface.
model gives reasonable predictions of the change in slip with fill level:
as the fill level increases from 10% to 20% to 30%, the slip decreases, as 5.3. Three-dimensional CFD model of rotary kiln
might be expected from the increased normal stress at the wall. How-
ever, the CFD model appears overly sensitive to changes in rotation A three-dimensional model of a rotary kiln using the DEM calibrated
speed. At 10% fill, the CFD simulations show a significant increase in rheology model was run to demonstrate that the model performs in a
slip as the speed is increased from 3 rpm to 6 rpm, whereas the DEM reasonable way in three-dimensions. The kiln is 15 m long with a dia-
simulations show decrease. Similarly, at 20% fill, the CFD simulations meter of 5 m rotating at a speed of 3 rpm. At the kiln inlet, solids en-
again show a significant increase in slip as the speed is increased from tered through a boundary patch covering the lower 4.6% of the area at
1.5 rpm to 3 rpm, whereas the DEM simulations show almost no change. a volume fraction of 1.0 and a velocity of 0.05 m s–1. Air entered at a
Further investigation of this effect is required. velocity of 2 m s–1 through the top 35% of the inlet face of the kiln. The
A complete compilation of results showing the bed behaviour and remaining area of the inlet face was treated as a “no slip” stationary
properties predicted by the CFD model are shown in Figs. A16–A33 of wall. At the kiln outlet, a constant pressure was assumed with Neumann
the Supplementary Data. In the plots, the outer coloured ring shows the conditions applied to the volume fractions and velocities. Initial con-
slip percentage between the bed and the kiln wall while the coloured ditions for the kiln were that the solids fill ratio was 24%, uniform
inner surface shows the value of the variable in the bed. For compar- along the kiln length, and that the bed and gas were stationary.
ison, plots from the DEM data for the various parameters are shown in The model was initially run for a real time of 21 s with an open end.
Figs. A1–A15 of the Supplementary Data. Plots of the bed surface coloured by axial velocity at various times are
Comparison for selected cases is given here. Direct side-by-side shown in Fig. 15. In order to simulate a more practically viable con-
comparison is shown in Fig. 13 for pressure and strain for the 10% fill figuration (given the relative short length of the kiln), a dam ring was
case, and, in Fig. 14, velocity predictions are compared for the 20% fill added at the discharge end of the kiln at 21 s (in simulation time), and
case. The results show reasonable agreement between the models, with the model was run for a further 28 s giving a total duration of 49 s.
the trends and distributions being in reasonable agreement. The main Fig. 16 shows the bed behaviour with the dam ring fitted. In practice,

253
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

1 pres 1 shear
9233 0
8573.5 -178.571
0.75 7914 -357.143
0.75
7254.5 -535.714
6595 -714.286
0.5 5935.5 -892.857
0.5
5276 -1071.43
4616.5 -1250
3957 -1428.57
0.25 0.25
3297.5 -1607.14
2638 -1785.71
1978.5 -1964.29
0 0
Y

Y
1319 -2142.86
659.5 -2321.43
0 -2500
-0.25 -0.25

-0.5 -0.5

-0.75 -0.75

-1 -1
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
X X

Fig. 13. Plots of Pressure (left) and shear stress (right) within the bed for rotary kiln with a fill of 10%. DEM results at the top and CFD results at the bottom. The flows are at near steady
state.

dam rings are often installed to increase retention time. processes, uniform heating of the bed is the aim. Transfer of heat
Fig. 17 shows the axial velocity field for the two cases. Without the through the bed occurs through several mechanisms. Advection of heat
dam ring, the axial velocity is higher than with the dam ring, though accompanies the movement of particles, and this heat transfer me-
the model did not reach steady state for this configuration, and the bed chanism is thus accounted for in both DEM and CFD models in a
was still emptying (i.e. evolving towards a low bed loading state) at 21 s straightforward way, provided the models have predicted the granular
of real time. After the dam ring is installed, the axial velocity reduces motion accurately. Particle-to-particle heat transfer occurs as a result of
with the flow being close to a steady state condition at time 49 s. conduction through particle contacts, conduction across the gas filled
Fig. 16(a) shows that the bed initially begins to build up after in- spaces between particles, and radiation between particle surfaces.
stallation of the dam, and after about 9 s (Fig. 16(b)), the bed surface is
approaching close to the steady state depth and shape.
For such a low aspect ratio (i.e. short) kiln, a discharge dam would 6.1. Heat transfer modelling using DEM
normally be needed because the truly steady state loading would be
quite small, and the retention time also correspondingly short, leaving As part of the hybrid CFD-DEM methodology presented in this
insufficient time for heat transfer to the charge and reaction. The dis- paper, DEM simulations were used to model as closely as possible heat
charge rings retard the particle exit increasing the residence time and transfer occurring within the bed microstructure via inter-particle and
thereby provides increased time for the heat transfer processes to occur. internal particle heat exchanges so as to provide closure relationships
Results presented in this section indicate that the three-dimensional for fine scale granular heat transfer for use in the macro-scale CFD
behaviour of a rotary kiln can be successfully modelled using this hy- model of the kiln. For a complete treatment of heat flow inside the bed,
brid technique. models for conduction, radiation and convection would be required.
Attention is focussed on conduction and radiation, the aspects for which
macro-scale (CFD) models require most input, given that CFD models
6. Heat transfer modelling in the granular bed routinely incorporate convective heat transfer. Convection (i.e. heat
transfer via the gas phase) has therefore not been included. It should be
Rotary kilns are typically operated at high temperature, with heat noted that Shi et al. (2008), using DEM simulations of heat transfer for a
transferred from a flame in the gas space to the granular bed surface, granular bed under kiln conditions, found that for particles with high
whence it is transported throughout the bed. In most industrial thermal conductivity (such as aluminium) particle-particle conduction

254
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

1 vx 1 vy
0.4 0.2
0.342857 0.171429
0.75 0.285714 0.142857
0.75
0.228571 0.114286
0.171429 0.0857143
0.5 0.114286 0.0571429
0.5
0.0571429 0.0285714
0 0
-0.0571429 -0.0285714
0.25 0.25
-0.114286 -0.0571429
-0.171429 -0.0857143
-0.228571 -0.114286
0 0
Y

Y
-0.285714 -0.142857
-0.342857 -0.171429
-0.4 -0.2
-0.25 -0.25

-0.5 -0.5

-0.75 -0.75

-1 -1
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
X X

Fig. 14. Plots of X component of velocity (left) and Y component of velocity (right) for rotary kiln with a fill of 20%. DEM results on the top line and CFD results on the bottom line. The
flows are at near steady state.

was dominant, whereas for particles with lower conductivity (such as 6.1.1. Conduction
glass beads) conduction through the interstitial gas (‘convection’) also Heat conduction in the DEM simulations considers both internal
played a significant role. However, radiation was neglected in the heat diffusion inside the particles and inter-particle heat exchange be-
computations of Shi et al. (2008). tween each contacting pair. The heat equation describing diffusion of
this thermal energy measured in terms of temperature T for a constant
thermal diffusivity, α, is:

Fig. 15. Change in bed surface as predicted by the 3D CFD model of the kiln, viewed from the exit with an open-ended kiln. The surface of the bed is coloured by its axial velocity. Times:
(a) 0.1 s; (b) 4.0 s; (c) 21.1 s.

255
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

Fig. 16. Change in bed surface as predicted by the 3D CFD model of the kiln, viewed from the exit with a dam ring at the kiln exit. The surface of the bed is coloured by its axial velocity.
Times: (a) 22.0 s; (b) 30.0 s; (c) 49.0 s.

approximation is not an unreasonable assumption except in the region

directly adjacent to the inter-particle contacts where significant spatial
variation related to the specific contacts can occur. To capture these
critical contact related aspects of the particle temperature field, a
thermal boundary layer model is associated with each such contact. Its
temperature is controlled by heat flux balances between heat transfer
into the volume of the particle (as represented by the 1-dimensional
overall temperature field within the particle) and heat transfer from the
other particle at the contact. The temperature difference between this
boundary layer region and the temperature of the outermost finite
difference shell provides an estimate of the local temperature gradient
within the particle created by the contact.
When two particles contact, their surfaces within the contact region
rapidly equilibrate to an intermediate temperature Tedge which gives
the same magnitude of local heat flux within both the contacting par-
ticles. This heat flux (given by Eq. (16)) is calculated using the local
temperature gradient at the surface of the particle at the contact point.
This is given by the temperature difference of Tedge from the boundary
layer temperature for that contact on that particle and the thickness of
the boundary layer region. The use of the boundary layer model and the
Fig. 17. A slice along the axis of the three-dimensional CFD kiln model showing predicted
longitudinal flow with (a) no dam at the exit at time 21.1 s; and (b) with a dam at the exit equilibrated surface temperature allows local temperature gradients to
at time 49 s. Red color indicates solids volume of 1.0 and blue indicates air; vectors in- be estimated in each particle at each particle contact. This is a very
dicate direction and magnitude of solid bed material. In the case with no dam, kiln is still important part of the model since these gradients are the dominant rate
emptying and the flow has not reached steady state; in the case with discharge dam, the controls for the thermal conduction in the granular media. If for ex-
flow is close to steady state. (For interpretation of the references to colour in this figure ample, only the outermost shell temperature in the 1-dimensional
legend, the reader is referred to the web version of this article.)
temperature solution was used then the temperature gradients at the
contact points would be significantly over-estimated leading to sig-
∂T nificant numerical diffusion of the heat.
= α∇2 T.
∂t (15) It is important to note that the contact area A used in the calculation
of the heat fluxes in Eq. (16) be the real physical contact area. In DEM
The amount of heat energy transferred across a contact area A, in a
in is common practice to decrease the spring stiffness so as to increase
time dt (typically this is the DEM timestep), for a given 1-D thermal
the timestep and accelerate the simulation. This is possible because the
gradient ∂T / ∂r is given by:
dynamical timescales for the collisions are much larger than the time-
∂T scales for stress wave propagation which does not need to be accurately
dQ = kA dt
∂r (16) predicted. If the contact springs are softened (as was the case in the
simulations here) then the geometric contact area calculated from the
where the conductivity k = αρc for density ρ and specific heat c.
DEM contacts are significantly enhanced by the softening. In this case
A finite difference approximation to the 1-D heat equation has been
the numerical estimates of the contact areas need to be rescaled to
used to model internal conduction within each particle in the DEM si-
correct physical ones using the ratio of the spring stiffness to the
mulation. Since the solution is an explicit one, a stability condition for
Young’s modulus of the material. If this is not done then the intra-
the computation timestep dt must be satisfied for a given spatial step dx
particle heat conduction will be significantly and unphysically too high.
(in the radial direction within the particle) such that
Particle temperatures and heat energies for the granular system are
dx 2 monitored to ensure that heat energy is conserved. Heat losses to the
dt ⩽ boundary are small for the timescales used to characterise the heat
2α (17)
transfer in the DEM model so there are not included.
The 1-dimensional heat conduction solution is used to track the
temperature within most of the particle. The use of such an

256
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

6.1.2. Radiation group of particles in a kiln the point contacts affect the heat transfer and
At the operating temperatures of rotary kilns, radiation is expected this needs to be considered. In principle, DEM simulations, such as
to be an important heat transfer mechanism. A radiation model for those described in Section 6.2, could be used to calculate the effective
radiative heat transfer between particles is used here. The heat flux q thermal diffusivity, α, which is related to the thermal conductivity, κ, by
due to radiative transfer between two surfaces h and e is given by: κ
α=
q = σε (Th4−Te4 ) (18) ρCp (20)

where σ is the Stefan-Boltzmann constant, ε is the emissivity, Th is the where Cp is the material specific heat. This thermal diffusivity would
surface temperature of the radiating particle, and Te is the temperature account for particle-particle conduction, particle-particle radiation and
of the surrounding environment. Heat transfer through radiative ex- small scale particle diffusion effects (such as jostling of particles, small
change is thus implemented by updating particle heat energies of par- particles slipping between large particles, etc). The determination of
ticles using the heat flux in Eq. (18). this diffusivity from DEM data is actually not straightforward, because
A number of assumptions have been made to simplify this pre- the diffusive effects are not easily separated from advective effects.
liminary model: shadowing of particles by smaller particles has not Consequently, as described below, we determine a best fit diffusivity for
been included (which is reasonable here because the particles used are the CFD model by comparing with the DEM results.
nearly mono-disperse); and radiative transfer for a given particle occurs CFD simulations were carried out for the same three cases analysed
only with a relatively small group of neighbours (which is reasonable by DEM simulation. A thermal diffusivity of 4 × 10−6 m2/s was found
since the bed is and remains densely packed). These neighbours can be to give the best match to the temperature profiles calculated for the
readily identified using the usual DEM search process for finding DEM model. Results for the case of a fill level of 20% and a rotation rate
neighbours. of 1.4 rpm are shown in Fig. 20. The locations of the red and blue re-
gions are determined by the bed flow part of the CFD model while the
6.2. Determination of thermal mixing using DEM simulations degree of diffusion (indicated by the orange to green to light blue areas)
is controlled by the heat transfer. The choice of the thermal diffusivity
Three two-dimensional DEM simulations of a rotary kiln were car- in the CFD model is made to ensure that the width of the thermal dif-
ried out to investigate the effects of fill level and rotation rate on the fusion layer and its rate of expansion is as similar as possible to that
predicted thermal mixing and heat transfer. Based on the results of shown by the DEM model. The plots in Fig. 20 show good agreement
preliminary simulations, we expect that the heat transfer from ad- between the predicted CFD and DEM temperature profiles in the bed.
vective mixing will dominate over particle-particle conduction since Only short term behaviour needs to be compared since the longer-term
particles do not remain in contact for long and the mixing is strong in mixing will be dominated by advection, which is already treated ef-
such systems. For typical kiln conditions (at temperatures of fectively the same way in both CFD and DEM models. In future, it may
400–800 °C), radiation should be the dominant heat transfer me- be possible to extract point-by-point values of the effective thermal
chanism after advection. Radiation and internal conduction were diffusivity from the DEM results, since it is expected that there will be
therefore included whilst inter-particle conduction was omitted. some spatial and directional variation.
The first case has a fill level of 7% and a rotation rate of 1.4 rpm.
Bed temperatures were initially set to 400 °C (673 °K) in the left half of 7. Discussion
the bed and 800 °C (1073 °K) in the right half. Temperature distribu-
tions for the simulation are shown at various times in Fig. 18 and Rotary kilns are used in several types of minerals processing op-
particles are coloured by average temperature (keeping in mind that the erations, and in related industries such as cement manufacture. Though
temperature varies within each particle according to the intra-particle the equipment is robust and the processes generally stable, issues such
conduction). Mixing proceeds as the hot and cold regions wrap around as ring accretion formation and minimization of energy consumption
each other with radiative heat transfer across the void spaces within the have often proven to be hard to tackle. In this paper, we have described
bed). By three minutes, the bed has thermally equilibrated. a hybrid hierarchical CFD-DEM modelling technique that can be used to
The second case has a fill level of 20% and a rotation rate of address the complexity of heat and mass transfer processes in rotary
1.4 rpm. Temperature distributions are shown at various times in kilns in a predictive way.
Fig. 19. Again recirculation of the hot and cold regions occurs rapidly The motivation for developing this new technique is that the con-
with the cold and hot regions wrapping around each other. Radiation siderable length of many rotary kilns, combined with the small particle
then transfers heat between exposed surfaces across each void space in size (ilmenite feed is typically less than 1 mm) and long process time-
the bed. The central region equalises first and by four minutes the full scales mean that practical design calculations using DEM to track every
bed has reached thermal equilibrium. The increase in bed volume particle in an industry scale kiln is still computationally intractable for
makes the central region more visible here than was the case for the 7% design purposes. As an example, fully resolved DEM simulation of a
fill. 50 m long kiln of inner diameter 3 m with a bed of 1 mm particles of
The third case used a 7% fill level but with a rotation rate of volume fill 10% would require 1011 particles. Even using softened si-
0.7 rpm. By four minutes, the bed has reached thermal equilibrium. mulation of collisions, many more than 108 timesteps would be needed
Since the rotation speed is lower here than for the 1.4 rpm case, reduced to reach thermal steady state, and this computation is not presently
mixing and thus a slower advection rate are expected. feasible even on a massively parallel machine. On the other hand, since
kiln internals do not involve much fine detail, a reasonably high re-
6.3. CFD simulations of thermal mixing in the rotary kiln solution 3D CFD simulation of such a bed could be carried out using a
million nodes, with timesteps much longer than that required in the
The CFD model of Section 5 was extended to include heat transfer DEM simulation. Such a CFD simulation is now quite routine with run
within the bed by solving the thermal energy equation within the times easily affordable for industrial application, if a smart technique,
CFX4.4 framework (AEA Technology, 2001) such as that of Schwarz and Lee (2007) is used to account for the long
thermal timescale. There would of course also be the additional over-
∂ρh
+ ∇ ·(ρ uh) = ∇ ·κ∇T head of the validation and calibration studies described in this paper,
∂t (19)
but the run-times for these are relatively short since they can mostly be
where h is the enthalpy, ρ is the density, u the velocity, κ the thermal carried out in 2D. Typical run-times for both DEM and CFD cases de-
conductivity. For a solid material κ is a material property but for a scribed in this paper are of order an hour on a single processor machine,

257
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

Fig. 18. DEM simulation of a kiln with 7% fill level and a rotation rate of 1.4 rpm. Particles are coloured by their average temperature (in units of K) and are shown at t = 0, 0.5, 1, 1.5, 2
and 3 min.

without any particular effort taken to optimize them. gaseous freeboard including combustion processes, and can also be
It should be remembered in this context that design simulations readily extended to incorporate heat transfer through the shell. On the
aimed at process, feedstock and equipment optimization typically re- other hand, DEM models are most naturally applied to modelling the
quire many iterations with different parameters, and also need to ac- details of particle dynamics in the granular bed. In the hybrid hier-
count for flame structure, radiative and convective heat transfer, solids archical approach, DEM simulations have been used to derive con-
reactions and so on. Extensively validated models, techniques and stitutive relationships for the granular bed, and these were then used to
software already exist for dealing with many of these phenomena validate a solids rheology model that is in turn, used in a CFD model of
within a CFD framework, but there are gaps in the characterization of an entire industrial-scale kiln, including radiation and convective heat
bed rheology and heat transfer that limit the application of the tech- transfer. This approach could also be extended to including solids re-
niques to a rotary kiln. actions. It should be emphasized that in the present hierarchical ap-
The hierarchical technique exploits the individual strengths of both proach, the CFD and DEM modelling are carried out sequentially, rather
DEM and CFD methods in a sequential manner in which DEM modelling than in the close-coupled concurrent technique commonly known as
is used to determine constitutive relationships which are then used in “CFD-DEM”, which is possible because of the large separation in the
CFD modelling. CFD models are naturally suited to dealing with the timescale of the particle scale processes and the kiln scale ones.

258
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

Fig. 19. DEM simulation of the pilot kiln with 20% fill level and rotation rate of 1.4 rpm. Particles are coloured by average temperature (in units of K) and are shown at t = 0, 0.5, 1, 2, 3
and 4 minutes.

The DEM simulations have identified important qualitative and PSD range, but is largely insensitive to macroscopic quantities of the
quantitative features of the effective bed viscosity and bed-shell friction bed such as fill level and rotation rate. This suggests that the details of
that need to be incorporated into continuum models suitable for ro- the bed microstructure are important for understanding the friction
tating kilns. The shear viscosity distribution is qualitatively similar for force between the bed and kiln wall. Shear viscosity and bed-shell
all cases studied in that its maximum values lie at the bottom of the bed friction are both sensitive to particle shape and further work in un-
close to the kiln wall. It varies subtly with some of the particle scale derstanding the effect of non-circular shape is desirable.
parameters varied. Models used in two-fluid CFD simulations to describe two-phase
The maximum viscosities have been found to be quite sensitive to all suspensions of particles (e.g. the kinetic theory of granular flow) are not
the parameters studied. The highest viscosity case was found to occur suitable for the granular bed. The viscous-plastic Drucker-Prager model
for a bed of non-circular particles with a large fill level and low rotation has been tested for suitability as the constitutive relation in the CFD
rate. The magnitude of the shear viscosity increases linearly with bed simulation stage, and with inclusion of a modified Coulombic friction
depth (and hence bed weight) over the range investigated, and also yielded a model which gave satisfactory agreement with the DEM re-
decreases with decrease in particle size. Bed-shell friction is found to be sults over the range of rotation speeds. Thus, the DEM simulations have
strongly dependent on particle size and shape, and the breadth of the pointed towards a bed rheology model that appears to have general

259
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

Fig. 20. CFD predicted temperature (left) profiles compared to the DEM temperature (right) profiles for a kiln rotating at 1.4 rpm and a 20% fill level at 10 s intervals.

applicability for granular beds in the avalanching regime. for a particular granular solid for an actual industrial situation.
It is likely that the constitutive relation thus would need to be ca- Parameters in the model will depend in a complex way on details of
librated empirically to determine the specific constitutive parameters particle size distribution, particle shape, particle roughness, and

260
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

potentially also cohesiveness. be successfully used for industrial simulation of large-scale rotary kilns.
An isothermal three-dimensional CFD model of a rotary kiln with
flow-through of solids was run using the DEM calibrated rheology Acknowledgements
model to demonstrate that the model performs in a reasonable way in
three-dimensions. The kiln was 15 m long with a diameter of 5 m ro- The authors thank the companies that sponsored the AMIRA man-
tating at a speed of 3 rpm, and the modelling was demonstrated both aged project, P581: Rotary Kiln Technology, which partially supported
with and without a dam-ring. The results indicated that the three-di- this work.
mensional solids movement in a rotary kiln can be successfully mod-
elled using this hybrid technique. Appendix A. Supplementary material
The second part of the paper was devoted to the heat transfer as-
pects of the modelling approach. Most of the heat transfer processes Supplementary data associated with this article can be found, in the
within an industrial rotary kiln can be simulated in a CFD framework online version, at http://dx.doi.org/10.1016/j.mineng.2018.01.035.
used well validated techniques. Most of the heat transfer within the
granular bed itself takes place as a result of advection of the solids, and References
provided the constitutive relationships for bed viscosity yield reliable
predictions of granular flow velocities, this dominant component of AEA Technology, 2001. CFX-4.4: Solver. AEA Technology, Harwell Laboratory,
heat transfer will be accurately simulated in the CFD framework. Gas- Oxfordshire, UK.
Boateng, A.A., Barr, P.V., 1996. A thermal model for the rotary kiln including heat
solids heat transfer can be also be adequately captured using well-es- transfer within the bed. Int. J. Heat Mass Transfer 39 (10), 2131–2147.
tablished correlations for local heat transfer. Remaining phenomena Boateng, A.A., Barr, P.V., 1997. Granular flow behaviour in the transverse plane of a
such as particle-particle conduction and diffusion of heat due to particle partially filled rotating cylinder. J. Fluid Mech. 330, 233–249.
Bui, R.T., Simard, G., Charette, A., Kocaefe, Y., Perron, J., 1995. Mathematical modelling
velocity fluctuation (e.g. ‘jostling’) can also be captured but, since they of the rotary coke calcining kiln. Can. J. Chem. Eng. 73, 534–545.
depend on particle-scale processes, may need to be calibrated. Chaudhuri, B., Muzzio, F.J., Tomassone, M.S., 2010. Experimentally validated compu-
2D simulations of heat transfer in a granular bed using DEM simu- tations of heat transfer in granular materials in rotary calciners. Powder Tech. 198,
6–15.
lations were used to check and calibrate the thermal components of the
Cleary, P.W., Metcalfe, G., Liffman, K., 1998. How well do discrete element granular flow
CFD simulations. The contribution of particle-scale phenomena to the models capture the essentials of mixing processes. Appl. Math. Modelling 22,
overall thermal mixing in the granular bed was modelling by means of a 995–1008.
Cleary, P.W., Sawley, M.L., 2002. DEM modelling of industrial granular flows: 3D case
thermal diffusivity fitted to provide the best match to the DEM thermal
studies and the effect of particle shape on hopper discharge. Appl. Math. Modelling
predictions. Acceptable agreement was then obtained between the DEM 26, 89–111.
and CFD simulations in terms of the evolution of the temperature dis- Cleary, P.W., 2004. Large scale industrial DEM modelling. Eng. Comput. 21, 169–204.
tribution within the bed, providing confidence that this procedure Cleary, P.W., 2008. The effect of particle shape on simple shear flows. Powder Tech. 179,
144–163.
could be successfully used for industrial simulation. Cleary, P.W., 2009. Industrial particle flow modelling using DEM. Eng. Comput. 26,
698–743.
8. Conclusion Eghlimi, A., Benito, R. Golab, C., 1999. The CFD investigation of flash dryer and rotating
kiln design. In: Proc. 2nd Inter. Conf. on CFD in the Minerals and Process Industries,
6–8 December, Melbourne, Australia, pp. 455–460.
Rotary kilns, as used in minerals processing operations and cement Feng, Y.Q., Swenser-Smith, T., Witt, P.J., Doblin, C., Lim, S., Schwarz, M.P., 2012. CFD
manufacture, involve complex multiphase heat and mass transfer pro- modelling of gas-solids flow in an internally circulating fluidized bed. Powder Tech.
219, 78–85.
cesses. Optimization of the equipment would ideally be guided by nu- Finnie, G.J., Kruyt, N.P., Ye, M., Zeilstra, C., Kuipers, J.A.M., 2005. Longitudinal and
merical modelling, but simulation models so far have tended to con- transverse mixing in rotary kilns: a discrete element method approach. Chem. Eng.
centrate on either the freeboard aspects such as flame structure and Sci. 60, 4083–4091.
Hanes, D.M., Inman, D.L., 1985. Experimental evaluation of a dynamic yield criterion for
heat transfer, or on the granular bed movement. A new generic hier- granular fluid flows. J. Geophys. Res. 90, 3670–3674.
archical approach has been proposed in this paper, sequentially com- Hassan, M., Schwarz, M.P., Feng, Y.Q., Rafique, M., Witt, P., Huilin, L., 2016. Numerical
bining 2D Discrete Element Method (DEM) simulations of slices of the investigation of solid circulation flux in an internally circulating fluidized bed with
different gas distributor designs. Powder Technol. 301, 1103–1111.
bed and 3D two-phase computational fluid dynamics (CFD) models of
Jonsén, P., Pålsson, B.I., Häggblad, H.Å., 2012. A novel method for full-body modelling of
the entire kiln, encompassing both bed and gaseous freeboard. The grinding charges in tumbling mills. Minerals Eng. 33, 2–12. http://dx.doi.org/10.
hierarchical technique exploits the individual strengths of both DEM 1016/j.mineng.2012.01.017.
and CFD methods in a sequential manner in which DEM modelling is Karlsson, T., 1995. A Finite Element Method for Simulation of Granular Flows, Tekniska
Högskolan I Luleå Research Report TULEA 1995:10. Luleå University of Technology,
used to determine constitutive relationships which are then used in CFD Sweden.
modelling. Karlsson, T., Klisinski, M., Runesson, K., 1998. Finite element simulation of granular
As a demonstration of the approach, 2D DEM simulations have been material flow in plane silos with complicated geometry. Powder Technol. 99, 29–39.
http://dx.doi.org/10.1016/S0032-5910(98)00087-4.
used to derive mean solids velocities in a bed for various rotation speeds Kritzinger, H.P., Kingsley, T.C., 2015. Modelling and optimization of a rotary kiln direct
and particle sizes, and these were then used to calibrate and validate reduction process. J. S. Afr. Inst. Min. Metall. 115, 419–424.
the solids rheology model used in the subsequent CFD model. A mod- Kwapinska, M., Saage, G., Tsotsas, E., 2006. Mixing of particles in rotary drums: a
comparison of discrete element simulations with experimental results and penetra-
ified Coulombic friction for the bed in the CFD model was shown to give tion models for thermal processes. Powder Technol. 161, 69–78.
satisfactory agreement with the DEM bed results over a range of rota- Küssel, U., Abel, D., Schumacher, M., Weng, M., 2009. Modeling of rotary kilns and ap-
tion speeds. In addition, the DEM simulations have identified important plication to limestone calcination. In: Proc. 7th Modelica Conference, Como, Italy,
Sep. 20–22, 2009, pp. 814–822, 10.3384/ecp09430084.
qualitative and quantitative features of the effective bed viscosity and Larsson, S., Gustafsson, G., Häggblad, H., Jonsén, P., 2017. Experimental and numerical
bed-shell friction that need to be incorporated into a continuum model. study of granular material flow using smoothed particle hydrodynamics (SPH):
The capability of the resultant CFD model was demonstrated by 3D Application in handling of potassium chloride (KCl). In: Computational Modelling
2017, Minerals Engineering International, Falmouth, June, 2017.
simulation of the bed of a large rotary kiln.
Liu, T.Y., Schwarz, M.P., 2009. CFD based multiscale modelling of bubble–particle col-
Similarly, thermal mixing simulations carried out using the 2D DEM lision efficiency in a turbulent flotation cell. Chem. Eng. Sci. 64, 5287–5301. http://
model for the granular bed were used to validate a CFD model of heat dx.doi.org/10.1016/j.ces.2009.09.014.
transfer. The contribution of particle-scale phenomena to the overall Luding, S., 2009a. Constitutive relations of dense granulates with friction and adhesion
from DEM. In: Oñate, E., Owen, D.R.J. (Eds.). Inter. Conf. on Particle-Based Methods
thermal mixing in the granular bed was modelling by means of a Particles 2009. CIMNE, Barcelona, 2009, ISBN: 978-84-96736-82-5, pp. 159–162.
thermal diffusivity fitted to provide the best match to the DEM thermal Luding, S., 2009b. From molecular dynamics and particle simulations towards con-
predictions. The acceptable agreement obtained between the DEM and stitutive relations for continuum theory. In: Koren, B., Vuik, K. (Eds.). Advanced
Computational Methods in Science and Engineering. Springer Series Lecture Notes in
CFD thermal simulations provides confidence that this procedure could

261
P.J. Witt et al. Minerals Engineering 119 (2018) 244–262

Computational Science and Engineering #71, Springer, Berlin, ISBN: 978-3-642- Schwarz, M.P., Liu, T., 2005. Multi-scale modelling of processes in the minerals industry.
03343-8, 2009, pp. 453–492. Computational Modelling 05, Cape Town, November, 2005.
Lun, C.K.K., Savage, S.B., Jeffrey, D.J., Chepurniy, N., 1984. Kinetic theories for granular Schwarz, M.P., Lee, J., 2007. Reactive CFD simulation of an FCC regenerator. Asia-Pac. J.
flow: inelastic particle in Couette flow and slightly inelastic particles in a general flow Chem. Eng. 2, 347–354. http://dx.doi.org/10.1002/apj.064.
field. J. Fluid Mech. 140, 223–256. Schwarz, M.P., Koh, P.T.L., Verrelli, D.I., Feng, Y.Q., 2015. Sequential multi-scale mod-
Mastorakos, E., Assias, A., Tsakiroglou, C.D., Goussis, D.A., Burganos, V.N., Payatakes, elling of mineral processing operations, with application to flotation cells. Miner.
A.C., 1999. CFD predictions for cement kilns including flame modelling, heat transfer Eng. 90, 2–16.
and clinker chemistry. Appl. Math. Modelling 23, 55–76. Shi, D., Vargas, W.L., McCarthy, J.J., 2008. Heat transfer in rotary kilns with interstitial
Metcalfe, G., Graham, L., Zhou, J., Liffman, K., 1999. Measurement of particle motions gases. Chem. Eng. Sci. 63, 4506–4516.
within tumbling granular flows. Chaos 9, 581–593. Syamlal, M., O’Brien, T.J., 1989. Computer simulation of bubbles in a fluidized bed. In:
Nedderman, R.M., 1992. Statics and Kinematics of Granular Materials. Cambridge Uni. AIChE Symp. Series No. 85, pp. 22–31.
Press, Cambridge. Thornton, G.J., Batterham, R.J., 1982, The transfer of heat in kiln, In: Proceedings of
Nicholson, T.A., 1995. Mathematical Modelling of the Ilmenite Reduction Process in Chemeca 82, 10th Aust Chem Eng Conf, Sydney, (Institution of Engineers, Aust),
Rotary Kilns, PhD Dissertation, University of Queensland. National Conference Publication No. 82/9, pp. 260–266.
Pereira, G.G., Sinnott, M.D., Cleary, P.W., Liffman, K., Metcalfe, G., Sutalo, I.D., 2011. Thurlby, J.A., 1988. A dynamic mathematical model of the complete grate/kiln iron-ore
Insights from simulations into mechanisms for density segregation of granular mix- pellet induration process. Met Trans B 19B, 103–112.
tures in rotating cylinders. Granular Matter 13, 53–74. Vun, S., Naser, J.A., Witt, P.J., 2010. Extension of the kinetic theory of granular flow to
Pereira, G.G., Cleary, P.W., 2013. Radial segregation of multi-component granular media include dense quasi-static stresses. Powder Technol. 204, 11–20. http://dx.doi.org/
in a rotating tumbler. Granular Matter 15, 705–724. 10.1016/j.powtec.2010.07.001.
Pereira, G.G., Tran, N., Cleary, P.W., 2014. Segregation in combined size and density Witt, P.J., Perry, J.H., Schwarz, M.P., 1998. A numerical model for predicting bubble
varying binary granular mixtures in a slowly rotating tumbler. Granular Matter 16, formation in a 3D fluidized bed. Appl. Math. Modelling 22, 1071–1080.
711–732. Witt, P.J. Mittoni, L.J., Schwarz, M.P., 1999. Application of CFD to industrial fluidised
Pereira, G.G., Cleary, P.W., 2017. Segregation due to particle shape in sheared binary bed processes. In: Werther, J.,(Ed.) Proceedings of the 6th International Conference
granular media. Granular Matter 19–23. http://dx.doi.org/10.1007/s10035-017- on Circulating Fluidized Beds, CFB-6, 22–27 August, 1999, Würzburg, Germany, pp.
0708-7i. 211–212.
Rahaman, M.F., Naser, J., Witt, P.J., 2003. An unequal granular temperature kinetic Zhang, Y., Campbell, C.S., 1992. The interface between fluid-like and solid-like behaviour
theory: description of granular flow with multiple particle classes. Powder Technol. in two-dimensional granular flows. J. Fluid Mech. 237, 541–568.
138, 82–92. http://dx.doi.org/10.1016/j.powtec.2003.08.050.

262