Numerical analysis of the thermal profile inside the wall of a rotary cement kiln

Análise numérica do perfil térmico na parede de um forno rotativo de cimento

Abstract

The present study aimed to evaluate the temperature profile along the inside of the wall of a clinker kiln from a cement industry. The problem was modeled by the equation of transient heat conduction in cylindrical coordinates, considering radial symmetry. Being the wall composed of different materials, even adopting constant physical properties, there is no analytical solution to the problem. The method of the lines was used, being the radial and axial directions discretized by finite differences and the resulting system of ordinary differential equations integrated in time until obtaining the temperature field in steady state. The obtained field was compatible with heat transfer fundamentals and presented a good fit in relation to industrial data. The main limitations of the modeling performed in this study include the fact that the gases and solids contained in the kiln have not been modeled, and the variation in thicknesses of the layers of the kiln wall has not been considered. The program developed in this study can be used to evaluate the performance of different refractories or to infer the refractory wear level from experimental kiln surface temperature profiles.

Keywords:
cement kiln; temperature profile; numerical analysis; method of lines

Resumo

O presente trabalho buscou avaliar o perfil de temperatura ao longo do interior da parede de um forno de clínquer de uma indústria de cimento. O problema foi modelado pela equação da condução transiente de calor em coordenadas cilíndricas, considerando-se simetria radial. Sendo a parede composta de diferentes materiais, mesmo adotando-se propriedades físicas constantes, não há solução analítica para o problema. Empregou-se o método das linhas, sendo as direções radial e axial discretizadas por diferenças finitas e o sistema de equações diferenciais ordinárias resultante integrado no tempo até a obtenção do campo de temperaturas em regime permanente. O campo obtido foi compatível com fundamentos de transferência de calor e apresentou bom ajuste em relação a dados industriais. As principais limitações da modelagem realizada incluem o fato de os gases e sólidos contidos no forno não terem sido modelados e a variação nas espessuras das camadas da parede do forno não ter sido considerada. O programa desenvolvido pode ser empregado na avaliação do desempenho de diferentes refratários ou na inferência do nível de desgaste dos refratários a partir de perfis experimentais de temperatura na superfície do forno.

Palavras-chave:
forno de cimento; perfil de temperatura; análise numérica; método das linhas

INTRODUCTION

Cement production is continually expanding, and one of the industry’s key challenges is dealing with high energy consumption and reducing the environmental impacts caused by the process ¹. Therefore, there has been a constant concern to understand, model, and optimize the cement production process in order to reduce its environmental impact ^2)-(5. Portland cement is the most widely used type of cement in the world and can be applied to the construction of houses, buildings, bridges, and other structures ⁶. In the process of producing this cement, one of the steps that most interferes with the final quality of the product is clinkerization, a process in which clinker is fired in a rotary kiln, undergoing various physical (phase changes) and chemical processes (endothermic and exothermic reactions) ^7),(8. Chemical reactions from this process vary as the temperature varies throughout the kiln ^9)-(11. Thus, it is of great importance for cement production that the temperature is within the proper range throughout the equipment, so as to allow the necessary reactions to occur to obtain a good quality clinker ¹².

In order to maintain a suitable temperature range inside the kiln and to protect the outside of the equipment, clinker kilns are lined with refractory material, usually magnesia-spinel bricks ¹³. These refractories play an important role in ensuring the energy efficiency of the kiln, thanks to their thermal insulation properties ^14),(15. Studying the thermal profile of the kiln wall can be useful to evaluate the performance of different refractories used. And since the temperature profile inside the kiln cannot always be measured directly due to technical limitations, numerical modeling is a good approach for understanding this profile. In addition, a wall temperature profile model would be useful to infer the refractory lining wear level by identifying deviations between the model’s predicted external surface temperature and the measured surface temperature. Several published studies ^16)-(19 have successfully used numerical methods to model the thermal behavior of various ceramic kilns. In this sense, the present study aimed to perform numerical modeling of the temperature profile inside the wall of a steady-state clinker kiln. In general, previous studies on clinker kilns proposed to predict the temperature on the shell by modeling the transformation processes inside the equipment and the energy transfer involved ²⁰. In this article, on the other hand, we propose an approach focused on the heat transfer phenomena through the wall, relating the temperature on the kiln shell surface to the temperature profile inside the equipment. Choosing this particular approach greatly reduces the computational cost of the calculations since the gases and solids inside the kiln are not comprised within the model. This choice also makes it possible for the algorithm to be coded in any programming language, without relying on any specific commercial simulation software.

MATERIALS AND METHODS

Description of the system: the equipment under study consisted of a clinker kiln of a cement industry located in Minas Gerais, Brazil, whose shape can be approximated by a 60 m long cylindrical tube with a 2 m radius. The kiln wall consisted of 3 layers: a layer of steel casing; a layer of refractory bricks for thermal insulation; and the coating (a layer of material adhered to the inner surface of the kiln) ²¹. The steel layer was composed of low-carbon steel. The refractory layer was considered to be made exclusively of a magnesia-spinel brick (Magkor series, RHI Magnesita), whose chemical composition is presented in Table A-I (^Appendix). The coating was assumed to have the same composition as the kiln output flow. The typical chemical composition of a clinker kiln output is shown in Table A-I (^Appendix). The coating physical properties were approximated by the solid bed properties of the kiln described in ²². Fig. 1 schematizes the cross-section view of the clinker kiln wall. The segments R₁, R₂, R₃, and R₄ measured 1.4, 1.5, 1.9, and 2.0 m, respectively. The thickness of each layer was considered constant along the kiln axis. In this study, only the kiln wall was enclosed in the system. The contents of the kiln (i.e., the materials and gases circulating inside it) were not considered.

Figure 1
Schematic of the cross-section of clinker kiln, showing the layers that make up the wall and their dimensions.
Figura 1
Esquema da seção transversal do forno de clínquer, mostrando as camadas que compõem a parede e suas dimensões.

Development of the partial differential equation: given the cylindrical geometry of the system, the starting point was the general heat equation in cylindrical coordinates ²³:

$\frac{1}{\alpha }\frac{\partial T}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)+\frac{1}{r^2}\left(\frac{{\partial }^{2}T}{\partial {\theta }^2}\right)+\frac{{\partial }^{2}T}{\partial z^2}+\frac{\dot{q}}{k}$ (A)

Since the kiln had rotational symmetry with respect to its central axis and was operated in constant rotation, it was considered that the temperature in the kiln wall did not change along the angular direction θ. Then, the derivative of the temperature in relation to θ was neglected. In addition, there was no generation or consumption of heat inside the walls, since there was no chemical reaction inside them. Thus, the generation term $\dot{q}$ was neglected. Finally, the term $\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)$ was expanded by applying the product rule for derivatives, and then, Eq. B was obtained:

$\frac{1}{r}\frac{\partial T}{\partial r}+\frac{{\partial }^{2}T}{\partial r^2}+\frac{{\partial }^{2}T}{\partial z^2}=\frac{1}{\alpha }\frac{\partial T}{\partial t}$ (B)

which is a partial differential equation (PDE) that relates the variation of temperature over time with its variation along the directions r and z. To solve this PDE, the method of lines was applied, in which the spatial dimensions were discretized, and the resulting ordinary differential equation (ODE) was integrated over time. The finite difference method was used to discretize the spatial dimensions, while the modified Euler method (predictor-corrector method) was used to integrate the ODE over time.

Construction and discretization of the domain: for the construction of the domain for the finite difference method, a longitudinal section of the clinker kiln was taken (Fig. 2). It is important to note that the kiln was oriented so that the flame was on the left side, near the origin of the coordinate system. The obtained rectangular domain was 60 m long and 0.6 m high. Fig. 3 shows an illustrative representation of the obtained domain. Thermal contact resistance between the layers was neglected since it was small in comparison to the thermal resistance of the layers. For the domain discretization, the number of nodes in which each dimension was divided was chosen and named npZ, npRs, npRr, and npRc. With these values, the ‘steps’ of the mesh were calculated, that is, the distance between two adjacent nodes. The mesh steps were named Δz, Δr_s, Δr_r, and Δr_c. The mesh steps were calculated by:

$\Delta z=\frac{L}{npz-1}$ (C)

$\Delta r_m=\frac{R_{1+1}-R_1}{npRm-1}$ (D)

where, in Eq. D, the subscript letter m corresponds to the material, which can be equal to s (steel), r (refractory), or c (coating), and the subscript l is l=1 for the coating, l=2 for the refractory bricks, and l=3 for the steel. The nodes obtained when discretizing the mesh were divided into 19 different groups according to their location in the mesh (Fig. 4). Each group of nodes had a different set of equations and boundary conditions related to its location in the domain. In Table A-II (^Appendix), descriptions of each group of nodes are presented, as well as the boundary conditions considered for each one.

Figure 2
Schematic drawing of the cylindrical kiln, showing the longitudinal section plane (yellow rectangle). The rectangular section indicated by the red arrow, which corresponded to the intersection between the section plane and the wall, was the section used to construct the mesh.
Figura 2
Desenho esquemático do forno cilíndrico, exibindo o plano de corte meridional (retângulo amarelo). O trecho retangular indicado pela seta vermelha, correspondente à intersecção entre o plano de corte e a parede, foi a seção utilizada para a construção da malha.

Figure 3
Illustration of the domain, showing the layers in which it divides, and their dimensions.
Figura 3
Ilustração do domínio, explicitando as camadas em que se divide e suas dimensões.

Figure 4
Representation of the mesh obtained by discretization of the domain, detailing the steps of the mesh, and the 19 groups of nodes with equations and boundary conditions in common. The number of nodes in the figure is merely illustrative and does not correspond to the actual amount used for resolution.
Figura 4
Representação da malha de pontos obtida por discretização do domínio, detalhando os passos da malha e os 19 grupos de pontos com equações e condições de contorno em comum. A quantidade de pontos na figura é apenas representativa, não correspondendo à quantidade real utilizada para resolução.

Application of finite difference equations: the formula chosen to express the first derivative present in the PDE was that of central finite differences ^24),(25:

$f\text{'}\left(x\right)\approx \frac{f\left(x+H\right)-f\left(x-H\right)}{2H}$ (E)

where f represents a generical function, f’ represents its derivative, and H is the ‘step’. By applying the Eq. E, the first-order derivative present in Eq. B can be represented as:

$\frac{\partial T}{\partial r}=\frac{T_{i+1,j}-T_{i-1,j}}{2\Delta r_m}$ (F)

where a new notation is introduced. The symbol T_i,j is used to represent the temperature at the node of the position i,j of the mesh, where i represents the position in the r direction, and j represents the position in the z-direction. The Eq. B also presents second-order derivatives and, for those, the formula of ‘symmetric’ finite differences was used ²⁶:

$f\text{'}\text{'}\left(x\right)\approx \frac{f\left(x+H\right)-2f\left(x\right)+f\left(x-H\right)}{H^2}$ (G)

By applying the formula from Eq. G, the second-order derivatives from Eq. B can be expressed as:

$\frac{{\partial }^{2}T}{\partial r^2}=\frac{T_{i+1,j}-2T_{i,j}+T_{i-1,j}}{\Delta r_m^2}$ (H)

$\frac{{\partial }^{2}T}{\partial z^2}=\frac{T_{i+1,j}-2T_{i,j}+T_{i-1,j}}{\Delta z^2}$ (I)

Substituting Eqs. F, H, and I into Eq. B gives:

$\frac{1}{r}=\left(\frac{T_{i+1,j}-T_{i-1,j}}{2\Delta r_m}\right)+\frac{T_{i+1,j}-2T_{i,j}+T_{i-1,j}}{2\Delta r_m^2}+\frac{T_{i,j+1}-2T_{i,j}+T_{i,j-1}}{2\Delta z^2}=\frac{1}{\alpha }\frac{\partial T}{\partial t}$ (J)

It is important to note that some constants in Eq. J vary depending on the wall layer to which they are being applied, as the values of ∆r_m and α vary for each material. Rearranging Eq. J gives Eq. K:

$A{T}_{i+1,j}+B{T}_{i-1,j}+C{T}_{i,j}+D\left(T_{i,j+1}+T_{i,j-1}\right)=\frac{\partial T}{\partial t}$ (K)

where $A=\alpha /\left(2r{\Delta }_m\right)+\alpha /\Delta r_m^2,B=-\alpha /\left(2r\Delta r_m\right)+\alpha /\Delta r_m^2,C=-2\alpha /\Delta r_m^2-2\alpha /\Delta z$, and $D=\alpha /\Delta z^2$. Eq. K, however, cannot be immediately applied to each one of the 19 node groups highlighted in Fig. 4. This is because, in some of these groups, there are problems with some of the variables. For example, the nodes in group 2 are at the edge of the mesh, having no nodes to their left, so the variable T_i,j-1 does not exist for these nodes. The strategy used to deal with these problems was the creation of ghost nodes. This technique is a mathematical tool which consists in assuming the existence of ghost nodes in the vicinity of the studied nodes, whose effect is exactly the effect generated by the boundary conditions. Thus, boundary conditions are indirectly applied to the equations, simulated by the ghost nodes, and this eliminates the mentioned problems presented by certain node groups. The following sections describe the mathematical manipulations to create ghost nodes for the various node groups.

Ghost nodes on the outer surface (groups 17, 18, and 19): in the nodes on the outer surface, the energy received by conduction through steel is delivered to the environment by convection and radiation. This means that, at steady state, the energy balance takes the form of:

${\dot{E}}_{in}={\dot{E}}_{out}$ (L)

where ${\dot{E}}_{in}$ is the rate of energy that enters the surface nodes by conduction, and ${\dot{E}}_{out}$ is the rate of energy that exits the surface nodes by convection and radiation. The term ${\dot{E}}_{in}$ can be replaced by the equation of heat conduction in the r direction. The term ${\dot{E}}_{out}$ can be replaced by the sum of two portions, one corresponding to convection heat transfer and the other to radiation. For convection heat transfer, Newton’s law of cooling can be used ²³. Radiation heat transfer can be described by the equation of the net radiant exchange between two surfaces, considering a small object within a large cavity ²³, since the kiln was relatively small in comparison to the place it was located and it was positioned far from any walls. After the described substitutions, Eq. M was obtained:

$-k_{s}S\frac{\partial T}{\partial r}=h_{ext}S\left(T_{i,j}-T_{ext}\right)+\epsilon \sigma S\left[{\left(T_{i,j}\right)}^4-\left(T_{ext}^4\right)\right]$ (M)

where k_s is the thermal conductivity of steel, S is the heat transfer area, h_ext is the convective heat transfer coefficient, T_i,j is the temperature at the surface node, T_ext is the temperature of the surroundings, ε is the emissivity of the surface, and σ is the Stefan-Boltzmann constant. The derivative ∂T/∂r present in Eq. M can be expressed by the central finite difference formula (Eq. F). The problem is that, since the outer surface nodes are on the edge of the mesh, there is no node ‘one step ahead’ of them in the r direction. Then we assumed a ghost node outside the mesh (Fig. 5), whose temperature T_f is unknown, and we proceed to apply the formula normally. Eq. N was obtained as a result:

$-k_{s}S=\frac{\left(-T_f+T_{i-1,j}\right)}{2\Delta r_s}{h}_{ext}S\left(T_{i,j}-T_{ext}\right)+\epsilon \sigma S\left[{\left(T_{i,j}\right)}^4-\left(T_{ext}^4\right)\right]$ (N)

Figure 5
Representation of the ghost node T_f on the steel surface. The black dots are real nodes, and the white dot is a ghost node.
Figura 5
Representação do ponto fictício T_f na superfície do aço. Os pontos pretos são reais e o branco é fictício.

It is now possible to isolate the value of T_f, obtaining the temperature of the ghost node as a function of other known values, as shown by Eq. O:

$T_f=-T_{i-1,j}-\frac{2\Delta r_s\left[h_{ext}\left(T_{i,j}-T_{ext}\right)+\epsilon \sigma \left(T_{i,j}^4-T_{ext}^4\right)\right]}{k_s}$ (O)

Ghost nodes on interfaces between materials (groups 5, 6, 7, 11, 12, and 13): on the interfaces between two materials, it was considered that the heat received by conduction in the r direction from the innermost layer was entirely transferred to the outermost layer. Thus, in steady state, the energy balance for each interface node takes the form:

${\dot{E}}_{in}={\dot{E}}_{out}$ (P)

Then, for the nodes on interfaces between two materials, called material X (outermost layer) and material Y (innermost layer), it is possible to calculate two different ghost nodes. The interface nodes can be assumed to belong to material Y, and a ghost node within the material X can be calculated by assuming it to be composed of Y. Or, conversely, the interface nodes can be assumed to be in material X, and a ghost node within the material Y can be calculated, assuming that it is composed of X. In both cases, the ghost temperature is calculated, such that the amount of heat transferred is the same as the ‘realistic’ situation, i.e., the situation where no ghost temperatures were used. The strategy of creating a ghost node simplifies the calculations since, assuming the interface consisting of a single material, the Eq. K can be applied with the thermophysical properties of a single material. Fig. 6 illustrates the ghost node strategy adopted at the interface.

Figure 6
Representation of the ghost nodes used for the interfaces between two materials in the cases in which all nodes are considered to belong to material Y (a) and material X (b). The black dots are real nodes, and the white dots are ghost nodes.
Figura 6
Representações dos pontos fictícios utilizados para os nós das interfaces entre dois materiais nos casos em que se consideram todos os pontos como pertencentes ao material Y (a) e material X (b). Os pontos pretos são reais e os brancos são fictícios.

Returning to the energy balance for the interface nodes (Eq. P), the energy coming from material Y that enters the node by conduction can be expressed through the equation of heat conduction in the r direction for material Y. The energy leaving the node in the direction of material X can also be expressed by the conduction equation in the r direction but considering the properties of material X. Thus, we have Eq. Q:

$-k_{y}S\frac{\partial T}{\partial r}=-k_{x}S\frac{\partial T}{\partial r}$ (Q)

where k_y is the thermal conductivity of material Y, k_x is the thermal conductivity of material X, and S is the heat transfer area. Expressing the derivatives of the previous equation in terms of central finite differences gives:

$-k_{y}S\frac{\left(T_{Yf}-T_{i-1,j}\right)}{2\Delta r_y}=-k_{x}S\frac{\left(T_{i+1,j}-T_{Xf}\right)}{2\Delta r_x}$ (R)

where T_Yf and T_Xf are the temperatures of the ghost nodes. By rearranging Eq. R, Eq. S is obtained:

$\frac{k_y}{2\Delta r_y}{T}_{Yf}+\frac{k_x}{2\Delta r_x}{T}_{Xf}=\frac{k_y}{2\Delta r_y}{T}_{i-1,j}+\frac{k_x}{2\Delta r_x}{T}_{i-1,j}$ (S)

The two situations represented in Fig. 6 should lead to the same result when applied to Eq. K, since the amount of heat transferred is always the same, leading to the same changes in temperature. Thus, Eq. T must be valid:

$$\begin{gathered} A_y{T}_{Yf}+B_y{T}_{i-1,j}+C_y{T}_{i,j}+D_y\left(T_{i,j+1}+T_{i,j-1}\right)= \\ {A}_x{T}_{i+1,j}+B_x{T}_{Xf}+C_x{T}_{i,j}+D_x\left(T_{i,j+1}+T_{i,j-1}\right) \end{gathered}$$ (T)

where A_x, B_x, C_x, and D_x are the terms of Eq. K applied to material X, while A_y, B_y, C_y, and D_y are the same terms, but applied to properties of material Y. Rearranging the previous equation gives:

$$\begin{gathered} A_y{T}_{Yf}-B_x{T}_{xf}=A_x{T}_{i+1,j}-B_y{T}_{i-1,j}+ \\ \left(C_x-C_y\right)T_{i,j}+\left(D_x-D_y\right)\left(T_{i,j+1}+T_{i,j-1}\right) \end{gathered}$$ (U)

Eqs. S and U form a system of two equations and two unknowns (T_Yf and T_Xf). The resolution of this system gives the value of the ghost temperatures. With these two values, the temperatures T_Yf or T_Xf can be substituted in Eq. K applied to the corresponding conductive material in place of T_i+1,j or T_i-1,j, respectively. In the present study, it was chosen to always use Eq. K applied to the material of the innermost layer for the interface nodes. Therefore, the ghost temperature T_Yf was used in place of T_i+1,j (the situation is shown in Fig. 6a).

Ghost nodes on the left surface (groups 2, 5, 8, 11, 14, and 17): at the adiabatic surface on the left, conduction in the z-direction is zero. This property can be written in terms of the heat conduction equation as:

$-kS\frac{\partial T}{\partial r}=0$ (V)

where the derivative ∂T/∂r can be rewritten in terms of central finite differences. However, the nodes on the left surface are at the edge of the mesh and do not have a ‘one step back’ node in the z-direction. To solve this, it was assumed the presence of a ghost node of unknown temperature T_f, present ‘one step back’ from each of these nodes (Fig. 7a). So, by rewriting Eq. V in terms of finite differences with a ghost node, we get Eq. W. From this equation, it was possible to find the temperature of the ghost node, as shown in Eq. X:

$-kS\frac{\left(T_{i,j+1}-T_f\right)}{2\Delta r}=0$ (W)

$T_f=T_{i,j+1}$ (X)

Figure 7
Representations of the ghost node T_f on: a) left surface; and b) right surface. The black dots are real nodes and the white dots are ghost nodes.
Figura 7
Representação do ponto fictício T_f na: a) superfície esquerda; e b) superfície direita. Os pontos pretos são reais e os brancos são fictícios.

Ghost nodes on the right surface (groups 4, 7, 10, 13, 16, and 19): for the adiabatic surface on the right, we proceeded analogously to the surface on the left, described in the previous topic. But this time, the ghost nodes were ‘one step ahead’ of the surface nodes in the z-direction (Fig. 7b). We have Eq. Y of finite-difference conduction. Isolating the ghost temperature value in this equation gives Eq. Z:

$-kS\frac{\left(T_f-T_{i,j-1}\right)}{2\Delta r}=0$ (Y)

$T_f=T_{i,j-1}$ (Z)

Calculation of convection heat transfer coefficient: to calculate the convective heat transfer coefficient on the outer surface of the kiln (h_ext), it was considered that the kiln could be modeled as a long cylinder subjected to an air flow. The rotation of the kiln (about 5 rpm) generated a relative movement between air and cylinder. Taking the cylinder as an inertial reference frame, it was as if there was a rotating flow of air around the cylinder. However, due to the lack of correlations in the literature that model rotary flow around a cylinder in the range of Reynolds number considered, it was approximated as a cross-flow in relation to the cylinder. The relative velocity (V) between air and kiln surface was calculated by:

$V=\text{rotation.kiln circumference}$ (AA)

Considering a rotation of 5 rpm and the kiln radius of 2 m (circumference=2.π.2=12.6 m), the relative velocity was ~1 m/s (62.8 m/min). To calculate the film temperature (23), the surface temperature of the cylinder was considered to be on average 250 °C (based on collected industrial data), and the free current temperature to be the temperature defined for the surroundings (25 °C). Thus, the film temperature was (250+25)/2≈140 °C. The properties of air at film temperature are summarized in Table I. From the data shown in Table I, the Reynolds number was calculated (Re_D=145719), as provided in Eq. AB. Hilpert’s empirical correlation ²³ was used to calculate the Nusselt number (Eq. AC). For 40000≤Re_D≤400000, we had C=0.027 and m=0.805. Substituting the relevant values in Eq. AC gave the value of h_ext as 2.9 W/(m².K).

$R{e}_D=\frac{VD}{\nu }$ (AB)

$\overline{N{u}_D}\equiv \frac{h_{ext}D}{k}=CR{e}_D^{m}P{r}^{1/3}$ (AC)

Table I
Atmospheric air properties at film temperature ²⁷.
Tabela I
Propriedades do ar atmosférico na temperatura de filme ²⁷.

Property	Value
Kinematic viscosity (ν)	2.745x10-5 m2/s
Thermal conductivity (k)	0.03374 W/(m.K)
Prandtl number (Pr)	0.7041

Algorithm implementation: in order to implement the computational problem-solving algorithm, Eq. K was applied to each one of the 19 sets of mesh nodes, using the ghost nodes when necessary and employing the properties of each respective material. Thus, 19 different equations were obtained, each for one of the specific sets of nodes. In order to find the steady state temperature profile, a possible method is to integrate the time differential equations using a numerical method from an initial value until a steady state is reached ²⁸. For this, then, an initial temperature matrix was needed, which would serve as initial values for the iterative procedure of the algorithm. Therefore, as in the previous study ²⁹, it was considered that the kiln was divided into different zones, each with a different temperature range. It was considered that the temperature varied linearly within these zones. The starting temperature of the algorithm for a node, as a function of its position in the kiln length (z-axis), is shown in Fig. 8. For simplicity, no variation in the r direction was considered for constructing the initial value matrix.

Figure 8
Graph of temperatures used as initial values for the algorithm as a function of kiln length, adapted from an available kiln internal temperature profile ³⁰.
Figura 8
Gráfico das temperaturas utilizadas como valores iniciais para o algoritmo em função do comprimento do forno, adaptado de perfil disponível de temperatura interna do forno ³⁰.

The 19 equations obtained were ODEs, and they were integrated over time by means of the modified Euler method (predictor-corrector method) as defined in ^{31), (32}. The algorithm was coded in Matlab language. The stopping criterion was a user-defined end time value (t_final). This meant that the algorithm response was a temperature profile corresponding to a specific time instant chosen by the user. In order to verify that the steady state was effectively reached, different values for the final time were used, and the results were compared to determine whether or not there was variation within a tolerance of 0.0001 °C. Step h of the modified Euler method was determined by:

$h=\frac{t_{final}}{nt}$ (AD)

where the number of steps nt was stipulated by the user. The values of some system and surrounding properties, as well as some physical constants used in the algorithm, are shown in Table II. Table III shows the method parameter values that had to be arbitrated for the execution of the algorithm. The description of each symbol displayed in these two tables can be found in section Nomenclature. For comparison purposes, industrial data on the temperature of the external kiln surface were obtained, collected by the industry that owned the studied clinker kiln. The temperatures on the kiln shell were measured by an optical pyrometer along a straight line parallel to the equipment axis.

Table II
Properties and thermophysical constants used in the algorithm.
Tabela II
Propriedades e constantes termofísicas utilizadas no algoritmo.

Property	Value	Source
ks	13.4 W/(m.K)	3
kr	4.0 W/(m.K)	22
kc	0.5 W/(m.K)	22
ρs	8238 kg/m3	23
ρr	3100 kg/m3	33
ρc	1046 kg/m3	22
cp,s	468 J/(kg.K)	23
cp,r	1130 J/(kg.K)	33
cp,c	800 J/(kg.K)	22
Text	298.15 K	*
σ	5.67x10-8 W/(m2.K4)	23
ε	0.54	34

• * defined by the authors.

Table III
Parameter values chosen by the user.
Tabela III
Valores de parâmetros arbitrados pelo usuário.

Parameter	Value
npRs	5
npRr	15
npRc	5
npZ	181
nt	20000
tfinal	400000

RESULTS AND DISCUSSION

After executing the algorithm, the temperature map shown in Fig. 9 was obtained as a result. The temperature distribution pattern displayed on the map was consistent with the principles of heat transfer. The innermost part of the wall (bottom of the figure) was at a higher temperature as it was in contact with hot material circulating inside the kiln. On the other hand, the outermost part of the wall (top of the figure) had the lowest temperatures. In addition, the nearest side of the flame (left side of the figure) was warmer than the far side, as expected.

Figure 9
Temperature map obtained as a result. The white dashed lines that divide the map horizontally mark the division between layers of different materials.
Figura 9
Mapa de temperatura obtido como resposta do algoritmo. As linhas tracejadas brancas que dividem o mapa horizontalmente marcam a divisão entre as camadas de diferentes materiais.

The experimental kiln shell temperature data were compared with the values provided by the algorithm for the wall surface nodes (Fig. 10a). It was seen that the modeled temperature seemed to adjust well to the actual temperatures since the readings of temperature are distributed around the curve of the model. The root mean square error (RMSE) between measured and predicted data was 47 ºC. In Fig. 10, it can be seen that the measured temperature values varied greatly along the length of the kiln, forming a sinuous distribution marked by peaks and valleys. In a published study on estimation of coating formation in a cement kiln ³⁵, a tortuous distribution was also obtained for the experimental temperatures measured on the kiln shell. The authors ³⁵ proposed that sharp variations in the measured temperature could be attributed to irregularities in the thickness of the coating layer. Kiln regions where a thicker coating layer has formed exhibit lower temperatures on the outer surface of the kiln due to the insulating property of the coating. Thus, since the model of the present study considered the thickness of the coating layer constant, it was expected that the present model could not predict the experimental temperature profile sinuosities.

Figure 10
Comparisons between kiln surface temperatures measured experimentally and predicted by the models for: a) this study; and b) reference study ²⁰.
Figura 10
Comparações entre temperaturas da superfície do forno medidas experimentalmente e preditas pelos modelos para: a) este estudo; e b) estudo de referência ²⁰.

Another published study on clinker furnace modeling also came to similar conclusions ²⁰. Comparing the temperature profile predicted by the model with the values measured along the outer surface of the kiln, the authors ²⁰ noted the presence of several fluctuations of up to ±100 K in the experimental temperatures, as shown in Fig. 10b. These fluctuations were attributed to irregularities of the coating layer and erosion in parts of the refractory. The authors appeared to be satisfied that the predicted temperature followed approximately the same general trend observed in the experimental values, claiming that the fit was as good as could be expected given the complexity of a real clinker kiln. According, to ²⁰ an exact comparison between predicted and experimental temperatures is not meaningful, since several characteristics of a real kiln were not considered in the model assumptions, such as the variation in the thickness of the steel and refractory layers, the variation in emissivity (due to shell coloration), and the variation in heat transfer coefficient (due to spacing of the blowers). Therefore, based on the results shown in Fig. 10 and the conclusions of previous works available in scientific literature, it can be said that the model of this study obtained a satisfactory fit to the temperature values measured on the kiln shell. Fig. 10a shows that the model from this study followed the general trend displayed by the experimental values. Oscillations in measured temperatures that were not fitted by the model may have been caused by coating layer irregularities, as they were not considered in the modeling.

CONCLUSIONS

The present study aimed to model the temperature profile inside a clinker kiln wall. The proposed model made use of the heat equation and numerical methods of differential equation resolution to obtain a temperature field inside the wall. The results obtained with the developed model proved to be coherent from the theoretical point of view, i.e., obeyed the principles of heat transfer. Moreover, when compared to industrial temperature measurements, the values predicted by the model were close to the experimental values, with a root-mean-square error (RMSE) of 47 ºC. The quality of the fit obtained by the model was similar to that of other models available in the scientific literature for the shell temperature of a clinker kiln, with deviations that can be explained by irregularities in the coating layer that were not included in the model. One of the main limitations of the model developed in this study is that, since the contents of the kiln were not included in the modeling, the algorithm relies on input for the temperature profile of the internal surface of the kiln, an input that was taken from the literature and may not exactly reflect the real condition of the studied kiln. In addition, simplifications were made for the kiln wall structure, such as the assumption of constant thickness for the wall layers and the consideration of only one type of refractory brick along the entire length of the kiln. A possible application for the model from this study is the comparison of the performances of different refractories by simulating temperature profiles employing different values for the properties of the refractory layer. The model of this study can also be used to infer the refractory wear level by identifying large positive deviations of the experimentally measured temperature from the predicted temperature of the model for the external surface of the clinker kiln. As the next steps for model refinement, it is suggested to consider the different types of refractory bricks along the kiln, with different properties, and to consider the variation of specific heat capacity with temperature.

NOMENCLATURE

cp,s	Specific heat capacity for steel at constant pressure	[J/(kg.K)]
cp,c	Specific heat capacity for coating at constant pressure	[J/(kg.K)]
cp,r	Specific heat capacity for refractory brick at constant pressure	[J/(kg.K)]
h	Modified Euler method time step	[s]
hext	Convective heat transfer coefficient on kiln surface	[W/(m2.K)]
k	Thermal conductivity	[W/(m.K)]
ks	Thermal conductivity for steel	[W/(m.K)]
kc	Thermal conductivity for coating	[W/(m.K)]
kr	Thermal conductivity for refractory brick	[W/(m.K)]
L	Kiln length	[m]
npRs	Number of nodes that steel layer thickness is divided	-
npRc	Number of nodes that coating layer thickness is divided	-
npRr	Number of nodes that refractory layer thickness is divided	-
npZ	Number of nodes that kiln length is divided	-
nt	Number of steps that total time interval is divided	-
NuD	Nusselt number	-
Pr	Prandtl number	-
$\dot{q}$	Heat generation rate	[W/m3]
r	Radial direction of cylindrical coordinates	-
ReD	Reynolds number	-
R1	Kiln radius from central axis to inner surface of coating layer	[m]
R2	Kiln radius from central axis to inner surface of refractory layer	[m]
R3	Kiln radius from central axis to inner surface of steel layer	[m]
R4	Kiln radius from central axis to outer surface	[m]
T	Temperature	[K]
t	Time	[s]
Text	Temperature of surrounding	[K]
tfinal	Final time	[s]
TXf	Ghost temperature for interface composed of X only	[K]
TYf	Ghost temperature for interface composed of Y only	[K]
z	Axial direction of cylindrical coordinates	-
α	Thermal diffusivity	[m2/s]
Δrs	Step in r direction for steel layer	[m]
Δrc	Step in r direction for coating layer	[m]
Δrr	Step in r direction for refractory layer	[m]
Δz	Step in z-direction	[m]
ε	Emissivity of surface material	-
θ	Angular direction of cylindrical coordinates	-
ν	Kinematic viscosity	[m2/s]
ρs	Steel specific mass	[kg/m3]
ρc	Coating specific mass	[kg/m3]
ρr	Refractory brick specific mass	[kg/m3]
σ	Stefan-Boltzmann constant	[W/(m2.K4)]

ACKNOWLEDGMENT

The authors acknowledge the financial support provided by CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior) through the Program “Demanda Social”.

REFERENCES

• ¹
  K.T. Kaddatz, M.G. Rasul, A. Rahman, Procedia Eng. 56 (2013) 413.
• ²
  D.C.Q. Rodrigues, A.P. Soares Jr., E.F. Costa Jr., A.O.S. Costa, Cerâmica 62, 362 (2016) 140.
• ³
  K.S. Mujumdar, K.V. Ganesh, S.B. Kulkarni, V.V. Ranade, Chem. Eng. Sci. 62 (2007) 2590.
• ⁴
  W. Barbosa, K.F. Portella, Cerâmica 65, 373 (2019) 54.
• ⁵
  L. Simão, N.J. Lóh, D. Hotza, F. Raupp-Pereira, J.A. Labrincha, O.R.K. Montedo, Cerâmica 64, 371 (2018) 311.
• ⁶
  T.O. Mason, F.M. Lea, “Cement”, Encyclopædia Britannica (2019) www.britannica.com
  » www.britannica.com
• ⁷
  D.C.Q. Rodrigues, A.P. Soares Jr., E.F. Costa Jr., A.O.S. Costa, Cerâmica 59, 350 (2013) 302.
• ⁸
  T.F. Anacleto, L.F. Turetta, E.F. Costa Jr., A.O.S. Costa, Cerâmica 64, 372 (2018) 507.
• ⁹
  H.M. Souza, A.P. Soares Jr., E.F. Costa Jr., A.O.S. Costa, Cerâmica 61, 357 (2015) 23.
• ¹⁰
  M.L.G. Renó, “Uso de técnicas de otimização robusta multi-objetivos na produção de cimento”, M.Sc. Diss., Univ. Fed. Itajubá (2007).
• ¹¹
  P.A. Alsop, Cement plant operations handbook: for dry process plants, Tradeship Publ. (2007).
• ¹²
  M.A. Aldieb, H.G. Ibrahim, in Proc. World Congr. Eng. Comput. Sci., San Francisco (2010) 699.
• ¹³
  J. Södje, S. Uhlendorf, H.-J. Klischat, Refract. Worldforum 5, 4 (2013) 53.
• ¹⁴
  D.G.M. Silva, W.L. Vasconcelos, Cerâmica 63, 367 (2017) 281.
• ¹⁵
  D.G.M. Silva, W.L. Vasconcelos, Cerâmica 64, 369 (2018) 30.
• ¹⁶
  T.S. Possamai, R. Oba, V.P. Nicolau, Appl. Therm. Eng. 48 (2012) 414.
• ¹⁷
  R. Oba, T.S. Possamai, V.P. Nicolau, Appl. Therm. Eng. 63 (2014) 59.
• ¹⁸
  A. Agrawal, P.S. Ghoshdastidar, Int. J. Heat Mass Transf. 106 (2017) 263.
• ¹⁹
  E. Morsch Filho, T.S. Possamai, V.P. Nicolau, R. Oba, Therm. Sci. Eng. Prog. 16 (2020) 100455.
• ²⁰
  E. Mastorakos, A. Massias, C.D. Tsakiroglou, D.A. Goussis, V.N. Burganos, A.C. Payatakes, Appl. Math. Model. 23 (1999) 55.
• ²¹
  K.S. Mujumdar, A. Arora, V.V. Ranade, Ind. Eng. Chem. Res. 45, 7 (2006) 2315.
• ²²
  K.S. Mujumdar, V.V. Ranade, Chem. Eng. Res. Des. 84, 3 (2006) 165.
• ²³
  F.P. Incropera, A.S. Lavine, T.L. Bergman, D.P. DeWitt, Fundamentals of heat and mass transfer, John Wiley Sons, Hoboken (2007).
• ²⁴
  J.W. Thomas, Numerical partial differential equations: finite difference methods, Springer, New York (1995).
• ²⁵
  Z. Li, Z. Qiao, T. Tang, Numerical solution of differential equations: introduction to finite difference and finite element methods, Cambridge Un. Press, Cambridge (2017).
• ²⁶
  D.M. Causon, C.G. Mingham, Introductory finite difference methods for PDEs, Bookboon, London (2010).
• ²⁷
  Y.A. Cengel, J.M. Cimbala, Fluid mechanics: fundamentals and applications, McGraw-Hill, New York (2013).
• ²⁸
  R.J. LeVeque, Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems, SIAM, Philadelphia (2007).
• ²⁹
  A. Atmaca, R. Yumrutaş, Appl. Therm. Eng. 66 (2014) 435.
• ³⁰
  K.E. Peray, The rotary cement kiln, Chem. Publ., New York (1986).
• ³¹
  G.S. Rao, Numerical analysis, New Age Int., New Delhi (2010).
• ³²
  N. Singh, Int. J. Adv. Res. Comput. Sci. 5, 3 (2014) 41.
• ³³
  Magnesita Refratários, “MAGKOR-B-RA-LP-JP”, Techn. data sheet (2010).
• ³⁴
  Engineering ToolBox, “Emissivity coefficients materials”, 34 Engineering ToolBox, “Emissivity coefficients materials”, www.engineeringtoolbox.com , ac. 09/02/2019.
  » www.engineeringtoolbox.com
• ³⁵
  S. Sadighi, M. Shirvani, A. Ahmad, Can. J. Chem. Eng. 89, 1 (2011) 116.
• ³⁶
  Magnesita Refractories, “Refractories for cement rotary kilns”, Hilden (2014).
• ³⁷
  D.C.Q. Rodrigues, A.P. Soares Junior, E.F. Costa Junior, A.O.S. Costa, An. Acad. Bras. Ciênc. 89, 4 (2017) 3123.

APPENDICES

Table A-I
Chemical composition of Magkor A1 refractory brick (wt%) ⁽³⁶⁾, and output stream of an industrial clinker kiln (wt% calculated from molar fraction data in ⁽³⁷⁾).
Tabela A-I
Composição química de tijolo refratário Magkor A1 (% mássica) ⁽³⁶⁾ e corrente de saída de forno de clínquer industrial (% mássica calculada a partir de dados de fração molar em ⁽³⁷⁾).

Compound	Magkor A1	Output stream
MgO	89.5	-
Al2O3	9.0	11.8
Fe2O3	0.4	3.0
CaO	0.8	65.3
SiO2	0.3	20.0

Table A-II
Mesh node groups with equations and boundary conditions in common.
Tabela A-II
Grupos de nós da malha com equações e condições de contorno em comum.

Group	Mesh location	Boundary conditions
1	Inner surface of the kiln (r=R1)	Constant temperature over time due to continuous heat supply from the flame
2	Lateral surface close to the flame (z=0) in the section composed only by coating (R1<r<R2)	Adiabatic surface (no heat flow perpendicular to the surface)
3	Inside the coating layer (0<z<L; R1<r<R2)	No boundary condition
4	Lateral surface away from the flame (z=L) in the section composed only by coating (R1<r<R2)	Adiabatic surface (no heat flow perpendicular to the surface)
5	Lateral surface close to the flame (z=0) at the interface between the coating and refractory (r=R2)	Adiabatic surface in the z-direction and heat received from the innermost layer is transferred by conduction to the outermost layer
6	Inside the interface between the coating and refractory (0<z<L; r=R2)	Heat received from the innermost layer by conduction in the r direction is transferred by conduction to the outermost layer
7	Lateral surface away from the flame (z=L) at the interface between the coating and refractory (r=R2)	Adiabatic surface in the z-direction and heat received from the innermost layer is transferred by conduction to the outermost layer
8	Lateral surface close to the flame (z=0) in the section composed only by refractory bricks (R2<r<R3)	Adiabatic surface (no heat flow perpendicular to the surface)
9	Inside the refractory layer (0<z<L; R2<r<R3)	No boundary condition
10	Lateral surface away from the flame (z=L) in the section composed only by refractory bricks (R2<r<R3)	Adiabatic surface (no heat flow perpendicular to the surface)
11	Lateral surface close to the flame (z=0) at the interface between refractory and steel (r=R3)	Adiabatic surface in the z-direction and heat received from the innermost layer is transferred by conduction to the outermost layer
12	Inside the interface between refractory and steel (0<z<L; r=R3)	Heat received from the innermost layer by conduction in the r direction is transferred by conduction to the outermost layer
13	Lateral surface away from the flame (z=L) at the interface between refractory and steel (r=R3)	Adiabatic surface in the z-direction and heat received from the innermost layer is transferred by conduction to the outermost layer
14	Lateral surface close to the flame (z=0) in the section composed only by steel (R3<r<R4)	Adiabatic surface (no heat flow perpendicular to the surface)
15	Inside the steel layer (0<z<L; R3<r<R4)	No boundary condition
16	Lateral surface away from the flame (z=L) in the section composed only by refractory bricks (R2<r<R3)	Adiabatic surface (no heat flow perpendicular to the surface)
17	Crossing between the lateral surface close to the flame (z=0) and the outer surface (r=R4)	Adiabatic surface in the z-direction; the heat received by conduction is dissipated by convection and radiation in the r direction
18	Outer surface of the kiln (0<z<L; r=R4)	The heat received by conduction is dissipated to the environment by convection and radiation
19	Crossing between the lateral surface close to the flame (z=L) and the outer surface (r=R4)	Adiabatic surface in the z-direction; the heat received by conduction is dissipated by convection and radiation in the r direction

Publication Dates

• Publication in this collection
  30 Oct 2020
• Date of issue
  Oct-Dec 2020

History

• Received
  20 Jan 2020
• Reviewed
  23 Apr 2020
• Reviewed
  10 July 2020
• Accepted
  14 July 2020