A parametric study of lycopodium dust flame

Dust flames are associated with two-phase combustion phenomena where flame characteristics depend on interactions between solid and gas phases. Since organic dust particles can be effectively utilized in energy production systems, investigation of this phenomenon is essential. In this study, an analytical model is presented to simulate the combustion process of moist organic dust. The flame structure is divided into three zones: preheat zone, reaction zone, and postflame zone. To determine the effects of moisture content and volatile evaporation, the preheat zone is also divided into four subzones: first heating subzone and drying subzone, second heating subzone, and volatile evaporation subzone. The results obtained from the presented model are in reasonable agreement with experimental data for lycopodium particles. An increase in moisture content causes a reduction in burning velocity owing to moisture evaporation resistance. Consequently, the effects of some important parameters, like volatilization temperature, volatilization Damköhler number and drying Damköhler number are investigated. In special cases, like high moisture content, low volatilization temperature, and high drying resistance, the second heating subzone is omitted.

of energy production [1][2][3][4]. Biomass has been identified as a major renewable resource [5]. Furthermore, power generation from biomass is a CO 2 -neutral method and will help us to preserve our environment for future generations. Biomass is presently estimated to contribute on the order of 10-14% of the world's energy supply [6].
Because they have a low energy density, it is not economical to transport biomass fuels to power plants. However, a small-scale biomass conversion system is suitable for local usage. The Stirling engine is one of the best available technologies for small-scale biomass power production; it can be driven by many kinds of heat sources and can be incorporated into biomass energy utilization research projects as a viable prime mover [7][8][9]. Stirling engines are fed up to 500 μm microscale biomass to enhance their efficiency.
Besides power generation issues, microscale organic dusts demonstrate an extraordinary combustive behavior. Dust explosion occurs in most organic and inorganic dusts. A dust explosion is initiated by the rapid combustion of flammable particulates suspended in the air. Any solid material that can burn in the air will have the same behavior, with a violence and speed that increase with the degree of subdivision of the material. Any solid material have been recognized as a threat to humans and their property for the last 150 years [10].
Generally, combustible particles ignite in two phases: a gaseous and a solid phase. Since most organic particles combust in a gaseous phase, the volatilization and combustion mechanisms of volatile particles needs to be investigated. So far, little research has been done in this area. For example, fundamental information, such as the structure and movement of a combustion zone in a dust particle cloud in a vertical duct, is still ambiguous [11,12]. In general, solid particles heat up until they reach the evaporation temperature; then volatiles come out of the particles, and finally these volatiles combust in a reaction zone.
Unlike dust combustion, flame propagation in a uniformly dispersed spray has been widely studied [13][14][15]. It has been shown for combustion in liquid sprays that if the overall equivalence ratio of the initial combustible mixture is greater than 0.7, the flame standoff distance is greater than the characteristic separation distance between the liquid droplets [11,16]. This assumption guarantees continuity conditions in the flame structure. This finding regarding spray combustion can be generalized for the combustion of organic dust particles. With respect to organic dust, it is accepted that for small n s (local number density of particles) in the reaction zone, fuel particles burn with an envelope diffusion flame surrounding each particle. In the present paper, it is assumed that the value of n s is large enough that the standoff distance of the envelope flame surrounding each particle is much larger than the characteristic separation distance between particles [16]. According to the preceding discussion, our analytical model is valid for ϕ u > 0. 7. In experiments, to measure the explosivity of other types of dust particles and to calibrate dust explosion testing, lycopodium has been used as a reference dust since it has good dispersability and flowability (ISO 1985) [10]. Researchers have conducted many studies on organic dust combustion both experimentally [10,17,18] and analytically [16,19,20].
Laminar flames of lycopodium in air in a dust concentration range of 125-190 g/m 3 has been investigated [21]. Based on the wall cooling effects in these experiments, the maximum achieved burning velocities is similar to those found by Kaesche-Krischer and Zehr [22]. The burning velocity was found to be approximately 0.25 m/s in both experiments.
The effect of gravity on lycopodium dust combustion has been investigated [23]. In the presence of gravity, dust burning velocity was higher compared to free-gravity conditions in an upward stream.
A thermal gravity analysis (TGA) test was used by Han et al. [10] to determine the combustive properties of lycopodium particles. They calculated the burning velocity, flame temperature, and flame length. In their observations, they found that a propagating lycopodium flame front is discontinuous and not smooth. The flame in the spatial area between independent flames or individual burning particles is not observed. The researchers reported that a lycopodium flame cannot move (propagate) continuously in comparison with premixed gaseous flames.
Proust [17] has reported the burning velocity and flame temperature for three different particles; lycopodium, starch, and sulfur. Proust reported that there appeared to be a significant discrepancy between the theoretical and experimental values of maximum flame temperatures. He thought that this difference was a consequence of severe heat losses, but in reality, the mixing with air of combustible vapors produced during the pyrolysis of the particles  could not proceed to completion at the onset of combustion; consequently, the oxidation in the reaction zone of the flame was not complete. Seshadri et al. [16] studied analytically the structure of premixed flames propagating in combustible systems containing uniformly distributed volatile fuel particles in an oxidizing gas mixture. Bidabadi et al. [19] presented an analytical model for lycopodium dust combustion considering the temperature difference between gas and particles. Bidabadi et al. [24] calculated the effect of radiation emitted from the preheat reaction zone to the preheat zone and reported that taking into consideration the radiation mechanism has improved the burning velocity. A brief overview of previous researches carried out on lycopodium dust combustion fields is presented in Table 1.
In the present study, the flame propagation mechanism and the structure of the combustion zone are analytically investigated in order to clarify the mechanisms of flame propagation through moist dust clouds. The moisture content changes the flame structure, burning velocity, and flame temperature. Our last model [25] was developed to investigate the effect of moisture content. Also in this article, the effects of the Damköhler numbers (Da vap and Da dry ) on the combustion phenomena of organic dust particles are studied. In fact, the Damköhler numbers (Da vap and Da dry ) represent the tendency of particles to evaporate. These parameters have a large effect on the combustive behavior of organic particles. For example, particles that tend to evaporate in a highly volatile way combust quickly. In the present model, it is presumed that the moisture content of the solid particles evaporates first, then the dried particles vaporize to yield a gaseous fuel (CH 4 ); finally, the produced gaseous fuel combusts in a thin reaction zone. In the next section, all of the required equations are introduced and then solved using some assumptions and simplifications. In the results section, some of the effective parameters on flame propagation are explained. Then the conditions in which the second heating subzone can be omitted and in which the presented model eventually turns into a three-subzone model is investigated.

Theoretical model
The Lewis number (Le) for organic dust materials is defined as follows: where λ, ρ, C, and D are the thermal conductivity of the gaseous mixture, mixture density, mixture-specific heat, and characteristic mass diffusivity, respectively. In this study, the Lewis number is assumed to be unity [16,25]. In this article, it is assumed that a reaction occurs in a thin zone O(1/Ze), whereas preheat and postflame zones are quite large. This assumption is based on a high Zeldovich number, which is defined as follows: where E, R, T f , and T u are the activation energy of the reaction, the universal gas constant, the flame temperature, and the fresh mixture temperature, respectively.
The Damköhler number has various definitions [26,27]. While it is traditionally defined as the ratio of characteristic fluid mechanical time to characteristic chemical time, the volatilization Damköhler number (Da vap ) used in this study is defined as the ratio of volatilization time (τ vap ) to characteristic reaction time (τ chem ) or the ratio of chemical reaction rate to volatilization rate: Another dimensionless number that appears in this study is the drying Damköhler number (Da dry ), which takes into account the effect of the moisture content. Moisture exists in biomass either as bound water, which is held chemically within the cell walls, or as free water, which is stored in the cell cavities. In the drying process, freshly cut or green biomass starts to dry, and initially, the free water evaporates. The fiber saturation point is reached when all the free water is gone, leaving only the bound water within the cell walls [28]. The volatilization resistance of water within biomass particles, which is calculated according to both free and bound water contents, is represented here as A one-step overall reaction is used as a combustion process model: , where the symbols F, O 2 , and P denote the fuel, oxygen, and product, respectively, and the quantities ν F , ν O 2 , and ν P denote their stoichiometric coefficients. The governing equations are introduced as follows.
• Solid fuel particle conservation: In the preceding equation, Y s , V , and ρ denote the mass fraction of organic dust particles, the burning velocity, and the density of the mixture, respectively. ρ u = 1.135 × 10 −3 g/cm 3 which is extracted from Seshadri et al. [16]. w dry and w vap are the particle drying rate and volatilization rate, which are expressed by where H , τ dry , τ vap , T , T dry , and T vap are the Heaviside function, constant characteristic time of drying and volatilization, mixture temperature, and threshold temperature of drying and volatilization, respectively. • Gaseous fuel conservation: where Y g and D u are respectively the mass fraction of gaseous fuel gained from the volatilization of organic dust particles and the binary diffusion coefficient of the gaseous limiting component. w chem is the rate of chemicalkinetic for lycopodium particles which is defines as below: In the preceding equation, B is the frequency factor. B = 3.5 × 10 6 [mol −1 s −1 ], and E = 96.2 [kJ/mol] which is extracted from Seshadri et al. [16]. • Energy conservation: where Q dry , Q vap , and Q are the drying heat, volatilization heat, and reaction heat per unit mass of consumed fuel particles, respectively. Q = 50,009 kJ/kg fuel, Q vap = 500.09 kJ/kg water, and λ u = 1.46 × 10 −3 J/(cm s K) which is extracted from Seshadri et al. [16]. Q dry = 2,257 kJ/kg water was taken from thermodynamic tables.
C is the heat capacity of a mixture: where ρ s , r p , n s , (C p ), and (C s ) are the particle density, particle radius, average number of particles per unit volume, heat capacity of the gas, and heat capacity of the mixture, respectively. ρ s = 1 g/cm 3 which is extracted from Seshadri et al. [16].

Nondimensionalization of governing equations
Dimensionless parameters are defined as follows: In the preceding equations, θ = 0 states the dimensionless unburned mixture temperature and θ = 1 the dimensionless flame temperature. In addition, T f and V u are respectively the flame temperature and burning velocity. Y FC is defined as Consequently, the dimensionless governing equations (5), (8), and (10) are obtained as follows: The definitions of q vap and q dry are as follows: Volatilization Damköhler Da vap and drying Damköhler Da dry numbers andŵ chem are expressed as follows: Also, θ dry , θ vap , and β are described as

Flame structure and evaluation of flame characteristics
To partition the flame structure and solve the governing equations, we should neglect some insignificant terms and use some simplifications in each zone. When Ze tends to positive infinity (Ze → +∞), the reaction rate is negligible everywhere except in a small zone located near Z = 0, where θ = 1. Based on the preceding discussion, the flame structure is divided into three regions: preheat zone: The preheat zone is a very important in analysis of flame structure, thus this zone itself is divided into four subzones for a better understanding [29]. These subzones are explained below.
The first heating subzone is a region where particles heat up to T = T dry , and there is no change in weight of the particles. In the drying subzone, particles' moisture evaporates in a temperature-increasing process and water vapor is produced. This water vapor diffuses to the previous and the next subzones. In the second heating subzone, the temperature of the particles increases up to volatilization temperature. Finally, in the volatile evaporation subzone, the particles vaporize to yield a gaseous fuel. The gaseous fuel also diffuses to other regions. Because of diffusion phenomena, there is no sharp change in the temperature and quantities of the gaseous fuel. The flame structure of organic dust particles is shown in Fig. 1.
To obtain the flame structure, an analytical method is used to solve the governing equations with existing boundary conditions and a matching condition.
(1) The required boundary conditions for the subzone I are as follows: At the border of the first heating and drying subzones, the boundary conditions are as follows: Equation (24) demonstrates the continuity in the domain, which means there is no change in the quantity of mass fraction of solid and gaseous fuel at the borders. In addition, the diffusivity of gaseous fuel causes no change in the derivative of gaseous fuel (mass fraction) concentration. At the interface of the drying and the second heating subzones, the boundary conditions are as follows: Finally, at the interface of the second heating and the volatile evaporation subzones, the boundary conditions are as follows: (2) At the interface of the volatile evaporation subzone and the reaction zone, the boundary conditions are as follows: (3) Convection and volatilization terms are considered to be negligible in the reaction zone in comparison with reaction and diffusion terms. The integration of reactive diffusive equations (15) and (16) in the reaction zone gives the following matching condition: As discussed previously, in the preheat zone the reaction term is neglected in favor of the convection and diffusion terms. q vap is also a negligible quantity (close to zero) in comparison with other terms of the energy conservation equation, which means that heat released from a reaction is greater than heat absorbed by particles for volatilization. Thus, q vap will not be considered in the energy equation. To solve the conservation equations analytically, q dry is also assumed to be negligible. To compensate for the effect of q dry , the moisture evaporation effect should be considered in the solution of the energy equation. Thus, the energy equation and its boundary conditions are written as follows: The dimensionless temperature is defined as To consider the effect of moisture content (M), a coefficient F[M, φ u ] is multiplied by the final solution of the problem where T f appears, and it is defined as follows: This coefficient is obtained by comparing the adiabatic temperatures of a gaseous fuel (C H 4 ), which is obtained from thermodynamics, to a gaseous fuel with water content. As shown in Fig. 2, it is clear that increasing the moisture content causes a decrease in the adiabatic flame temperature of the gaseous fuel and F[M, φ u ]. ϕ u is defined in the following equations [16]: For wet particles the definition of ϕ u and ϕ g are written as a function of moisture content (M): In the preceding equations, Y Fu and Y FC are defined as follows:

First heating subzone (−∞ < Z < Z d )
In this subzone (−∞ < Z < Z d ), the particles only heat up to T = 100 • C and start to dry. Consequently, there is no change in the mass fraction of the particles: where Y FC is as defined in Eq. (13). No gaseous fuel is produced in this subzone. To evaluate the mass fraction of the gaseous fuel, the diffusion and convection terms are considered while other terms are assumed to be negligible. Thus: The boundary conditions force B 2 to be zero, and the quantity of B 1 is obtained by applying the boundary condition between the first heating and the drying subzones. This boundary condition is also coupled with other matching conditions.

Drying subzone
In this subzone, moisture evaporation causes the particles to lose their weight. The mass fraction of the particles can be determined using the following equation: Now this equation is used to calculate Z h2 : Also, Z d is calculated using Eq. (31) by considering the fact that the drying temperature of the material is known, For gaseous fuel, the mass fraction distribution is the same as in the previous subzone:

Second heating subzone
This subzone is similar to the first heating subzone; boundary temperatures are the only difference between them. In addition, there is no change in the weight of the particles and no gaseous fuel is produced. However, gaseous fuel diffuses from the next subzone (volatile evaporation subzone). Thus, Z vap is calculated using Eq. (31) and the pyrolysis temperature of the material (T vap ).

Volatile evaporation subzone
This subzone is the most important region in the flame structure. In this subzone, gaseous fuel is produced and then diffuses to other regions. The mass fraction of the particles is determined to as dy s dZ = − y s Da vap , The amount of unburned solid mass is calculated using the following expression: Finally, the gaseous fuel in this subzone is evaluated using the following equation: The constants in the gaseous fuel distribution in all subzones are dependent on each other and can be determined by solving eight equations simultaneously using existing boundary conditions.

Postflame zone (0 + < Z < +∞)
Before analyzing the reaction zone, the postflame zone is investigated in depth: In the postflame zone, the available mass fraction of solid particles is approximately equal to the quantity of y s at the end of the volatilization zone (i.e., Z = 0 − ).

Reaction zone
In this region, the convective and volatilization terms, in the conservation equations, are too small in comparison with the diffusive and reactive terms. To analyze the structure of this zone, the following expressions are introduced: where η, y * , and t are expanded using the expansion parameter ε = 1/Ze. This expansion parameter is presumed to be small. Substituting the expansion relations into the gaseous fuel and energy conservation equations (15) and (16), they yield the following expressions, respectively: In the preceding equations, b = y g F /ε, in which y g F is the amount of gaseous fuel at the beginning of the reaction zone.
Combining Eqs. (48) and (49) and using some mathematical operations, it is possible to obtain the burning velocity, which is the main purpose of this article: In the preceding equations, φ u is the equivalence ratio based on fuel available in the particles in the ambient reactant stream. Both expressions on the right-hand side of the matching equation are too small (according to Eq. 46) and can be neglected: Now, replacing the left-hand side with correlations leads to a relation between Da vap , Ze, Le, and θ vap : where F 1 , I 1 , and I 2 are obtained from Eq. (45): (52)

Results and discussion
To predict flame characteristics, such as burning velocity and mass fractions of the solid and gas, the implicit expression for the burning velocity in Eq. (50) should be solved simultaneously with Eq. (51). The iterative process starts with a known φ u and a guess of the temperature. The temperature is plugged into Eq. (50), and the burning  velocity is obtained from this equation. Then both the guessed temperature and the obtained burning velocity are inserted into Eq. (51). If the result is less than 0.001, then both the temperature and the burning velocity are correct. Otherwise, we must make another guess about the temperature and repeat the iterative process. To evaluate the accuracy of the presented model, the burning velocity and the flame temperature are compared with flame velocity and flame temperature using the calculations of previous researchers. Other researchers calculated the burning velocity of lycopodium particles as a function of mass particle concentration using a variety of methods.
To compare our results with the experimental data, some initial data (physical characteristics), such as r p , T vap , Ze, Le, and Da vap of the lycopodium particle, are essential. These parameters are obtained from previous works (T vap ≈ 180 • C, from a Derivative Thermogravimetry (DTG) analysis of a Thermogravimetry (TG) diagram, which Han measured for lycopodium dust particles [10], Ze ≈ 5.5 [16], and Le = 1 [16]). As shown in Fig. 3, the radius of a lycopodium particle is approximately r p = 15.64 μm. This measurement is completely compatible with other results ( Table 2). As is shown in Fig. 3, it is reasonable to assume that these particles are mono-size in lycopodium clusters, as reported in Table 2. Using Eq. (51) it is possible to determine Da vap . Since other researchers have calculated the flame velocity and flame temperature for dry lycopodium particles, F[M, φ u ] = 1 is used in Eq. (51). Consequently, Da vap for dry lycopodium particles would be 0.7.
As shown in Fig. 4a, our presented analytical model results for burning velocity as a function of particle concentration are compared with other researchers' experimental results. Also, the flame temperature obtained for the presented model is compared with that of previous works (Fig. 4b). According to Fig. 4, the evolution of the burning velocity and flame temperature as a function of mass particle concentration is in good agreement with experimental data.
There are no experimental data for determining Da dry . However, since the interaction between volatile fuel and the matrix of solid fuel particles is stronger than the interaction between the moisture content and matrix of solid fuel particles, it is possible to assume Da dry < Da vap . Based on the physical characteristics of lycopodium, it is assumed that those of other organic dust particles are similar and on the same order as lycopodium dust. As discussed earlier, the overall equivalence ratio of the initial combustible mixture must be larger than 0.7; thus, we present the results for ϕ u > 1 to ensure flame continuity.
According to Fig. 5a, an increase in moisture content causes a reduction in flame temperature. This is due to the latent heat of water. In fact, if the moisture content increases, more heat is consumed to produce the same amount of gaseous fuel. It is worth mentioning that in gas flame analysis, two types of resistance exist against flame propagation: heat transfer resistance and chemical reaction resistance. In dry dust flames, volatilization resistance is also added to those resistances. In wet dust, moisture evaporation resistance also exists before flame propagation. In fact, the interactions of these resistances in each flame (fuel) control the flame characteristics, like the burning velocity and flame temperature. It is clear that mass transfer resistance exists in all types of the aforementioned flames.
According to Fig. 5b, the burning velocity decreases when the amount of moisture content surges from 0 to 30% at r p = 15.5 μm, Da dry = 0.3, Da vap = 0.7, and Le = 1. As the moisture content in the particle increases, the moisture evaporation resistance also increases. Therefore, flame propagation slows down, which leads to a decrease in the burning velocity. Also, an increase in the equivalence ratio ϕ u (particle number density) reduces the volatilization resistance, and consequently, the burning velocity increases.
As shown in Fig. 5c, it is concluded that an increase in the moisture content causes an increase in the mass fraction of solid particles. This can be illustrated by the following equation:  In the first heating and the drying subzones, the presence of moisture yields a higher mass fraction of solid particles (Y s ). But in the second heating and the volatile evaporation subzones, there is no moisture in the particles. In other words, by omitting H 2 O from the numerator of the preceding equation, the denominator of the equation does not change. Thus, for the second heating and the volatile evaporation subzones, the mass fraction (Y s ) of the wetter particles becomes less than that of the drier particles. Since the start points of the drying and volatilization temperatures are considered to be constant, their dimensionless distances from the reference point (Z = 0) are almost constant. But the end point of the drying process (which is the start point of the second heating subzone) is greatly affected by the moisture content. By increasing the moisture content, the drying subzone becomes bigger.
Volatilization temperature is one of the essential parameters in the behavior of organic dusts. As shown in Fig. 6a, b, a decrease in the volatilization temperature causes an increase in the burning velocity and mass fraction of the produced gaseous fuel. In fact, if the volatilization temperature decreases, the initiation point of volatilization happens sooner (Fig. 6c).
Clearly any increase in the moisture evaporation resistance (Da dry ) will result in a decrease in the burning velocity. Furthermore, it will cause the drying process subzone to become bigger and, consequently, will affect the start point of the second heating zone, as shown in Fig. 7a. An increase in Da vap will also result in a decrease in the rate of volatilization, which means a higher volatilization resistance for a constant reaction rate. Consequently, the burning velocity and mass fraction of the gaseous fuel decrease, and a larger amount of mass fraction (Y s ) is left unburned, as shown in Fig. 7b. Finally, the conditions in which the presented model is transformed into a three-subzone model are discussed. From Fig. 8 it is clear that the moisture content, volatilization temperature, and Da dry are effective parameters regarding the structure of the preheat zone. It is also obvious that an increase in the moisture content or Da dry or a decrease in the volatilization temperature will cause the proposed model to transform into a three-subzone model.

Conclusion
An analytical model was presented to determine the effects of moisture content on flame characteristics and to investigate the effective parameters for organic dust combustion. To analyze the presented model, it is assumed that, first, the particles vaporize to yield a familiar chemical compound (methane) before starting to combust in the reaction zone. The flame structure is divided into three zones: preheat zone, reaction zone, and postflame zone. The preheat zone is further subdivided into four subzones: the first heating subzone and the drying subzone, the second heating subzone, and the volatile evaporation subzone. This division makes it easier to find the counteraction among various parameters and to view separately the effects of each parameter on the flame structure.
Then the governing equations are nondimensionalized, and these equations with their boundary and matching conditions are simultaneously solved. Consequently, the burning velocity and temperature profiles obtained from this model are compared with previously published experimental data. The comparison demonstrates that the proposed model accurately predicts the behavior of burning velocity and flame temperature.
The results make it clear that an increase in the moisture content causes the flame velocity and flame temperature to decrease. Other physical characteristics of organic particles also affect the flame structure. Volatilization Damköhler number (Da vap ) and drying Damköhler number (Da dry ) are the determining factors in the combustion of dust particles. It is also seen that an increase in Da vap causes a reduction in the mass fraction of available gaseous fuel and a subsequent decrease in the burning velocity. An increase in the volatilization temperature leads to a decrease in the gaseous fuel mass fraction. In addition, the onset of the volatilization process and the maximum of the diagram move toward the reaction zone.
The present research has shown that an increase in the moisture content or Da dry or a decrease in the volatilization temperature will cause a four-subzone model to switch to a three-subzone model. All the topics discussed in connection with the proposed model are valid for a three-subzone model. However, to avoid repetition, the results of the secondary model are not presented.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.