Coal and Biomass Co-Combustion: CFD Prediction of Velocity Field for Multi-Channel Burner in Cement Rotary Kiln ^†

Abstract

This paper represents the medialization of alternative fuels co-combustion, in a cement rotary kiln, established on the commercial computational fluid dynamic (CFD) software ANSYS FLUENT. The focus is placed on the key issues in the flow field, mainly on how they are affected by turbulence models and co-processing conditions. Real data, from a Moroccan cement plant, are used for model input. The simulation results have shown a potential effect of the physics model on turbulent and gas-solid flow prediction. The CFD results can be taken as a guideline for improved co-processing burner design and reduce the effect of using alternative fuels.

1. Introduction

Alternative fuels (AF) have different chemical and physical properties compared to fossil fuels [1]. Thus, AF particles differ from coal in shape and size, their aerodynamics, and their heating up and combustion mechanism. Co-firing coal/OP (olive pomace) modeling in a cement kiln burner is an intricate process that includes gas and particle phases as well as a turbulence effect on the combustion mechanism and heat transfer.

Although much work [2,3,4,5,6] has been done on AF co-combustion, few researches on CFD (computational fluid dynamic) modeling of co-firing AF under rotary cement kiln burner conditions are available in the literature. Even though the efficiency of different K-epsilon models has been improved in recent years on predicting PF combustion, most improvements have been achieved by minimizing the calculation cost. Nonetheless, it is possible to further improve efficiency by comparing these models to choose a well-predicting model for novel alternative fuel. With this goal, this work seeks to develop a comparison between the three K-epsilon models (standard (STD), RNG (re-normalization group), and RKE (realizable K-epsilon)) combined with the eddy dissipation model, in a 2D modeling case of co-firing coal with olive pomace, under real rotary cement kiln burner conditions. The object is to explore how the flow field is affected by the turbulence models and how the co-processing environment can change the turbulence shape in the rotary kiln. The goal is to have a model that could well predict the turbulence behavior in the kiln environment, could take into account the difference in physical properties of AF versus coal, and could be able to give an overview of the kiln process parameters when changing the co-firing rate.

The remainder of the paper is organized as follows: The proposed mathematical models are presented in Section 2. The computational methodology is discussed in Section 3. Section 4 shows the results of different cases. Finally, Section 5 concludes with a summary.

2. Mathematical Models

The co-combustion in rotary cement kiln burner startup by carried combustibles particles, from the pulverizer to the furnace, using pressurized air. The flow in the rotary kiln is a gas-solid turbulent flow with chemical reactions. The Eulerian-Lagrangian approach is usually used and recommended for modeling dilute gas-solid flows in pulverized fuels (PF) combustion applications. In the Eulerian-Lagrangian approach, the motion of the continuous phase is modeled using the Eulerian framework, and the motion of dispersed-phase particles is modeled using the Lagrangian framework [7]. In this work, the focus was placed on the comparison of K-epsilon model varieties coupled with turbulence-chemistry interaction model (FR/ED (finite rate/eddy dissipation), i.e., how each model impacts the flow field).

2.1. Gas-Phase Governing Equations

The time-averaged steady-state Navier–Stokes equations, as well as the mass and energy conservation equations, are solved. The governing equations for the conservations of mass, momentum, energy, and species are given as:

$$\nabla \cdot \left(\rho \overrightarrow{u}\right)=S_m,$$

(1)

$$\nabla \left(\rho \overrightarrow{u}\overrightarrow{u}\right)=-\nabla P+\nabla \cdot \left(\stackrel{=}{\tau }\right)+\rho \overrightarrow{g}+\overrightarrow{F},$$

(2)

$$\nabla \cdot \left(\overrightarrow{u}\left(\rho E+P\right)\right)=\nabla \cdot \left(k_{eff}\nabla T-{\displaystyle \sum _i{h}_i{\overrightarrow{J}}_i}+\left(\stackrel{=}{\tau }\cdot \overrightarrow{u}\right)\right)+S_h,$$

(3)

$$\nabla \cdot \left(\rho {\overrightarrow{u}}_i{Y}_i\right)=-\nabla \cdot {\overrightarrow{J}}_i+R_i+S_i,$$

(4)

The source S_m in Equation (1) is the mass added to the continuous phase due to the combustion of particles. Equation (2) is the steady-state momentum equation, where P is static pressure, $\stackrel{=}{\tau }$ is the shear forces, and it is called a stress tensor. The third term on the right-hand side (RHS) represents the gravitational body force and the last term is the drag force source term. The first three terms on the RHS of the energy equation in Equation (3) represent the energy transfer due to conduction, species diffusion, and viscous dissipation, respectively. $S_h$ is a source term including the heat of chemical reaction and radiation. Radiation heat transfer is calculated using the P1 model, and the absorption coefficient is evaluated by the weighted sum of the gray gas model (WSGG) [8,9]. Equation (4) is the species mass fraction conservation equation for the $i^{th}$ species. $R_i$ is the chemical reaction production rate of species $i$ (this rate will be calculated by the FR/ED turbulence-chemistry interaction models) and $S_i$ is the source term added from the dispersed phase. $\overrightarrow{J_i}$ is the turbulent diffusion flux of species $i$ [10].

2.2. Particle-Phase Governing Equations

Alternative fuels are bigger, lighter in density, and are less spherical compared to coal, which affects the motion of particles in the furnace. In the present work, for simplification, we considered AF as a spherical particle; also, the AF and coal particles were modeled separately with two discrete phases following a Rosin-Rammler size distribution [11]. The particles were tracked in Lagrangian frame reference using a stochastic model [12] as shown in Equation (5).

$$\frac{d{u}_p}{dt}=\frac{3\rho C_D}{4d_p{\rho }_p}\left|{\overline{u}}_p-\overline{u}\right|\left(\overline{u}-{\overline{u}}_p\right)+\overline{g}\frac{{\rho }_p-\rho }{{\rho }_p},$$

(5)

where $\overline{g}$ is the gravitational force, and $C_D$ is the drag coefficient, which is an empirical function of $\mathrm{Re}$ as described by Morsi and Alexander [13]. Here $\mathrm{Re}$ is the relative Reynolds number, and it’s written as,

$$\mathrm{Re}=\frac{\rho d_p\left|u_p-u\right|}{\mu },$$

(6)

In this work, the combustion of solid particle conversion is treated as heating, devolatilizing, and char burning process. For computing volatilization, the single kinetic rate model [14] was used. The kinetic-diffusion limited model [15] was used for char combustion modeling. For both volatilization and char combustion, only convection and heat released by reactions were computed.

3. Numerical Computation and Strategy

3.1. Geometry and Grid Description

The rotary kiln used in this paper was assimilated to a real kiln with 46 m in length and 3.8 m in diameter, and it was specially equipped with a multichannel burner (see Figure 1). For simplification, the calculation was started from the fuel injection (i.e., from the burner nozzle end). Therefore, the domain of calculation considered was a 2D axisymmetric rectangle with a multi-entry for fuels and airs. Also, the burner entry geometry was modified for modeling facilities (for more details, see [16]).

Based on the geometric and calculation purposes, a structured non-uniform mesh was used, and a fine grid with 14,040 elements was selected, based on the mesh independence. To optimize convergence, the mesh had a higher density near to the burner and on the axial direction, and was coarse at the rest of the domain.

3.2. Numerical Method and Boundary Condition

The simulation of the 2D domain was implemented on the commercial CFD tool, ANSYS Fluent, based on the finite volume method. A pressure-based coupled algorithm (Coupled) was used to solve the Lagrangian-Eulerian model for gas and particle equations. To achieve numerical stability, the under-relaxation factors were carefully modified for each k-epsilon model. The limited residuals values were taken as 10⁻³ for continuity, 10⁻⁵ for energy and radiation, and 10⁻⁴ for all the rest of the parameters.

To be consistent with the real kiln operating condition, the inlet for temperature and velocity was similar to the real case data (see

Table 1. Operating condition data.

Boundaries	Velocity m/s	Temperature °K
Secondary air inlet	5	1126
Jet air inlet	218	500
Fuels inlet	Transport air 30 Coal and AF 30	500
Swirl inlet	Axial 212 Tangential 45	500
Central air inlet	200	500

). A velocity inlet was specified at all the inlet boundaries with direction normal to boundary except for the swirl inlet that had two velocity components inlet. All the walls were treated as adiabatic and no-slip walls. The outlet was set as a pressure outlet condition. The K and epsilon inlet conditions were computed from the turbulent intensity and hydraulic diameter.

3.3. Case Study

As mentioned earlier, this paper discusses the effect of the different K-epsilon turbulence models on the velocity field prediction in a cement rotary kiln under co-combustion parameters. To this aim, six different cases were investigated and are summarized in

Table 2. Investigated cases details.

Case	Turbulence Model	Turbulence Chemistry Interaction Model	Discrete Phase Model (DPM)
1	STD	OFF	OFF
2		ON	ON
3	RNG	OFF	OFF
4		ON	ON
5	RKE	OFF	OFF
6		ON	ON

.

4. Results and Discussion

To investigate the effect of the turbulence model on the flow field, we started with an initial test without the DPM model and species reaction for the three k-epsilon models. Figure 2 shows the axial velocity streamlines without the DPM model and species reaction simulations. One can observe that the STD and RNG models predict the overall flow field similarly, however, the velocity field in RKE simulation extends further into the rotary kiln than STD and RNG, before taking a uniform distribution (see Figure 3).

The three models differ in the capture of the recirculation zone length. The RNG has a wide inner bubble recirculation zone compared to the STD and RKE at around 6.8–10.6 m. The RKE simulation gives a small inner bubble recirculation zone but a longer recirculation circle at around 9.7–14 m. The recirculation zone given by the STD is around 4.5 m. The difference comes from the epsilon equation forms and the solution method.

From Figure 4, for the STD and RNG cases, it can be seen that when the pulverized particles mixe into the flow, the uniform velocity field is deteriorated from 13.24 to 9.21 m/s for STD and from 10.8 to 8.58 m/s for RNG. Also, the recirculation zone is smaller than the cases without DPM as shown in Figure 5; this is due to the interaction of particle-flow.

Generally, particles interact with the flow through the drag forces and thermodynamics effect, causing a reduction in the gas turbulence and an increase in the dissipation rate of that energy which explain the decay of velocity; for more details, see [17]. For the RKE model case, we remark that the uniform velocity is higher in the case with particle injection than the DPM-off case (7.22 m/s vs. 2.54 m/s).

From the results illustrated, it is clear that the RKE model in both test cases gives the lower velocity compared to the RNG and STD models, which may be due to the dissipative nature of the RKE model. Another important point is that the STD and RNG give the lower reverse velocities, as shown in Figure 5, ensuring a speedy mixture of the hot gases with cold burner stream and the stabilization of the flame in the kiln compared to the RKE cases, which impacts the species and temperature distribution that will be discussed in another paper. Current numerical simulations of AF co-combustion are compared with the numerical results obtained by [18,19] who have investigated the simulation of oxy-coal combustion with RANS (Reynolds averaged Navier stokes) and LES (Large eddy simulation) models. The results obtained for the standard and RNG models show the same trend.

5. Conclusions

The effect of the turbulence model (i.e., STD, RNG, and RKE K-epsilon models) on the prediction of the velocity field of the real rotary kiln case was investigated.

Six different cases were simulated and predicted results were verified with available literature study. Results of the size and the position of the recirculation zone closer to the wall showed the sensitivity of the velocity field to the particle-flow interaction, especially in the RKE cases. A longer recirculation circle at around 9.7–14 m was remarked in RKE case and the recirculation zone given by the STD was around 4.5 m. Nevertheless, the STD and RNG models gave a good estimation of the flow field and recirculation zone size and length, while the RKE model performed best.