A quantitative theory for heterogeneous combustion of nonvolatile metal particles in the diffusion-limited regime
Daoguan Ning Andreas Dreizler


































ning@rsm.tu-darmstadt.de


Reactive Flows and Diagnostics, Department of Mechanical Engineering, Technical Univeristy of Darmstadt, Otto-Berndt-Str. 3, 64287 Darmstadt, Germany


The paper presents an analytical theory quantitatively describing the heterogeneous combustion of nonvolatile (metal) particles in the diffusion-limited regime. It is assumed that the particle is suspended in an unconfined, isobaric, quiescent gaseous mixture and the chemisorption of the oxygen takes place evenly on the particle surface. The exact solution of the particle burn time is derived from the conservation equations of the gas-phase described in a spherical coordinate system with the utilization of constant thermophysical properties, evaluated at a reference film layer. This solution inherently takes the Stefan flow into account. The approximate expression of the time-dependent particle temperature is solved from the conservation of the particle enthalpy by neglecting the higher order terms in the Taylor expansion of the product of the transient particle density and diameter squared. Coupling the solutions for the burn time and time-dependent particle temperature provides quantitative results when initial and boundary conditions are specified. The theory is employed to predict the burn time and temperature of micro-sized iron particles, which are then compared with measurements, as the first validation case.
The theoretical burn time agrees with the experiments almost perfectly at both low and high oxygen levels. The calculated 
particle temperature matches the measurements fairly well at relatively low oxygen mole fractions,  whereas the theory overpredict the particle peak temperature due to the neglect of evaporation and the possible transition of the combustion regime.



Analytical solution Burn time Temperature evolutionHeterogeneous combustion Nonvolatile particle Stefan flow

















§ NOVELTY AND SIGNIFICANCE STATEMENT

For the first time, we present a comprehensive and quantitative analytical theory elucidating the heterogeneous combustion of nonvolatile (metal) particles in the diffusion-limited regime.
This novel theoretical model exhibits a remarkable capacity for quantitative prediction, obviating the need for supplementary information from numerical simulations or experimental data.
The derivation process of analytical solutions for burn time and time-dependent particle temperature from conservation equations is elaborated, offering transparency and insight into the model's foundations.
To demonstrate the practical utility of the theory, we apply it to analyze the combustion of iron particles, providing valuable mathematical perspectives on the underlying processes. The model's predictions for burn time and temperature align closely with experimental results, offering a partial validation of the theory within the realm of its applicable assumptions. This pioneering work contributes a robust and versatile analytical framework, advancing our understanding of diffusion-limited combustion phenomena of nonvolatile particles.



§ CREDIT AUTHORSHIP CONTRIBUTION STATEMENT

D. Ning: Conceptualization, Investigation, Methodology, Writing – original draft, review & editing.
A. Dreizler: Funding acquisition, Project administration, Writing – review & editing.



§ INTRODUCTION

The combustion of nonvolatile metal particles or powders presents a promising strategy for achieving carbon-free conversion of clean energy on demand. Unlike the vapor-phase combustion observed in volatile droplets or the heterogeneous combustion of coal and biomass particles, where both mass and diameter typically decrease during combustion, nonvolatile metal particles exhibit a distinctly burning behavior. As the nonvolatile burning process progresses, the condense-phase reaction product accumulates within the particle, leading to an increase in both particle mass and diameter.
This unconventional combustion feature deviates significantly from the established theories developed for vapor-phase or heterogeneous combustion of droplets and volatile particles, such as the classical d^2-law <cit.>, which originated from the vaporization of liquid droplets. When applied to the scenario of nonvolatile particle combustion, these theories lack a rigorous physical foundation.
Recognizing the deficiency in understanding the combustion of nonvolatile particles, recent years have witnessed a surge in research efforts. Iron particles, in particular, have garnered substantial attention, with numerous studies aimed at unraveling the intricacies of their combustion behavior. Other metals, like aluminium and silicon, are also expected to burn heterognerously at elevated pressures and reduced oxygen mole fractions, where the boiling points of the metals exceed their combustion temperatures.

In experimental measurements, multiple concepts have been proposed to investigate the fundamentals of single iron particle combustion. Wright et al. <cit.> measured the burn time of single iron particles by injecting single particles, ranging in diameter from 30 to 60µm, into preheated O_2/Ar mixtures at 1000K and atmospheric pressure. It is oberved that the particle burn time (defined as the time to peak luminosity) reduces as the oxygen concentration increases. Later, Ning et al. <cit.> further quantified the burn time of micron-sized iron particles, ignited by a focused laser beam and subsequently burning at room temperature. It is inferred that the combustion process from ignition to the peak temperature is limited by external oxygen diffusion, supported by the inverse proportionality between the burn time and the oxygen mole fraction in the ambient. Using the same particle generation and ignition apparatus, the evolution of the surface temperature and diameter during iron particle combustion were also quantified by Ning et al. <cit.>. The measurement revealed that particle diameter increases almost linearly with elevating particle temperature, reaching a plateau during the cooling phase. This observation suggests that a relatively rapid oxidation rate is sustained throughout the particle burn time, with further oxidation occurring at a significantly slower rate after the particle peaks in temperature. In addition to laser ignition accomplished at low gas temperature, the combustion of iron particles in hot environments provided by the exhaust of lean-premixed gaseous flames <cit.> or a drop-tube furnace <cit.> has been more extensively investigated. It is commonly observed that as the oxygen mole fraction increases, the burn time shortens, and simultaneously, the maximum particle temperature rises until a moderate oxygen level. These trends align with typical characteristics of particle combustion in the diffusion-dominated regime (but unnecessarily the diffusion-limited regime).


In numerical simulations, investigations into single iron particle combustion have been conducted from various perspectives, including particle-resolved modeling <cit.>, the point-particle approximation <cit.>, and molecular-scale surface chemisorption <cit.>. Based on the assumption of the diffusion-limited combustion regime, Thijs et al. <cit.> performed quantitative modeling of burn time and temperature for laser-ignited iron particles, corresponding to experimental setups <cit.>.  The simulated burn time aligns well with measurements across a wide range of oxygen levels. However, calculated particle peak temperatures only match experimental results at relatively low oxygen mole fractions. At high oxygen levels, the numerical simulation tends to overestimate the experimental outcome. To account for this discrepancy, Fujinawa et al. <cit.> speculated that at relatively high oxygen levels, the rate-limiting mechanism for iron particle combustion may transition to the diffusion of ions in the liquid oxide layer.
Employing reactive molecular dynamic simulations, Thijs et al. <cit.> demonstrated that the continuum model, constructed for the diffusion-limited regime, effectively describes the combustion behavior of iron particles with dimensions on the order of tens of microns up to the peak temperature not exceeding approximately 2500K. This finding indicates that conservation equations with the continuum assumption remain applicable for the theoretical analysis of nonvolatile particle combustion, provided that the particle size is significantly larger than the mean free path, and the mass accommodation coefficient of the oxygen is close to unity.

In theoretical frameworks, analytical solutions offer a comprehensive overview of general problems, complementing experiments and numerical simulations to enhance our understanding. Unfortunately, the quest for justified analytical models with predictive capability and quantitative accuracy for nonvolatile particle combustion remains elusive, impeded by the complexity of solving the nonlinear differential equations governing the problem.
Bidabadi and his collaborators  <cit.> attempted theoretical analysis of single iron particle combustion using several mathematical approaches. However, these theoretical models were formulated on the basis of self-contradictory assumptions: 1) the oxidation rate, described by an Arrhenius-like expression, was assumed to be proportional to the particle surface area throughout the entire combustion process, a justification applicable only to the kinetic-limited regime; 2) the burn time was evaluated using an analytical solution derived for the lifetime of vapor-phase droplet combustion in the diffusion-limited regime <cit.>. Consequently, the results obtained from these models lack a rigorous physical justification. 
Accounting for the variation of particle diameter during oxidation, Hazenberg <cit.> derived an analytical solution for the burn time of nonvolatile particles in the diffusion-limited combustion regime. While the expression explicitly describes the dependence of burn time on other parameters, 
for the mass diffusivity of oxygen, it needs estimation from one-to-one numerical simulations or using measured particle temperatures from experiments <cit.>. Moreover, the model becomes less rigorous with increasing oxygen levels due to the neglect of Stefan flow. Therefore, there is a need to further develop a comprehensive analytical theory for quantitatively predicting the combustion characteristics of nonvolatile particles in the diffusion-limited regime. This study aims to uncover the underlying objectives in this regard.

The paper is structured as follows. In <ref>, the interfacial heat and mass transfer rates are derived from the conservation equations of the gas phase considering relevant boundary conditions. Subsequently, the burn time of nonvolatile particle combustion in the diffusion-limited regime is formulated explicitly. 
The transient solution for the particle temperature is analytically solved using Taylor expansion and appropriate approximations. The coupling of the solutions for the burn time and time-averaged temperature of the particle yields quantitative results for defined boundary and initial conditions.
In <ref>, the theory is applied to analyze the combustion of iron particles, considering various effects of particle and gas-phase properties. The theoretical results are then compared with experiments to validate the proposed theory. Finally, the main conclusions drawn from this study are summarized in <ref>.









§ MATHEMATICAL FORMULATION AND SOLUTIONS


The problem being addressed is illustrated in <ref>. An isolated, spherical metal particle is suspended in a quiescent, isotropic, unconfined gaseous mixture. The ignition of the pure metal particle is simplified by initiating it at its ignition temperature. The subsequent nonvolatile combustion process is assumed to be limited by the external diffusion of the oxygen (i.e., the diffusion-limited regime), and the structure of the particle remains spherically symmetric, facilitating even chemisorption of the oxygen on the particle surface. 
During combustion, the heat generated by exothermic oxidation is transported outward, while the oxygen is transported inward. Meanwhile, the mass of the particle increases due to the incorporation of the oxygen, resulting in a growth of the particle size. As the particle temperature may exceed the melting points of the metal and oxide, the solid particle may transform into a liquid droplet consisting of both liquid metal and liquid oxide. Nevertheless, it will still be referred to as a particle in this context. For the sake of simplicity, the phase transitions of the particle are neglected, adopting constant specific heat. The temperature, T_g, and composition, Y_i, of the gas vary between the particle surface and the far-field as r →∞, where r is the radial coordinate. The combustion time and temperature history of the particle will be solved analytically.





 §.§ Interfacial heat and mass transfer rates

 During the combustion of non-volatile fuel particles, gas-phase reactions are absent and the combustion product is formed on the particle surface. Therefore, employing the quasi-steady-state assumption, the gas phase conservation equations for mass, the mass fraction of the oxygen, and energy are written in a spherical coordinate system as

    ∂/∂ r(r^2 ρ_gu) = 0,


    ∂/∂ r(r^2 ρ_guY_O) - ∂/∂ r(r^2  ρ_g𝒟_O∂ Y_O/∂ r) = 0,


    ∂/∂ r(r^2 ρ_gu  h_g) - ∂/∂ r(r^2 λ_g∂ T_g/∂ r) - ∂/∂ r( ∑_i=1^Nr^2 ρ_g𝒟_i h_i  ∂ Y_i/∂ r)= 0,

where ρ, u, Y, 𝒟, h, and λ represent the density, velocity, mass fraction, mass diffusion coefficient, specific enthalpy, and conductivity, respectively.  The subscripts, g, O, and i denote the gas phase, oxygen, and individual species, respectively. N is the number of species in the gas phase.
The boundary conditions at the particle surface are

    [left=r=r_p:]align
            ρ_g 𝒟_O ∂Y_O/∂r|_s -ρ_g u_s Y_O,s = ṁ_O/4πr_p^2,
        

            Y_O = Y_O,s,
        

            λ_g ∂T_g/∂r|_s =  h_T (T_g,∞ - T_p) = Q̇_conv/4πr_p^2, 
        

            T_g  = T_p,
        
    

with ṁ_O being the mass consumption rate of the oxygen, h_T the heat transfer coefficient, and Q̇_conv the convective heat transfer rate between the particle and the gas phase. The boundary conditions at the far-field read 

    [left=r→∞: ]align
        Y_O = Y_O,∞,

        T_g  = T_g, ∞.
    

To derive the interfacial heat and mass transfer rates, we will assume that the thermo-physical properties remain constant. This assumption is found to be satisfactory when the properties are evaluated at a virtual film layer with the reference conditions:

    T_f := T_p + A_f(T_g,∞-T_p),

    and

    Y_O,f := Y_O,s + A_f(Y_O,∞-Y_O,s),

with A_f being the weighting coefficient, commonly taken as 1/3 <cit.>. 

We will first derive the interficial mass transfer rate of the oxygen. Since only the oxygen 
has a net mass flux,
integrating (<ref>) from r_p to ∞ gives

    r^2 ρ_gu = -ṁ_O/4π,

where the negative sign indicates that ṁ_O is inward.
Substituting (<ref>) into (<ref>) gives


    ∂/∂ r(-ṁ_O/4πY_O - r^2  ρ_g𝒟_O∂ Y_O/∂ r) =0.

Integrating <ref> from r_p to r, and applying the boundary condition (<ref>) yields

    4π r^2 ρ_g𝒟_Od Y_O/d r = ṁ_O(1- Y_O).

Integrating <ref> from r to ∞ yields

    ln(1-Y_O,∞/1-Y_O(r)) = - 1/rṁ_O/4πρ_g𝒟_O.

Applying the boundary condition (<ref>) to the above equation (<ref>) gives the first expression for ṁ_O:
 
    ṁ_O = -4π r_pρ_g𝒟_Oln(1+B_M),

 with B_M being the Spalding mass transfer number of the oxygen, defined analogously to that of vaporized fuel <cit.> as
 

    B_M := Y_O,∞-Y_O,s/Y_O,s-1.

We introduce the Sherwood number, Sh,  defined as

    Sh := 2ln( 1+B_M)/B_M,

which inherently takes the Stefan flow effect into account. Except for the atmosphere with (nearly) pure oxygen <cit.>, in the diffusion-limited combustion regime:

    Y_O,s≈ 0  ⇒ B_M≈ -Y_O,∞.

Substituting (<ref>) and (<ref>) into (<ref>) 
yields

    ṁ_O    = Shπ d_pρ_g𝒟_O  Y_O,∞
             = -2 π d_pρ_g𝒟_Oln(1-Y_O,∞).

The explicit expression of the interficial heat transfer rate, Q̇_conv, can be derived by following a similar procedure as for ṁ_O.
We consider the facts that the net mass flux of each inert  species is zero:

    r^2 ρ_gu  Y_i  -  r^2 ρ_g𝒟_i∂ Y_i/∂ r = 0, i∈ [1,N] & i≠ i_O


We consider the facts that the net mass flux of all inert species is zero:

    r^2 ρ_gu∑_i=1, i≠O^N Y_i  -  ∑_i=1, i≠O^Nr^2 ρ_g𝒟_i∂ Y_i/∂ r = 0,


and the specific enthalpy of the gaseous mixture encompasses those of individual species:
 
    h_g = ∑_i=1^N h_i  Y_i = h_O Y_O +  ∑_i=1, i≠ i_O^N  h_i Y_i.

Substituting <ref> into <ref> yields

    ∂/∂ r( - ṁ_O/4π h_O - λ_g r^2 ∂ T_g/∂ r) = 0.

Integrating <ref> from r_p to r and applying the boundary condition (<ref>) gives

    4 π r^2 λ_gd T_g/d r = ṁ_Oc_p,O( T_p -  T_g +  Q̇_conv/ṁ_O c_p,O),

where the relation h = c_pT has been used. Integrating <ref> from r to ∞ gives T_g as a function of the radial distance:

    ln(T_p-T_g,∞  + Q̇_conv/ṁ_O c_p,O/T_p-T_g(r) + Q̇_conv/ṁ_O c_p,O) = -1/rṁ_Oc_p,O/4 πλ_g

At r=r_p, <ref> provides the second expression for ṁ_O:

    ṁ_O = - 4π r_pλ_g/c_p,Oln(1+B_T),

where B_T is the Spalding heat transfer number in our case, defined as

    B_T := ṁ_Oc_p,O/Q̇_conv(T_p - T_g,∞).

Substituting <ref> into <ref> yields


    Q̇_conv = 4 π r_pλ_gln(1+B_T)/B_T  (T_g,∞ - T_p).

We introduce the Nusselt number, Nu, defined as

    Nu:=2ln(1+B_T)/B_T.

Substituting <ref> into <ref> gives

    Q̇_conv = Nuπ d_pλ_g (T_g,∞ - T_p)

Finally, the two expressions for ṁ_O, i.e., <ref>, relate B_T and B_M:

    1+ B_T=(1+B_M)^1/Le_Oc_p,O/c_p,g,

where Le_O = λ_g/(ρ_g c_p,g𝒟_O) is the Lewis number of the oxygen at the film layer. In the follows, the mass diffusivity of the oxygen and the thermal conductivity and specific heat of the gaseous mixture  evaluated at the film layer will be denoted as  (ρ_g𝒟_O)_f, λ_f, and c_p,f, respectively.
As aforementioned, in <ref>, ρ_gD_O,m and λ_g should be evaluated at the film layer:


    {   (ρ_gD_O,m)_f = ρ_g (T_f, Y_i,f) · D_O,m(T_f, Y_i,f)
       λ_g,f = λ_g(T_f,Y_i,f)
        .







 §.§ Analytical solution for particle burn time

During non-volatile combustion, the particle mass increases by consuming the oxygen,  which is described as

    dm_p/dt = π/6d(ρ_p d_p^3)/dt = ṁ_O,

with m_p and ρ_p being the time-dependent mass and density of the particle, respectively. As derived in <ref> (see Supplementary Material), ρ_p is a function of the density of the metal, ρ_m, and that of the metal oxide, ρ_mo:

    ρ_p =ρ_1 + ρ_2(d_0/d_p)^3,

with

    ρ_1 := ρ_mo/s+(1-s)ϱ,

and

    ρ_2 := (ρ_m-ρ_mo)s/s+(1-s)ϱ,

where s is the stoichiometric oxide-to-oxygen mass ratio, and ϱ is the oxide-to-metal density ratio.
Substituting <ref> into <ref> yields

    d(ρ_1d_p^3 + ρ_2d_0^3 )/dt = -12 d_p(ρ_g𝒟_O)_fln( 1-Y_O,∞).

Because ρ_2d_0^3 is constant, <ref> simplifies as

    ρ_1 d_pd( d_p)/dt = -4 (ρ_g𝒟_O)_fln( 1-Y_O,∞).

Integrating <ref> from t_0 (t_0=0: d_p = d_0) to t (t_0⩽ t ⩽ t_b) gives the particle diameter evolution:

    d_p^2 = d_0^2 -8 ( ρ_g𝒟_O)_fln( 1-Y_O,∞)/ρ_1t.

where d_0 is the initial particle diameter.
Since the particle size grows during combustion, we define an expansion coefficient:

    K := -8 ( ρ_g𝒟_O)_fln( 1-Y_O,∞)/ρ_1.

Consequently, <ref> becomes

    d_p^2 = d_0^2 + Kt.

<Ref> resembles the size evolution of a shrinking droplet during vaporization (i.e., d^2 = d_0^2 - β t with β the evaporation coefficient <cit.>). 
Substituting the final condition of combustion  (i.e., t=t_b: d_p = d_1 and d_1 is the final particle diameter) into <ref> gives the analytical solution for the burn time:

    t_b = ρ_1( d_0^2-d_1^2 )/8 ( ρ_g𝒟_O)_fln( 1-Y_O,∞)
        = ρ_1( 1-ϵ^2)d_0^2/8 (ρ_g𝒟_O)_fln( 1-Y_O,∞),

with ϵ:=d_1/d_0 defining the diameter expansion ratio.
<Ref> indicates a d_0^2-dependence of the burn time in the diffusion-limited regime.
Besides, using the time-dependent density  and diameter of the particle, given in <ref> and <ref> respectively, the transient increase of the particle mass is determined as

    m_p = π/6[ρ_1 (d_0^2 + Kt)^3/2 + ρ_2 d_0^3  ].

Substituting <ref> into <ref> and normalize the result by the initial particle mass gives the fraction of remaining unburned metal in the particle during combustion:

    m_m/m_0 = (1-s)[ρ_1/ρ_m(1 + K/d_0^2t)^3/2 + ρ_2/ρ_m]  + s.





 §.§ Analytical solution for transient particle temperature

For some metals (e.g., Fe),  a majority of the combustion process is in the liquid phase, where the specific heats of the metal and the oxide are close and independent of temperature <cit.>.  Therefore, the specific heat of the particle does not change significantly during combustion.
Assuming a constant specific heat of the particle, the  governing equation of the particle temperature is expressed as

    c_p,pd(m_pT_p)/dt = 
        c_p,pm_pdT_p/dt + c_p,pT_pdm_p/dt
        =
        Q̇_conv + Q̇_chem + Q̇_rad,

with c_p,p being the specific heat of the particle, 
m_p = ρ_pπ d_p^3/6 the particle mass, Q̇_conv the convective heat transfer rate, Q̇_chem the chemical heat release rate, and Q̇_rad the radiative heat loss rate. Since only nonvolatile combustion is considered in this work, the heat loss via evaporation is neglected. Although evaporation may become important when the particle temperature approaches the boiling point of the metal or the oxide, it is beyond the scope of strictly-defined nonvolatile combustion and makes the analytical solution seemingly impossible to complete. Consequently, the results derived here may be limited to scenarios where the maximum particle temperature is considerably below the boiling points.

After applying the chain rule, <ref> writes

    c_p,pm_pdT_p/dt + c_p,pT_pdm_p/dt  = Q̇_conv + Q̇_chem + Q̇_rad + Q̇_vap.


The particle mass reads

    m_p = ρ_pπ d_p^3 /6,

where ρ_p and d_p are the time-dependent particle density and diameter, respectively.

Combining <ref> gives the rate of particle mass change in the diffusion-limited combustion regime: 

    dm_p/dt  = ṁ_O =  Sh (ρ_g𝒟_O)_fπ d_p Y_O,∞,

Rewriting <ref>, the convective heat transfer rate reads

    Q̇_conv = Nuλ_fπ d_p (T_g,∞-T_p).


The chemical heat release rate is written as

    Q̇_chem = ṁ_O_2ϕΔ h,

where ϕ is the stiochiometric fuel-to-oxygen mass ratio, and Δ h is the combustion enthalpy per gram fuel. The radiative heat transfer rate is described by the Stefan-Boltzmann law:

    Q̇_rad = - π d_p^2 εσ(T_p^4 - T_0^4),



where ε is the emissivity of the particle surface and σ is the Stefan-Boltzmann constant.
Substituting <ref> into <ref> and rearranging it yields:

    c_p,pρ_p d_p^2/6dT_p/dt   = Nuλ_fT_g,∞ + Sh (ρ_g𝒟_O)_f Y_O, ∞ϕΔ h 
        - ( Nuλ_f + Sh (ρ_g𝒟_O)_f Y_O, ∞ c_p,p) T_p
        - d_pεσ(T_p^4 - T_0^4).

With  ρ_p = ρ_1 + ρ_2(d_0/d_p)^3 (<ref>) and d_p^2 = d_0^2 + Kt (<ref>), the following equation is obtained:

    ρ_p d_p^2 = ρ_1 d_p^2  + ρ_2 d_0^3/d_p = ρ_1 (d_0^2 + Kt) + ρ_2 d_0^3/√(d_0^2 + Kt).

In order to solve <ref> analytically, <ref> needs to be approximated. To this end, we consider the Taylor expansion at t=0 for the second term on the right-hand side of <ref>, which expresses as

    d_0^3/√(d_0^2 + Kt) = d_0^2/√(1 + (K/d_0^2)t) =  ∑_n=0^∞( K/d_0^2 t )^n (2n-1)!!/n!(-2)^n d_0^2.

Recalling <ref>, one can find the following inequality valid for many metals with examples given in <ref>, where thermal expansion is neglected:

    |K/d_0^2 t | = d_p^2/d_0^2 - 1 ≤d_1^2/d_0^2 - 1 = (ρ_mℳ_mo/ρ_moℳ_m1/x)^2/3 - 1 < 1, ∀ t ∈ [0,t_b],

in which ℳ_m and ℳ_mo are the molar mass of the metal and the metal oxide, respectively, x is the numbers of metal atom in the chemical formula of the metal oxide, M_xO_y. Therefore, higher order terms in <ref> become smaller.
By neglecting 𝒪(t^2) 
[It is possible to find a solution when the second order is included but the improvement is negligibly small and the expression is cumbersome. Therefore, we do not present it in the paper.], 
<ref> is approximated as

    d_0^3/√(d_0^2 + Kt)≈  d_0^2  - K/2t, ∀ t ∈ [0,t_b].


Therefore, we can neglect 𝒪(t^3) and <ref> is approximated as

    d_0^3/√(d_0^2 + Kt)≈  d_0^2  - K/2t + 3K^2/8t^2, ∀ t ∈ [0,t_b].


Then, <ref> is also approximated as

    ρ_p d_p^2    ≈ρ_1 (d_0^2 + Kt) + ρ_2 ( d_0^2  - K/2t)
        = (ρ_1 + ρ_2) d_0^2 + (ρ_1 - ρ_2/2)Kt

The maximum truncation errors in <ref> introduced by the neglect of higher order terms in the Taylor series, <ref>, reach at t_b, which do not exceed 5% as provided in <ref> for several sample metals.
Substituting <ref> into <ref> and reaggranging it yields

    (1 + 2ρ_1 -ρ_2/2(ρ_1 +ρ_2)d_0^2K t )dT_p/dt   ≈6Nuλ_fT_g,∞ + Sh (ρ_g𝒟_O)_f Y_O, ∞ϕΔ h/c_p,p(ρ_1 +ρ_2)d_0^2
         - 6Nuλ_f + Sh (ρ_g𝒟_O)_f Y_O, ∞ c_p, p/c_p,p(ρ_1 +ρ_2)d_0^2 T_p
        - 6 d_pεσ(T_p^4 - T_0^4)/c_p,p(ρ_1 +ρ_2)d_0^2

By defining constant parameters as follows:

    𝔅 := 2ρ_1 -ρ_2/2(ρ_1 +ρ_2)d_0^2K,



    𝔇 := 6Nuλ_fT_g,∞ + Sh (ρ_g𝒟_O)_f Y_O, ∞ϕΔ h/c_p,p(ρ_1 +ρ_2)d_0^2,



    𝔈 := 6Nuλ_f + Sh (ρ_g𝒟_O)_f Y_O, ∞ c_p, 𝓅/c_p,𝓅(ρ_1 +ρ_2)d_0^2,

and a non-linear function:

    𝐅 (T_p) := 6 d_pεσ(T_p^4 - T_0^4)/c_p,p(ρ_1 +ρ_2)d_0^2≈6εσ(T_p^4 - T_0^4)/c_p,p(ρ_1 +ρ_2)d_0

in which the approximation is obtained by replacing d_p with d_0. Alternatively, d_p can also be replaced by d_1. The difference in the solved maximum particle temperatures that will be discussed in <ref> is only around 1% between these two approximations, which can be regarded as negligible. The exact solution is bounded by those solved using these two approximations.
Substituting <ref> into
<ref> simplifies it as

    (1+𝔅t)dT_p/dt≈𝔇-𝔈T_p - 𝐅 (T_p).

Rearranging  <ref> and integrating both sides yields

    ∫dT_p/𝔇-𝔈T_p - 𝐅 (T_p)≈∫dt/1+𝔅t

To derive an explicit solution for T_p(t) from <ref>, we first neglect radiation, i.e., 𝐅 (T_p) = 0. Then the general solution of <ref> expresses as

    T_p(t) ≈ k_1(1 + 𝔅t)^-𝔈/𝔅 + 𝔇/𝔈,

with k_1 the constant of integration, which can be determined as

    k_1 = T_p,ign - 𝔇/𝔈,

using  the initial condition:

    t = 0: T_p(0) = T_p,ign.

Finally, the explicit solution of <ref> reads

    T_p(t) ≈(T_p,ign - 𝔇/𝔈)(1 + 𝔅t)^-𝔈/𝔅 + 𝔇/𝔈.












When radiation is considered, 
an explicit solution cannot be found but T_p can be expressed implicitly as

    t ≈1/𝔅[ exp( ∫_T_p,ign^T_p(t)𝔅dz/𝔇-𝔈z - 𝐅 (z)) - 1 ].

The integral in <ref> can be solved analytically using partial fraction decomposition, yielding

    ∫_T_p,ign^T_p(t)𝔅dz/𝔇-𝔈z - 𝐅 (z) = 𝔅∑_ω_iln(T_p,ign - ω_i ) - ln(T_p - ω_i )/4𝔉ω_i + 𝔈,

where ω_i is the complex roots of the quartic equation 𝔇-𝔈ω - 𝐅 (ω) = 0 and 𝔉  is a constant defined as

    𝔉:=6εσ/c_p,p(ρ_1 +ρ_2)d_0







 §.§ Coupling of solutions and computational procedure


The analytical solutions for particle burn time and temperature evolution have been derived separately. However they are conjugated via the mass diffusivity of the oxygen at the film layer, (ρ_g𝒟_O)_f. The mass diffusivity is unknown a priori since it is evaluated at the film layer temperature, T_f, whose determination needs the particle temperature as shown in <ref>. Since the particle temperature is transient, so is the film layer temperature. 

Considering the fact that (ρ_g𝒟_O)_f∝ T_f^0.7  for diatomic gases, which is close to be a linear correlation, the following approximate relation is obtained:

    (ρ_g 𝒟_O)_f(1-A_f)T_p + A_f T_g,∞

This approximation implies that, on the one hand, the variation of (ρ_g𝒟_O)_f during particle combustion is not significant owing to a constant far-field gas temperature. On the other hand, due to the linear proportionality, an averaged (ρ_g𝒟_O)_f could be evaluated at a time-averaged film layer temperature:


    <T_f> = (1-A_f)<T_p> + A_f T_g,∞

with <T_p> the time-averaged particle temperature:

    <T_p> = 1/t_b∫_0^t_bT_p(t)dt



To obtain quantitative results, 

the time-average temperature and gas composition at the film layer are determined iteratively using the procedure illustrated in <ref>.
The thermal physical and transport properties of the gaseous mixture and species at the film layer are evaluated using Cantera <cit.>.








§ APPLICATION TO IRON PARTICLE COMBUSTION



In this section, the theory developed above is applied to analyze the burning properties of micron-sized iron particles, which serves as a typical example of non-volatile, heterogeneous combustion at low to moderate oxygen concentrations. This choice is justified by the much higher boiling points of iron and iron oxides (>3000K at 1atm) compared to the peak temperature of the burning particles in air at room temperature (∼2500K <cit.>). Additionally, systematic experiments providing quantitative data have been conducted, offering an opportunity to validate the theory. 



 §.§ Iron-oriented assumptions and parameters


The basic assumptions, which were adopted to reasonably simplify the practical problem, leading to the successful derivation of the theory, such as constant thermophysical properties and homogeneous distributions of temperature inside the particle, remain unchanged. 
The main assumptions related to iron particle combustion are three-fold.
First, the ignition temperature is assumed to be 1100K, as suggested by the calculations using a state-of-the-art model <cit.>.
Second, the oxidation reaction is assumed to be Fe + O_2 →FeO, and combustion enthalpy per gram iron is determined by the NASA thermodynamic database. 
Third, the reaction stops upon completely conversion to FeO, and further oxidation is neglected.
The physical parameters used for performing the quantitative analysis of iron particle combustion are listed in <ref>.






 §.§ Theoretical results and discussion

In this section, sample results on the combustion characteristics of iron particles under different boundary and initial conditions will be presented. The mechanisms underling the theoretical results are explained. 







  §.§.§ Effects of particle size and surface radiation


<Ref> illustrates the quantitative solutions for the temperature and the fraction of unburnt iron of a  40 µm and a 60µm iron particles during combustion in the 21%O2/79%N2  at 300K, both with and without considering radiation. 
When the radiative heat loss is neglected, the two particles reach the same peak temperature of 2578K at approximately 14ms and 32ms, respectively. By the same time, iron has burnt out, as indicated by the disappearance of unburnt iron mass in <ref>(a). 
The independence of the particle peak temperature on the diameter arises because d_0 cancels out from the analytical solution for the particle peak temperature, as derived in <ref>. This solution is obtained by substituting the particle burn time solution, <ref>, into the radiation-excluded expression of particle temperature history, <ref>.

    T_p, max   ≈(T_p,ign - 𝔇/𝔈)
         ( 1 + (2ρ_1-ρ_2)(1-ϵ^2)/2(ρ_1+ρ_2))^-𝔈/𝔅
          + 𝔇/𝔈.

As time continues to elapse, the particle temperature still increases until reaching a plateau, corresponding to the steady-state solution for the particle temperature. This plateau represents the time-independent term, 𝔇/𝔈, in <ref>. The steady-state solution is achieved when the rates of the conductive heat loss and the chemical heat release are balanced.
However, this solution is nonphysical because it yields a negative value for the unburnt mass of iron. 
When radiation is considered, the dependence of the particle peak temperature on the particle diameter emerges, indicating that the larger particle has a lower peak temperature. The maximum temperatures of the small and large particles are 2456K and 2402K, about 60K and 90K lower than their non-radiative counterparts, respectively. Additionally, <ref> shows that the effect of radiative heat loss on particle burn times is negligible. The overall temporal evolution of the particle temperature, including both the physical and nonphysical solutions, clearly depicts that the maximum particle temperature is constrained by the burnout of fuel rather than heat losses.




The effect of radiation is more profound on the nonphysical steady-state solution of the particle temperature, making the maximum particle temperature is more close the the steady-state solution, 2601K and 2508K, for the small and the big particles, respectively. This will help us to explain the physical reason for the size dependence of the maximum particle temperature in the follows.
Since we do not have an explicit solution for the particle peak temperature when radiation is considered, an intuitive explanation for the size dependence of the particle peak temperature could be made approximately based on the the steady-state solution. In the steady-state, both the convective heat loss rate Q̇_conv and the chemical heat release rate Q̇_chem scale with the particle diameter as shown by <ref>, respectively. Whereas, the radiative heat loss scales with the particle diameter squared. Consequently, the steady-state particle temperature (i.e., T_p,ss) is related to the radiative heat loss as

    T_p, ss∝Q̇_rad/d_p = - π d_pεσ(T_p, ss^4 - T_0^4), ∀ T_p,ss>T_0.

The above expression suggests that the bigger the particle, the lower the steady-state particle temperature. 


The time history of the particle temperature and the fraction of unburnt iron shown in <ref> are re-plotted versus the time coordinated normalized by the initial particle diameter squared (i.e., t/d_0^2) in <ref>. 
Curves for different particle diameters overlap when radiation is neglected. As the initial particle diameter, d_0, drops out in 𝔇/𝔈 and 𝔈/𝔅, the effect of the initial particle diameter on the particle temperature, described by <ref>, appears only in the term:
 
    𝔅t = 2ρ_1 -ρ_2/2(ρ_1 +ρ_2)K ·t/d_0^2.

As the first factor of the right hand side of <ref> is independent of d_0, the particle temperature history scales with 1/d_0^2 when radiation is neglected. This scaling behavior is also valid for the fraction of unburnt iron, as described by <ref>. Even when accounting for radiation, the curves for 40 and 60µm coincide very well. For particle peak temperatures below approximately ∼2500K, radiation appears to have a weak impact on controlling the temperature evolution of burning iron particles, particularly for particles with relatively small diameters (e.g., <100µm). This scaling is consistent with the expanding d^2-law derived in <ref>. 

 


 
In other words, when the temperature trajectories of burning particles with different initial diameters, given by <ref>, 
In other words, when the temperature trajectories of burning particles with different initial diameters, given by <ref>, are plotted versus the time coordinated normalized by the initial particle diameter squared (i.e., t/d_0^2), they overlap as shown in fig. XXX. Besides, the last term of <ref> is the steady-state solution for the maximum particle temperature without considering radiation, which is due to the thermal balance between the conductive heat loss and heat release. 




  §.§.§ Effects of gas-phase pressure and temperature

<Ref> depicts the burn time and temperature history of a 50µm iron particle under varying gas-phase pressure, P_∞, and temperature, T_∞.
Increasing the ambient pressure from 1atm to 10atm exhibits no discernible impact on the diffusion-limited combustion process of the particle.   The independence of particle burn time from ambient pressure is attributed to the negligible influence of pressure on the mass diffusivity of the oxygen at the film layer, denoted as (ρ_g𝒟_O)_f, a critical factor controlling the duration of particle combustion in the diffusion-limited regime. The gas density, ρ_g, scales linearly with pressure and the binary diffusion coefficient of oxygen, 𝒟_O, is inversely proportional to pressure according to Chapman-Enskog theory <cit.>.

    This phenomenon can be explained through the relationship between gas density, ρ_g, and pressure, where the former is directly proportional to the latter. Meanwhile, the mass diffusion coefficient of the oxygen in a binary mixture (e.g., O2/N2), denoted as 𝒟_O, follows an inverse proportionality to pressure, as outlined in Chapman-Enskog’s theory <cit.>:

    𝒟_12 = 3/16√(2π k^3_BT^3/ℳ_12)/Pπσ^2_12Ω,

where  k_B denotes the Boltzmann constant, ℳ_12 = ℳ_1ℳ_2/(ℳ_1+ℳ_2) represents the reduced molecular mass for the species pair (subscripts 1 and 2 denote the two kinds of molecules in the gaseous mixture),
σ_12 is the reduced collision diameter, and Ω is the collision
integral. 

Additionally, the thermal conductivity of ideal gases remains independent of pressure, except at extremely low pressures where the mean free path approaches or surpasses the dimensions of fuel particles (Knudsen number Kn  1). Consequently, the governing equation for particle temperature (<ref>) and its solutions (<ref>), depicting the time history of the particle temperature, are unaffected by ambient gas pressure.



In contrast to pressure, the burn time of the particle shortens with increasing the gas-phase temperature, ending at a higher maximum temperature. <Ref>(a) summarizes the peak particle temperature, T_p,max, and burn time normalized by the initial particle diameter squared, t_b/d_0^2, as a function of the gas-phase temperature for 10 and 100µm iron particles burning in the 14%O2/86%N2 mixture. The calculation is conducted for particle peak temperatures below 2500K, ensuring negligible evaporation. As shown in <Ref>(a), the peak temperature shows an almost linear correlation with the gas-phase temperature. Given that the role of radiation is significantly weaker in determining the maximum particle temperature compared to convection, the linear dependence of the particle peak temperature on the far-field gas temperature can be derived by rearranging <ref>, yielding 





    T_p, max   ≈ (1-𝔊)Nu Le_Oc_p,f/Nu Le_Oc_p,f + Sh Y_O,∞c_p,pT_g,∞
       + (1-𝔊)Y_O,∞ϕΔ h/Nu Le_Oc_p,f + ShY_O,∞c_p,p
        + T_p,ign𝔊 ,

with the T_g,∞-independent constant 𝔊 defined as

    𝔊 := ( 1 + (2ρ_1-ρ_2)(1-ϵ^2)/2(ρ_1+ρ_2))^3/4NuLe_Oc_p,f + ShY_O,∞c_p,p/(2-ρ_2/ρ_1 )ln(1-Y_O,∞),


where Le_O and c_p,g are the Lewis number of the oxygen and the specific heat of the gaseous mixture evaluated at the film layer, respectively. In <ref>, all the other parameters are independent on the far-field gas temperature, T_g,∞, and thus the linearity between T_p,max and T_g,∞ is determined by the first term on the right hand side of the equation. 

In <ref>(b), it is observed that as the far-field gas temperature increases, the particle burn time decreases in proportion to T_f^-0.7. This behavior is attributed to the fact that the mass diffusivity of the oxygen at the film layer is directly proportional to T_f^0.7, which increases with the rise in ambient gas temperature. Additionally, the larger particle exhibits a larger value of t_b/d_0^2 at the same gas temperature. This can be attributed to increased heat loss through surface radiation. The heightened surface radiation leads to a reduction in the time-averaged particle temperature and, consequently, the time-averaged film-layer temperature, <T_f>. This reduction in <T_f> contributes to a decrease in the mass diffusivity of the oxygen, resulting in larger values of t_b/d_0^2.




  §.§.§ Effect of initial particle temperature


Laser heating is a commonly employed experimental technique for initiating combustion in metal particles. Despite the advantages of being able to flexibly adjust the gas composition and temperature in laser-ignition experiments, there is a tendency for the particle to rapidly reach a high-temperature state in the diffusion regime without undergoing the typical ignition process. Since oxidation during laser heating is considered negligible, the resulting laser-driven high-temperature state can be treated as an initial particle temperature for theoretical analysis. This section examines the impact of the initial particle temperature on the combustion characteristics of iron particles.



<Ref> shows the calculated time histories for different initial temperature of a 50µm iron particle burning in a mixture of 21%O2/79%N2 at 300K and atmospherical pressure. Despite a temperature difference of 400K between neighboring cases at the onset of combustion, the maximum particle temperatures at burnout are very close, with a difference of about 35K. This weak effect of the initial temperature on the maximum particle temperature is attributed to the combustion progress closely approaching a steady state, where the initial condition has been largely relaxed.
Moreover, the significant elevations of the initial particle temperature only induce very minor reductions (i.e., <10%) in the particle burn time. In the diffusion-limited regime, the insensitivity of the burn time to the initial particle temperature is due to the controlling parameter of particle burn time being the mass diffusivity of the oxygen ρ_g𝒟_O, evaluated using a time-averaged film layer temperature <T_f>. According to <ref>, the role of the initial particle temperature in determining ρ_g𝒟_O is significantly weakened by the two-fold averaging operations. 
For example, the calculated <T_f> in the cases of T_p,0=1100K and 1900K shown in <ref> are 1423K and 1638K, respectively. This means that an 73% variation in T_p,0 causes only a 15% change in <T_f>. Since ρ_g𝒟_O∝<T_f>^0.7, the difference in the mass diffusivity of the oxygen further reduces to 9.7%, which finally determines the discrepancies between the particle burn times.
The weak coupling between the burn time and temperature of the particle is a key feature of heterogeneous combustion in the diffusion-limited regime.  




 §.§ Comparison with experiments

To validate the predictive capability of the theory, the theoretical burn times and temperatures of micron-sized iron particles are compared with a series of measurements reported in Refs. <cit.>. These measurements were conducted for isolated iron particles burning in O2/N2 mixtures at room temperature. As the particle temperature right after laser ignition could not be directly measured for most cases in the experiments <cit.>, an initial partial temperature of 1500K is assumed, identical to that used in Ref. <cit.>.



<Ref>(a) reveals a comparison between the theoretical and experimental burn times of laser-ignited iron particles within the size range of 25–55µm. 
The calculated burn times, based on the current theoretical model, exhibit almost perfect agreement with the measurements for initial particle diameters larger than 38µm. Whereas for smaller particles, the theory slightly underestimates the burn time. This high level of agreement validates the quantitative accuracy of the current theoretical model in predicting particle burn times, not only at normal but also at elevated oxygen levels. 
Additionally, <ref>(b) compares the particle peak temperatures predicted by the theory with those approximately measured using a spectrometer <cit.>. In the experiment, the arithmetic mean diameters of the adopted particle size distributions are 39µm and 54µm, respectively. For the theoretical calculation, the mean diameters are used. The results show that at oxygen mole fractions below approximately 26%, the theory reasonably captures the maximum particle temperature, although the measurements slightly surpass the calculations that include radiation. This discrepancy is likely attributed to the overestimation of the experimental approach, as the maximum particle temperatures are deduced from the measured gray-body emission spectra.  As explained in Ref. <cit.>, these spectra are prone to be dominated by the emissions of relatively hot particles appearing within the time period of data acquisition due to the strong dependence of thermal radiation on temperature. 
The maximum particle temperatures calculated without considering radiation are marginally higher than their radiation-included counterparts, especially at relatively low oxygen mole fractions. However, this discrepancy becomes more significant with increasing oxygen levels because the particle peak temperature rises, making radiative heat loss more prominent. 
Furthermore, the current theoretical results closely resemble particle-resolved direct numerical simulations for a 50µm particle <cit.>, which also incorporates the surface evaporation of liquid iron and oxide. When the mole fraction of oxygen exceeds 26%, both the numerical simulation and the present theory remarkably overpredict the particle peak temperature. As stated previously, evaporation is neglected in the theory, leading to the overpredication of the particle temperature at relatively high oxygen levels. In addition, this overestimation can also be attributed to the transition of the rate-limiting mechanism from external diffusion to either surface chemisorption <cit.> or internal diffusion of ions across the liquid oxide layer <cit.>, which are neither considered in the theory nor the numerical simulation of Thijs et. al <cit.>.  
Given the excellent agreement between the measured burn time and the calculation based on the assumption of the external-diffusion-limited combustion regime, it is plausible that the other two potential rate-limiting mechanisms may only come into play for a brief period near the end of the liquid-phase combustion of iron particles, specifically in proximity to the particle peak temperature. For the majority of the iron particle combustion process, it is likely that external diffusion of the oxygen predominantly dictates the dynamics, making the inner structure of the particle less influential for predicting burn time.


Finally, the weak coupling between burn time and particle temperature suggests that while burn time is a necessary parameter, it alone is not sufficient to validate particle combustion modes. Time-resolved particle temperature emerges as a crucial experimental information for validating theoretical and numerical models of particle combustion.  <Ref> shows a comparison between the time histories calculated using the current theory and those measured using two-color pyrometry <cit.>. 

    In contrast to previous experiments <cit.>, a higher laser power was used to ignite the particles, resulting in an initial heating of the particles to a relatively high temperature, which was measured directly at a relatively low oxygen level (i.e., 14%). As the oxygen level increases, the span between the particle peak temperature and intial temperature exceeds the the dynamic range of the two-color pyrometer. In the experiment, the initial temperature was not covered.

For the theoretical calculation, the initial temperature and mean particle size are taken from the experiment.  It worth noting that the initial particle temperature was only directly measured at the relatively low oxygen mole fraction of 14%. At the oxygen mole fraction of 17%, the same initial particle temperature is assumed in the theory because the particles were heated up with an identical laser power of 40W.
As illustrated in <ref>(a), the time history of the theoretical particle temperature generally aligns well with the measurement. Approaching the maximum particle temperature, the theory slightly underpredicts the measurement. The comparison ends at the peak temperature, since further oxidation beyond stoichiometric FeO during subsequent reactive cooling is not considered.  As at higher oxygen mole fraction of 17±0.5%, where the particle temperature right after ignition was not measured, the initial event of the measured temperature  trajectory cannot represent the start of the combustion process. To establish an absolute time coordinate for the comparison with the theory, in <ref>(b) the measured temperature profile is shifted  to match with the burn time extrapolated from a fitted correlation between burn time and the initial particle diameter <cit.>. 
The good agreement confirms that the current theoretical model accurately predicts the time-dependent temperature of burning iron particles, providing a robust validation for scenarios where the particle peak time is much lower than the boiling point of iron, leading to a negligible impact of evaporation.


    the oxygen levels in the ambient are not above the air condition.








§ CONCLUSIONS AND OUTLOOK


This work introduces a quantitative analytical theory for the diffusion-limited combustion of nonvolatile (metal) particles. 
Analytical solutions for the burn time and time-dependent particle temperature are derived from the conservation equations. A simple and effective approach is proposed to couple the solutions, which demonstrates a remarkable predictive capability without relying on numerical or experimental inputs. The application of the model to analyze iron particle combustion provides valuable insights into the underlying processes.

First of all, the peak temperature of the particle is constrained by the burnout of iron rather than the balance between the heat generation and heat loss of the particle. When radiation is excluded, the time-dependent parameters of the particle, including temperature, diameter, and the fraction of unburnt metal, follows a scaling of 1/d_0^2 along the time coordinate. Taking the surface radiation into account do not change the scaling significantly, whereas it introduces a negative size-dependence of the particle temperature. 

Furthermore, the gas-phase pressure, if not very low, plays no role in the diffusion-limited combustion process of nonvolatile particles when the evaporation is negligible. In contrast, increasing the far-field gas temperature accelerates the combustion process, which results in an elevated particle peak temperature following a nearly linear trend with respect to the gas temperature. As the the mass diffusivity of the oxygen is insensitive to the initial particle temperature, the burn time in the diffusion limited regime only slightly changes even for large variations of the initial particle temperature. By the time when the particle peak temperature is reached,  the initial temperature has been largely relaxed. Consequently, the peak temperature is weakly influence by the initial condition.

Finally, the theoretical predictions for the burn time and temperature of isolated iron particle are compared with their experimental counterparts.  
The  quantitatively accurate alignment between the theoretical and experimental results validates the model within the scope of its applicable assumptions.
The presented theory not only stands as a valuable tool for predicting combustion characteristics of iron particles but also serves as a foundation for future studies exploring the intricacies of nonvolatile particle combustion in diverse scenarios.




§ DECLARATION OF COMPETING INTEREST

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to
influence the work reported in this paper.



§ ACKNOWLEDGMENTS

This work is founded by the Hessian Ministry of Higher Education, Research, Science and the Arts under the cluster project Clean Circles. 














model1-num-names






















§ SIZE-DEPENDENT PARTICLE DENSITY 


An universal expression of the oxidation reaction of metal with oxygen reads

    O_2 (g) + 2x/yM (c) = 2/yM_xO_y (c).

Therefore, a stoichiometric oxide-to-oxygen mass ratio  can be defined as

    s:=2/yℳ_mo/ℳ_O_2,

with ℳ_mo and ℳ_O_2 the molar mass of the metal oxide (M_xO_y) and that of O_2.
The mass of the particle encompassing those of total oxygen and metal elements reads

    m_p = m_M + m_O,

where subscripts, M and O, denote the metal and oxygen elements, respectively. The total mass of metal elements is essentially the initial mass of the metal particle:

    m_M = m_0 = π/6ρ_md_0^3,

with d_0 the initial particle diameter and ρ_m the density of the metal. Because of the conservation of O elements, the total mass of O elements in the particle equals that of consumed O_2:

    m_O = m_O_2.

According to <ref>, the mass of metal oxide is:

    m_mo = sm_O_2 = sm_O=  s(m_p-m_M) =  s(m_p-m_0),

The remaining mass of the metal is

    m_m = m_p - m_mo = (1-s)m_p + sm_0.

Assuming the volumetric shrinkage due to mixing between liquid metal and liquid metal oxide (if they are miscible) is negligible, the volume of the particle is

    m_p/ρ_p = V_p = V_mo + V_m = m_mo/ρ_mo + m_m/ρ_m.

Substituting <ref> into <ref> yields


    ρ_p   = ρ_mρ_mo/sρ_m+(1-s)ρ_mo + (ρ_m-ρ_mo)sρ_m/sρ_m+(1-s)ρ_moV_0/V_p
       = ρ_mo/s+(1-s)ρ_mo/ρ_m + (ρ_m-ρ_mo)s/s+(1-s)ρ_mo/ρ_mV_0/V_p

Now let's define a dimensionless parameter, oxide-to-metal density ratio:

    ϱ := ρ_mo/ρ_m.

Introducing another dimensionless parameter, diameter expansion ratio of the particle:

    ϵ := d_1/d_0 = √(m_1/ρ_mo/m_0/ρ_m) = √(ℳ_mo/x/ℳ_m1/ϱ),

with d_1 the particle diameter at burnout, ℳ_m the molar mass of the metal, ϱ also reads

    ϱ = 1/ϵ^3·ℳ_mo/xℳ_m = (d_0/d_1)^3ℳ_mo/xℳ_m.

<Ref> allows us to determine ϱ by measuring the diameter expansion ration, ϵ. Therefore, estimating the densities of liquid metals and liquid oxides becomes unnecessary if they are not known a priori. However, the oxidation state at burnout is still needed, which may be clear for some metals, such as FeO for iron at approximately the end of diffusion-limited combustion.
Then substituting <ref> into <ref> gives:

    ρ_p = ρ_mo/s+(1-s)ϱ + (ρ_m-ρ_mo)s/s+(1-s)ϱV_0/V_p

Assuming that the particle remains spherical during oxidation, the particle volume then expresses as:

    v_p=π/6d_p^3.

Therefore, the particle density as a function of diameter reads

    ρ_p   = ρ_mo/s+(1-s)ϱ + (ρ_m-ρ_mo)s/s+(1-s)ϱ(d_0/d_p)^3 
       =ρ_1 + ρ_2(d_0/d_p)^3,

with

    ρ_1 := ρ_mo/s+(1-s)ϱ,

and

    ρ_2 := (ρ_m-ρ_mo)s/s+(1-s)ϱ.
 
Although <ref> has been obtained previously by Hazenberg <cit.>, the alternative derivation procedure elaborated here is new and more straightforward.



§ COMPOSITION OF MULTI-SPECIES GAS AT THE FILM LAYER


 When only two gaseous species (i.e., oxygen and inert gas) exist, the gas composition in the film layer can be easily determined using <ref>. When several inert species present, the determination of the gas composition in the film is a bit complicated (we only consider the situation where only one kind of oxidizer exists). To accomplish this task, we first define a mean molar mass of inert species:

    ℳ_I := ∑_j ≠ j_O^NX_j ℳ_j/∑_j ≠ j_O^NX_j,

where N is the number of species and j_O is the index of the oxygen. As the molar ratios between inert species are constant, ℳ_I is constant too, provided that the diffusivities of all the inert species are identical. At the film layer, the mass fraction of the oxygen can be formulated with ℳ_I:

    Y_O,f = X_O,fℳ_O/X_O,fℳ_O + (1-X_O,f)ℳ_I.

Rearranging <ref> yields the molar fraction of the oxygen at the film layer:

    X_O,f = Y_O,fℳ_I/Y_O,fℳ_I+(1-Y_O,f)ℳ_O.

The molar fractions of inert species are than can be determined as

    X_j,f = (1-X_O,f)X_j,∞/∑_j ≠ j_O^NX_j,∞, j∈ [1,N] & j≠ j_O.

Finally, λ_f =λ(T_f, X_f), (ρ_g𝒟_O)_f = ρ_g(T_f, X_f) 𝒟_O(T_f, X_f), c_p,f = c_p,g(T_f, X_f), and c_p,O = c_p,O(T_f) are determined. X_f is the array of the molar fractions of all gaseous species at the film layer.

