Improvement of heat- and mass transfer modeling for single iron particles combustion using resolved simulations

ABSTRACT

In this work, we use a boundary layer resolved model to improve a Lagrangian point particle model to simulate the combustion of single iron particles. By resolving the full boundary layer, mass and heat transfer are accurately modeled, including Stefan flow. Therefore, the model is suitable to improve point particle models. This work focuses on the first stage of iron combustion, which lasts up to the maximum temperature. Temperature- and composition-dependent properties are used and phase transitions from solid to liquid and liquid to gas are taken into account. The Nusselt and Sherwood correlations are investigated in conditions typical for iron particle combustion. It is found that the $1/2$-film temperature is the best film rule to use to model heat- and mass transfer for iron particle combustion. The boundary layer resolved model is used to validate the point particle models. Then, the model is systematically elaborated by including a temperature-dependent particle density, slip velocity and Stefan flow. The individual and combined effect of these phenomena on the burn duration are investigated. Including all these effects decreases the time to maximum temperature by around $25\mathrm{\%}$. Furthermore, it is shown that if one neglects physical phenomena like slip and Stefan flow, but uses the $1/3$-film rule instead of the $1/2$-film rule, errors cancel and still reasonable agreement is obtained with experiments.

KEYWORDS:

• Metal fuel
• iron
• Nusselt and Sherwood number
• boundary layer resolved
• point particle

Introduction

Micron-sized iron is considered as a promising metal fuel since it is inherently carbon-free, recyclable, compact and also cheap and widely available (Bergthorson et al. 2015). In turn, the iron-oxide particles can be captured and reduced back into iron particles using renewable energy sources, thereby enabling a carbon-free energy cycle.

For the development of practical iron fuel burning setups, it is needed to get a better understanding of the fundamental characteristics of iron combustion. Multiple experimental studies (Poletaev and Khlebnikova 2020; Sun, Dobashi, Hirano 2000; Tang et al. 2009; Wright, Higgins, Goroshin 2016) have contributed to the understanding of iron dust combustion. Ning et al. (2021a, 2021b, 2022) performed single particle combustion experiments where the particles are ignited by a laser. The combustion process was divided into two stages, where the first stage lasts up to the maximum temperature. The burning process duration of the first stage is sensitive to the oxygen concentration and that the maximum particle temperature is a function of the oxygen concentration, saturating at elevated oxygen concentrations. The second stage, after the peak temperature, has a negligible effect on oxygen concentration. In this last stage of combustion, an intensity jump was observed, which was also observed during the combustion of other metal particles (Dreizin, Suslov, Trunov 1993). Li et al. (2022) investigated the ignition and combustion process of single micron-sized iron particles in the hot gas flow of burned methane-oxygen-nitrogen mixture, and draw similar conclusions as Ning et al. (2021a).

In our previous work (Thijs et al. 2022), a boundary-layer resolved transient model for iron particle combustion was presented. By resolving the full boundary layer, mass and heat transfer are accurately modeled, including Stefan flow. However, the development of accurate point particle models is of importance to simulate industrial burners with many particles, where the boundary layer around each particle cannot be resolved. The development of point particle models gained more attention in the last couple of years. Soo et al. (2018) developed a numerical model for the combustion of a single metal particle, which was extended by Hazenberg and van Oijen (2021) by taking into account more realistic parameters and combustion products were not assumed to vanish. Schiemann, Fischer, and Bergthorson (2017) developed a numerical model which assumes the oxidation from $\mathrm{F}\mathrm{e}$ to $\mathrm{F}{\mathrm{e}}_3{\mathrm{O}}_4$, but used a predefined maximum temperature based on the dissociation temperature of $\mathrm{F}{\mathrm{e}}_3{\mathrm{O}}_4$. Mi, Fujinawa, and Bergthorson (2022) used solid-phase oxidation kinetics described by a parabolic rate to improve their iron combustion model in the solid phase, but they only focused on the ignition phase. The above described models make assumptions on the effect of slip, Stefan flow, or evaporation.

To enable modeling of single particle combustion without resolving the full boundary layers, reliable heat and mass transfer coefficients must be used. To model the heat transfer in the boundary layer, a common correlation used for the Nusselt number was proposed by Ranz and Marshall (1952)

(1) $\mathrm{N}\mathrm{u}=2+0.6\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{1/2}$(1)

This correlation has been obtained by measuring the temperature of water droplets in an air atmosphere for ambient temperatures and is valid up to $\mathrm{R}\mathrm{e}=200$. The maximum temperature difference tested between the water droplet and the air stream was $169\mathrm{K}$ and provides therefore a small variation in the thermophysical properties of air across the boundary layer.

Richter and Nikrityuk (2012) performed numerical calculations of heat and fluid flow past spherical particles and non-spherical particles of various shapes. They developed a new correlation which is dependent on the sphericity $\varphi$ and is valid from $\mathrm{R}\mathrm{e}=10$ up to $\mathrm{R}\mathrm{e}=250$. For a spherical particle, $\varphi =1$, their correlation becomes

(2) $\mathrm{N}\mathrm{u}=1.76+0.55\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{1/2}+0.014\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{2/3}.$(2)

They considered a temperature difference between the particles and fluid of $100\mathrm{K}$, but used all constant material properties. Note that this correlation does not approach the theoretical limit of $\mathrm{N}\mathrm{u}=2$ if $\mathrm{R}\mathrm{e}\to 0$.

The reference temperatures which are used to determine the transport properties differ in the literature. Ranz and Marshall (1952) determined their transport properties with the $1/2$-film rule, which is defined as the average between the far-field temperature and the particle temperature. Whitaker (1972) used the gas temperature as a reference temperature, while Hubbard, Denny, and Mills (1975) stated that for engineering calculations, the most appropriate reference state for the evaluation of properties is the $1/3$ rule. The latter is also done in the earlier work regarding metal fuels (Hazenberg and van Oijen 2021; Mi, Fujinawa, Bergthorson 2022; Ravi, de Goey, van Oijen 2022).

Several authors (Fiszdon 1979; Whitaker 1972) investigated the effect of large temperature differences on the Nusselt correlations. Ellendt et al. (2018) concluded that the Ranz-Marshall correlation can be applied to high-temperature differences in the laminar regime, with the properties calculated using the $1/2$-film rule. They performed numerical simulations with Comsol Multiphysics 5.2, but did not show any results regarding their derived Nusselt numbers as a function of the Reynolds number. Instead, they only stated that their data fit the well-established Ranz and Marshall correlation very well ($R^2=0.99$). More recently, Jayawickrama et al. (2021) investigated the Nusselt number with Stefan flow under large temperature differences. They proposed to use a volume averaged temperature, dependent on the Stefan velocity, Reynolds number and temperature.

In this work, we will use our boundary-layer resolved model (Thijs et al. 2022) to further improve the point particle models for the first stage of iron particle combustion. The Nusselt and Sherwood correlations are reviewed in flame-like conditions and the choice of reference temperature on the Nusselt correlations is investigated. Then, the point particle model is compared to the boundary-layer resolved model and the latter is used to investigate the effect on the time to maximum temperature when neglecting a temperature-dependent density, slip or Stefan flow.

The paper is organized as follows. Section 2 describes the methodology of the boundary-layer resolved model and the Lagrangian point particle model. Section 3 presents the results and discusses the above-mentioned points. Finally, conclusions are provided in Section 4.

Methodology

In this section, the methodology for simulating the combustion of single iron particles is described. First, the boundary layer resolved model is introduced, which is used as a validation tool for the Lagrangian point particle model.

Boundary layer resolved model

Figure 1 shows the computational domain used for the boundary layer resolved model. A 2D axi-symmetric model for the combustion of a single iron particle with diameter $d_{\mathrm{p}}$ was developed, of which more details can be found in Thijs et al. (2022).

Figure 1. Boundary layer resolved geometry (Thijs et al. 2022).

For the boundary layer resolved model the continuum approach will be used for modeling the mass and heat transfer. The governing equations for mass, species, momentum and energy in the gas phase ${\mathrm{\Omega }}_{\mathrm{g}}$ read:

(3) $\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left(\rho \vec{u}\right)=0,$(3)

(4) $\rho \frac{\mathrm{\partial }{Y}_i}{\mathrm{\partial }t}+\rho \left(\vec{u}\cdot \mathrm{\nabla }\right)Y_i=-\mathrm{\nabla }\cdot j_i+{\omega }_{\mathrm{g},i},$(4)

(5) $\rho \frac{\mathrm{\partial }\vec{u}}{\mathrm{\partial }t}+\rho \left(\vec{u}\cdot \mathrm{\nabla }\right)\vec{u}=\mathrm{\nabla }\cdot \tau ,$(5)

(6) $\rho c_{\mathrm{p}}\frac{\mathrm{\partial }T}{\mathrm{\partial }t}+\rho c_{\mathrm{p}}\vec{u}\cdot \mathrm{\nabla }T-\mathrm{\nabla }\cdot \left(k\mathrm{\nabla }T\right)=Q,$(6)

Here, $\rho$, $c_{\mathrm{p}}$, $T$, $k$, $\mu$, $\vec{u}$ and $p$ are the density, specific heat capacity at constant pressure, temperature, thermal conductivity, dynamic viscosity, velocity and pressure of the gas phase, respectively. The equations are solved on a body-fitted mesh with the “Comsol Multiphysics 5.5” package.

While in Thijs et al. (2022) a kinetic model was used in the gas phase, this is not considered in this work. Therefore, only the species ${\mathrm{O}}_2$, ${\mathrm{N}}_2$, $\mathrm{F}\mathrm{e}$ and $\mathrm{F}\mathrm{e}\mathrm{O}$ are considered for the gas phase. The density of the gas mixture is based on the ideal gas law, the transport properties are based on kinetic gas theory and thermodynamic properties of these species are described by the 7-coefficient NASA polynomials. For the $\mathrm{F}\mathrm{e}$-containing species, data from the NASA database (McBride 2002), Goos et al. (2005). and Rumminger et al. (1999) are used. For the remaining species, the values from Konnov (2019) are used.

We only consider the conversion of $\mathrm{F}\mathrm{e}$ to $\mathrm{F}\mathrm{e}\mathrm{O}$, which is the relevant oxide up to the maximum particle temperature (Thijs et al. 2022). We assume that mass and species transport is infinitely fast within the condensed phase ${\mathrm{\Omega }}_{\mathrm{p}}$ (Sun, Dobashi, Hirano 2000; Thijs et al. 2022) and only the energy equation has to be solved, which is given by

(7) ${\rho }_{\mathrm{c}}{c}_{\mathrm{p},\mathrm{c}}\frac{\mathrm{\partial }T}{\mathrm{\partial }t}-\mathrm{\nabla }\cdot \left(k_{\mathrm{c}}\mathrm{\nabla }T\right)=Q,$(7)

with ${\rho }_{\mathrm{c}}$, $c_{\mathrm{p},\mathrm{c}}$ and $k_{\mathrm{c}}$ the density, specific heat at constant pressure and thermal conductivity of the condensed phase. The apparent heat capacity method (COMSOL 2019) is used to model the phase change of the iron particle. Temperature dependent densities are used for $\mathrm{F}\mathrm{e}$ (Hixson, Winkler, Hodgdon 1990; Watanabe et al. 1980) and FeO (Lee and Gaskell 1974; Saxena et al. 1993). A temperature dependent thermal conductivity is used and based on the values reported by Akiyama et al. (1992) and Ho, Powell, and Liley (1972). Thermodynamic properties are calculated with the NASA-7 polynomials with the coefficients obtained from the NASA database (McBride 2002).

Evaporation of the liquid iron droplet is taken into account via the saturated vapor pressure. The concentration of gaseous species on the particle surface can be written as

(8) $c_i{|}_{{\mathrm{\Gamma }}_{\mathrm{p}}}=\left(\frac{M_i}{R_{\mathrm{u}}T}\right)\sum _{j=1}^2{X}_j{p}_{\mathrm{v},i,j},$(8)

where $i$ and $j$ is either $\mathrm{F}\mathrm{e}$ and $\mathrm{F}\mathrm{e}\mathrm{O}$ and $X_i{p}_{\mathrm{v},i,j}$ is the molar fraction in the condensed phase times the vapor pressure, stemming from Raoult’s law. The vapor pressure is a function of temperature and is based on thermochemical equilibrium calculations (Ning et al. 2021b). These calculations were performed with Cantera for pure $\mathrm{F}\mathrm{e}(\mathrm{l})$ and $\mathrm{F}\mathrm{e}\mathrm{O}(\mathrm{l})$ at $2000$ – $3500\mathrm{K}$ and $1\mathrm{a}\mathrm{t}\mathrm{m}$.

Lagrangian particle model

For the point particle model we extend the model of Hazenberg and van Oijen (2021). The rate of change of $m_{\mathrm{p}}$ and $m_{\mathrm{F}\mathrm{e}}$ are related to the rate of oxidation and the rate of evaporation via

(9) $\frac{\mathrm{d}{m}_{\mathrm{p}}}{\mathrm{d}t}={\dot{m}}_{{\mathrm{O}}_2}\mathrm{D}{\mathrm{a}}^{\ast }-\underset{i}{\hspace{0.17em}\sum }\frac{\mathrm{d}{m}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{p},i}}{\mathrm{d}t},$(9)

(10) $\frac{\mathrm{d}{m}_{\mathrm{F}\mathrm{e}}}{\mathrm{d}t}=-\frac{1}{s}{\dot{m}}_{{\mathrm{O}}_2}\mathrm{D}{\mathrm{a}}^{\ast }-\frac{\mathrm{d}{m}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{p},\mathrm{F}\mathrm{e}}}{\mathrm{d}t},$(10)

with $m_{\mathrm{p}}$ the mass of the particle, ${\dot{m}}_{{\mathrm{O}}_2}$ the mass flux of oxygen to the particle, $\mathrm{D}{\mathrm{a}}^{\ast }$ the normalized Damköhler number (Hazenberg and van Oijen 2021) and $s$ the stoichiometric mass ratio. Here, $s=\frac{1}{2}{M}_{{\mathrm{O}}_2}/M_{\mathrm{F}\mathrm{e}}=0.29$.

The rate of change of the particle enthalpy is described by

(11) $\frac{\mathrm{d}{H}_{\mathrm{p}}}{\mathrm{d}t}=q+q_{\mathrm{r}\mathrm{a}\mathrm{d}}+{\dot{m}}_{{\mathrm{O}}_2}\mathrm{D}{\mathrm{a}}^{\ast }{h}_{{\mathrm{O}}_2}-\sum _i{h}_{i,v}\frac{\mathrm{d}{m}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{p},i}}{\mathrm{d}t}$(11)

with $q$ the convective heat flux, $q_{\mathrm{r}\mathrm{a}\mathrm{d}}$ the radiative heat flux, $h_{{\mathrm{O}}_2}$ the mass-specific enthalpy of the consumed oxygen evaluated at the particle temperature $T_{\mathrm{p}}$ and $h_{i,v}$ the mass-specific enthalpy of the evaporated species.

The mass flux and heat flux can be described as

(12) ${\dot{m}}_{{\mathrm{O}}_2}=2\pi \mathrm{S}\mathrm{h}{r}_{\mathrm{p}}{\rho }_{\mathrm{O}2,\mathrm{f}}{D}_{\mathrm{f}}{X}_{{\mathrm{O}}_2,\mathrm{g}},$(12)

(13) $q=2\pi \mathrm{N}\mathrm{u}{r}_{\mathrm{p}}{k}_{\mathrm{f}}\left(T_p-T_{\mathrm{g}}\right),$(13)

with ${\rho }_{\mathrm{O}2,\mathrm{f}}$, $D_{\mathrm{f}}$ and $k_{\mathrm{f}}$ the density of oxygen, the mixture-averaged diffusion coefficient and the thermal conductivity derived at the film temperature and composition, $T_{\mathrm{g}}$ the gas temperature, $X_{\mathrm{O}2,\mathrm{g}}$ the mole fraction of oxygen in the gas phase and $\mathrm{N}\mathrm{u}$ and $\mathrm{S}\mathrm{h}$ the new Nusselt and Sherwood numbers. Here, ${\rho }_{\mathrm{O}2}$ is calculated via the ideal gas law

(14) ${\rho }_{\mathrm{O}2,\mathrm{f}}=\frac{M_{\mathrm{O}2}p}{R_{\mathrm{u}}{T}_{\mathrm{f}}},$(14)

with $M_{\mathrm{O}2}$ the molar mass of oxygen, $p$ the pressure and $R_{\mathrm{u}}$ the universal gas constant.

In Section 3.1.2 we derive new Nusselt and Sherwood correlations which are accurate for low-Reynolds numbers and large temperature differences. We derive the transport properties at a reference temperature $T_{\mathrm{f}}$ and composition $X_{{\mathrm{O}}_2,\mathrm{f}}$ (Hubbard, Denny, Mills 1975)

(15) $T_{\mathrm{f}}=T_{\mathrm{p}}+A_{\mathrm{f}}\left(T_{\mathrm{g}}-T_{\mathrm{p}}\right),$(15)

(16) $X_{{\mathrm{O}}_2,\mathrm{f}}=X_{{\mathrm{O}}_2,\mathrm{p}}+A_{\mathrm{f}}\left(X_{{\mathrm{O}}_2,\mathrm{g}}-X_{{\mathrm{O}}_2,\mathrm{p}}\right),$(16)

In this work we will use the $1/2$-film rule $\left(A_{\mathrm{f}}=1/2\right)$ as the reference temperature and composition. The Nusselt and Sherwood correlations are

(17) $\mathrm{N}\mathrm{u}=2+\left(0.02\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{1/2}+0.33\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{2/3}\right)\mathrm{\Phi },$(17)

(18) $\mathrm{S}\mathrm{h}=2+\left(0.02\mathrm{S}{\mathrm{c}}^{1/3}\mathrm{R}{\mathrm{e}}^{1/2}+0.33\mathrm{S}{\mathrm{c}}^{1/3}\mathrm{R}{\mathrm{e}}^{2/3}\right)\mathrm{\Phi },$(18)

with $\mathrm{P}\mathrm{r}=\frac{\mu c_p}{k}$ the Prandtl number and $\mathrm{S}\mathrm{c}=\frac{\nu }{D}$ the Schmidt number. A factor $\mathrm{\Phi }$ is introduced to account for the effect on the transport properties due to the large temperature differences

(19) $\mathrm{\Phi }=1-a\left[1-{\left(\frac{T_{\mathrm{p}}}{T_{\mathrm{g}}}\right)}^b\right],$(19)

with $a=1.1$ and $b=0.2$.

To model Stefan flow induced by the oxygen diffusion toward the particle, theory from droplet condensation/evaporation is used. The Nusselt and Sherwood numbers are corrected for Stefan flow via (Bird, Lightfoot, Stewart 2002; Both, Mira, Lehmkuhl 2022)

(20) $\mathrm{N}{\mathrm{u}}_{\mathrm{S}\mathrm{t}\mathrm{e}\mathrm{f}\mathrm{a}\mathrm{n}}=\mathrm{N}\mathrm{u}ln\left(1+B_{\mathrm{T}}\right)/B_{\mathrm{T}},$(20)

(21) $\mathrm{S}{\mathrm{h}}_{\mathrm{S}\mathrm{t}\mathrm{e}\mathrm{f}\mathrm{a}\mathrm{n}}=\mathrm{S}\mathrm{h}ln\left(1+B_{\mathrm{m}}\right)/B_{\mathrm{m}},$(21)

with the mass Spalding number defined as (Bird, Lightfoot, Stewart 2002)

(22) $B_{\mathrm{m}}=\frac{Y_{{\mathrm{O}}_2,\mathrm{g}}-Y_{{\mathrm{O}}_2,\mathrm{p}}}{Y_{{\mathrm{O}}_2,\mathrm{p}}-1}.$(22)

Analogous to the mass Spalding number, a Spalding heat transfer number can be derived as (Bird, Lightfoot, Stewart 2002)

(23) $1+B_{\mathrm{T}}={\left(1+B_{\mathrm{m}}\right)}^{{\varphi }_m},$(23)

with ${\varphi }_m=\frac{c_{\mathrm{p},\mathrm{o}\mathrm{x}}}{c_{\mathrm{p}}}\frac{1}{\mathrm{L}\mathrm{e}}\frac{\mathrm{S}\mathrm{h}}{\mathrm{N}\mathrm{u}}$ and $\mathrm{L}\mathrm{e}=\frac{\mathrm{S}\mathrm{c}}{\mathrm{P}\mathrm{r}}$ the Lewis number.

The rate of evaporation can be described with and without evaporation induced Stefan flow. Since the mass flux of oxygen is much larger than the evaporate mass flux (Thijs et al. 2022), we can neglect the evaporation induced Stefan flow. Therefore, the evaporation rate becomes

(24) $\frac{\mathrm{d}{m}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{p},i}}{\mathrm{d}t}={\rho }_i{A}_{\mathrm{p}}{k}_{d,evap}{X}_i,$(24)

with $i$ either $\mathrm{F}\mathrm{e}$ or $\mathrm{F}\mathrm{e}\mathrm{O}$, $k_{d,evap}$ the evaporated mass transfer coefficient and $X_i$ the molar fraction of the $\mathrm{F}\mathrm{e}$ or $\mathrm{F}\mathrm{e}\mathrm{O}$ vapor at the particle surface. The vapor molar fraction is determined according to Equation 8(8) $c_i{|}_{{\mathrm{\Gamma }}_{\mathrm{p}}}=\left(\frac{M_i}{R_{\mathrm{u}}T}\right)\sum _{j=1}^2{X}_j{p}_{\mathrm{v},i,j},$(8) .

The above described equations are solved within MATLAB, were we used the ode15s function for solving the differential equations.

Results

Boundary layer resolved model

The boundary layer resolved model is used to investigate the influence of several mentioned physical phenomena in order to derive modeling strategies for the Lagrangian point particle model. The Nusselt and Sherwood correlations are reviewed in flame-like conditions, with a focus on the reference temperature used to determine transport properties. Then, the 0D model is compared to the boundary-layer resolved model and the latter is used to investigate the effect of different modeling strategies.

In the analysis below, the effect of Stefan flow on the Nusselt and Sherwood correlations will be neglected. Both, Mira, and Lehmkuhl (2022) investigated different methods subjected to flame-like conditions to incorporate Stefan flow effects in the heat mass transfer. They stated that the Bird’s correction, as described in Section 2, together with the Abramzon–Sirignano model, are the most relevant models to model the effect of Stefan flow.

Validation

First, simulations have been performed to investigate the effect of the domain size and the grid resolution on the obtained heat transfer coefficients. It is known that the effect of domain size becomes more important when the Reynolds number is decreased (Tajfirooz et al. 2021). we found that a domain of $L_{\mathrm{r}}=40d_{\mathrm{p}}$, $L_{\mathrm{z},\mathrm{u}}=30d_{\mathrm{p}}$, $L_{\mathrm{z},\mathrm{d}}=100d_{\mathrm{p}}$ with 120000 elements is sufficient to obtain an undisturbed flow when $\mathrm{R}\mathrm{e}>1$, while we use $L_{\mathrm{r}}=200d_{\mathrm{p}}$, $L_{\mathrm{z},\mathrm{u}}=100d_{\mathrm{p}}$, $L_{\mathrm{z},\mathrm{d}}=200d_{\mathrm{p}}$ with around 720000 elements if $\mathrm{R}\mathrm{e}<1$.

In order to verify the level of confidence of the used model, a validation of the program for different flow regimes is performed. Note that in Thijs et al. (2022) we validated the physical model used for the combustion of iron particles. Here, we will use the boundary-layer resolved model to investigate the heat transfer coefficients, and compare it to existing well-established correlations. In the model, the particle is kept stationary while the gas flows across it. Therefore, the particle Reynolds number can be determined via

(25) $\mathrm{R}{\mathrm{e}}_{\mathrm{f}}=\frac{u_0{d}_{\mathrm{p}}}{{\nu }_{\mathrm{f}}},$(25)

with $u_0$ the relative velocity, ${\nu }_{\mathrm{f}}$ the kinematic fluid viscosity determined at the reference temperature.

The numerically obtained Nusselt numbers are compared to the correlations proposed by Ranz and Marshall (1952) (Equation (1)(1) $\mathrm{N}\mathrm{u}=2+0.6\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{1/2}$(1) ) and Richter and Nikrityuk (2012) (Eq 2)(2) $\mathrm{N}\mathrm{u}=1.76+0.55\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{1/2}+0.014\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{2/3}.$(2) . Constant material properties corresponding to air at $300\mathrm{K}$ are used and we consider a temperature difference of $100\mathrm{K}$ between the particle and the gas. The results are shown in Figure 2. Good agreement with the Richter correlation can be seen.

Figure 2. Comparison of numerically obtained Nusselt numbers and Nusselt numbers obtained with correlations from literature. Note, constant material properties at $300\mathrm{K}$ are used.

Nusselt and Sherwood correlations at large temperature differences

To enable modeling of single particle combustion without resolving the full boundary layers, reliable heat and mass transfer coefficients must be used. In the previous Section we showed that the boundary-layer resolved model shows good agreement with existing correlations in a large range of Reynolds number. However, particles of $d_{\mathrm{p}}=50$ m or smaller are practically relevant for iron particle combustion, the particle Reynolds number will be in the lower Reynolds number regime, with $\mathrm{R}\mathrm{e}<10$.

We perform steady simulations to investigate the heat transfer in the boundary layer for different Reynolds numbers. Since the mass transfer coefficient is analogous to the heat transfer coefficient, the analysis shown below also holds for the Sherwood number (Bird, Lightfoot, Stewart 2002). We consider different cases with $T_{\mathrm{p}}=400$, $1500$ and $2500\mathrm{K}$ and $T_{\mathrm{g}}=300\mathrm{K}$, which we will refer to as a cold flow. Also the opposite case is investigated, with $T_{\mathrm{p}}<T_{\mathrm{g}}$, to investigate the effect of a cold particle in a hot flow. The latter case is important before the thermal runaway happens for the iron particle, while the first case is of importance after thermal runaway.

The average heat transfer coefficient for the particle is calculated as

(26) $h_{\mathrm{c}}=\frac{{\int }_A{q}_n\mathrm{d}A}{A\left(T_{\mathrm{p}}-T_{\mathrm{g}}\right)}$(26)

with $q_n$ the normal heat flux at the particle surface. Subsequently, the Nusselt number is calculated as

(27) $\mathrm{N}{\mathrm{u}}_{\mathrm{f}}=\frac{h_{\mathrm{c}}{d}_{\mathrm{p}}}{k_{\mathrm{f}}},$(27)

with $k_{\mathrm{f}}$ the thermal conductivity determined at the reference temperature and composition.

Figure 3 shows the Nusselt number correlations obtained with the boundary-layer resolved model in the case of a hot flow (curves with markers) and a cold flow (curves without markers) past an iron particle with diameter $d_{\mathrm{p}}=50$ m. The Reynolds number and Nusselt number are determined at the same reference temperature, which is either the $1/2$-film or the $1/3$-film temperature rule. Three different temperature differences are used, namely $10$ K (solid line), $1200$ K (dotted line) and $2200$ K (dashed line). For $A_{\mathrm{f}}=1/2$ it can be seen that the Nusselt correlations are sensitive to the large temperature differences in the higher Reynolds number regime, but are close to the theoretical limit of $\mathrm{N}\mathrm{u}=2$ in the lower Reynolds number regime. The opposite can be seen if $A_{\mathrm{f}}=1/3$. This means that the $1/3$-film rule would be a wise choice if the particle Reynolds number is large. Then, the Nusselt and Sherwood numbers are independent of the temperature differences. However, if micron sized particles are considered, the Reynolds number is generally small. Therefore, we conclude that the $1/2$-film rule is the best rule to use when modeling the combustion of single iron particles, since it is independent of temperature and approaches $\mathrm{N}\mathrm{u}=2$ in the lower Reynolds number regime.

Figure 3. Nusselt correlation for different reference temperatures. Line represent either a hot gas flow (with markers) or a cold gas flow (without markers) past a iron particle with a particle diameter $d_{\mathrm{p}}=50$ m. Temperature differences between the gas and the particle are: 10 K (solid line), 1200 K (dotted line) and 2200 K (dashed line).

Why the 1/2-file rule is more accurate in the low-Reynolds number regime can be demonstrated by investigating the 1D temperature profile in the boundary layer. When neglecting convection, the total heat conduction rate of a sphere can be written as (Liu et al. 2006)

(28) $q=4\pi r_{\mathrm{p}}{\int }_{T_{\mathrm{g}}}^{T_{\mathrm{p}}}k\mathrm{d}T,$(28)

with the thermal conductivity $k$ a function of temperature and written as

(29) $k=a{T}^n,$(29)

with constants $a$ and $n$. The average thermal conductivity in the range from $T_{\mathrm{p}}$ to $T_{\mathrm{g}}$

(30) $k_{\mathrm{f}}=\frac{1}{T_{\mathrm{p}}-T_{\mathrm{g}}}{\int }_{T_{\mathrm{g}}}^{T_{\mathrm{p}}}a{T}^n\mathrm{d}T=\frac{a}{n+1}\frac{T_{\mathrm{p}}^{(n+1)}-T_{\mathrm{g}}^{(n+1)}}{T_{\mathrm{p}}-T_{\mathrm{g}}}.$(30)

For air, the thermal conductivity scales almost linear with temperature ($n=0.7$). If we substitute $n=1$, and rewrite the equation we find

(31) $k_{\mathrm{f}}=\frac{a}{2}\left(T_{\mathrm{p}}+T_{\mathrm{g}}\right)=k\left(\frac{1}{2}\left(T_{\mathrm{p}}+T_{\mathrm{g}}\right)\right),$(31)

which shows that the conduction of the film is equal to the conductivity at the $1/2$-film temperature. Substituting Equation (31)(31) $k_{\mathrm{f}}=\frac{a}{2}\left(T_{\mathrm{p}}+T_{\mathrm{g}}\right)=k\left(\frac{1}{2}\left(T_{\mathrm{p}}+T_{\mathrm{g}}\right)\right),$(31) into Equation (28)(28) $q=4\pi r_{\mathrm{p}}{\int }_{T_{\mathrm{g}}}^{T_{\mathrm{p}}}k\mathrm{d}T,$(28) results in an expression of the heat flux, which is dependent on the average thermal conductivity.

(32) $q=4\pi r_{\mathrm{p}}{k}_{\mathrm{f}}\left(T_{\mathrm{p}}-T_{\mathrm{g}}\right)$(32)

From this we conclude that when neglecting convection, using the 1/2 rule will result in an accurate heat flux.

In Figure 2 it was shown that the correlation proposed by Richter and Nikrityuk (2012) gives good agreement with the boundary layer resolved model for a large range of Reynolds numbers. However, this correlation does not approach the theoretical limit of $\mathrm{N}\mathrm{u}=2$ if $\mathrm{R}\mathrm{e}\to 0$. In order to come up with a Nusselt correlation which is accurate in the regime $0.01<\mathrm{R}\mathrm{e}<10$ and applicable to high temperature differences, a similar approach as Wittig, Nikrityuk, and Richter (2017) and Richter and Nikrityuk (2012) is used. Steady simulations are performed where a small temperature difference of 5 K is applied to a particle of 50 µm. By using small temperature difference the impact of the reference temperature is almost negligible. The following equation is then fitted through the obtained Nusselt number

(33) $\mathrm{N}\mathrm{u}=2+a\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{1/2}+b\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{2/3},$(33)

with $a$ and $b$ the fitting parameters, which are found equal to $0.02$ and $0.33$, respectively. Figure 4 shows the results of the boundary-layer resolved model, the new proposed correlations and the correlations of Ranz and Marshall (1952) and Richter and Nikrityuk (2012). It can be seen that the new correlation is more accurate in this low Reynolds regime $0.01<\mathrm{R}\mathrm{e}<10$ than the existing correlations of Ranz and Marshall (1952) and Richter and Nikrityuk (2012).

Figure 4. Boundary-layer resolved obtained Nusselt numbers and different correlations.

To account for the changing thermophysical parameters in the boundary layer in the higher Reynolds number regime ($1<\mathrm{R}\mathrm{e}<10$), Equation (33)(33) $\mathrm{N}\mathrm{u}=2+a\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{1/2}+b\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{2/3},$(33) is corrected with a factor $\mathrm{\Phi }$ which is dependent on the temperature difference. A correction factor $\mathrm{\Phi }$, which is a function of the ratio $T_{\mathrm{p}}/T_{\mathrm{g}}$, can be fitted with Equation (19)(19) $\mathrm{\Phi }=1-a\left[1-{\left(\frac{T_{\mathrm{p}}}{T_{\mathrm{g}}}\right)}^b\right],$(19) , with $a=1.1$ and $b=0.2$.

Lagrangian point particle model

In this section, we present the results of a single iron particle burning in an ${\mathrm{O}}_2$-${\mathrm{N}}_2$ and ${\mathrm{O}}_2$-$\mathrm{A}\mathrm{r}$ atmosphere. The initial conditions are chosen such that the laser-ignited experiments as performed by Ning et al. (2021a) are mimicked. First, the Lagrangian particle model is compared to the boundary layer resolved model. Then, the Lagrangian particle model is used to investigate the effect of temperature-dependent particle density, slip velocity and Stefan flow on the time to maximum temperature. If the size of the particle becomes too small, modeling the heat and mass transfer using the continuum approach becomes invalid (Liu et al. 2006). Therefore, we will focus on particles with a diameter of $d_{\mathrm{p}}=50$ µm to ensure that we are still in, or at least close to the continuum regime (Sundaram, Puri, Yang 2016).

Comparison between boundary layer resolved and Lagrangian particle model

Temperature profiles obtained with the boundary layer resolved model are compared to the Lagrangian particle models. We consider a 50 µm particle burning in an ${\mathrm{O}}_2$-${\mathrm{N}}_2$ and ${\mathrm{O}}_2$-$\mathrm{A}\mathrm{r}$ atmosphere, with $X_{{\mathrm{O}}_2}=0.21$. The models include temperature-dependent particle density, slip velocity, Stefan flow and evaporation. Figure 5 shows this comparison of the temperature profiles between the boundary layer model and the particle model. The solid lines are obtained with the Nusselt and Sherwood correlations as proposed by Equations (17(17) $\mathrm{N}\mathrm{u}=2+\left(0.02\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{1/2}+0.33\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{2/3}\right)\mathrm{\Phi },$(17) , 18(18) $\mathrm{S}\mathrm{h}=2+\left(0.02\mathrm{S}{\mathrm{c}}^{1/3}\mathrm{R}{\mathrm{e}}^{1/2}+0.33\mathrm{S}{\mathrm{c}}^{1/3}\mathrm{R}{\mathrm{e}}^{2/3}\right)\mathrm{\Phi },$(18) ). Due to the high initial particle temperature, the particle burns in the diffusion-limited regime and then reaches a maximum temperature of around $2500\mathrm{K}$ when $\mathrm{F}\mathrm{e}$ is fully oxidized into $\mathrm{F}\mathrm{e}\mathrm{O}$. Two plateaus are visible in the temperature profile. The first plateau is due to the transition of $\mathrm{F}\mathrm{e}(\mathrm{s})$ into $\mathrm{F}\mathrm{e}(\mathrm{l})$ and the second plateau due to the transition of $\mathrm{F}\mathrm{e}\mathrm{O}(\mathrm{l})$ into $\mathrm{F}\mathrm{e}\mathrm{O}(\mathrm{s})$. Good agreement between the boundary layer resolved model and the Lagrangian particle model is obtained, when using the $1/2$-film rule and the new proposed Nusselt and Sherwood correlations. Note that agreement is good for both the inert gases (argon or nitrogen), since the Nusselt and Sherwood correlations scales with the dimensionless Prandtl and Schmidt numbers.

Figure 5. Temperature profiles for a $50$ µm particle burning in (a) Nitrogen and (b) Argon with $X_{{\mathrm{O}}_2}=0.21$. Temperature profiles are compared between the boundary layer resolved model (markers) and the Lagrangian particle model (line). The solid lines are obtained with the Nusselt and Sherwood correlations as proposed by Equations (17)(17) $\mathrm{N}\mathrm{u}=2+\left(0.02\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{1/2}+0.33\mathrm{P}{\mathrm{r}}^{1/3}\mathrm{R}{\mathrm{e}}^{2/3}\right)\mathrm{\Phi },$(17) and (18(18) $\mathrm{S}\mathrm{h}=2+\left(0.02\mathrm{S}{\mathrm{c}}^{1/3}\mathrm{R}{\mathrm{e}}^{1/2}+0.33\mathrm{S}{\mathrm{c}}^{1/3}\mathrm{R}{\mathrm{e}}^{2/3}\right)\mathrm{\Phi },$(18) ), while the dashed line is obtained with the Ranz- Marshall correlations and the 1/3-rule.

In Section 3.1.2 we showed that the $1/2$-film rule is a wise choice. However, when modeling the combustion of micron-sized particles, the $1/3$-film rule is often used (Hazenberg and van Oijen 2021; Mi, Fujinawa, Bergthorson 2022; Ravi, de Goey, van Oijen 2022). The dashed lines in Figure 5 are obtained with the $1/3$-film rule and Ranz-Marshall correlations. It can be seen that one underestimates the time to maximum temperature when using these settings.

In Table 1 the time to maximum temperature, the maximum temperature and the amount of iron evaporation as function of $X_{{\mathrm{O}}_2}$ for a $50$ µm particle obtained with the different models are stated. For the Lagrangian particle model, we distinguish again between the new proposed method (new correlation and $A_{\mathrm{f}}=1/2$) and the previous used methods (Ranz-Marshall and $A_{\mathrm{f}}=1/3$). We see that $t_{\mathrm{m}\mathrm{a}\mathrm{x}}$ decreases with an increasing oxygen concentration, while $T_{\mathrm{m}\mathrm{a}\mathrm{x}}$ increases. Since the particle is in a diffusion-limited regime up to the maximum temperature, an increasing oxygen concentration will enhance the oxygen flux toward the particle and therefore results in a lower $t_{\mathrm{m}\mathrm{a}\mathrm{x}}$. At the same time, the particle temperature will increase with increasing oxygen concentration since the heat release rate is larger than the heat loss rate. Excellent agreement between the boundary-layer resolved model and the Lagrangian point particle model with the new proposed settings can be seen. When using the Ranz-Marshall correlation and $1/3$-film rule one consistently underestimates the time to maximum temperature, while it obtains good agreement with the maximum temperature. Therefore, it is of importance to use the $1/2$-film rule to model the mass and heat transport in a Lagrangian particle model.

Table 1. Comparison between the boundary-layer resolved model (BLR) and the Lagrangian particle model (0D) for a $50$ µm particle burning in ${\mathrm{O}}_2$ – ${\mathrm{N}}_2$. We distinguish between the new proposed method (new correlation and $A_{\mathrm{f}}=1/2$) and the previous used methods (Ranz-Marshall and $A_{\mathrm{f}}=1/3$).

$X_{{\mathrm{O}}_2}$	Model	$t_{\mathrm{m}\mathrm{a}\mathrm{x}}$ [ms]	$T_{\mathrm{m}\mathrm{a}\mathrm{x}}$ [K]	$\mathrm{F}\mathrm{e}$ evap [%]
$0.13$	BLR	$38.9$	$1898$	$2.6\cdot {10}^{-3}$
	$1/2$ film	$38.1$	$1883$	$2.5\cdot {10}^{-3}$
	$1/3$ film	$31.8$	$1875$	$2.4\cdot {10}^{-3}$
$0.17$	BLR	$27.9$	$2251$	$5.2\cdot {10}^{-2}$
	$1/2$ film	$27.6$	$2231$	$4.8\cdot {10}^{-2}$
	$1/3$ film	$22.9$	$2233$	$4.8\cdot {10}^{-3}$
$0.21$	BLR	$21.3$	$2524$	$0.32$
	$1/2$ film	$21.2$	$2508$	$0.32$
	$1/3$ film	$17.6$	$2518$	$0.34$
$0.26$	BLR	$16.2$	$2784$	$1.24$
	$1/2$ film	$16.0$	$2769$	$1.30$
	$1/3$ film	$13.2$	$2783$	$1.41$
$0.31$	BLR	$13.1$	$2967$	$2.85$
	$1/2$ film	$12.5$	$2954$	$2.9$
	$1/3$ film	$10.3$	$2967$	$3.15$

Effect of modeling strategies

Next, the effect of different modeling strategies on the time to maximum temperature of a $50$ µm iron particle in an ${\mathrm{O}}_2$-${\mathrm{N}}_2$ is investigated. We will systematically elaborate the model by including a temperature-dependent particle density, slip velocity and Stefan flow. In Table 2 the 5 different cases are listed.

Table 2. Summary of test cases used to investigate the effect on $t_{\mathrm{m}\mathrm{a}\mathrm{x}}$.

Case	$A_{\mathrm{f}}$	$\rho =\mathrm{f}(T)$	Slip velocity	Stefan flow
1	1/2	no	no	no
2	1/2	yes	no	no
3	1/2	yes	yes	no
4	1/2	yes	yes	yes
5	1/3	no	no	no

Figure 6 shows the time to maximum temperature as a function of the oxygen concentration for the five different cases. It can be seen that adding more physical phenomena improves the agreement with the experimentally obtained results of Ning et al. (2021a, 2022). In Thijs et al. (2022) we also compared the burn times of $d_{\mathrm{p}}=40$ m particles, where we obtained excellent agreement. Here we will only focus on the results of $d_{\mathrm{p}}=50$ µm. We will use $X_{{\mathrm{O}}_2}=0.21$ as the reference composition to estimate the effect on $t_{\mathrm{m}\mathrm{a}\mathrm{x}}$.

Figure 6. Effect of different modeling strategies on $t_{\mathrm{m}\mathrm{a}\mathrm{x}}$, the cases are summarized in Table 2.

Using a temperature-dependent particle density instead of a constant particle density and neglecting slip and Stefan flow, $t_{\mathrm{m}\mathrm{a}\mathrm{x}}$ will decrease by around $7.5\mathrm{\%}$. The density of $\mathrm{F}\mathrm{e}$ and $\mathrm{F}\mathrm{e}\mathrm{O}$ decrease as a function of temperature, which means that the particle volume increases more than in the case of a constant density. A larger particle means a larger reacting surface area resulting in a smaller $t_{\mathrm{m}\mathrm{a}\mathrm{x}}$.

Considering the effect of slip velocity on heat and mass transfer will decrease $t_{\mathrm{m}\mathrm{a}\mathrm{x}}$ by another $6\mathrm{\%}$. Slip will enhance the mass and heat transfer, resulting in larger $\mathrm{N}\mathrm{u}$ and $\mathrm{S}\mathrm{h}$ numbers. Since the terminal velocity of particles scales with the particle diameter, the effect of slip will be less when smaller particles are considered and larger when larger particles are considered.

Finally, in case 4 it is shown that modeling the oxygen induced Stefan flow will decrease $t_{\mathrm{m}\mathrm{a}\mathrm{x}}$ by another $12\mathrm{\%}$. Similar as slip, Stefan flow will enhance the mass transfer toward and heat transfer away from the particle. Note that the effect of Stefan flow increases when the oxygen concentration increases. For $X_{{\mathrm{O}}_2}=0.13$, $t_{\mathrm{m}\mathrm{a}\mathrm{x}}$ decreases by $7.5\mathrm{\%}$ with respect to case 3, while it decreases by $20\mathrm{\%}$ with respect to case 3 for $X_{{\mathrm{O}}_2}=0.36$.

Ravi, de Goey, and van Oijen (2022) showed that by using the $1/3$-film rule, good agreement is obtained with the experiments of Ning et al. (2021a, 2022). However, Ravi, de Goey, and van Oijen (2022) neglected the effect of slip and Stefan flow ($\mathrm{N}\mathrm{u}=2$), and did not take into account a temperature-dependent particle density. In Figure 7 we show the effect of the film rule on the temperature profile for a $50$ µm particle burning $X_{{\mathrm{O}}_2}=0.21$. Here, the results are shown for cases 1, 4 and 5. When comparing case 1 and case 5, it can be seen that using $A_{\mathrm{f}}=1/3$ will decrease the time to maximum temperature with around $13\mathrm{\%}$. Furthermore, it can be seen that when neglecting physical phenomena like slip and Stefan flow, but using $A_{\mathrm{f}}=1/3$, one will still obtain reasonable agreement with the experiments of Ning et al. (2021a, 2022). However, different errors cancel each other. Note that the available experimental data is limited, and more independent experiments of single iron particle combustion in normal air are needed to further validate the model.

Figure 7. Effect of reference temperature on the temperature profile for a $50$ µm particle burning $X_{{\mathrm{O}}_2}$. The results are shown for case 1, 4 and 5, the cases are summarized in Table 2.

Conclusions and future directions

In this paper, a boundary layer resolved model is used to improve a 0D Lagrangian particle model for the first stage of diffusion-limited iron particle combustion. First, the boundary layer resolved model is used to derive new Nusselt and Sherwood correlations which apply to iron-particle combustion conditions. We added an extra term that accounts for the temperature difference between the particle and the gas, while we derive transport properties with the $1/2$-film temperature.

The boundary layer resolved model is used to validate the Lagrangian point particle model. Good agreement between the boundary layer resolved model and the Lagrangian particle model is obtained, when using the new proposed Nusselt and Sherwood correlations. However, when using the Ranz Marshall correlation and the $1/3$-film rule, one underestimates the time to maximum temperature of single metal particles.

With the point particle model, the effect of different modeling strategies on the time to maximum temperature of a $50$ µm laser-ignited iron particle is investigated. We systematically elaborated the model by including a temperature-dependent particle density, slip velocity and Stefan flow. Including all these effects decreases the time to maximum temperature by around $25\mathrm{\%}$, and results in a good agreement with experiments. It is also shown that if one neglects physical phenomena like slip and Stefan flow, but uses the $1/3$-film rule instead of the $1/2$-film rule, still reasonable agreement is obtained with experiments, but different errors cancel each other.

In future research, the point particle model can be further improved. Other stages of combustion which are not diffusion limited could be investigated. The model could be elaborated with a two-film model to take into account the oxidation of the evaporated iron in the boundary layer. Furthermore, the Knudsen effect can be included to enable investigating the transition regime for smaller particles.

Acknowledgments

We would like to thank X. Mi for his valuable feedback. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under Grant Agreement no. 884916. and Opzuid (stimulus / European Regional Development Fund) Grant agreement No. PROJ-02594.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the H2020 European Research Council [884916]; Opzuid [PROJ-02594].