SpringerOpen Newsletter

Receive periodic news and updates relating to SpringerOpen.

Open Access Highly Accessed Research

Counteracting ring formation in rotary kilns

M Pisaroni1, R Sadi2 and D Lahaye1*

Author Affiliations

1 Scientific Computing Group, Delft Institute of Applied Mathematics, Faculty of Electrical Engineering, Mathematics and Computer Science, Delft University of Technology, Delft, The Netherlands

2 Almatis B.V., Theemsweg 30, Botlek, Rotterdam, 3197 KM, The Netherlands

For all author emails, please log on.

Journal of Mathematics in Industry 2012, 2:3  doi:10.1186/2190-5983-2-3

The electronic version of this article is the complete one and can be found online at: http://www.mathematicsinindustry.com/content/2/1/3

Received:27 April 2012
Accepted:9 October 2012
Published:24 October 2012

© 2012 Pisaroni et al.; licensee Springer

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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 kiln; computational fluid dynamics; ring formation

1 Introduction

1.1 Rotary kilns

Rotary kilns [1] are long cylindrical industrial furnaces that for their operation are slightly tilted as shown in Figure 1. 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.[2-4].

thumbnailFigure 1. General layout of a direct fired, counter-current fed rotary kiln.

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 calcium-aluminate 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 <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M1">View MathML</a> degrees Celsius. This cement is made by fusing a mixture of a 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.

1.2 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 [5]. 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 2 we motivate the development of a computational fluid dynamics model for the kiln. In Section 3 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 4 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 5 finally we report on an experimental validation at the plant that confirms the effectiveness of the change in operating conditions that we propose.

2 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. 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 2 for one of the raw materials to be processed.

thumbnailFigure 2. Phase diagram of one of the materials processed in the kiln.

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 3. 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.

thumbnailFigure 3. Geometry model of the kiln.

3 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.

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 3 and Figure 4. 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.

thumbnailFigure 4. Mesh model of the kiln.

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.

3.1 The geometry

The geometry model of the kiln is shown in Figure 3. Figure 3(a) and Figure 3(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.

3.2 Grid generation

The mesh model is shown in Figure 4. Figure 4(a), Figure 4(b) and Figure 4(c) give an exterior view of the mesh in the air inlet region, an interior view of the mesh and a detailed view of the mesh on the burner, respectively. We employed a polyhedral mesh of 2.8 million cells with local refinement in the critical inlets and burner regions. The main difficulties in meshing this geometry were found in balancing the accuracy required capture 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 10) 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 12 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 converge more robustly with respect to change in their parameters. A more detailed analysis of polyhedral meshes can be found in [6].

3.3 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.[7,8].

Species are defined by their mass fraction defined as

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M2">View MathML</a>


where N is the number of species in the reacting mixture, <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M3">View MathML</a> 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

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M4">View MathML</a>


where <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M5">View MathML</a> is the density of the gas and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M6">View MathML</a> its three dimensional velocity field, respectively. The conservation of species for <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M7">View MathML</a> can then be written as

NS: Conservation of species

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M8">View MathML</a>


where <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M9">View MathML</a> the ith component of the diffusion velocity <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M10">View MathML</a> of species and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M11">View MathML</a> the chemical reaction rate of species . The conservation of momentum for the gas can for <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M12">View MathML</a> be expressed as:

NS: Conservation of momentum

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M13">View MathML</a>


where p denotes pressure and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M14">View MathML</a> and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M15">View MathML</a> 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 <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M16">View MathML</a> as independent variable. To introduce this quantify, we will first denote the enthalpy of species as <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M17">View MathML</a>. This quantify is the sum of the of a sensible and chemical part, i.e.,

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M18">View MathML</a>


where the last term represents the enthalpy of formation of the species at a particular reference temperature <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M19">View MathML</a>. 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

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M20">View MathML</a>


We will denote the heat diffusion coefficient and the temperature by λ and T, respectively. The energy flux <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M21">View MathML</a> 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.,

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M22">View MathML</a>


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 <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M23">View MathML</a> and heat release due to combustion denoted as <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M24">View MathML</a>, where

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M25">View MathML</a>


The work done by the gas on the species can be expressed as <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M26">View MathML</a>. With all these quantities introduced, the conservation of energy in terms of <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M16">View MathML</a> can be expressed as:

NS: Conservation of energy

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M28">View MathML</a>


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 flame-generated 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 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.

3.3.1 RANS model

In constant density flows, Reynolds averaging consist in splitting any quantity ξ in mean and fluctuating component (<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M29">View MathML</a>). In variable density flow Favre [9] mass-weighted averages are usually preferred, i.e., <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M30">View MathML</a>. Any quantity f can therefore be split into:

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M31">View MathML</a>

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

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M32">View MathML</a>


the equation for conservation of species for <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M7">View MathML</a>

RANS: Conservation of species

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M34">View MathML</a>


the equation for conservation of momentum for <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M12">View MathML</a>

RANS: Conservation of momentum

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M36">View MathML</a>


and finally the equation for conservation of momentum

RANS: Conservation of energy

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M37">View MathML</a>


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: <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M38">View MathML</a>.

• Species chemical reaction rates: <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M39">View MathML</a>.

3.3.2 Turbulence model

Using the turbulence viscosity assumption by Boussinesq [10], the Reynolds stresses can be represented as

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M40">View MathML</a>


where <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M41">View MathML</a> is the turbulent dynamic viscosity and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M42">View MathML</a> the Kronecker delta. The turbulent kinetic energy k in turn can be expressed as

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M43">View MathML</a>


Modeling the turbulent viscosity <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M44">View MathML</a> 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[11]. Heat release effects on the Reynolds stresses are not explicitly taken into account in this approach and the turbulent viscosity is modeled as

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M45">View MathML</a>


where ε is the rate of energy dissipation. In this model the critical coefficient <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M46">View MathML</a> 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 (14) and the eddy viscosity definition in Equation (16) it is possible to obtain the following expression for the normal Reynolds stress <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M47">View MathML</a> in an incompressible strained mean flow U

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M48">View MathML</a>


where <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M49">View MathML</a>. It can be shown that <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M50">View MathML</a>, which by definition is a positive quantity, become negative (non-realizable), when the strain is large enough to satisfy

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M51">View MathML</a>


The easiest way to ensure the realizability is to make <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M46">View MathML</a> in (16) variable [12].

The critical coefficient <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M46">View MathML</a> 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

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M54">View MathML</a>





where <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M57">View MathML</a> is the mean rate of rotation tensor viewed in a rotating reference frame with the angular velocity <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M58">View MathML</a>. The parameters <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M59">View MathML</a> and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M60">View MathML</a> in (19) can be computed as





The turbulent kinetic energy k and its dissipation rate ε in Equation (16) are described by the following two balance equations



where <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M67">View MathML</a> is the production term of turbulent kinetic energy due to the mean velocity gradients, <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M68">View MathML</a> the production of turbulent kinetic energy due to buoyancy, <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M69">View MathML</a> the dilatation dissipation term that accounts for the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M70">View MathML</a> and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M71">View MathML</a> user defined source terms for turbulent kinetic energy and dissipation, and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M72">View MathML</a> and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M73">View MathML</a> the turbulent Prandtl numbers for k and ε, respectively. <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M74">View MathML</a>, <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M75">View MathML</a> and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M76">View MathML</a> 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 [11] was developed to repair this deficiency and addresses as such an issue that is of primary importance in our study.

3.3.3 Combustion model

The averaged equation for conservation of species (11) can be rewritten in compact form for <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M7">View MathML</a> as

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M78">View MathML</a>


where <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M79">View MathML</a> is the mass diffusion flux of species . The previous equation is solved in a CFD code for <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M80">View MathML</a> 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 Nth mass fraction is determined as one minus the sum of the <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M80">View MathML</a> solved mass fractions. To minimize numerical error, the Nth species should be selected as that species with the overall largest mass fraction.

In turbulent flows the mass diffusion flux is computed as

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M82','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M82">View MathML</a>


where <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M83','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M83">View MathML</a> is the turbulent Schmidt number and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M84','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M84">View MathML</a> is the molecular diffusivity of species .

The species chemical reaction rate unclosed term <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M11">View MathML</a> must be modeled with a combustion model. A combustion model describes the two-way interaction between properties of the turbulent flow produced by the flame and the chemical reactions. It serves to compute the reaction state space and the quantities it influences, namely the fluid density, viscosity, and temperature. It accounts for the processes that occur at length and time scales that we cannot resolve on a grid either in space or time due to limitations in computational resources. The choice of combustion model is decided by knowing the Damkohler number, defined as <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M86">View MathML</a>, where <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M87">View MathML</a> is the mixing time scale and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M88">View MathML</a> is the reaction time scale. When the Damkohler number is very large, as in the case of the kiln, the reaction rate is controlled by the turbulent mixing that brings reactants together at the molecular scale. In this limit, the Standard Eddy Break Up (EBU) [13] model is fairly accurate because it assumes that the reaction occurs instantaneously upon micromixing.

The EBU combustion model tracks individual mean species concentrations on the grid through transport equations. The reaction rates used in these equations are calculated as functions of the mean species concentrations, turbulence characteristics and, depending on the specific model used, temperature. A mean enthalpy equation is solved in addition to the species transport equations. The mean temperature, density and viscosity are then calculated knowing the mean enthalpy and species concentrations. In the EBU used, the individual species in the global reaction are assumed to be transported at different rates according to their own governing equations.

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 reaction of the form

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M89','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M89">View MathML</a>


where F stands for fuel, O oxidiser and P products of the reaction, the reaction rate for reaction is assumed to be

<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M90','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M90">View MathML</a>


where <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M91">View MathML</a>, <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M92','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M92">View MathML</a>, v is the molar stoichiometric coefficient for species j in reaction , M is molecular weight of species. Equation (30) 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 (<a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M93">View MathML</a>). <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M94','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M94">View MathML</a>, <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M95','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M95">View MathML</a>, <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M96">View MathML</a> are respectively the mean concentrations of fuel, oxidizer, and products. <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M97','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M97">View MathML</a> and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M98','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M98">View MathML</a> are model constants with typical values of 0.5 and 4.0 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 6 species and 4 reactions to account for a fuel that is a mixture of different alkanes. This mixture consists for 95% of CH4 and for 5% of C2H6, C3H8 and C4H10.

The above models are discretized by a finite volume technique using second order upwinding for the convective terms [14-16]. 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 [17].

3.4 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 [18,19].

3.5 Software implementation

Simulations were performed using the STAR-CCM+ software suite [20] on a ten-nodes Linux cluster having Intel Duo and Quad Core processors at a clock speed between 2.20 GHz and 3.16 GHz running a Slackware 13 64-bit distribution. Iterating the three-dimensional combustion model to equilibrium state required between <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M99','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M99">View MathML</a> and <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M100','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M100">View MathML</a> non-linear iterations and between three and three and a half days of computation time.

4 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 1 and differ in the amount of preheated air that flows through the secondary air inlet shown in Figure 3. In the first and second line in Table 1 the ratio of volume air to volume of fuel 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 5. This figure shows that with in the new operating conditions the oxigen is transported further down in the kiln.

thumbnailFigure 5. Effect of change in operating conditions on the streamlines of methane and oxigen in the kiln.

Table 1. Standard and new configuration of operating conditions of 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 the temperature and radiative heat distribution on the inside wall of the kiln. These will be elaborated separately in the two subsections.

4.1 Effect on temperature

In Figure 6, Figure 7 and Figure 8 we show how the change in operating conditions affects the computed temperature profile. Figure 6 shows the temperature distribution on an axial section through the air inlet. Figure 7 shows the temperature evaluated at the centroid of each cell of the surface mesh on the inside wall versus the centroid’s axial distance. 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 7 are therefore thicker for small values of the abscis. Figure 8 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 0.01. The vertical and horizontal color bars indicate the temperature of the flame and the wall, respectively.

thumbnailFigure 6. Effect of change in operating conditions on the temperature on an axial section of the kiln through the air inlet.

thumbnailFigure 7. Temperaturevs.position (blue standard operating conditions, red new operating conditions).

thumbnailFigure 8. Effect of change in operating conditions on the temperature on the inside wall of the kiln.

Figure 6 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 7 shows that this coolant reduces the wall peak temperature drops by 3.5% from <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M101','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M101">View MathML</a> degrees Celcius to <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M102','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M102">View MathML</a> 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 8 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.

4.2 Effect on radiative heat transfer

In Figure 9, Figure 10 and Figure 11 and we show how the change in operating conditions affects the computed radiative heat transfer on the inside wall of the kiln. Figure 9 and Figure 10 are the equivalent for radiative heat transfer of Figure 7 and Figure 8 for the temperature, respectively. Figure 11 finally gives an view on the radiative heat transfer distribution from the inside of the kiln.

thumbnailFigure 9. Outside view on the effect of change in operating conditions on radiative heat transfer on the inside wall.

thumbnailFigure 10. Incident radiationvs.position (blue standard operating conditions, red new conditions).

thumbnailFigure 11. Inside view on the effect of change in operating conditions on radiative heat transfer on the inside wall.

Figure 9 confirms that as expected changing to the new operating conditions results in a lowering of the peak in radiative heat transfer. Figure 10 shows that peak in the radiative heat transfer drops by 13.5% from <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M103','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M103">View MathML</a> W/m2 to <a onClick="popup('http://www.mathematicsinindustry.com/content/2/1/3/mathml/M104','MathML',630,470);return false;" target="_blank" href="http://www.mathematicsinindustry.com/content/2/1/3/mathml/M104">View MathML</a> W/m2 and that the region of high values is reduced in size. Figure 11 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.

5 Experimental validation at the plant

The change in operating conditions we propose was experimentally validated at the plant for the first time on August 28th, 2011. 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 12. Four hours after the decision was taken, it was observed that the ring stopped growing in size as shown in Figure 12(a). After twenty-four hours, the ring started to break into lumps due to the vibration of the drive gears as shown in Figure 12(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 12(c). Careful analysis has revealed that with the new operating conditions, the final cement product still meets the stringent quality specifications. After August 28th, 2011 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.

thumbnailFigure 12. Experimentally observed time evolution of the ring after the new operating conditions have been imposed.

6 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 reduction 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 28th, 2011, the kiln has been operation without unscheduled shut-downs, resulting in hugely important cost savings.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

MP developed and validated the CFD model. RS provided supervision from within the industrial partner. DL provided supervision from within the academic partner.


We thank Marco Talice for his advice worth more than platina and the CD-ADAPCO London office for the support in using their software.


  1. Boateng AA: Rotary Kilns: Transport Phenomena and Transport Processes. Butterworth-Heinemann, Stoneham; 2006. OpenURL

  2. Mastorakos E, Massias A, Tsakiroglou CD, Goussis DA, Burganos VN, Payatakes AC: CFD predictions for cement kilns including flame modelling, heat transfer and clinker chemistry.

    Appl Math Model 1999, 23(1):55-76. Publisher Full Text OpenURL

  3. Marias F, Roustan H, Pichat A: Modelling of a rotary kiln for the pyrolysis of aluminium waste.

    Chem Eng Sci 2005, 60(16):4609-4622. Publisher Full Text OpenURL

  4. Giannopoulos D, Kolaitis DI, Togkalidou A, Skevis G, Founti MA: Quantification of emissions from the co-incineration of cutting oil emulsions in cement plants: part i: NOx, CO and VOC.

    Fuel 2007, 86(7-8):1144-1152. Publisher Full Text OpenURL

  5. Tran HN, Barham D: An overview of ring formation in lime kilns.

    Tappi J 1991, 74(1):131-136. OpenURL

  6. Peric M: Flow simulation using control volumes of arbitrary polyhedral shape. ERCOFTAC Bulletin September 2004, 62.

  7. Williams FA: Combustion Theory. 2nd edition. Westview Press, New York; 1994. OpenURL

  8. Kuo KK: Principles of Combustion. 2nd edition. Wiley-Interscience, New York; 2005. OpenURL

  9. Poinsot T, Veynante D: Theoretical and Numerical Combustion. 2nd edition. Edwards, Ann Arbor; 2005. OpenURL

  10. Hinze JO: Turbulence. McGraw-Hill, New York; 1975. OpenURL

  11. Shih TH, Liou WW, Shabbir A, Zhu J: A new k-ϵ Eddy-viscosity model for high Reynolds number turbulent flows - model development and validation.

    Comput Fluids 1995, 24(3):227-238. Publisher Full Text OpenURL

  12. Reynolds WC:

    Fundamentals of Turbulence for Turbulence Modeling and Simulation. 1987. [Lecture Notes for Von Karman Institute]

    Agard Report No. 755.


  13. Magnussen BF, Mjertager BH: On mathematical modeling of turbulent combustion. In 16th Symp. (Int.) on Combustion. The Combustion Institute, Pittsburgh; 1976:719-727. OpenURL

  14. Demirdzic I, Peric M: Finite volume method for prediction of fluid flow in arbitrarily shaped domains with moving boundaries.

    Int J Numer Methods Fluids 1990, 10(7):771-790. Publisher Full Text OpenURL

  15. Mathur SR, Murthy JY: A pressure-based method for ustructured meshes.

    Numer Heat Transf, Part B, Fundam 1997, 31(2):195-215. Publisher Full Text OpenURL

  16. Murthy JY, Mathur SR: Radiative heat transfer in axisymmetric geometries using an unstructured finite-volume method.

    Numer Heat Transf, Part B, Fundam 1998, 33(4):397-416. Publisher Full Text OpenURL

  17. Trottenberg U, Oosterlee C, Schüller A: Multigrid. Academic Press, San Diego; 2001. OpenURL

  18. Modest MF: Radiative Heat Transfer. McGraw-Hill, New York; 1993. OpenURL

  19. Siegel R, Howell JR: Thermal Radiation Heat Transfer. 3rd edition. Hemisphere, Washington; 1992. OpenURL

  20. Star-CCM+Version 6.06.011 user guide; 2011.