Counteracting ring formation in rotary kilns

Avoiding the formation of rings in rotary kilns is an issue of primary concern to the cement production industry. We developed a numerical combustion model that revealed that in our case study rings are typically formed in zones of maximal radiative heat transfer. This local overheating causes an overproduction of the liquid phase of the granular material, which tends to stick to the oven’s wall and to form rings. To counteract for this phenomenon, we propose to increase the amount of secondary air injected to cool the oven. Experimental validation at the plant has confirmed that our solution is indeed effective. For the first time in years, the kiln has been operating without unscheduled shut-downs, resulting in hugely important cost savings.


Rotary kilns
Rotary kilns [] are long cylindrical industrial furnaces that for their operation are slightly tilted as shown in Figure . They are used in a wide range of material processing industries. The material to be processed is fed into the upper end of the cylinder, leaving a considerable amount of freeboard or empty space. As the kiln rotates about its axis, the material bed gradually moves down towards the lower end, and undergoes a certain amount of stirring and mixing. The material leaves the oven at the lower end to be further processed. The heating of the kiln serves to drive the specific bed reactions, which, for either kinetic and thermodynamic reasons, require high temperatures. The CFD simulation of these devices has been considered in e.g. [-].
We will look into direct fired kilns in which the energy necessary to heat the material to the level required for the intended reactions is generated by the combustion of hydrocarbon fuels continuously fed to a burner placed in the freeboard. This energy is subsequently transferred by heat exchange between the gas phase and the material bed. The lateral surface of the kiln is covered by an isolating lining of refractory material. Hot gases flow in a direction opposite to that of the material bed through the kiln. The heat transfer between the freeboard and the bed is a rather complex phenomenon as it occurs along various paths determined by the physics of radiation. Various transport processes therefore constitute the working principle of the ovens considered.
We study a rotary kiln used by Almatis B.V. in Rotterdam for a production of calciumaluminate cement, a white, high purity hydraulic bonding agent providing controlled setting times and strength required in today's high performance refractories operating at temperatures up to , degrees Celsius. This cement is made by fusing a mixture of a http://www.mathematicsinindustry.com/content/2/1/3 calcium-bearing limestone and an aluminum-bearing material. The following two considerations will be important in our approach to the modeling of the kiln. The first is that the material bed occupies no more than five percent of the total volume of the kiln. The second is that the combustion of fuel by the burner is aided by inflow of preheated air through a secondary air inlet.

Ring formation
As material slides and tumbles slowly through a heated rotary kiln, a thin layer of dust invariably forms on its inner lateral surface. Some zones of the kiln may be particularly prone to particle accumulation in such a way that the combination of particular thermal and flow conditions results in the formation of cylindrical deposits of solid crystalline structures referred to as rings. As such rings grow thicker, they form a dam in the freeboard, hindering the flow of material and flue gasses through the kiln.
In the kiln we consider, we observed in particular so-called front-end/mid-kiln rings []. They are formed close to the burner and are presumably caused by the direct impingement of the burner flame on the lateral surface causing local overheating. These are the most common and most troublesome type of rings as they cannot be reached from the kiln's exterior and are therefore impossible to remove while the kiln is in operation. In severe cases, ring dams grow rapidly and cause the unscheduled shutdown of the kiln. In the last three years we registered on average at least one ring formation per month, of which seventy percent caused the unforeseen shutdown of the kiln for on average three days. Each of these kiln outages causes very important production and turnover losses. This paper is structured as follows. In Section  we motivate the development of a computational fluid dynamics model for the kiln. In Section  we give details of this model's operational setup and geometry, the polyhedral mesh discretizing the spatial coordinates and the partial differential equations governing the gas flow in the oven, the combustion of hydrocarbon fuel in the burner, the reactions of the chemical species and the radiative heat transfer in the freeboard. In Section  we present numerical results that convincingly show that our model allows to localize and counteract the ring formation by increasing the amount of air through the secondary air inlet. In Section  finally we report on an experimental validation at the plant that confirms the effectiveness of the change in operating conditions that we propose.

Localizing and counteracting ring formation
Experience at the plant has shown that neither trial and error nor simplified mathematical modeling approaches provide sufficient insight to eliminate the problem of ring formation. http://www.mathematicsinindustry.com/content/2/1/3

Figure 2 Phase diagram of one of the materials processed in the kiln.
Measuring the temperature inside the kiln, using thermocouples for instance, has proven to be difficult due to the harsh operating conditions. We therefore decided to develop a sufficiently detailed mathematical model enabling to link the ring formation to the temperature, radiative heat and gas flow profiles in the kiln. This model will be described in more details in the next section and has revealed that rings are typically formed in zones of maximal radiative heat transfer. In a phase-diagram this local overheating can be linked to the overproduction of the amount of liquid phase of the material. Such a phase diagram is shown in Figure  for one of the raw materials to be processed.
The insight linking the ring formation to local overheating proved to be instrumental in altering the operating conditions of the kiln in such a way to avoid ring formation. Our solution avoids drastic changes in the geometry of the kiln, in the configuration of the burner and in the properties of the fuel or the lining. We instead alter the fuel/air composition by altering the amount of preheated air injected through the secondary air inlet shown in Figure . This allows us to eliminate peaks in temperature and radiative heat transfer profiles, to limit the amount of liquid phase production and therefore avoid ring formation. The change in operating conditions we propose does maintain the temperature distribution required by the stringent quality specifications on the end product.

Computational fluid dynamics model
The detailed mathematical model we developed is a multi-physics model that takes the following phenomena into account: the reactive gas flow and temperature, chemical species and radiative heat transfer distribution in the kiln, the turbulent non-premixed combustion of hydro-carbon gasses in the burner, the insulating properties of the lining, the rotary motion of the kiln and the forced convection on the outside surface.
The material bed only occupies a small fraction of the volume of the kiln and has a negligible limited impact of temperature distribution. We therefore do not take the material bed into account in our model and simulate an empty kiln. http://www.mathematicsinindustry.com/content/2/1/3

Figure 3 Geometry model of the kiln.
In the kiln hot gases are generated by a flame projected from a burner-pipe placed inside the freeboard. The setup of the burner is shown in Figure  and Figure . The burner injects fuel in axial and radial direction and is cooled by an amount of air forced through the cooling slots. The flow of cooling air through the rectangular air inlet breaks the circular symmetry of the configuration and requires a model in three spatial dimensions to be resolved.
The most important physical phenomenon that takes place in this burner region is the turbulent non-premixed combustion of the fuel injected from the burner with the secondary air. Combustion, even without turbulence, is an intrinsically complex process involving a large range of chemical time and length scales. Some of the chemical phenomena controlling flames take place in short times over thin layers and are associated with very large mass fractions, temperature and density gradients. The full description of chemical mechanisms in laminar flames may require hundreds of species and thousands of reactions leading to considerable numerical difficulties. Turbulence itself is probably the most complex phenomenon in non-reacting fluid mechanics. Various time and length scales are involved and the description of turbulence remains to date an open questions. The modeling of the kiln therefore requires resorting to a set of assumptions that are described in the remainder of this section.

The geometry
The geometry model of the kiln is shown in Figure . Figure (a) and Figure (b) give an exterior view on the complete kiln and a more detailed interior view of the burner region, respectively. The challenging aspect of this geometry is that the inlets of the burner are a factor thousand times smaller than the axial length of the kiln, imposing challenges in the mesh generation process.

Grid generation
The mesh model is shown in Figure . Figure  the flow around small features in the burner with the overall computational cost. Polyhedral meshes provide a balanced solution in complex mesh generation problems of this kind.
Tetrahedra are the simplest type of volume elements. As their faces are plane segments, both face and volume centroid locations are well defined. A disadvantage is that tetrahedra cannot be stretched too much. To achieve a reasonable accuracy a much larger number of control volumes is needed than if structured meshes are used. Furthermore, as tetrahedral control volumes have only four neighbors, and computing gradients at cell centers using standard approximations can be problematic.
Polyhedra offer the same automatic meshing benefits as tetrahedra while overcoming their disadvantages. A major advantage of polyhedral cells is that they have many neighbors (typically of order ) allowing gradients to be much better approximated. Obviously more neighbors implies more storage and computational operations per cell, but this is more than compensated by a higher accuracy. Polyhedral cells are also less sensitive to stretching than tetrahedra. A polyhedron with  faces for instance has six optimal directions which, together with the larger number of neighbors, leads to a more accurate solution with a lower cell count. Comparisons in many practical tests have verified that with polyhedral meshes, one needs about four times fewer cells, half the amount of memory and a tenth to a fifth of computing time compared to tetrahedral meshes to reach solutions of the same accuracy. In addition, solvers on polyhedral meshes were found to http://www.mathematicsinindustry.com/content/2/1/3 converge more robustly with respect to change in their parameters. A more detailed analysis of polyhedral meshes can be found in [].

Governing reacting flow equations
In this section we present the conservation equations for reacting flows we used. The equations are derived from the Navier-Stokes (NS) equations by adding terms that account for reacting flows. The reacting gas is a non-isothermal mixture of multiple species which must be tracked individually. As heat capacities change significantly with temperature and composition, the transport coefficients require specific attention. In this subsection we will describe the Navier-Stokes and Reynolds-Averaged Navier-Stokes flow model, the non-realizable K-Epsilon turbulence model and the Standard Eddy Break Up combustion model. A more detailed derivation of these equations can be found in e.g.
Species are defined by their mass fraction defined as where N is the number of species in the reacting mixture, m the mass of species in a volume V and m the mass of gas in the volume, respectively. The conservation of mass can then be written as NS: Conservation of mass where ρ = m/V is the density of the gas and u i its three dimensional velocity field, respectively. The conservation of species for =  : N can then be written as NS: Conservation of species where V ,i the ith component of the diffusion velocity V of species andω the chemical reaction rate of species . The conservation of momentum for the gas can for j =  :  be expressed as: where p denotes pressure and τ ij and f ,j are the components of the Reynolds stress tensor and of the volume force acting on species , respectively. We will express the conservation of energy using the sensible enthalpy of the mixture h s as independent variable. To introduce this quantify, we will first denote the enthalpy of species as h . This quantify is the sum of the of a sensible and chemical part, i.e., where the last term represents the enthalpy of formation of the species at a particular reference temperature T  . The enthalpy of the mixture is then defined as a weighted average that again can be decomposed into a sensible and chemical part, as follows We will denote the heat diffusion coefficient and the temperature by λ and T, respectively. The energy flux q i is the sum a heat diffusion term derived from Fourier's Law and a term associated with the diffusion of species with different enthalpies, i.e., In the energy conservation equation the following three source terms will play a role: the heat source due to radiative heat flux denoted as Q, the viscous heating term denoted as , where = τ ij ∂u i ∂x j and heat release due to combustion denoted asω T , wherė The work done by the gas on the species can be expressed as ρ N = Y f ,i V ,i . With all these quantities introduced, the conservation of energy in terms of h s can be expressed as: NS: Conservation of energy Turbulent combustion results from the two-way interaction between chemistry and turbulence. When a flame interacts with a turbulent flow, the combustion modifies the turbulence in two ways. The heat released induces high flow accelerations through the flame front and the temperature changes generate large changes in kinematic viscosity. These phenomena may either generate or damp turbulence and are referred to as flamegenerated turbulence and relaminarization due to combustion, respectively. The turbulence conversely modifies the flame structure. This may either enhance the chemical reactions or completely inhibit it, leading to flame quenching. Compared to premixed flames, turbulent non-premixed flames exhibit some specific features that have to be taken into account. Non-premixed flames do not propagate as they localized on the fuel-oxidizer interface. This property is useful for safety purposes but it also has consequences on the chemistry-turbulence interaction. Without propagation speed, a non-premixed flame is http://www.mathematicsinindustry.com/content/2/1/3 unable to impose its own dynamics on the flow field and is therefore more sensitive to turbulence.
The description of the turbulent non-premixed combustion processes in a computational fluid dynamics model may be achieved using three levels of accuracy in the computations. Either a Reynolds Averaged Navier Stokes (RANS), a Large Eddy Simulations (LES) or a Direct Numerical Simulations (DNS) model can be adopted. In current engineering practice, the RANS model is extensively used because it is less demanding in terms of resources. Its validity however is limited by the closure models describing turbulence and combustion and the need for some form of callibration. Considering the complexities and the dimensions of our kiln, using the RANS model is the only feasible choice.

.. RANS model
In constant density flows, Reynolds averaging consist in splitting any quantity ξ in mean and fluctuating component (ξ = ξ + ξ ). In variable density flow Favre [] mass-weighted averages are usually preferred, i.e., f = ρf ρ . Any quantity f can therefore be split into: The RANS equations derived from the reacting Navier-Stokes equation given above are then given by the equation for conservation of mass RANS: Conservation of mass the equation for conservation of species for =  : N RANS: Conservation of species the equation for conservation of momentum for j =  :  RANS: Conservation of momentum and finally the equation for conservation of momentum RANS: Conservation of energy The averaging procedure introduces unclosed quantities that have to be modeled. Without entering in the details we list here the two main unclosed terms that will be described in the next sections: • Reynolds stresses: ρ u i u j .

.. Turbulence model
Using the turbulence viscosity assumption by Boussinesq [], the Reynolds stresses can be represented as where μ t = ρν t is the turbulent dynamic viscosity and δ ij the Kronecker delta. The turbulent kinetic energy k in turn can be expressed as Modeling the turbulent viscosity μ t is the central problem in turbulence computations. Many approaches exist. In this work we use a classical turbulence model developed for non-reacting flows, namely the Realizable K-Epsilon model []. Heat release effects on the Reynolds stresses are not explicitly taken into account in this approach and the turbulent viscosity is modeled as where ε is the rate of energy dissipation. In this model the critical coefficient C μ is a function of mean flow and turbulence properties, rather than assumed to be constant as in the standard model. This allows to satisfy certain mathematical constraints on the normal stresses consistent with the physics of turbulence and is referred to as realizability.

From the Boussinesq relationship in Equation () and the eddy viscosity definition in Equation () it is possible to obtain the following expression for the normal Reynolds stress u  in an incompressible strained mean flow
where ν t = μ t ρ . It can be shown that u  , which by definition is a positive quantity, become negative (non-realizable), when the strain is large enough to satisfy The easiest way to ensure the realizability is to make C μ in () variable []. http://www.mathematicsinindustry.com/content/2/1/3 The critical coefficient C μ can be expressed as a function of the mean strain and rotation rates, the angular velocity of the system rotation and the turbulence fields as follows where ij is the mean rate of rotation tensor viewed in a rotating reference frame with the angular velocity ω k . The parameters A  and A s in () can be computed as The turbulent kinetic energy k and its dissipation rate ε in Equation () are described by the following two balance equations where P k is the production term of turbulent kinetic energy due to the mean velocity gradients, P b the production of turbulent kinetic energy due to buoyancy, Y M the dilatation dissipation term that accounts for the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, S k and S ε user defined source terms for turbulent kinetic energy and dissipation, and σ k and σ ε the turbulent Prandtl numbers for k and ε, respectively. C ε , C ε and C ε are model constants. Another weakness of traditional K-Epsilon turbulence models is their modeling of the dissipation rate ε. Indeed, the well-known spreading (or dispersion) rate anomaly refers to the fact that traditional models do reasonably well in predicting the spreading rate of a planar jet but perform unexpectedly poor for rounds jets. This weakness can be traced back to a deficiency in traditional ε-equations. The realizable model proposed by Shih [] was developed to repair this deficiency and addresses as such an issue that is of primary importance in our study. http://www.mathematicsinindustry.com/content/2/1/3

.. Combustion model
The averaged equation for conservation of species () can be rewritten in compact form for =  : N as where J is the mass diffusion flux of species . The previous equation is solved in a CFD code for N - species where N is the total number of fluid phase chemical species present in the system. Since the mass fraction of the species must sum to unity, the N th mass fraction is determined as one minus the sum of the N - solved mass fractions. To minimize numerical error, the N th species should be selected as that species with the overall largest mass fraction.
In turbulent flows the mass diffusion flux is computed as where Sc t is the turbulent Schmidt number and D is the molecular diffusivity of species .
The species chemical reaction rate unclosed termω must be modeled with a combus- The reaction rate is modeled through an expression that takes the turbulent micromixing process into account. This is done through dimensional arguments. Thus, for a reac- (   ) http://www.mathematicsinindustry.com/content/2/1/3 where F stands for fuel, O oxidiser and P products of the reaction, the reaction rate for reaction is assumed to bė , v is the molar stoichiometric coefficient for species j in reaction , M is molecular weight of species. Equation () essentially states that the integrated micromixing rate is proportional to the mean (macroscopic) concentration of the limiting reactant divided by the time scale of the large eddies ( k ε = τ mix ). Y F , Y O , Y P are respectively the mean concentrations of fuel, oxidizer, and products. A ebu and B ebu are model constants with typical values of . and . respectively. The values of these constants are fitted according to experimental results and they are suitable for most cases of general interest.
In our simulations we used a reduced combustion mechanism with  species and  reactions to account for a fuel that is a mixture of different alkanes. This mixture consists for % of CH  and for % of C  H  , C  H  and C  H  .
The above models are discretized by a finite volume technique using second order upwinding for the convective terms [-]. The flow equations are solved in a segregated approach in which the SIMPLE algorithm realizes the velocity-pressure coupling. The energy equation is solved for the chemical thermal enthalpy using again a segregated approach. The temperature is computed from the enthalpy according to the equation of state. At each outer non-linear iteration the resulting linear systems are solved using an algebraic multigrid preconditioner for a suitable Krylov subspace acceleration [].

Additional information
In situations in which the media separating hot walls is transparent for thermal radiation as in the case of dry air, radiation can only occur as a surface phenomenon. In the our case however, the gas in the freeboard of the kiln will absorb, emit and scatter the thermal radiation intensity emitted from the hot walls of the kiln. This process is governed by the radiative transfer equation (RTE) that is implemented in a Participating Media Radiation Model. The model is discretized in solid angle by the discrete ordinate method described in detail in [, ].

Software implementation
Simulations were performed using the STAR-CCM+ software suite [] on a ten-nodes Linux cluster having Intel Duo and Quad Core processors at a clock speed between . GHz and . GHz running a Slackware  -bit distribution. Iterating the threedimensional combustion model to equilibrium state required between , and , non-linear iterations and between three and three and a half days of computation time.

Computational results
In this section we discuss how our combustion model validates the difference between two operating conditions of the oven. These operating conditions are listed in Table  and differ in the amount of preheated air that flows through the secondary air inlet shown in Figure . In the first and second line in Table  the ratio of volume air to volume of fuel http://www.mathematicsinindustry.com/content/2/1/3  equals nine and twelve, respectively. The ratio of nine corresponds to the kiln's operation conditions that were commonly used before the start of this work and will be referred to as the standard operating conditions. The ratio of twelve correspond to the new operation conditions that we propose. In switching from the standard to the new configuration, the flux of fuel injected in axial and radial direction and the flux of air through the cooling slot is kept the same. The effect of this switch in operation conditions on the mass of oxygen and methane is shown in Figure . This figure shows that with in the new operating conditions the oxigen is transported further down in the kiln. The objective of this section is to show how according to the combustion model described in the previous section the new configuration is less prone to the development of rings than the standard configuration. The two key arguments in this demonstration are http://www.mathematicsinindustry.com/content/2/1/3  the temperature and radiative heat distribution on the inside wall of the kiln. These will be elaborated separately in the two subsections.

Effect on temperature
In Figure Figure  shows the temperature distribution on an axial section through the air inlet. Figure  shows the temperature evaluated at the centroid of each cell of the surface mesh on the inside wall versus the centroid's axial distance. http://www.mathematicsinindustry.com/content/2/1/3

Figure 8 Effect of change in operating conditions on the temperature on the inside wall of the kiln.
This wall forms the interface between the hot gas and the innermost part of the insulating lining. For positions closer to low end of the kiln the temperature values shown vary more with the circumferential angle due to the higher temperature gradients caused by the injected air. The lines in Figure  are therefore thicker for small values of the abscis. Figure  finally shows the temperature distribution on the inside wall of the kiln. In this figure the flame propagating from the burner is represented by an iso-surface that encloses a region in which the methane concentration is larger than or equal to .. The vertical and horizontal color bars indicate the temperature of the flame and the wall, respectively. Figure  clearly illustrates the effect of switching to the new operating conditions as it shows in the new conditions, the relatively cold preheated air is transported further into the kiln acting as a coolant on the wall. Figure  shows that this coolant reduces the wall peak temperature drops by .% from , degrees Celcius to , degrees Celcius. It also shows that switching to the new operating conditions does not significantly alter the global distribution characteristics of the temperature which is of paramount importance in the material production process. Figure  is particularly interesting as it shows that in both configurations the peak in temperature is situated in a zone at four and halve to seven meters from the burner. This zone coincides with the zone in which rings are typically formed. It is therefore plausible that a reduction in peak temperature by switching to the new operating conditions will result in a reduction of the amount of liquid phase of the material being formed, and therefore impede the formation of rings. This hypothesis will be confirmed by looking into the computed radiative heat transfer distribution as we will do in the next subsection.

Effect on radiative heat transfer
In Figure Figure  confirms that as expected changing to the new operating conditions results in a lowering of the peak in radiative heat transfer. Figure  shows that peak in the radiative heat transfer drops by .% from ,, W/m  to ,, W/m  and that the region of high values is reduced in size. Figure  gives another representation of this fact. By switching to the new operating conditions, the material in the kiln is less likely to absorb an excessive amount of heat, effectively limiting the amount of liquid phase of the material and therefore the formation of rings. http://www.mathematicsinindustry.com/content/2/1/3 Figure 11 Inside view on the effect of change in operating conditions on radiative heat transfer on the inside wall.

Experimental validation at the plant
The change in operating conditions we propose was experimentally validated at the plant for the first time on August th, . That day the formation of a massive ring was reported. After assessing the situation, it was decided to change to the new operating conditions. The effect of this decision over time is shown in Figure . Four hours after the decision was taken, it was observed that the ring stopped growing in size as shown in Figure (a). After twenty-four hours, the ring started to break into lumps due to the vibration of the drive gears as shown in Figure (b). These lumps are transported through the kiln and leave it at the lower end by the rotary motion. After forty hours, the kiln is clean and in stable operating conditions as shown in Figure (c). Careful analysis has revealed that with the new operating conditions, the final cement product still meets the stringent quality specifications. After August th,  the benefits of the changing to the new operating conditions has been repeatedly observed when maintenance of other circumstances required switching back to the standard operating conditions.

Conclusions
We developed a numerical model allowing to access the effectiveness of measures implemented to counteract the formation of rings in a rotary cement kiln in use by Almatis B.V. in Rotterdam. In this three-dimensional combustion model, the gas flow, the temperature profile, radiative heat distribution and the concentration of hydro-carbon species in the kiln is taken into account. Simulations show that deluting the air-fuel mixture with air reduces peaks in radiative heat transfer in zones critical to ring formation. This reduc-http://www.mathematicsinindustry.com/content/2/1/3 tion results in turn in less heat being absorbed by the granular material bed, effectively reducing the amount of material liquid phase prone to sticking to the kiln's surface and to forming rings. The validity of our model has been experimental observed at the Almatis plant in Rotterdam. Since August th, , the kiln has been operation without unscheduled shut-downs, resulting in hugely important cost savings. http://www.mathematicsinindustry.com/content/2/1/3