Prediction of Pulverized Coal Combustion Behavior around Tuyere by Using LES and Extended CPD Model

Abstract

Influence of pulverized coal into the blast furnace is an invaluable technology for reducing ironmaking costs. In order to increase the rate of pulverized coal injection, it is necessary to investigate the combustion behavior of pulverized coal around the tuyere, and many numerical simulations as well as experimental studies have examined this issue. In this paper, the combustion behavior of individual pulverized coal particles was investigated by using the Extended CPD model, which is a new devolatilization model, and the Large Eddy Simulation (LES) model. The volume fraction of CO+CO₂, which is the indicator of how much carbon included in pulverized coal is changed to the gas phase, was reasonably predicted by the new simulation model. The results of this study suggested that the devolatilization processes of pulverized coals with sizes under 40 µm varies widely between individual particles.

1. Introduction

Pulverized coal injection (PCI) into the blast furnace has become a key technology in the blast furnace process from the viewpoint of prolonging the life of aging coke ovens and improving the cost-effectiveness of blast furnace operation by reducing the coke ratio, and injection of large amounts of pulverized coal into blast furnaces is a global trend, including both Japan and other countries.^1,2) When the pulverized coal ratio is 100 kg/t-pig, the influence on pressure loss in the furnace is small.^3,4) However, coke powder in the tuyere increases as the amount of pulverized coal injected is increased to 150 kg/t-pig and 200 kg/t-pig, and this deteriorates the permeability of the lower part of the blast furnace.⁵⁾ Because improvement of the combustibility of pulverized coal is an effective countermeasure for deterioration of furnace permeability due to powdering of coke,⁶⁾ a lance capable of achieving high combustion efficiency has been developed,²⁾ and mixing different types of pulverized coal is performed to improve combustibility.⁷⁾ In those studies, numerical simulation of pulverized coal combustion simultaneously with pulverized coal combustion experiments has been employed as a useful tool for grasping the behavior of pulverized coal combustion.

The challenge of numerical simulation of the combustion field in pulverized coal injection into the blast furnace has a long history, and studies by one-dimensional (1D) simulation were done in the 1980s.⁸⁾ Simulations expanded to 2D in 1990s, and efforts were made to elucidate the flow, heat transfer and combustion reaction around the tuyeres, including the raceway.^9,10) Recently, 3D simulation of the combustion field in blast furnace pulverized coal injection has been reported.^11,12,13) However, in those studies, the time averaged Navier-Stokes simulation, or so-called Reynolds averaged Navier-Stokes (RANS) simulation, is used. Although the calculation cost of a RANS simulation in which any scale eddy is modeled by a turbulence model such as the k-ε model is reasonable, the RANS simulation method is lacking in reproducibility of combustion flows which are essentially unsteady phenomena, such as flame fluctuation, and it is difficult to understand the structure of the combustion flow and detailed combustion characteristics of individual pulverized coal particles. In addition, a simple model such as a one-step model¹⁴⁾ or two-step model¹⁵⁾ is often applied to modelling of devolatilization accompanying coal pyrolysis, which is one of the important factors determining the characteristics of the combustion flows. These models are considered to have low applicability to rapid heating conditions such as pulverized coal injection into blast furnaces because they employ hypothetical chemical species to match the calorific value and elemental composition of coal, and ignore the particle heating rate dependence of the devolatilization rate, the yield of volatile matter and the volatile composition, which are very important physical quantities.

On the other hand, the Large Eddy Simulation (LES) method, which can accurately reproduce turbulence flow fields, has recently been applied to the simulation of the pulverized coal combustion flow, and in addition to the temperature distribution and the major chemical species, prediction technology for more detailed combustion characteristics such as the NO concentration and particle dispersion behavior has been developed.^16,17,18,19)

Devolatilization models such as the chemical percolation devolatilization (CPD) model²⁰⁾ and FLACHCHAIN model,²¹⁾ which predict the devolatilization rate and volatile matter composition according to the thermal history of pulverized coal particles based on a chemical network structure model constructed from information such as the ¹³C NMR analysis values of coal, ultimate analysis values and proximate analysis values, have been proposed, and attempts to apply the proposed models have focused mainly on simulation of pulverized coal boilers for thermal power generation. However, these models are insufficient in terms of affinity with numerical simulation because polycyclic aromatic compounds such as the tar content are replaced with virtual chemical species.²²⁾ Very recently, Umemoto et al.²³⁾ proposed the extended CPD (Ex-CPD) model, which concretely describes the devolatilization process of polycyclic aromatic compounds in the initial pyrolysis stage by adding rapid pyrolysis experimental data obtained by a Curie point pyrolizer to the CPD model, but there have been no reports which verified its applicability to high-speed and high-temperature conditions such as pulverized coal injection into blast furnaces.

Pulverized coal (PC) injected from lances enters the raceway through the tuyere. Since the depth of the raceway is about 1 m to 1.5 m,²⁵⁾ assuming the velocity of the hot blast around the tuyere is 150 m/s, the residence time of the pulverized coal in the raceway can be estimated to be about 10 ms. In order to accurately grasp the combustion flow at such a minute time scale, it is necessary to accurately reproduce unsteadiness, including the ignition behavior of the pulverized coal near the lance exit and the dispersion of the coal particles. The LES method is suitable in this case. In the Ex-CPD model proposed by Umemoto et al., the gas composition of the volatiles released from the initial pyrolysis process can be concretely considered by using CO, CO₂, CH₄ and H₂O as the composition of light gas and benzene, naphthalene and phenanthrene as the composition of tar, which is unable to be determined in the conventional CPD model. Therefore, under blast furnace PC injection conditions, where the accuracy of the initial combustion behavior of the pulverized coal is important, the Ex-CPD model can reproduce the pulverized coal combustion phenomena with higher accuracy than the conventional CPD model.

In order to verify the applicability of the LES and Ex-CPD models to the injection of pulverized coal into blast furnaces, we first compared the results of numerical simulations using the LES and Ex-CPD models under the simplest pulverized coal injection conditions and the results of combustion experiments using a coke-packed bed type small combustion furnace. The devolatilization behavior of each pulverized coal particle was also examined based on the simulation results.

2. Numerical Simulation

2.1. Governing Equations of Gas Phase and Analytical Method

In this study, we used the finite volume code for unstructured grid systems FrontFlowRed-Comb (NuFD/FFR extended by Kyushu University, CRIEPI, Kyoto University and NuFD). Turbulence is also considered by LES. LES takes advantage of the difference in eddy characteristics due to the eddy scale to directly calculate larger eddies depending on the system as the GS (Grid scale) component, while using small eddies (SGS (Subgrid scale)), which show universal behavior.²⁴⁾ The unsteady flow can be reproduced accurately by considering the SGS component in the turbulent flow model. This is made possible by applying a filter of approximately the same size as the grids to the governing equation of the gas phase. When a phenomenon with density change such as a turbulent combustion field is treated by numerical simulation, Favre averaging is used. Favre averaging is defined as:   

$$\tilde{\varphi }\equiv \frac{\overline{\rho \varphi }}{\overline{\rho }}$$(1)

  

$$\varphi =\tilde{\varphi }+\overline{\varphi }$$(2)

where, ̄ is the spatial filtering operation and  ͂ is the Favre averaging operation. The LES equation of the gas phase is composed of the Favre averaged continuity equation, Navier-Stokes equation, energy conservation equation and the conservation equations of each species.²⁶⁾   

$$\frac{\partial \overline{\rho }}{\partial t}+\frac{\partial }{\partial x_j}\left(\overline{\rho }{\tilde{u}}_j\right)=S_m$$(3)

  

$$\frac{\partial \left(\overline{\rho }{\tilde{u}}_j\right)}{\partial t}+\frac{\partial }{\partial x_j}\left(\rho {\tilde{u}}_i{\tilde{u}}_j+\overline{p}{\delta }_{ij}-{\tilde{\tau }}_{ij}-{\tau }_{u_i{u}_j}\right)-\overline{\rho }{g}_i={\tilde{S}}_{ui}$$(4)

  

$$\begin{array}{l}{\displaystyle \frac{\partial \left(\overline{\rho }\tilde{h}\right)}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_j}}[\overline{\rho }{\tilde{u}}_j\tilde{h}-\overline{\rho }a{\displaystyle \frac{\partial \tilde{h}}{\partial x_j}}+\sum _{k=1}^n{\tilde{h}}_k\left(\overline{\rho }a-\overline{\rho }{D}_k\right){\displaystyle \frac{\partial {\tilde{Y}}_k}{\partial x_j}}-{\tau }_{u_{j}h}]\\ =S_h\end{array}$$(5)

  

$$\frac{\partial \left(\overline{\rho }{\tilde{Y}}_k\right)}{\partial t}+\frac{\partial }{\partial x_j}\left(\overline{\rho }{\tilde{u}}_j{\tilde{Y}}_k-\overline{\rho }{D}_k\frac{\partial {\tilde{Y}}_k}{\partial x_j}-{\tau }_{u_{j}h}{Y}_k\right)=\overline{{\dot{\omega }}_k}+{\overline{S}}_{Y_k}$$(6)

The terms ${\tau }_{u_i{u}_j}$, ${\tau }_{u_{i}h}$, ${\tau }_{u_i{Y}_k}$ in Eqs. (4), (5), (6) are the SGS terms that emerge by the filtering and are modeled by the Dynamic Smagorinsky model.^27,28) The SMAC method²⁹⁾ is employed for pressure and velocity correction. The discretized scheme of the convective term in Eq. (4) is a 9:1 blend of the second order central differencing scheme and the first order upwind differencing scheme. The first order Euler implicit method is used for the time integration algorithm. The time step is dt = 1.6 μs, and physical quantities averaged between 7500 and 5000 steps from the start of the simulation are calculated.

2.2. Governing Equations of Particles and Analytical Method

Pulverized coal particles are tracked individually in a Lagrangian manner. The PSI-Cell method³⁰⁾ was used for the interaction between the gas phase and the PC particles. In the PSI-Cell method, the gas phase mass, momentum, substance and enthalpy are exchanged at the mass point where particles are present. This makes it possible to give the sum of the influence of the particles existing in the cell as a source term in the governing equations concerning preservation of the gas phase mass, momentum, substance and enthalpy. The preservation formula of the particle position (x_p,i), velocity (u_p,i), temperature (T_p), mass (m_p) is composed of the following equations.   

$$\frac{d{x}_{p,i}}{dt}=u_{p,i}$$(7)

  

$$m_p\frac{d{u}_{p,i}}{dt}=\frac{1}{2}{C}_d{A}_p{Re}_p\left|u_i-u_{p,i}\right|\left(u_i-u_{p,i}\right)$$(8)

  

$$C_d=\frac{24\left(1+0.15{Re}_p{}^{0.687}\right)}{{Re}_p}$$(9)

  

$${Re}_p=\frac{D_p\left|u_i-u_{p,i}\right|}{\nu }$$(10)

  

$$\frac{d{T}_p}{dt}=-\frac{1}{c_{p,p}{m}_p}\left[A_p{q}_p{}^{″}-A_p{\epsilon }_p\sigma \left({\mathrm{\Theta }}_R^4-T_p^4\right)-Q_p-L_v\frac{d{m}_p}{dt}\right]$$(11)

  

$$\frac{d{m}_p}{dt}=\frac{dV}{dt}+\frac{dC}{dt}$$(12)

where, ${\mathrm{\Theta }}_R^4$ is (I/4σ)^1/4. The radiation intensity I was obtained by the discrete ordinate method.³¹⁾ dV/dt and dC/dt are the release rate of volatile matter in pulverized coal and the reaction rate of fixed carbon, respectively. The influence of turbulence on particle motion is considered by the Dynamic random walk SGS model.³²⁾ The Dynamic random walk SGS model is a model for representing the SGS component (u_i′) of the gas phase velocity ($u_i={\overline{u}}_i+u_i\prime$) in Eq. (8). In the Dynamic random walk SGS model, u_i′ is given by the normal distribution of average “0” and standard deviation “$\sqrt{2/3k_s}$”. Since the Dynamic random walk SGS model is a model that considers the SGS component, high precision of particle dispersion can be expected. In previous studies^16,17,18,19) and others, the authors verified the accuracy of the particle motion model by comparison with experiment results and attempted to elucidate the phenomena of particle dispersion behavior. Furthermore, Zhang et al.³³⁾ investigated the influence of particle size on particle dispersion behavior, and their simulation using the same particle motion model as in this study showed that the dispersion behavior of each particle size was in good agreement with the experimental results.

2.3. Pulverized Coal Combustion Models

This section describes the simulation models for pulverized coal combustion.

2.3.1. Devolatilization Model

Pulverized coal particles injected from the lance into the blowpipe are rapidly heated by heat transfer and radiation from the hot blast and furnace wall, and release volatile matter mainly containing hydrocarbons. In this study, we use the Ex-CPD model and one-step model. The original CPD model is one of the network type primary pyrolysis models that consider the chemical structure of coal. The chemical structure of coal is assumed to be a structure consisting of an aromatic nucleus cluster, which is the smallest unit constituting coal, and a bond linking it. The initial coal structure is determined from only ¹³C NMR, ultimate analysis and proximate analysis. In the CPD model, the yields of char, light gas and tar can be obtained by calculating the number of bonds that are cleaved and decrease by pyrolysis from the kinetics. However, because this model considers only two types of bonds, namely, degradable bonds and non-degradable bonds, and only one type of functional group, the composition of the light gases and tar produced is unknown. The Ex-CPD model newly considers two degradable bonds, one non-degradable bond and four functional groups, and the proportion of 1 to 3 rings present in the coal is decided from the results of a rapid thermal cracking test and the average molecular weight of the aromatic nucleus cluster obtained from the NMR data. The light gases released in the pyrolysis process expressed as CO, CO₂, CH₄, and H₂O as well as tar released in the pyrolysis process are expressed concretely as benzene, naphthalene and phenanthrene. As a result, the Ex-CPD model has high applicability to numerical analysis. On the other hand, in the one-step model, devolatilization is simplified as follows.   

$$Coal\to VM+Char$$(13)

The volatile matter release rate is given by the following equations.   

$$\frac{dV}{dt}=k_v\left(V^{*}-V\right)$$(14)

  

$$k_v=A_{v}exp\left(-\frac{E_v}{R{T}_p}\right)$$(15)

where, V* is the mass fraction of the volatile matter contained in the pulverized coal, and V is the mass fraction of the volatile matter released in the gas phase at the present time in V*. Also, A_v is the frequency factor, E_v is the activation energy, R is the gas constant and T_p is the particle temperature. It is known that the mass fraction of volatile matter (V*) in pulverized coal varies with the heating history. V* and volatile content V_pro obtained by the proximate analysis are expressed by the following formula using the Q factor (Q_fac).   

$$V^{*}=Q_{fac}\times V_{pro}$$(16)

In this study, the particle diameter of pulverized coal particles in the release process of volatile matter is constant.

2.3.2. Char Gasification

Pulverized coal particles become char mainly containing carbon when the release of volatile matter is completed. CO and H₂ are generated by the gasification reaction of fixed carbon (FC) in char. The gasification reaction of char assumed one exothermic reaction (Eq. (17)) and two endothermic reactions (Eqs. (18) and (19)). The gasification reaction by O₂ is faster than the gasification reaction by CO₂ and H₂O, but in pulverized coal injection into the blast furnace, the amount of oxygen in the hot blast is less than the amount of oxygen required for complete combustion of the pulverized coal. Under such conditions, gasification by CO₂ and gasification by H₂O also have a large influence on the combustion field.   

$$C+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.O_2\stackrel{k_1}{⟶}CO$$(17)

  

$$C+C{O}_2\stackrel{k_2}{⟶}2CO$$(18)

  

$$C+H_{2}O\stackrel{k_3}{⟶}CO+H_2$$(19)

The random pore model proposed by Bhatia et al.³⁴⁾ was utilized in calculating the char gasification rate with a gasifying agent.   

$$\frac{d\gamma }{dt}=A_{0i}{P}_{Ai}^{n}exp\left(-E/R{T}_p\right)\left(1-\gamma \right)\sqrt{1-{\mathrm{\Psi }}_{i}ln\left(1-\gamma \right)}$$(20)

  

$$\gamma =1-M_{C_i}/M_{C_0}$$(21)

where, Ψ_i is the pore structure coefficient, P_Ai is the partial pressure of the gasifying agent, γ is the conversion rate of char, $M_{C_i}$ is the fixed carbon amount in the simulation and $M_{C_0}$ is the initial fixed carbon amount. In this study, we used the reaction rate parameter of DT coal,²³⁾ whose composition and element ratio are close. In this study, the particle diameter of the pulverized coal particles in the gasification process of char was assumed to be constant.

2.4. Gaseous Phase Reactions

Gaseous combustion between the volatile matter and hot blast is calculated by the following two-step global reactions.   

$$C{H}_4+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.O_2\stackrel{k_4}{⟶}CO+2H_2$$(22)

  

$$H_2+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.O_2\stackrel{k_5}{⟶}2H_{2}O$$(23)

  

$$CO+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.O_2\stackrel{k_6}{⟶}C{O}_2$$(24)

  

$$CO+H_{2}O\stackrel{k_7}{⟷}C{O}_2+H_2$$(25)

  

$$C_6{H}_6+3O_2\stackrel{k_8}{⟶}6CO+3H_2$$(26)

  

$$C_{10}{H}_8+5O_2\stackrel{k_9}{⟶}10CO+4H_2$$(27)

  

$$C_{14}{H}_{10}+7O_2\stackrel{k_{10}}{⟶}14CO+5H_2$$(28)

  

$$C_{16}{H}_{10}+8O_2\stackrel{k_{11}}{⟶}16CO+5H_2$$(29)

The scale similarity filtered reaction rate model (SSFRRM), which can consider the influence of the SGS component of LES in the turbulent combustion model,³⁴⁾ was adopted. The SSFRRM model is a model which was developed for combustion simulations using LES. The reaction rate filtered by the spatial filter necessary for the combustion simulation by LES is directly calculated from the GS component of the temperature and the chemical species concentration by the scale similarity rule. In the SSFRRM model, ${\overline{\dot{\omega }}}_F$ is expressed as the follow equation.   

$$\begin{array}{l}{\overline{\dot{\omega }}}_F=\overline{\dot{\omega }\left(\rho ,Y_k,T\right)}\\ \hspace{1em} =\overline{\dot{\omega }\left(\overline{\rho },{\tilde{Y}}_k,\tilde{T}\right)+K_1\left(\overline{\overline{\dot{\omega }\left(\overline{\rho },{\tilde{Y}}_k,\tilde{T}\right)}}-\overline{\overline{\dot{\omega }\left(\overline{\overline{\rho }},{\tilde{\tilde{Y}}}_k,\tilde{\tilde{T}}\right)}}\right)}\end{array}$$(30)

where, K₁ is a model constant and K₁ = 1.0, which is the recommended value, is adopted. $\dot{\omega }$(ρ, Y_k, T) is expressed by the following Arrhenius type equation.   

$$\dot{\omega }\left(\rho ,Y_k,T\right)=A_k{T}_g{}^{n}exp\left(-\frac{E_k}{RT}\right)\cdot {\left[F\right]}^d{\left[O\right]}^e$$(31)

The values of the parameters in Eq. (31) are listed in Table 1. The parameter of the aromatic used in this study is the same as that of benzene.

Table 1. Kinetic parameters for gas reactions.

	Ak	n	Ek [J/mol]	Reference
k4	3.0×108	0	1.26×105	32)
k5	6.8×1015	−1	1.68×105	32)
k6	2.2×1012	0	1.67×105	33)
k7	Functional Chemistry			34)
k9	2.4×1011	0	1.26×105	33)
k10	2.4×1011	0	1.26×105	33)
k11	2.4×1011	0	1.26×105	33)
k12	2.4×1011	0	1.26×105	33)

Equation (25) is a reaction called the water-gas-shift reaction. In the gasification reaction of pulverized coal, H₂O and CO₂ play important roles, as shown in the Eqs. (18) and (19). Therefore, it is important to accurately ascertain the amounts of H₂O and CO₂ contained in the gaseous phase. However, general-purpose values of the parameters in the overall reaction of the water-gas-shift reaction have not yet been obtained. Therefore, in this study, the functional chemistry model proposed by Watanabe et al.³⁸⁾ was adopted in order to obtain an appropriate rate constant of water-gas-shift reactions. In the functional chemistry model, a table of reaction rate parameters is prepared in advance by carrying out multiple elementary reaction analyses under the gas temperature and chemical species concentration of the assumed flow field. In simulation of a pulverized coal combustion field, reaction rate parameters corresponding to the temperature and chemical species concentration are determined for each calculation grid from the table. CHEMKIN-PRO³⁹⁾ for elementary reaction analysis was used. In CHEMKIN-PRO, the dependence between chemical reactions can be clarified by solving thousands of reaction combinations (reaction mechanisms) and by simulating the composition change and the temperature change with the progress of the reaction. For the detailed reaction mechanism, GRI-Mech 3.0⁴⁰⁾ consisting of 53 chemical species and 325 elementary reactions was used. A complete mixing reactor model was used as the reaction model. By adopting this method, it became possible to reproduce the water-gas-shift reaction with a low computational load and high precision.

3. Experimental and Numerical Setup

3.1. Experimental Equipment and Numerical Setup

The experimental apparatus used in this study is shown in Fig. 1. The equipment consists of a coke packed bed, tuyere and blowpipe. The coke packed bed part has a width of 400 mm, depth of 1000 mm and height of 1400 mm. The inner diameter of the blowpipe portion is Φ 90 mm, the tuyere diameter is Φ 65 mm and the blowpipe has eight observation windows. Cooling N₂ gas was injected from the observation windows at a rate of 8 Nm³/h. High temperature hot blast was produced by an LPG burner and then mixed with oxygen, adjusted to the composition in Table 2, and sent to the blowpipe section. From the lance, pulverized coal carried by the N₂ carrier gas is blown into and mixed with the hot blast in the blowpipe. The injected pulverized coal is heated rapidly by heat transfer and radiation from the coke in the raceway and the furnace wall, flame and hot blast. The pulverized coal moves through the blowpipe while burning and enters the raceway from the tuyere. To measure the gas composition, sampling probes were inserted at positions 37 mm and 300 mm from the lance tip, and the gas was sampled in the radial direction (horizontal direction).

Fig. 1.

Experimental apparatus.

Table 2. Hot blast conditions.

Volume [Nm3/h]	Temperature [K]	Gas composition [vol%]
269.4	1723	O2	N2	CO2	H2O
		32.1	50.5	7.46	9.95

The analysis geometry and boundary conditions of the simulation in this study are shown in Fig. 2. The coordinate system in the subsequent description follows Fig. 2. The blowpipe, lance, observation window and tuyere reproduced the shape of the experimental apparatus. Since the evaluation point of the pulverized coal combustion behavior in this paper is inside the blowpipe, the raceway has been simplified. The number of calculation grids was about 3.2 million, and they were arranged so that the lattice density becomes higher near the lance tip. A temperature boundary condition of 1923 K was set on the wall of the downstream side to consider the effect of the coke in the raceway.

Fig. 2.

Boundary conditions.

3.2. Coal Properties

Table 3 shows the properties of the coal used in the study. The object coal in this study has a ratio (fuel ratio) of volatile matter (VM) to fixed carbon (FC) of 1.61. The particle size distribution of the coal is shown in Fig. 3. In this simulation, the cumulative mass fraction (η) measured in the experiment is approximated by the following Rosin-Rammler distribution equation. For the simulation, pulverized coal particles are charged at every time step according to the distribution of Eq. (32) in the range of 5 to 200 μm. The density of the pulverized coal was 1000 kg/m³.   

$$\eta =100\times \left[1-exp\left(-{\left(\frac{D_p}{93.1}\right)}^{1.08}\right)\right]$$(32)

Table 3. Coal properties.

Volatile matter [wt%]	33.5
Fixed carbon [wt%]	51.8
Ash [wt%]	12.4
C [wt%, d.a.f]	82.6
H [wt%, d.a.f]	5.24
N [wt%, d.a.f]	2.19
S [wt%, d.a.f]	0.49
O [wt%, d.a.f]	9.44
LHV [kJ/kg]	28670

Fig. 3.

Particle diameter distribution.

3.3. Devolatilization Model Parameters

In this study, numerical simulations were carried out using the Ex-CPD model (Case 1) and the 1-step model (Case 2). The average molecular weight of the aromatic nucleus cluster necessary for setting up the Ex-CPD model was obtained from the result of ¹³C NMR. In addition, the parameters of the DT coal,²³⁾ whose composition and element ratio are close to those of one ring, two rings and three rings in the tar component, were used. Table 4 shows the 4 devolatilization model parameters used in the one-step model. For Case 2, the following equation was used to estimate the devolatilization parameters.⁴¹⁾ The Q factor was 1.2.   

$$A_v={10}^{22.1-0.067\times C\%}$$(33)

  

$$E_v=23.7\times {10}^{-3}exp\left(\frac{C\%}{6.35}\right)+232\mathrm{\hspace{0.17em} }907$$(34)

Table 4. Devolatilization model parameters.

	Model	Av [1/s]	Ev [J/mol]	Qfac	Ref
Case 1	Ex-CPD	–	–	–
Case 2	1-step	1.675×1017	2.35×105	1.2	37)

4. Results

Figure 4 shows the volume fraction of CO + CO₂ and N₂ in the z direction at the position 300 mm from the lance tip, together with the experimental results and the simulation results for each case. Here, r/R is a dimensionless radius, which is normalized by dividing the measurement position by the blowpipe radius. r/R = 0 represents the center portion of the blowpipe, and r/R = 1 represents the wall surface. In both Case 1 and Case 2, the volume fraction of CO + CO₂, which is an indicator of the conversion rate from the carbon contained in the pulverized coal to the gas phase, and the volume fraction of N₂, which is not involved in the reaction, were in good agreement with the experimental results in the radial direction.

Fig. 4.

Volume fraction of CO+CO₂ at 300 mm from lance tip.

Figure 5 shows the volume fraction of CO + CO₂ on the central axis of the blowpipe against the distance from the lance tip for each case and the experiment. Figure 6 shows the volume fraction of CO + CO₂ in the xy plane at z = 0. From Fig. 5, both Case 1 and Case 2 show good agreement with the experimental results at the measurement position. However, the volume fraction of CO + CO₂ shows different values in Case 1 and Case 2 except at the measurement points of the experiment. From Figs. 5 and 6, it can be seen that Case 1 using the Ex-CPD model has a higher initial devolatilization rate. In order to compare the devolatilization rate parameters between the Ex-CPD model and the one-step model, the A_v, E_v and Q factors in the one-step model were estimated from the simulation results of the Ex-CPD model. The results are shown in Table 5. The method is given in Appendix 1.

Fig. 5.

Volume fraction in axial direction.

Fig. 6.

Volume fraction of CO+CO₂.

Table 5. Estimated devolatilization model parameters.

	Model	Av [1/s]	Ev [J/mol]	Qfac	Ref
Case 1	Ex-CPD	4.639×106	5.17×105	1.17
Case 2	1-step	1.675×1017	2.35×105	1.2	37)

Figure 7 shows K_v in Eq. (15) for the particle temperature using the obtained A_v and E_v. In comparison with Case 1 using the Ex-CPD model, Case 2 using the one-step model has higher sensitivity to K_v, that is, the devolatilization rate, than the particle temperature. In general, it is difficult to make general use of the parameters of simple models such as one-step models. Therefore, the model developer must decide the parameter with which the combustion field to be analyzed can be reproduced each time. In this study, Case 2 using the one-step model agreed with the experimental results as well as Case 1 using the Ex-CPD model. However, as shown in Fig. 7, when the particle temperature changes, K_v of Case 2 changes by about 15 digits. Thus, unless the application range of the model parameters is carefully assessed before the parameters are used, there is a risk of immediately over-evaluating or under-evaluating. Furthermore, it is difficult to determine the parameters of the one-step model in advance when changing the coal type, reducing atmosphere, oxidizing atmosphere, CO₂ partial pressure and H₂O partial pressure. On the other hand, since Ex-CPD is based on the chemical structure, the devolatilization parameters can be clearly determined chemically. Therefore, the Ex-CPD model can be applied to various combustion fields. In the future, the applicability of the Ex-CPD model to blast furnace pulverized coal injection conditions will be verified by comparing the experimental results under various conditions with the simulation results using the Ex-CPD model.

Fig. 7.

Devolatilization parameter (K_v) for each case.

5. Discussion

A major feature of the Ex-CPD model is that the devolatilization rate and Q factor are determined depending on the particle heating rate. Therefore, the devolatilization behavior of individual particles in Case 1 was examined. Figure 8 shows the relationship between K_v and the temperature rise rate for each particle size assuming the particle temperature of 1000 K. The heating rate of the particles was determined by the least squares method from the times at which the particle temperature reached 600 K, 800 K, 1000 K, and 1200 K. From Fig. 8, K_v and the particle heating rate are positively correlated. In addition, K_v of the small particles has a large range due to the large dispersion of the temperature heating rate of each particle. In the following, particles with a diameter of not more than 40 μm are defined as small-diameter particles.

Fig. 8.

Effect of coal particle heating ratio on devolatilization parameter (K_v).

Figure 9 shows the velocity vector of the plane perpendicular to the y-axis at the lance tip. N₂, which is the carrier gas injected from the center of the lance, forms a flow perpendicular to the lance exit surface, while the flow at the lance wall is drawn into the flow of surrounding hot blast with large momentum by the Coanda effect. The volume fraction of N₂ and the position of the pulverized coal particles in the cross section are shown in Fig. 10. The distribution of the pulverized coal particles in the yz plane at 100 mm from the lance tip is shown in Fig. 11. The small-diameter particles injected from the center of the lance are heated while surrounded by the N₂ carrier gas, and the rise in temperature becomes slow. On the other hand, the small-diameter particles discharged from near the lance wall are drawn into the hot blast together with the N₂ carrier gas, and the temperature rises immediately. As shown in Fig. 11, the positions of the small-diameter particles vary in the z-direction, and this variation is considered to be the cause of the variation in the heating rate of the small particles. In contrast, as shown in Fig. 10, large particles move in the direction perpendicular to the lance exit surface due to the large inertial force of the particle itself, separate from the N₂ carrier gas, and are then heated by the hot blast. As shown in Fig. 11, there is no scattering of the particle positions of the larger diameter particles. Therefore, the heating rate does not change greatly from particle to particle, and the devolatilization rate also differs little.

Fig. 9.

Gas velocity vectors around lance tip.

Fig. 10.

N₂ volume fraction distribution and particle distribution with particle diameter.

Fig. 11.

Particle distribution at 100 mm from lance tip.

The simulation using the Ex-CPD model suggested that the combustion behavior of each particle in the pulverized coal combustion flow was different for each particle size, and was also different for each particle in the case of small-diameter particles. At present, a lance for increasing the oxygen concentration around pulverized coal particles by introducing oxygen from an outer pipe using the pulverized coal injection lance as a double tube,⁴²⁾ and a lance for rapidly heating pulverized coal by blowing a readily-combustible gas reducing agent instead of oxygen into the outer tube¹³⁾ have been studied and put into practical use. In these complicated pulverized coal combustion fields, the thermal history of the particles seems to differ greatly depending on the combustion conditions. Because the oxygen concentration around the pulverized coal particles and the temperature of the ambient gas are greatly different for each particle, it is anticipated that the thermal history of each particle will also be different. It is difficult to determine the general devolatilization parameters of various pulverized coal combustion flows by a one-step model without considering the effect of the thermal history. However, the Ex-CPD model, which can determine the devolatilization rate variably for each thermal history of particles, may be applicable to combustion flows where the thermal history of each particle is significantly different simply by using parameters which are specific to pulverized coal and independent of its thermal history. In the future, we will compare the experimental results with the Ex-CPD model under the various conditions described above in order to verify the applicability of the Ex-CPD model to pulverized coal injection into blast furnaces.

6. Conclusion

In this study, in order to verify the applicability of the LES and Ex-CPD models to pulverized coal injection into a blast furnace, we first compared the results of a numerical simulation using a combined LES and Ex-CPD model under the simplest pulverized coal injection condition and the results of a combustion experiment using a coke-packed bed type small combustion furnace. As a result, the following findings were obtained.

(1) It was confirmed that the numerical simulation using the LES and Ex-CPD model could reproduce the experimental results for the sum of the volume fraction distribution of CO₂ and CO, which shows the gasification rate of carbon in pulverized coal, and the volume fraction distribution of N₂, which is not involved in the reaction.

(2) The simulation using the Ex-CPD model showed that the devolatilization rate of pulverized coal differs depending on the particle size, and also differs from particle to particle in the case of pulverized coal particles smaller than 40 μm.

(3) The Ex-CPD model can be used for various blast furnace conditions without changing parameters.

Symbols

A_p: surface area of the particle [m²]

A_v: pre-exponential factor for the volatile matter evolution rate equation [1/s]

A_k: frequency factor

A₀: pre-exponential factor for char gasification reaction [1/s]

a: thermal diffusivity [m²/s]

C: char mass [kg]

C_d: drag coefficient [-]

c_p,p: specific heat of the particle [J/(kg·K)]

D: diffusion coefficient of the chemical species [m²/s]

D_p: particle diameter [m]

E: activation energy of char gasification reaction [J/mol]

E_k: activation energy for Eq. (31) [J/mol]

E_v: activation energy for the volatile matter evolution rate equation [J/mol]

F: fuel concentration [mol/m³]

g: gravitational acceleration [m/s²]

h: enthalpy [J/kg]

I: radiation intensity [J/(m²·s)]

K1: model constant [-]

k_s: kinetic energy of the SGS component [m²/s²]

L_v: latent heat of vaporization [J/kg]

M_Ci: amount of fixed carbon in the simulation [kg]

M_C0: initial fixed carbon amount [kg]

m_p: particle mass [kg]

O: oxidizer concentration [mol/m³]

P_Ai: partial pressure of gasifier [Pa]

Q_p: particle source term [J/s]

Q_fac: Q factor [-]

q_p″: heat flow rate at the particle surface [J/(m²·s)]

R: gas constant of the ideal gas [J/(mol·K)

Re: Reynolds number [-]

Re_p: particle Reynolds number [-]

S: source term

T_g: gas temperature [K]

T_p: particle temperature [K]

t: time [s]

u: velocity [m/s]

u_p: particle velocity [m/s]

V: mass of volatile matter evolved from a coal particle [kg]

V^*: mass of volatile matter in coal particle in the simulation [kg]

x: spatial coordinates

x_p: particle position coordinates

Y: mass fraction of the species [-]

Greeks

γ: conversion of char [-]

δ: Kronecker Delta

ε_p: emissivity [-]

η: cumulative mass fraction [-]

Θ_R: radiation temperature [K]

μ: viscosity coefficient [Pa·s]

ρ: density [kg/m³]

σ: Stefan-Boltzmann constant [J/(m²·s·k⁴)]

τ: stress tensor [Pa]

ϕ: physical quantity

Ψ_i: structure parameter in random pore model [-].

$\dot{\omega }$: reaction rate [mol/(s·m³)]

Operators

$\overline{•}$: spatial filter

$\widetilde{•}$: Favre filter

Appendix 1

In this analysis, the devolatilization and the surface reaction of char occurred simultaneously. To investigate only the devolatilization process, a simulation for Case 1 without the surface reaction formula for char shown in Eqs. (17), (18) and (19) was conducted, and the time variation of 80 particles of pulverized coal was monitored. Figure A1 shows an example of the relationship between time and the reciprocal of the particle temperature. The relationship between the time of the sampled particles and the reciprocal of the particle temperature is approximated by a convex downward quadratic function or a linear function with a monotonic decrease. First, the relationship is approximated by a quadratic function, and if α₂ < 0 in Eq. (A1), it is approximated by the linear function with the monotonic decrease (Eq. (A2)). As a result, 40 particles are calculated by Eq. (A1) and 40 particles are calculated by Eq. (A2). The coefficients in Eqs. (A1) and (A2) are determined by the least squares method.   

$$\frac{1}{T_p}={\alpha }_0+{\alpha }_{1}t+{\alpha }_2{t}^2\mathrm{\hspace{0.17em} } \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{\hspace{0.17em} } {\alpha }_0<0,\mathrm{\hspace{0.17em} } {\alpha }_1<0,\mathrm{\hspace{0.17em} } {\alpha }_2>0$$(A1)

  

$$\frac{1}{T_p}={\alpha }_0+{\alpha }_{1}t\mathrm{\hspace{0.17em} } \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{\hspace{0.17em} } {\alpha }_0<0,\mathrm{\hspace{0.17em} } {\alpha }_1<0$$(A2)

Fig. A1.

Time history of reciprocal particle temperature.

The average value of the correlation coefficient of the 80 pulverized coal particles obtained from the approximation and the simulation result is 0.996.

When the relationship between time and the reciprocal of the particle temperature is approximated by the quadratic function, the devolatilization volatile release rate in Eq. (14) and Eq. (15) lead to Eq. (A3), and the analytical solution is given by Eq. (A4) or Eq. (A5).   

$$\begin{array}{l}{\displaystyle \frac{1}{\left(V^{*}-V\right)}}dV=\\ A_{v}exp\left\{{\displaystyle \frac{E_v}{R}}\left({\displaystyle \frac{{\alpha }_1{}^2}{4{\alpha }_2}}-{\alpha }_0\right)\right\}exp\left\{-{\displaystyle \frac{E_v{\alpha }_2}{R}}{\left(t+{\displaystyle \frac{{\alpha }_1}{2{\alpha }_2}}\right)}^2\right\}dt\end{array}$$(A3)

  

$$\begin{array}{l}-log\left(V^{*}-V\right)=\\ A_{v}exp\left\{{\displaystyle \frac{E_v}{R}}\left({\displaystyle \frac{{\alpha }_1{}^2}{4{\alpha }_2}}-{\alpha }_0\right)\right\}\sqrt{{\displaystyle \frac{\pi R}{4E_v{\alpha }_2}}}erf\left\{\sqrt{{\displaystyle \frac{E_v{\alpha }_2}{R}}}\left(t+{\displaystyle \frac{{\alpha }_1}{2{\alpha }_2}}\right)\right\}\\ +C_0\end{array}$$(A4)

  

$$\begin{array}{l}{V}^{*}-V=\\ C_{1}exp\left[-A_{v}exp\left\{{\displaystyle \frac{E_v}{R}}\left({\displaystyle \frac{{\alpha }_1{}^2}{4{\alpha }_2}}-{\alpha }_0\right)\right\}\sqrt{{\displaystyle \frac{\pi R}{4E_v{\alpha }_2}}}erf\left\{\sqrt{{\displaystyle \frac{E_v{\alpha }_2}{R}}}\left(t+{\displaystyle \frac{{\alpha }_1}{2{\alpha }_2}}\right)\right\}\right]\end{array}$$(A5)

where, C₁ = exp(−C₀) is the integral constant and erf is the error function. When the relationship between the reciprocal of the particle temperature and time is approximated by the linear function, the analytical solution is given by Eq. (A6) or Eq. (A7).   

$$log\left(V^{*}-V\right)=\frac{A_{v}R}{E_v{\alpha }_1}exp\left(-\frac{E_v\left({\alpha }_0+{\alpha }_{1}t\right)}{R}\right)+C_0$$(A6)

  

$$V^{*}-V=C_{1}exp\left[\frac{A_{v}R}{E_v{\alpha }_1}exp\left(-\frac{E_v\left({\alpha }_0+{\alpha }_{1}t\right)}{R}\right)\right]$$(A7)

where, C₁ = exp(−C₀) is the integral constant. A_v and E_v are obtained as follows. First, the integral constant C₁ is determined. Let t₀ be the time step at which devolatilization starts and t_end be the time step at which devolatilization ends, with n time steps between t₀ ≦ t ≦ t_end. Taking the sum of the two sides of Eq. (A5) for each time step, C₁ is Eq. (A8), and taking the sum of the two sides of Eq. (A7) for each time step, C₁ is Eq. (A9).   

$$\begin{array}{l}{C}_1=\\ {\displaystyle \frac{n{V}^{*}-\sum _{k=1}^n{V}_k}{exp\left[-A_{v}exp\left\{{\displaystyle \frac{E_v}{R}}\left({\displaystyle \frac{{\alpha }_1{}^2}{4{\alpha }_2}}-{\alpha }_0\right)\right\}\sqrt{{\displaystyle \frac{\pi R}{4E_v{\alpha }_2}}}\right]\sum _{k=1}^{n}exp\left[erf\left\{\sqrt{{\displaystyle \frac{E_v{\alpha }_2}{R}}}\left(t_k+{\displaystyle \frac{{\alpha }_1}{2{\alpha }_2}}\right)\right\}\right]}}\end{array}$$(A8)

  

$$C_1={\displaystyle \frac{n{V}^{*}-\sum _{k=1}^n{V}_k}{exp\left(-{\displaystyle \frac{A_{v}R}{E_v{\alpha }_1}}\right)\sum _{k=1}^{n}exp\left\{exp\left(-{\displaystyle \frac{E_v\left({\alpha }_0+{\alpha }_1{t}_k\right)}{R}}\right)\right\}}}$$(A9)

Next, A_v is calculated. If A_vk is constant between t_k≦t≦t_k+1, Eq. (A10) can be derived by taking the difference between the time steps in Eq. (A4), and Eq. (A6) can be derived by taking the difference between the time steps in Eq. (A11).   

$${\displaystyle \begin{array}{l}{A}_{vk}=\\ {\displaystyle \frac{log\left\{{\displaystyle \raisebox{1ex}{$\left(V*-V_k\right)$}\!\left/ \!\raisebox{-1ex}{$\left(V*-V_{k+1}\right)$}\right.}\right\}}{A_{v}exp\left\{{\displaystyle \frac{E_v}{R}}\left({\displaystyle \frac{{\alpha }_1{}^2}{4{\alpha }_2}}-{\alpha }_0\right)\right\}\sqrt{{\displaystyle \frac{\pi R}{4E_v{\alpha }_2}}}\left[erf\left\{\sqrt{{\displaystyle \frac{E_v{\alpha }_2}{R}}}\left(t_{k+1}+{\displaystyle \frac{{\alpha }_1}{2{\alpha }_2}}\right)\right\}-erf\left\{\sqrt{{\displaystyle \frac{E_v{\alpha }_2}{R}}}\left(t_k+{\displaystyle \frac{{\alpha }_1}{2{\alpha }_2}}\right)\right\}\right]}}\end{array}}$$(A10)

  

$$A_{vk}={\displaystyle \frac{log\left\{{\displaystyle \raisebox{1ex}{$\left(V*-V_k\right)$}\!\left/ \!\raisebox{-1ex}{$\left(V*-V_{k+1}\right)$}\right.}\right\}}{{\displaystyle \frac{R}{E_v{\alpha }_1}}\left\{exp\left(-{\displaystyle \frac{E_v\left({\alpha }_0+{\alpha }_1{t}_{k+1}\right)}{R}}\right)-exp\left(-{\displaystyle \frac{E_v\left({\alpha }_0+{\alpha }_1{t}_k\right)}{R}}\right)\right\}}}$$(A11)

A_v is determined by averaging n-1 of A_vk obtained from Eqs. (A10) and (A11).   

$$Av={\displaystyle \frac{\sum _{k=1}^{n-1}A{v}_k}{n-1}}$$(A12)

Finally, E_v is determined such that the amount of devolatilization (V′₀) is 0 at time t₀. Figure A2 shows an example of the time variation of the amount of devolatilization obtained above and the amount devolatilization in Case 1. The correlation coefficient between the time variation of the amount of devolatilization in Case 1 and the time variation of the amount of devolatilization obtained by the above method is 0.97. As shown by Eq. (16), the Q factor is given by the ratio of the total amount of devolatilization (V*) as shown in Fig. A2 and the amount of devolatilization (V_pro) obtained by the proximate analysis value.

Fig. A2.

Time history of released volatile matter.

References

• 1)   T.  Ariyama,  M.  Sato,  R.  Murai,  K.  Miyagawa,  K.  Nozawa and  T.  Kamijo: Tetsu-to-Hagané, 81 (1995), 1114 (in Japanese).
• 2)   M.  Sato,  R.  Murai,  T.  Ariyama,  A.  Maki,  A.  Shimomura and  K.  Mori: Tetsu-to-Hagané, 84 (1998), 37 (in Japanese).
• 3)   M.  Deguchi,  K.  Shibata and  M.  Fukuoka: Tetsu-to-Hagané, 76 (1990), 691 (in Japanese).
• 4)   Y.  Iwanaga: Tetsu-to-Hagané, 77 (1991), 71.
• 5)   K.  Yamaguchi,  T.  Uno,  T.  Yamamoto,  H.  Ueno,  N.  Kon-no and  S.  Matsuzaki: Tetsu-to-Hagané, 82 (1996), 641 (in Japanese).
• 6)   M.  Sato,  R.  Murai,  T.  Ariyama,  A.  Maki,  A.  Shimomura and  K.  Mori: Tetsu-to-Hagané, 85 (1999), 717 (in Japanese).
• 7)   K.  Kunitomo,  T.  Orimoto,  T.  Nishimura,  M.  Naito and  J.  Yagi: Tetsu-to-Hagané, 90 (2004), 190 (in Japanese).
• 8)   J.  He,  M.  Kuwabara and  I.  Muchi: Tetsu-to-Hagané, 72 (1986), 1847 (in Japanese).
• 9)   H.  Aoki,  H.  Nogami,  H.  Tsuge,  T.  Miura and  T.  Furukawa: ISIJ Int., 33 (1993), 646.
• 10)   K.  Takeda and  F. C.  Lockwood: ISIJ Int., 37 (1997), 432.
• 11)   R.  Murai,  M.  Sato and  T.  Ariyama: ISIJ Int., 44 (2004), 2168.
• 12)   Y.  Shen,  B.  Guo,  A.  Yu,  D.  Maldonado,  P.  Austin and  P.  Zulli: ISIJ Int., 48 (2008), 777.
• 13)   A.  Murao,  Y.  Kashihara,  K.  Takahashi,  N.  Oyama,  H.  Matsuno and  M.  Sato: Tetsu-to-Hagané, 101 (2015), 653.
• 14)   S.  Badzioch and  P. G. W.  Hawksley: Ind. Eng. Chem. Process Des. Dev., 9 (1970), 521.
• 15)   H.  Kobayashi,  J. B.  Howard and  A. F.  Sarofim: Proc. Combust. Inst., 16 (1977), 411.
• 16)   H.  Watanabe,  K.  Tanno,  Y.  Baba,  R.  Kurose and  S.  Komori: Turbulence, Heat and Mass Transfer 6, Begell House, Danbury, (2009), 1027.
• 17)   R.  Kurose,  H.  Watanabe and  H.  Makino: KONA Powder Part. J., 27 (2009), 144.
• 18)   M.  Muto,  H.  Watanabe,  R.  Kurose,  S.  Komori,  S.  Balusamy and  S.  Hochgreb: Fuel, 142 (2015), 152.
• 19)   S.  Ahn,  K.  Tanno and  H.  Watanabe: Appl. Therm. Eng., 124 (2017), 1194.
• 20)   T.  Fletcher,  A.  Kerstein,  R.  Pugmire and  D.  Grant: Energy Fuels, 4 (1990), 54.
• 21)   S.  Niksa: Combust. Flame, 100 (1995), 384.
• 22)   O.  Stein,  G.  Olenik,  A.  Kronenburg,  F.  Cavallo Marincola,  B. M.  Franchetti,  A. M.  Kempf,  M.  Ghiani,  M.  Vascellari and  C.  Hasse: Flow Turbul. Combust., 90 (2013), 859.
• 23)   S.  Umemoto,  S.  Kajitani,  K.  Miura,  H.  Watanabe and  M.  Kawase: Fuel Process. Technol., 159 (2017), 256.
• 24)   S.  Kida and  S.  Yanase: Turbulent Dynamics, Asakura Shoten, Tokyo, (1999), 93 (in Japanese).
• 25)   Y.  Matsui,  Y.  Yamaguchi,  M.  Sawayama,  S.  Kitano,  N.  Nagai and  T.  Imai: Tetsu-to-Hagané, 92 (2006), 919 (in Japanese).
• 26)   P.  Moin,  K.  Squires,  W.  Cabot and  S.  Lee: Phys. Fluids, 3 (1991), No. 11, 2746.
• 27)   M.  Germano,  U.  Piomelli,  P.  Moin and  W. H.  Cabot: Phys. Fluids A, 3 (1991), 1760.
• 28)   D. K.  Lilly: Phys. Fluids A, 4 (1992), 633.
• 29)   A. A.  Amsden and  F. H.  Harlow: J. Comput. Phys., 6 (1970), 322.
• 30)   C. T.  Crowe,  M. P.  Sharma and  D. E.  Stock: J. Fluids Eng., 99 (1977), 325.
• 31)   W. A.  Fiveland: J. Thermophys., 2 (1988), 309.
• 32)   K.  Lei,  N.  Taniguchi and  T.  Kobayashi: JSME Int. J. Ser. B, 45 (2002), 164.
• 33)   W.  Zhang,  K.  Tainaka,  S.  Ahn,  H.  Watanabe and  T.  Kitagawa: Adv. Powder Technol., 29 (2018), 2322.
• 34)   S.  Bhatia and  D.  Perlmutter: AIChE J., 26 (1980), 379.
• 35)   P.  DesJardin and  S.  Frankel: Phys. Fluids, 10 (1998), No. 9, 2298.
• 36)   W. P.  Jones and  R. P.  Lindstedt: Combust. Flame, 73 (1988), 233.
• 37)   C. K.  Westbrook and  F. L.  Dryer: Combust. Sci. Technol., 27 (1981), 31.
• 38)   H.  Watanabe,  K.  Tanno,  H.  Umetsu and  S.  Umemoto: Fuel, 142 (2015), 250.
• 39)  CHEMKIN-PRO, Release 15113, Reaction Design, Canonsburg, PA, (2013).
• 40)   G. P.  Smith,  D. M.  Golden,  M.  Frenklach,  N. W.  Moriarty,  B.  Eiteneer,  M.  Goldenberg,  C. T.  Bowman,  R. K.  Hanson,  S.  Song,  W. C.  Gardiner, Jr.,  V. V.  Lissianski and  Z.  Qin: Gri-Meck 3.0, http://combustion.berkeley.edu/gri-mech/, (accessed 2019-10-02).
• 41)   T.  Maffei,  A.  Frassoldati,  A.  Cuoci,  E.  Ranzi and  T.  Faravelli: Proc. Combust. Inst., 34 (2013), 2401.
• 42)   Z.  Zhou,  Q.  Xue,  C.  Li,  G.  Wang,  X.  She and  J.  Wang: Appl. Therm. Eng., 123 (2017), 1096.