Downloaded from http://rsta.royalsocietypublishing.org/ on January 15, 2017

Phil. Trans. R. Soc. A (2011) 369, 2574–2583 doi:10.1098/rsta.2011.0046

Direct numerical simulation of particulate flows with heat transfer in a rotating cylindrical cavity BY ROBIN SCHMIDT*

AND

PETR A. NIKRITYUK

CIC VIRTUHCON, Technische Universität Bergakademie Freiberg, Reiche Zeche/Fuchsmühlenweg 9, 09599 Freiberg, Germany 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. Keywords: heat transfer; particles; rotation; Euler–Lagrange

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 *Author for correspondence ([email protected]). One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

2574

This journal is © 2011 The Royal Society

Downloaded from http://rsta.royalsocietypublishing.org/ on January 15, 2017

DNS of particulate flows

2575

Tw

W

g

R0 T0

d

Figure 1. Scheme of set-up.

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 (dp = 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. Phil. Trans. R. Soc. A (2011)

Downloaded from http://rsta.royalsocietypublishing.org/ on January 15, 2017

2576

R. Schmidt and P. A. Nikrityuk Table 1. Nomenclature. roman symbols cp g n p t T u m r d

heat capacity (J kg−1 K−1 ) gravitational constant (m s−2 ) surface normal (–) pressure (N m−2 ) time (s) temperature (K) velocity vector (m s−1 ) mass (kg) radius (m) diameter (m)

Greek symbols e l m r U

solid volume fraction (–) heat conductivity (W K−1 m−1 ) dynamic viscosity (kg m−1 s−1 ) density (kg m−3 ) angular velocity (rad s−1 )

subscripts i g p w 0

particle number gas particle wall initial

Table 2. Physical properties of air and coal.

r (kg m−3 ) cp (J kg K−1 ) l (W m K−1 ) m (Pa s) l/rcp (m2 s−1 )

air

coal

ratio coal/air

1.225 1006.43 0.0242 1.8 × 10−5 1.96 × 10−5

2000 1200 0.17

1633 1.2 7

7.1 × 10−8

1/276.06

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. Phil. Trans. R. Soc. A (2011)

Downloaded from http://rsta.royalsocietypublishing.org/ on January 15, 2017

DNS of particulate flows

2577

(a) Fluid flow and heat transfer Taking into account the assumptions made above, the conservation equations for the fluid phase take the form V · u = 0, vu + rV · (uu) = −Vp + mV2 u vt vT rcp + rcp V · (T u) = V · (lVT ). vt r

and

(2.1) (2.2) (2.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: u = Ui × (x − x p ).

(2.4)

It 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 ls . The temperature Tp,i is calculated for each particle as follows: Npart   DTij vT vTp,i = −lg dA − ls Sij , (2.5) r p cp p vt vn Dx j=0 where 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 Dx is the grid cell length. To simulate the contact resistance at the contact area Sij between two particles we used ls 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: mi

dU i = G i + F i + F collis i dt

Phil. Trans. R. Soc. A (2011)

and

mi

dUi = T Mi + T collis Mi . dt

(2.6)

Downloaded from http://rsta.royalsocietypublishing.org/ on January 15, 2017

2578

R. Schmidt and P. A. Nikrityuk

Here, mi is the mass of the particle, U i and Ui 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: 



Fi =

−pn dA + m Si

and

(V × u) × n dA Si

 T Mi =

(r − r i ) × [(V × u) × n]dA,

(2.7)

Si

where 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]: U l = U l −

x J ij ml

and

Ul = Ul −

xrl k × J ij , I1

l = 1, 2.

(2.8)

Here x = 2((m1 m2 )/(m1 + m2 )), and Ii is the mass moment of inertia, in the case of cylinders Ii = 0.5 mi ri2 . 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 J ij = h1 (k · h)k + h2 (k × h) × k,

(2.9)

where k is the collision direction or the direction between the centres of mass, h is the relative velocity at the contact point and h1 and h2 are the dissipation parameters, which are calculated as follows: h1 =

a+1 2

and

h2 =

b+1 . 2(1 + (w2 /k1 ) + (w1 /k2 ))

(2.10)

Here, wi is the mass fraction of each particle (mi /(m1 + m2 )) and k is the coefficient of the moment of inertia I , which is 1/2 for cylindrical particles. The coefficients a and b used in equation (2.10) are the restitution coefficient in the normal direction and the surface roughness, respectively. The coefficient a takes values between zero and unity. If a = 1 or a = 0, the collisions are elastic or inelastic, respectively. The surface roughness can be changed within the limits of −1 ≤ b ≤ 1. An increase in b corresponds to the increase in the surface roughness and, as a result, tangential velocity decreases. If b becomes greater than 0 the direction of the tangential velocity is reversed. For further information see Rao & Nott [4]. In this work, we used a = 0.9 and b = 0. To avoid multiple Phil. Trans. R. Soc. A (2011)

Downloaded from http://rsta.royalsocietypublishing.org/ on January 15, 2017

DNS of particulate flows

2579

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 [A]u x = Q ux ,

(3.1)

where [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 a was zero and b = −0.5. The surrounding medium was air. The comparison between experimental and numerical data is presented in figure 2a, which shows the cylinder positions at different times. Figure 2b,c depicts the comparison of the time histories of the mean centre of mass position in the x- and ydirections. The convective heat transfer was validated against the experiment of Kuehn & Goldstein [6]. Figure 2d 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 × 104 . To study the dynamics of the particle heating, we used the particle volume-averaged and gas volume-averaged temperatures, Phil. Trans. R. Soc. A (2011)

Downloaded from http://rsta.royalsocietypublishing.org/ on January 15, 2017

2580

R. Schmidt and P. A. Nikrityuk

(a)

(b) 0.35

0s

(c) 0.125

0.1 s

0.100

0.25

y/L

x/L

0.30

0.20

0.050

0.15 0.10

0.075 0.025

0.1

0

0.2 0.3 t (s)

0.4

0.5

0

0.1

0.2 0.3 t (s)

0.4

0.5

(d)

0.3 s

370

0.5 s

350 340

0.94

358 348

x (m)

T (K)

360

34

8

346 344

340 338 330 326 328

330

0

0.06

y (m)

0.07

0.094

0.08 x (m)

0.09

Figure 2. Code validations: (a) the comparison of the cylinder positions between the experiment of Zhang et al. [5] and the simulations; (b) and (c) time evolutions of the mean centre of mass positions in the x-direction and y-direction, respectively (open circles with solid lines, experiment; solid lines, present work); (d) calculated and experimentally [6] defined temperature profiles along the vertical symmetry axis of two coaxial cylinders (filled circles, experiment; open squares, present work, grid 50 × 50; dashed line, present work, grid 100 × 100).

0

0.015 T pmean

T gmean

0.020 0.6 0.5 0.4 0.6 0.3 0.4 0.2 0.2 0.1 0

0.010 0.005

1 5

2

3

10 15 revolutions

20

25

0

5

10 15 revolutions

20

25

g

Figure 3. Time histories of the gas volume-averaged temperature, Tmean : (a) calculated using the DNS model with 30 particles for different kiln rotation speeds and the particle volume-averaged p temperature, Tmean (black line, 30 r.p.m.; red line, 15 r.p.m.); (b) calculated using different models for 30 and 100 particles in the kiln with U = 30 r.p.m (2 s per rotation; red solid line: DNS, dp = 14 mm; red dashed line: linear, dp = 14 mm; black solid line: DNS, dp = 28 mm; black dashed line: linear, dp = 28 mm).

which were calculated as follows: x0 y0 ((T − T0 )/(Tw − T0 ))3 dx dy p , Tmean = 0 0 Vpart and g Tmean

x0 y0 =

0

0

Vpart =

((T − T0 )/(Tw − T0 )) (1 − 3) dx dy , Vgas

Phil. Trans. R. Soc. A (2011)

 x0  y0 3 dx dy 0

(4.1)

0

Vgas = pR02 − Vpart . (4.2)

Downloaded from http://rsta.royalsocietypublishing.org/ on January 15, 2017

2581

DNS of particulate flows (a)

0

13.5

4.5

18

(b)

9 revolutions

24.5 revolutions

1.00 0.95 0.90 0.85 0.80 0.75 0.70 0.65 0.60 0.55 0.50 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0

0

4.5

13.5

18

9 revolutions

24.5 revolutions

1.50× 10–02 1.43× 10–02 1.36× 10–02 1.29× 10–02 1.22× 10–02 1.15× 10–02 1.08× 10–02 1.01× 10–02 9.40× 10–03 8.70× 10–03 8.00 × 10–03 7.30 × 10–03 6.60 × 10–03 5.90 × 10–03 5.20 × 10–03 4.50 × 10–03 3.80 × 10–03 3.10 × 10–03 2.40 × 10–03 1.70 × 10–03 1.00 × 10–03

p

Figure 4. Snapshots of the non-dimensional temperature Tmean distribution calculated using the DNS model at different times (rp = 14 mm and U = 30 r.p.m).

12.15

12.18

12.41 revolutions

12.98

13.08

13.11 revolutions

Figure 5. Time series of temperature counterplots for the contact period of a particle with the wall (t = 24.3–26.225 s, rp = 14 mm and U = 30 r.p.m.).

p g The physical meaning of Tmean and Tmean is the mean temperature of all particles g is and of the gas phase in the rotary kiln, respectively. The time evolution of Tmean shown in figure 3a. It can be seen that the gas phase is heated rapidly within one g reaches its maximum, it to two rotations of the kiln. After the temperature Tmean stays approximately constant at a nearly equal level. The temperature oscillations are caused by the vortices induced by the particle movement. The heating history p is depicted in figure 3b. It can be seen that of the particles in the form of Tmean 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 p calculated by the DNS through the contact points. The comparison of Tmean 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

Phil. Trans. R. Soc. A (2011)

Downloaded from http://rsta.royalsocietypublishing.org/ on January 15, 2017

2582

R. Schmidt and P. A. Nikrityuk

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 ls 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 − T0 )/(Tw − T0 )) distribution in the kiln at different times in figure 4a. In order to explore the temperature distribution inside the particles, figure 4b 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. 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.

References 1 Chaudhuri, B., Muzzio, F. & Tomassone, M. 2006 Modeling of heat transfer in granular flow in rotating vessels. Chem. Eng. Sci. 61, 6348–6360. (doi:10.1016/j.ces.2006.05.034) 2 Feng, Z. & Michaelides, E. 2009 Heat transfer in particulate flows with direct numerical simulation (DNS). Int. J. Heat Mass Transfer 52, 777–786. (doi:10.1016/j.ijheatmasstransfer. 2008.07.023) Phil. Trans. R. Soc. A (2011)

Downloaded from http://rsta.royalsocietypublishing.org/ on January 15, 2017

DNS of particulate flows

2583

3 Shi, D., Vargas, W. & McCarthy, J. 2008 Heat transfer in rotary kilns with interstitial gases. Chem. Eng. Sci. 63, 4506–4516. (doi:10.1016/j.ces.2008.06.006) 4 Rao, K. & Nott, P. 2008 An introduction to granular flow. Cambridge, UK: Cambridge University Press. (doi:10.1017/CBO9780511611513) 5 Zhang, S., Kuwabara, S., Suzuki, T., Kawano, Y., Morita, K. & Fukuda, K. 2009 Simulation of solid–fluid mixture flow using moving particle methods. J. Comput. Phys. 228, 2252–2565. (doi:10.1016/j.jcp.2008.12.005) 6 Kuehn, T. & Goldstein, R. 1978 An experimental study of natural convection heat transfer in concentric and eccentric horizontal cylindrical annuli. J. Heat Transfer 100, 635–640. (doi:10.1115/1.3450869) 7 Ferziger, J. & Peri´c, M. 2002 Computational methods for fluid dynamics, 3rd edn. Berlin, Germany: Springer. 8 Zienkiewicz, O. 1971 The finite element method in engineering science. Columbus, OH: McGrawHill. 9 Ananiev, S., Nikrityuk, P. & Meyer, B. 2009 Simulation of the heat and flow dynamics in a one-roll mill using the fictitious boundary method. In Turbulence, heat and mass transfer, vol. 6 (eds K. Hanjalic, Y. Nagano & S. Jakirlic), pp. 1063–1066. New York, NY: Begell House.

Phil. Trans. R. Soc. A (2011)