IMPLEMENTATION OF A MATHEMATICAL MODELLING OF A ROTARY CEMENT KILNS

Rotary cement kiln is the main equipment in the cement industry that has complex dynamic behavior, where any changes will affect the quality of the product and the consumed energy. A one-dimensional model of rotary kiln is needed to understand kiln’s behavior and improve kiln operating and design to achieve the optimum condition of product quality and energy required. In this study, the onedimensional mathematical model of a dry rotary cement kiln with pulverized coal combustion is developed. This model consists of a set of nonlinear ordinary differential equations and nonlinear algebraic equations that describe material and energy balance equations. The model has been solved numerically by using Matlab R2015a, and it has been validated by comparing the result with published experimental data. Based on the result, the steady-state simulation shows that the behavior of the model developed is appropriate with the results presented in the literature. It can be concluded that the model is accurate (error < 6%) to describe the profile of temperature and bed composition along with the kiln. It can be used to obtain a better understanding of kiln’s behavior and improve the kiln operating and design to achieve the optimum condition.


INTRODUCTION
Rotary kilns are commonly unit that used in chemical, metallurgical, and pharmaceutical process industries. Rotary kilns are used for heating, reacting, drying processes of solid material. In many cases, it can be used to combine those processes, such as the uses of rotary kiln in the cement industry [1] . Rotary cement kiln is the main equipment in the cement industry that can be used to produce a clinker from raw material through a complex reaction within the kiln by using a high-temperature condition. Process along the kiln consists of complex reaction and heat transfer between inside wall, freeboard gas, and solid material. Besides, rotary cement kilns are not only the main equipment but also the unit with the highest energy required. This unit requires a continuous heat supply by the combustion of fuel to support the heat of reactions needed [2] . These reasons make the process in the rotary kiln have a dynamic behavior, where any changes will affect the quality of the product and the consumed energy. So, the performance of the whole process in the industry will depend on the performance of the rotary cement kilns [3] .
One of the factors that can be used to make sure that the product has achieved the qualification desired is bed temperature. This parameter will determine that the reactions along the kiln occur as expected. But, due to the rotary kilns are operated in the high-temperature condition, these units are generally not equipped with a temperature sensor. It makes the profile of temperature along the kiln is unknown. Therefore, the model that can be used to show the behavior of each component along the kiln is needed. Due to these reasons, it prompted the researchers to develop an appropriate model to simulate the rotary kiln both by the one-dimensional model and three-dimensional model using CFD. The model can be used for obtaining a better understanding of kiln's behavior and improving the kiln operating and design to achieve the optimum condition of product quality and energy required.
The researchers have tried to develop a rotary kiln model that considers a heat transfer phenomenon and the kinetic reaction of clinker formation along with the kiln. The first model has studied by Spang. Spang has developed a dynamic model with a simple flame model as a heat supply of the kiln. This result shows that the profile of temperature and the composition along the kiln obtained is matched with the actual profile of kiln qualitatively. Still, as quantitatively, this model is not accurate enough, especially the wall temperature [4] . In 1989, Barr et al. [5] have studied the model of rotary kiln. Their mathematical model developed has been compared with the pilot plant of kiln [6] . Then, Bar's model has been developed by Li et al. [1] , especially for evaluate the heat transfer coefficient of the covered wall and exposed bed in the rotary kiln. Both Barr's model and Li's model only focus on the heat transfer phenomena along with the kiln and negligible the chemical reactions occurred [1] . In 2006, Mujumdar and Ranade [7] developed a one-dimensional model 1D of rotary kiln. Their result shows that the model is entirely accurate to predict the composition of the outlet clinker. But the maximum bed temperature achieved is under the range of actual temperature. In 2017, Csernyei and Straatman [2] develops a numerical model that considers heat transfer on the outside wall of the kiln due to the installation of a cooling fan on the outer wall. Due to this model uses a different model of clinker formation with Mujumdar's model, Csernyei's model only success in improving the bed temperature profile. Still, the composition model has a significant error [2] .
The present study aims to obtain the one-dimensional mathematical model with an accurate model to show the profile temperature and composition along with the kiln. This study will combine the model developed by Mujumdar as a clinker formation model and Csernyei as a heat transfer model.

Mathematical Model
A schematic of rotary cement kiln is shown in the Figure. 1 . Rotary cement kiln operated at 2°-5°tilt to the horizontal with the rotation speed of approximately 1-5 rpm [2] . The feed of cement kiln, which is a mixture of limestone, silica, alumina, and iron oxide called a raw meal. Based on Figure 1 . Raw meal enters the elevated end and flows down by gravity. In another direction, the hot gas flows as a heating medium. When the solid bed flows down, the heat transfer that occurs along the kiln makes the temperature of solid bed reach the temperature target. The reaction of clinker formation, both as endothermic and exothermic occurs. Rotary cement kiln is divided into four stages. They are preheating, calcining, burning, and cooling [8] . Description of the model and solution methodology used is described in this chapter.
The mathematical model developed in this studied is a complete one-dimensional dry kiln model that includes mass and energy balance for three phases within the rotary kiln: freeboard gas (pulverized coal included), solid bed, and wall. The model will follow these assumptions: The model developed for a steady-state process. Specific heat, heat transfer coefficient, heat reaction, solid, and gas velocity along the kiln are constant. The mixing effect is neglected. Both the freeboard gas phase and solid bed phase follow the plug flow model. Chemical reactions that occur are simplified to be five main reactions, which described in Table 1 . The Arrhenius equation determines the rate of reactions [7] .

Bed Height Model
It requires considering the appropriate model to describe the axial flow of solid particles along with the kiln in developing a model for simulating the unit operation of the rotary cement kiln. The appropriate model of axial motion is essential. It is due to axial movement having an impact on the calculation of residence time of solid particle and the bed height calculation, where the two parameters impact the accuracy of the model that was developed.
The calculation of Bed height along the kiln depends on some aspects. They are volumetric flow rate of the solid particle (Φv), angle of inclination of the kiln (γ, radian), repose angle of the solid particle (β, radian), Inner radius of the kiln (R, m) and speed of kiln rotation (n, rps). These parameters are represented in Figure 2 . One of the models of bed height calculation that commonly used is the Kramers model. The Kramers model is represented by Equation 1. (1)

General Material Balance
The general rate equation of any compound in a solid phase is given in Equation 2. The amount of that compound flows into the region is equal to the amount generated or consumed by a chemical reaction in the region.
Where is the mass fraction of compound, is the cross-sectional area of solid bed [ 2], is the velocity of the solid bed [ ∕ ], is the bulk density of the bed [ ∕ 3], is the specific reaction rate [ ∕ 3. ] and the term of is the molecular weight of compound [ ∕ ]. By summing the material balance of each component, the overall mass balance can be expressed by Equation 3. Based on the assumptions above, the reaction rate of each component is shown in Equation 4. According to the reactions in Table 1 , the compounds that include in this model are , 2 2 ), tri-calcium aluminate ( 3 i.e., 3 . 2 3 ) and ferrite phase ( 4 i.e., 4 . 2 3 . 2 3 ). Therefore, there is 13 set of a nonlinear ordinary differential equation of the mass balance equation.

Heat Transfer Model
The specific energy balance equation of the cement rotary kiln is shown in Equation 5, 6, 7, and 8. The term coefficient heat transfer of each phase was obtained by the previous study [9] . To estimate the temperatures of the inner wall and outer shell, a steady-state heat balance across kiln walls are expressed by the following equations.
The equation was solved by the Newton Raphson method for the nonlinear equation. and are the heat transfer rates due to radiation and conduction, respectively, between the internal kiln wall and the bed. and are the heat transfer rates due to radiation and Conduction, respectively, between the internal kiln wall and the freeboard gas.
is the heat transfer due to conduction on kiln shell. and are the losses from the outer shell due to convection and radiation.
is the heat transfer due to conduction on kiln shell. and are the losses from the outer shell due to convection and radiation.

Combustion Model
The fuel that commonly uses the combustion process within the kiln is pulverized coal. In this study, the model of combustion follows the single reaction bellow (Equation 9).
Based on the reaction between the solid particle of pulverized coal and oxygen, the combustion rate model of rotary kiln can be modeled by a shrinking core model (SCM). SCM is a model that can be used to describe the reaction between solid particles and fluid. SCM divided into two types, that is: 1. SCM for a constant spherical particle 2. SCM for shrinking spherical particle For the first type, this type will be used for modeling the reaction that produces a solid phase. Whereas, the second type will be used for modeling of reaction without solid phase product occurs [10] . Based on the single reaction of combustion (Equation 9), the combustion reaction within the kiln does not produce any solid particle, so the model that be used is the second type. Based on The mechanism of second type SCM is shown in 3 . Based on Figure 3 , the mechanism of second type SCM consist of three steps, there are : 1. Diffusion of reactant A from main body gas to the solid surface through gas film 2. The reaction between reactant A and particle B on the solid surface 3. Diffusion of gas product from solid surface to the main body through gas film Based on the SCM model, the reaction rate of a single particle of pulverized coal can be calculated using Equation 10. The rate reaction in equation 10 is the rate reaction based on the surface of each particle. In this type, there are two resistances, surface reaction and gas film may play a role. Based on the literature, for coal combustion, the resistance that controls all of the reaction is a surface reaction. The coal combustion has high porosity, so it does not offer any resistance to gases leaving the reaction zone, and the gas film resistance can be ignored [7] . Equation 10 can be simplified as Equation 11.
The total reaction rate of pulverized coal combustion based on the gas volume through the kiln can be calculated by Equation 12. Then, the material balance of pulverized coal can be calculated by Equation 13.

Solution Methodology
In this study, this simulation consists of three sub-models; there are bed height variation model, bed model, and freeboard gas model. These models were solved separately by using Matlab 2015a. The solution procedure is described in Figure 4 . Based on Figure 4 , the first model that has to be solved is the bed height variation model. The bed height model used is Kramers's   . Based on the Kramers model, the variation of the bed height in rotary cement kilns depends on the angle of kiln tilt, the inner radius of the kiln, speed of kiln rotation, and volumetric of the solid particle on the kiln. Due to this model is an ordinary differential equation, so function ODE 45 is chosen as the method to solve this model. Then, the result of this calculation was used to calculate The speed of the particle in the rotary kiln, heat transfer area, and the cross-sectional area of the bed phase and freeboard gas phase were used to solve the bed model and freeboard gas model in the next step. Using the solution of the first step, the bed model was explained by assuming a guessed temperature profile of freeboard gas along with the kiln. The temperature profile of freeboard gas was predicted by setting the gas temperature in the three-point along with the kiln. They are inlet temperature, outlet temperature, and the maximum temperature of freeboard gas. The maximum flame temperature was variated at the range of 1900 K to 2400 K. This is because the model is consists of the simultaneous ordinary differential equation for bed temperature model and bed composition and system of nonlinear algebraic equation for wall temperature calculation. Therefore, to solve this model, the combination of the ODE 45 and Newton Raphson method can be used. The result obtained from this step is the profile of the bed temperature and profile of chemical composition along with the kiln. The solution of bed model was used to solve the last model, the freeboard gas model. The new freeboard gas profile obtained in this step was compared with the guessed gas temperature profile. The maximum allowable error was set at 1%.

Validation
Before simulating rotary kiln, the proposed model is tested with the steady-state operation to reproduce the result obtained by Mujumdar and Ranade using the data in Table 2 and Table 3 .

Simulation and Model Validation of Rotary Kiln
Rotary cement kiln is generally operated at very high operating temperatures, why it makes the data sampling inside the kiln is restricted during the kiln operations. The computational model can be used to show the profile of the temperature of the bed, the temperature of freeboard gas, and the composition of the clinker along with the kin. For this purpose, the computational model was used to simulate industrial cement kilns' behavior and optimize the operating conditions. Before carrying out the

FIGURE 5
Typical mass fraction profiles in the bed region.

FIGURE 6
Typical temperature profile along with the kiln.
optimization study of rotary kiln, the proposed model was compared with the steady-state operation shown in Table 2 and Table 3 . The result of the validation process is shown in Table 4 . The model was able to predict the performance of kiln satisfactorily. Two aspects make the present model more accurate. For the first, this model considers the variation of bed height along with the kiln, which most of the models for cement kiln assume a constant bed height. Developing an accurate model for the variety of the bed height model is essential. It controls the residence time, axial velocity of solid, heat transfer area along with the kiln. All parameters have a direct effect on the reaction of the clinker formation and heat transfer process along with the kiln. Besides that, assuming a constant bed height will fix the residence time of the solid particle in the kiln for any variation operating condition, so it will not be able to analyze the influence of the kiln's operating parameter. The second aspect is this model consist of a mathematical model. It describes the freeboard gas temperature profile based on mass and energy balance. The previous model just predicts the gas temperature profile based on the interpolation of three-point data (inlet gas temperature, maximum temperature, and outlet temperature). This point is obtained by trial and error. Because the profile is based on the interpolation data, the gas temperature profile looks like a linear graph.
In contrast, the freeboard gas temperature profile should be a nonlinear graph based on the experimental measurement. So, it makes the profile of freeboard gas is not accurate enough to predict the temperature along with the kiln. Temperature is one of the critical parameters of the kiln operation, so if the model did not accurately predict the temperature profile, the bed composition is not accurate.
Typical profiles of bed height, mass fraction, and temperature along with the kiln obtained by simulation are shown in Figure  5 , 7 . For those figures, it demonstrated that abscissa of 0 denotes the position of solid particle enters the rotary kiln. On the other hand, the abscissa of 1 indicates the location of the burner end. Figure 7 shows the profile of bed height along with the rotary cement kiln. It shows that the height of solid on the kiln decrease along with the kiln. It is due to the solid particle become liquid phase when the position of the particle is close to the burner end, where at this position, the reaction of clinker formation occurs. Figure. 5 shows the profile of the composition of the solid particle along with the kiln. After the raw meal feed to the kiln, calcination reaction and 2 formation occur. The next reaction follows this reaction. They are the reaction of 3 and C4AF formation. The solid particle becomes melt after it reaches the burning zone. The burning area is the hottest region along with the kiln so that the solid temperature at this zone is the highest temperature along with the kiln. In Figure. 6, it is shown  that the maximum temperature of the solid particle occurs when the particle reaches 80% of kiln length. At the burning zone, the amount of melt increase until the temperature of the bed crosses the solidus temperature and then reaches the maximum temperature. After that, the temperature decreases and reaches the cooling zone.
The Influence of Operating Parameter to The Amount of Energy Required The influence of rotational speed and kiln tilt to the kiln performance is determined by energy required by coal combustion per unit of weight of clinker coming out of the kiln ( ) and the composition of 3 on the product. The dimensional of kiln and input composition is the same as the data in Table 2 and Table 3 , but the rotational speed and kiln tilt are varied. The rotational speed is varied in range 3 -8 RPM and the kiln tilt varies in the range of 2°-5°. To compare the kiln performance of each variable based on the energy required, the product composition should be the same (error 1%). Thus the mass flow rate of the coal is adjusted to achieve a 50% mass of C3S in the product. In addition, to analyze the influence of all variables on the composition of 3 on the product, the mass flow rate of the pulverized coal is adjusted at 1.25 kg/s. Then the composition of C3S on the product is obtained from the result of the simulation. Then the result will be analyzed by using Design Experiment software to check the influence of each parameter.
A design expert will model the effect of all variables (kiln tilt and speed of rotation) to the energy required and composition of  (5) forward stepwise regression. The best model optimization by RSM is the model with the smallest error [11] . This section will discuss the effect of operating parameters (such as the speed of rotation and kiln tilt) to the amount of energy required. The significant factor analysis was evaluated by identifying the significant factor in the regression model developed by Central Composite Design (CCD). The response of energy required for each parameter was fitted to four types of models by using Design-Expert software. They are linear, two-factor interaction (2FI), quadratic, and cubic polynomials. The best model should have a high coefficient of determination R-square ( 2 ), low standard deviation (SD), low PRESS (predicted residual sum of the square), and the model was not aliased [12] . The sequential model sum of squares is carried out to find a precision type of regression model for the next step (optimization process). The results are shown in Table 5 . The comparative result showed that the cubic model has the lowest SD and the highest 2 . Although the cubic model has the lowest SD and the highest 2 , it was found that the cubic model to be aliased. It lacks runs to support a full cubic model in case of a CCD. While the quadratic model has met all of the criteria (high 2 , low SD, low PRESS, and not aliased). In general, the higher the R-square, the better the model fit the data. The value of R2 is 0.9994, which indicates the quadratic model, explains all the variability of the response data around its mean. Besides, the quadratic model has a lower SD of 3.76, which revealed that the predicted values would be closer to the actual value. Therefore, Design-Expert software suggested the quadratic model as the best model.
The significance of the quadratic model and the interaction of parameters on the responses are determined by statistical analysis of variance (ANOVA) in Table 5 . The importance of each parameter is checked by the probability of error value (p-value) [12] . Model terms with a p-value of less than 0.05 are considered significant. Based on the ANOVA result, it can be concluded that both speed-of rotation and kiln tilt have a substantial impact on the amount of energy required (due to the p-value of each parameter < 0.05). Thus, the model of the amount of energy needed is written by Equation. 14. In Equation 14, the positive sign before the terms specifies the synergistic effect, whereas negative sign specifies the antagonistic effect. As the rotational speed and kiln tilt increase, the amount of energy requires an increase (see Figure. 8). It is due to the inclination slope, and rotational speed is inversely proportional to the residence time. As the rotational speed and kiln tilt decrease, the residence time of bed solid increases. Since solid spends more time in the kiln, re-solidification of the melt is more in the kiln that low rotated and low angle of slope. Moreover, the change of speed rotation and kiln tilt can also change the area of bed exposed to the freeboard gas. As the increase of kiln rotational speed and angle of tilt, the area of bed exposed to the freeboard gas decrease (see Figure. 9). Therefore, at higher RPM and angle of the tilt, the heat transfer between freeboard gas and the solid bed isn't maximized and needs more supply energy from the coal combustion to achieve the target of specification of product.

CONCLUSION
The present study shows that combining a model of reaction clinker formation by Mujumdar and a heat transfer model that considers the installation of a cooling fan on the outside kiln wall can produce the appropriate model with the maximum error less than 6% (less than another model). So, it can be concluded that the present model is accurate enough, and it can be used for obtaining a better understanding of kiln's behavior and improving the kiln operating and design to achieve the optimum condition. Besides, by doing the simulation of rotary cement kiln by variation of speed rotation and angle of inclination, it can be concluded that speed of rotation and angle of kiln tilt significantly influenced both the amount of energy required. Both two parameters have a synergistic effect on the energy required.