## Abstract

The purpose of this work was the direct numerical simulation of heat and fluid flow by granular mixing in a horizontal rotating kiln. To model particle behaviour and the heat and fluid flow in the drum, we solve the mass conservation, momentum and energy conservation equations directly on a fixed Eulerian grid for the whole domain including particles. At the same time the particle dynamics and their collisions are solved on a Lagrangian grid for each particle. To calculate the heat transfer inside the particles we use two models: the first is the direct solution of the energy conservation equation in the Lagrangian and Eulerian space, and the second is our so-called linear model that assumes homogeneous distribution of the temperature inside each particle. Numerical simulations showed that, if the thermal diffusivity of the gas phase significantly exceeds the same parameter of the particles, the linear model overpredicts the heating rate of the particles. The influence of the particle size and the angular velocity of the drum on the heating rates of particles is studied and discussed.

## 1. Introduction

Heat transfer in particulate flows is an essential component in many industrial applications, such as gasification in fixed-bed and fluidized-bed reactors or heterogeneous catalysis in chemical reactors. However, the heat transfer in granular mixing is still not well understood [1,2], which makes it difficult to predict the optimal input parameters to obtain better effectiveness of the mixing process in such applications. Recently, with a significant increase in computer power, and especially with the further development of multi-core processors, molecular dynamic or discrete element modelling (DEM) became the state of the art in the simulation of granular flows. In spite of significant progress in the development of DEM to simulate flow, mixing and heat transport in particulate flows, basically, no gas flow is taken into account when modelling granular flows in rotating vessels (e.g. [1]). The heat transfer between and inside the particles is modelled using linear models with semi-empirical relations for the effective heat transfer coefficient. Recently, Shi *et al*. [3] adopted a coupled DEM–computational fluid dynamics technique to simulate conductive and convective heat transfers in a rotary kiln. The gas flow was governed by the Navier–Stokes equation written in a Euler–Euler formulation, and the particle dynamics was treated using DEM. The temperature evolution of the particles was calculated using a standard linear model in which the effective heat transfer coefficient was calculated from a Nusselt number relation. All models cited above use the assumption that the temperature within each particle is spatially uniform. Motivated by this fact, we present a direct numerical simulation (DNS) model of the fluid flow and the heat transfer in particulate flows which is capable of modelling transient heat transfer in both gas and solid phases, respectively. As an application of our method, we consider granular mixing in a system representing rotary drums and kilns similar to the study of Chaudhuri *et al*. [1] and Shi *et al*. [3].

## 2. Problem and model formulation

The drum system under consideration consisted of 30 coal particles of 0.028 m diameter in a cylindrical vessel of diameter 0.238 m. The cylinder was filled with air and the coal particles and rotated with an angular velocity of 30 r.p.m. The wall temperature was set to 400 K. The initial temperature of the air and coal particles was 300 K. The drum started rotating from rest. In order to study the influence of the number of particles and the rotational velocity of the kiln on the system behaviour, we increased the number of particles to 100 (*d*_{p}=0.014 m) and decreased the angular velocity to 15 r.p.m. The principal scheme of the system is shown in figure 1. Definitions of the symbols used in figure 1 are given in table 1. The material properties of air and the coal particles are summarized in table 2.

To model the particle behaviour and the heat and fluid flow in the drum, we solved the mass conservation, momentum and energy conservation equations directly on a fixed Eulerian grid for the whole domain including particles. At the same time the particle dynamics and their collisions were solved on a Lagrangian grid for each particle. Particle–fluid interaction was treated using both grids. To proceed with the governing equations, the following basic assumptions were made: the fluid was incompressible and Newtonian; the particles were cylindrical in shape; the particle collisions were inelastic, where the surface was rough, and were modelled with a hard sphere model; the properties of the fluid and solid were temperature independent.

### (a) Fluid flow and heat transfer

Taking into account the assumptions made above, the conservation equations for the fluid phase take the form 2.1 2.2 and2.3

The gas flow and the particles are coupled by setting up a no-slip boundary condition on the particle interfaces.

The heat transfer in the gas and in the particles is calculated directly using equation (2.3). However, the heat transfer modelling inside each particle is done in Lagrange space. Thus, for the calculation of the temperature inside the particles the velocity ** u** in equation (2.3) has to be calculated as follows:
2.4It should be noted that the rotational velocity guarantees the effect of temperature tracking by particle rotation. To treat the heat transfer by the translational motion, we used a simplified solid body tracking technique, where the temperature field inside the particles was aligned with the ‘solid’ cells. Thus, the temperature field in each particle was tracked together with the volume fraction of the solid.

For comparison with the direct modelling of the temperature field inside the solid particles, we used our so-called linear model. This model assumes that the temperature field is uniform within the particles. This is equal to an infinite thermal conductivity *λ*_{s}. The temperature *T*_{p,i} is calculated for each particle as follows:
2.5where the integral on the right-hand side (r.h.s.) of equation (2.5) defines the total heat flux between the particle and the gas phase; the second term on the r.h.s. accounts for the conductive heat transport between the particles and between the particle and the wall (*j*=0); and where Δ*x* is the grid cell length. To simulate the contact resistance at the contact area *S*_{ij} between two particles we used *λ*_{s} from table 2.

### (b) Particle motion and collisions

The particle motion is modelled using the Lagrange approach, solving the impulse conservation equations for every particle. In particular, the particle translation and rotation are governed by the following equations:
2.6Here, *m*_{i} is the mass of the particle, *U*_{i} and *Ω*_{i} are the translation and rotation velocity, respectively, and *G*_{i}, *F*_{i} and *T*_{i} are the force of gravity, the hydrodynamic force and moment, respectively. The bulk flow and the particles are coupled via the hydrodynamic force *F*_{i} and the torque *T*_{Mi} is calculated directly for each particle using the following equations:
and
2.7where ** r** denotes the coordinates of the surface and

*r*_{i}the coordinates of the centre of mass of particle

*i*. The forces

*F*_{i}

^{collis}and

*T*_{i}

^{collis}in equation (2.6) are responsible for the influence of the interparticle collisions on the particle trajectories. In this work, these forces are modelled implicitly by calculating the translational and rotational velocities of the particle after the bilateral collision as follows [4]: 2.8Here

*ξ*=2((

*m*

_{1}m

_{2})/(

*m*

_{1}+

*m*

_{2})), and

*I*

_{i}is the mass moment of inertia, in the case of cylinders . If a collision with the wall occurs the wall is handled as particle

*j*with infinite mass. The collision impulse

*J*_{ij}between particle

*i*and

*j*is given by 2.9where

**is the collision direction or the direction between the centres of mass,**

*k***is the relative velocity at the contact point and**

*h**η*

_{1}and

*η*

_{2}are the dissipation parameters, which are calculated as follows: 2.10Here,

*w*

_{i}is the mass fraction of each particle (

*m*

_{i}/(

*m*

_{1}+

*m*

_{2})) and

*k*is the coefficient of the moment of inertia

*I*, which is 1/2 for cylindrical particles. The coefficients

*α*and

*β*used in equation (2.10) are the restitution coefficient in the normal direction and the surface roughness, respectively. The coefficient

*α*takes values between zero and unity. If

*α*=1 or

*α*=0, the collisions are elastic or inelastic, respectively. The surface roughness can be changed within the limits of −1≤

*β*≤1. An increase in

*β*corresponds to the increase in the surface roughness and, as a result, tangential velocity decreases. If

*β*becomes greater than 0 the direction of the tangential velocity is reversed. For further information see Rao & Nott [4]. In this work, we used

*α*=0.9 and

*β*=0. To avoid multiple collisions we used a small time step. However, when multiple collisions were detected, we made up to 100 internal collision iterations in the time step, to solve the situation in a sequence of bilateral collisions.

## 3. Numerics and code validation

The transport equation has been discretized by a finite volume-based method. The semi-implicit method for pressure-linked equations (SIMPLE) algorithm with a collocated–variable arrangement was used to calculate the pressure and the velocities [7].

The set-up of the no-slip and Dierichlet boundary conditions on the particle surface was done using the original implicit fictitious boundary method, which was introduced by Zienkiewicz [8]. The finite volume formulation was developed by Ananiev *et al*. [9]. The main idea of this method is the direct implementation of the boundary conditions at the interface cells through the modification of the coefficient matrix. To illustrate the method, we consider a five-point difference scheme, e.g. for equation (2.2), given by
3.1where [*A*] is the sparse matrix of known coefficients, *u*_{x} is the column vector of unknowns and *Q*_{ux} is a r.h.s. vector. To fix the values in the immersed region, the values of [*A*] were set at unity and the corresponding values were subtracted from *Q*_{ux}.

The validation of the code in the sense of the particle collisions was done in comparison with the recent experiment of Zhang *et al*. [5]. Following the experiment [5], we considered the collapse of a stack of 33 aluminium cylinders owing to gravity. A uniform grid with 1000 × 308 control volumes (cv) was used. The time step was equal to 10^{−5} s. The restitution coefficient *α* was zero and *β*=−0.5. The surrounding medium was air. The comparison between experimental and numerical data is presented in figure 2*a*, which shows the cylinder positions at different times. Figure 2*b,c* depicts the comparison of the time histories of the mean centre of mass position in the *x*- and *y*-directions. The convective heat transfer was validated against the experiment of Kuehn & Goldstein [6]. Figure 2*d* shows the temperature profile on the symmetry line in the top part of the system. Good agreement can be seen for both validation cases.

## 4. Results and discussion

In this section, we discuss the results of the DNSs of coal particles heated in a rotary kiln. The simulations were performed using a grid with 400×400 control volumes in the *x*- and *y*-directions. The time step was set to 1×10^{−4} s. The simulations were stopped after the modelling at 45 *s*. Thus, the total number of time steps corresponds to 45×10^{4}. To study the dynamics of the particle heating, we used the particle volume-averaged and gas volume-averaged temperatures, which were calculated as follows:
4.1and
4.2The physical meaning of *T*^{p}_{mean} and *T*^{g}_{mean} is the mean temperature of all particles and of the gas phase in the rotary kiln, respectively. The time evolution of *T*^{g}_{mean} is shown in figure 3*a*. It can be seen that the gas phase is heated rapidly within one to two rotations of the kiln. After the temperature *T*^{g}_{mean} reaches its maximum, it stays approximately constant at a nearly equal level. The temperature oscillations are caused by the vortices induced by the particle movement. The heating history of the particles in the form of *T*^{p}_{mean} is depicted in figure 3*b*. It can be seen that the particles are heated more slowly than the gas phase. This can be explained, first, by the *thermal diffusivity* of air being 276.06 times higher than that of coal. Second, the heat transfer between the particles and the kiln walls is limited through the contact points. The comparison of *T*^{p}_{mean} calculated by the DNS and that calculated by the linear model shows that independently of the total number of particles in the kiln the linear model overpredicts the heating rate of the particles. The differences between the models are explained by the fact that the linear model is equivalent to a DNS model with infinite thermal conductivity of the coal particles. Additional simulations demonstrated that the DNS model produces results very close to those of the *linear* model when the thermal conductivity of the coal particles has a 100 times higher value than the original values of *λ*_{s} given in the table 2.

To demonstrate the heat transfer between the kiln walls, gas phase and particles, we plot the non-dimensional temperature ((*T*−*T*_{0})/(*T*_{w}−*T*_{0})) distribution in the kiln at different times in figure 4*a*. In order to explore the temperature distribution inside the particles, figure 4*b* presents the same snapshots with refined scaling. It can be seen that in comparison with the high temperatures of the gas phase in the upper part of the kiln the interstitial gas between the particles stays cold or is quickly cooled down by the particles. The temperature is not homogeneously distributed within the particles, owing to their low *thermal diffusivity*.

In particular, because of kiln rotation and the friction force between the particles and the kiln walls the particles move in a so-called cascading or ‘avalanche’ regime. Thus, the particles are heated by two basic mechanisms. The first corresponds to the heat transfer through direct contact between the particles and the kiln wall. The second is the convective heat transfer between the particles and the gas in the upper part of the particulate bed.

Figure 5 shows a series of snapshots of the temperature distribution in the kiln and in a particle (zoomed view) during one revolution of the kiln. It can be seen that the particle is heated at the contact point with the wall. When the particle arrives at the top point of the granular bed it starts rolling down over the other particles. Thus, inside the kiln we have counterclockwise rotation of the granular bed and clockwise rotation of the gas phase in the upper part of the kiln.

## 5. Conclusion

The mixing and heating of coal particles in the rotating kiln were investigated numerically using DNSs. The model developed is capable of modelling transient heat transfer in the gas phase and inside the particles. The particle collisions were taken into account. The analysis of the simulation results showed that in our case the heat transfer between particles and the wall plays a significant role in the heating of particles. The use of a linear model, where the heat transfer inside the particles is neglected, overpredicts the heating rate of the particles significantly and cannot be used in this case. A detailed analysis of the temperature distribution inside the particles showed a non-homogeneous distribution due to the low *thermal diffusivity* of the particles.

## Acknowledgements

The authors appreciate the financial support of the Government of Saxony and the Federal Ministry of Education and Science of the Federal Republic of Germany as a part of the Centre of Innovation Competence VIRTUHCON.

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society