2.1. Experimental set-up and procedures

Figure 3 is the schematic diagram for our experimental set-up. The chamber is made of four vertical aluminium framing posts, with glass panels covering the front and rear sides. The top and bottom sides of the chamber remain open to ensure proper airflow. The surrounding environment is maintained at 25 °C under air-conditioned conditions, with relative humidity 50 %.

Figure 3.

Schematic diagram of the experimental set-up used in the present droplet combustion study.

To support the droplet, we use a low thermal conductive ceramic fibre to minimise heat transfer from the fibre to the droplet (Lai & Pan Reference Lai and Pan2023). The fibre has thermal conductivity 0.12 W m^–1. Two fibre diameters are employed: 35 μm and 100 μm. The fibre is suspended from the top of the chamber, with one end attached on a sticky tape across the diagonal extrusions. To prevent the droplet sliding down along the fibre, we tie a small knot (∼140 μm in size) at the end of the fibre. Before each experiment, a droplet of the desired liquid is carefully placed onto the knot using a syringe needle. The starting diameter of the droplet is 0.5–1.0 mm.

After igniting the droplet, we monitor and record its burning and shrinkage using a high-speed CCD camera (Phantom v7.3) positioned in front of the chamber. The video shooting is undertaken at 1000 frames per second, with exposure time 1 ms. We further add an LED light source (LEDD1B, 14.4W) behind the rear side of the chamber. This is to ensure that the recording is made in a backlighting manner. Once recorded, the video is stored and converted into a sequence of images with resolution 800 × 600 pixels. The total number of the images is 500–3000, depending on the droplet lifetime. These images are then post-processed using MATLAB for the subsequent data acquisition and analysis.

In the experiment, we select a range of common pure liquid fuels as well as practically used fuel blends. For each fuel, we run its burning 50 times. The initial diameter D ₀ of a droplet is chosen at a redefined starting moment t = 0 after the early-time heat-up period during which the droplet undergoes a simple heating prior to combustion. To further avoid possible near-contact interferences from the fibre, we use only the data from D = D ₀ down to D = 10 % D ₀ larger than the fibre diameter. We additionally carry out the corresponding droplet evaporation experiments (10 times) at elevated temperature 70 °C as a control case, confirming that the D ²-law still holds and is not influenced by the presence of the fibre.

Since our goal is to determine the shrinkage exponent n from the measured instantaneous diameter D of a droplet, it is crucial to ensure accurate measurement of D. The methodology for this measurement is described next.

2.2. Using a volume-equivalent sphere to determine the effective droplet diameter

The presence of the fibre causes a droplet to deviate from a perfect spherical shape (see figure 2 a). Accurately determining the instantaneous droplet diameter D from images is critical, as it directly impacts the measured value of n in (1.2). Taking the D ²-law as an example, one may have an impression from (1.1) that D ² represents the area of the droplet, and hence take D to be the diameter of the projection area of the droplet (Shang et al. Reference Shang, Yang, Xuan, He and Cao2020; Wang et al. Reference Wang, Huang, Qiao, Ju and Sun2020b ). However, this approach can introduce a significant error in determining n, as it neglects the fact that the droplet shrinkage occurs in three dimensions rather than two. Furthermore, if n deviates from 2, then determining D would depend on the specific shrinkage kinetics. This creates a paradox, as n is initially unknown and is precisely what we aim to determine.

A rational approach to measuring D lies in the fact that the shrinkage of a burning droplet results from the consumption of its volume at a rate $\dot{V}$ that is accompanied by the energy required for vapourisation, ρ _L $\dot{V}$ $\Delta H$ $ _{v\textit{ap}} $ , where $\Delta H$ $ _{v\textit{ap}} $ is the latent heat of vapourisation, and ρ _L is the density of the liquid fuel. Since it is the volume change that drives the phase transition, it is essential to use the droplet’s volume to determine its effective diameter, especially when the droplet deforms due to the fibre.

In order to more accurately measure the shrinkage exponent n, we adopt the volume-equivalent-sphere approach (Nomura et al. Reference Nomura, Ujiie, Rath, Sato and Kono1996; Han et al. Reference Han, Yang, Zhao, Fu, Ma and Song2016), in which the effective diameter D of a deformed droplet is determined by treating the droplet as a sphere having the same volume as the actual droplet. To achieve this, we assume that the droplet profile is axisymmetric, with the interface position r(z) measured from the axis z of symmetry. We perform a volume integral to obtain the volume for this axisymmetric body, representing the droplet’s volume after subtracting the volume of the knot (with diameter d _knot ) at the end of the fibre. By equating the net volume to that of the equivalent sphere, we can determine the effective diameter D of the droplet according to

(2.1) \begin{equation} V=\pi \int r^{2}\,{\rm d}z-\frac{\pi }{6}d_{knot}^{3}=\frac{\pi }{6}D^{3}. \end{equation}

It is important to note that subtracting the knot’s volume is necessary because it guarantees D to be identically zero when the droplet has completely vapourised at the end of its burning.

It has been shown that this volume-equivalent-sphere approach is able to effectively address the droplet deformation issue, and successfully recovers the D ²-law (Mandal & Bakshi Reference Mandal and Bakshi2012; Chen et al. Reference Chen, Yang, Yang and Wei2024). In the study of Chen et al. (Reference Chen, Yang, Yang and Wei2024), the exponent n is determined by the dynamic slope method. Since this method involves the droplet lifetime t _life , which can be influenced by the presence of the fibre, it remains unclear whether this influence significantly affects the measured values of n in droplet combustion – where fibre heating effects may become more pronounced at elevated temperatures – even though such effects do not appear to impact the D ²-law for droplet evaporation. To address this concern, in addition to using a thinner fibre, we also determine the value of n by directly fitting the data to (1.2), as described next.

2.3. Determining the shrinkage exponent through direct data fitting

To fit the measured D–t data to (1.2), we start with the alternative form

(2.2) \begin{equation} \left(D/D_{0}\right)^{n}=1-K't, \end{equation}

where K ^′ = K/D ₀ ⁿ is the modified burning rate constant, which is also unknown. Therefore, we have two unknown parameters, n and K´, that need to be determined.

The purpose of directly fitting data to (2.2) is to minimise the influence of the fibre as much as possible. This is achieved by selecting the data sufficiently far from the end of the burning. In other words, instead of fitting the entire data set all the way to t _life , we select the data points only up to D = D _f at a specific moment t _f not close to t _life . At this moment, (2.2) must satisfy

(2.3) \begin{equation} \left(D_{\kern-1pt f}/D_{0}\right)^{n}=1-K't_{\kern-1pt f}. \end{equation}

The ratio between (2.2) and (2.3) yields

(2.4) \begin{equation} D/D_{\kern-1pt f}=\left[1+G\left(t_{\kern-1pt f}-t\right)\right]^{s}, \end{equation}

where the parameters are now transformed to

(2.5a,b) \begin{equation} s=1/n,\quad G=K'/(1-K't_{\kern-1pt f}). \end{equation}

Now we define Y = ln(D/D _f ) and T = t _f – t. Instead of directly fitting the data to (1.2), we perform a regression of the measured experimental data Y ^exp = ln (D/D _f ) to the model

(2.6) \begin{equation} Y=s\ln (1+G T), \end{equation}

with the new parameters s and G replacing n and K ^′.

The advantage of the recast model above is that the selected data range does not have to start from the initial point D = D ₀ at t = 0 – it can start at an arbitrary time, and end at any later moment. This flexibility allows for the exclusion of both the early heat-up period and the late-time period during which the fibre effects become more pronounced. Additionally, this approach offers the benefit of performing a consistency check to verify whether different data ranges still conform to the same kinetic law.

The regression is performed using the least square method for minimising the sum of squared errors, E, between the experimental value Y _i ^exp = ln [D(T _i )/D _f ] and the model value Y _i = Y(T _i ) at T _i (i = 1 … N) across all points in the selected data set:

(2.7) \begin{equation} E=\sum _{i=1}^{N}\left(Y_{i}^{\textit{exp} }-Y_{i}\right)^{2}. \end{equation}

The error minimisation is carried out under ∂E/∂s = 0 and ∂E/∂G = 0, giving

(2.8a,b) \begin{equation} s=\frac{\sum _{i=1}^{N}Y_{i}^{\textit{exp} }\ln \left(1+GT_{i}\right)}{\sum _{i=1}^{N}\left[\ln \left(1+GT_{i}\right)\right]^{2}},\quad s=\frac{\sum _{i=1}^{N}{Y_{i}^{\textit{exp} }}/({1+GT_{i}})}{\sum _{i=1}^{N}{\ln (1+GT_{i})}/({1+GT_{i}})}. \end{equation}

Combining the above equations gives an equation for G. Once G is found by numerically solving the combined equation, the value of n can be easily determined from s = 1/n by substituting the value of G into (2.8a ) or (2.8b ).

Using the volume-equivalent-sphere approach described in § 2.2 and the direct fitting method outlined above, we can later demonstrate that evaporation of a fibre-suspended droplet still follows the classical D ²-law very well, thereby validating the reliability of our measurement technique in capturing the overall shrinkage kinetics. In contrast, for droplet combustion, the measured n values across different fuels will be shown to consistently fall in a narrow range significantly greater than 2 with little variation, largely independent of fibre diameter or fibre support method. This consistency strongly indicates that the observed deviation from the D ²-law must arise from droplet combustion, rather than from the fibre’s presence or associated deformation effects.

4.1. Droplet vapourisation due to flame-driven buoyant convection

As demonstrated in the preceding section, we have shown that flame-driven droplet combustion exhibits a significant departure from the D ²-law, with the extent of deviation being influenced by fuel volatility. In this section, we present a comprehensive theory to explain these findings.

First and foremost, why n deviates from 2 has to be traced back to the method used to determine n, which is derived from a heat balance over a vapourising fuel droplet of radius R. Assuming that the droplet is rapidly mixed internally, the temperature within it can be considered uniform, as if it were an infinitely conductive sphere (Aggarwal, Tong & Sirignano Reference Aggarwal, Tong and Sirignano1984). Under this assumption, the heat required to vapourise the droplet (with density ρ _L and volume V = 4πR ³/3) is primarily supplied from the surrounding gas, driven by the temperature difference ΔT between the flame and the droplet surface (of area A = 4πR ²), as illustrated in figure 9. The corresponding energy balance can be expressed as

(4.1) \begin{equation} \rho _{L}\dot{V}{\Delta} H_{\textit{vap}}=-h\,{\Delta} \textit{TA}. \end{equation}

Figure 9.

Schematic figures for different local flow and boundary layer structures around a burning droplet to account for variations of n observed in the experiments. (a) Slipping flow with little soot particle contamination on the droplet surface, giving n = 8/3 ≈ 2.67. (b) Shear flow with severe soot particle contamination on the droplet surface, which results in n = 35/13 ≈ 2.69. (c) Straining flow resulting from strong airflow impingent arising from vigorous flame burning, yielding n = 2.6.

Here, the vapourisation energy provided by the latent heat ΔH _vap on the left-hand side is associated with the mass loss of the droplet ρ _L $\dot{V}$ = 4πR ² $\dot{R}$ ρ _L (< 0). The heat transfer coefficient h can be evaluated from the Nusselt number, which depends on the detailed thermal boundary layer structure around the droplet:

(4.2) \begin{equation} \textit{Nu}\equiv \textit{hR}/k, \end{equation}

where k is the thermal conductivity of the gas phase. With (4.2), (4.1) can be rewritten as

(4.3) \begin{equation} \rho _{L}{\Delta} H_{\textit{vap}}\frac{{\rm d}R^{2}}{{\rm d}t}=-2k\, \textit{Nu}\,{\Delta} T. \end{equation}

If the heat transfer is purely by conduction, i.e. Nu = 1, then (4.3) reduces to

(4.4a) \begin{equation} \frac{{\rm d}R^{2}}{{\rm d}t}=-2k{\Delta} T/{\Delta} H_{\textit{vap}}\rho _{L}, \end{equation}

which recovers the classical D ²-law given in (1.1). The corresponding vapourisation rate constant is

(4.4b) \begin{equation} K=8\alpha \left(\rho /\rho _{L}\right) \left(C_{p}{\Delta} T/{\Delta} H_{\textit{vap}}\right), \end{equation}

where α = k/ρC _p is the thermal diffusivity of the gas phase (with gas density ρ and heat capacity C _p ). The term C _p ΔT/ΔH $ _{v\textit{ap}} $ is known as the Spalding heat transfer number (Abramzon & Sirignano Reference Abramzon and Sirignano1989).

If there is a convection in the gas phase that is much stronger than conduction when the Péclet number ${\textit{Pe}}= \textit{UR}/\alpha\gg1$ , then h can be enhanced by the flow speed U according to

(4.5) \begin{equation} {\textit{Nu}}\sim {\textit{Pe}}^{\phi }, \end{equation}

where φ is the exponent, depending on detailed flow structures. If U is a given imposed flow speed and ΔT is fixed, as in forced convection, then (4.5) implies Nu ∝ R ^φ , leading to n = 2 – φ in (4.3). This results in a negative departure from the D ²-law (Fang et al. Reference Fang, Tseng, Yang and Wei2026). For example, under strong forced convection driven by a constant airflow where φ $=$ 1/2, the corresponding exponent is n $=$ 1.5 (Starinskaya et al. Reference Starinskaya, Miskiv, Nazarov, Terekhov, Terekhov, Rybdylova and Sazhin2021).

In droplet combustion under gravity, however, ΔT will not only set up a buoyant convection to affect Nu in (4.5), but also vary with the droplet radius R that determines the fuel consumption and degree of combustion. Such convection is set up by the red-hot flame of width 2W, causing an upward buoyant flow with velocity (Deen Reference Deen1998)

(4.6) \begin{equation} U\sim \left(g\beta {\Delta} T W\right)^{1/2}. \end{equation}

This velocity scale comes from balancing the driving buoyancy force ρβ $\Delta T$ g W ³ to the inertial force ρU ² W ² on the fluid within the lower envelope of the flame where burning is most intense. Here, β ( $=$ 1/T _∞) is the volume expansion coefficient of the air (of density ρ _∞ at temperature T _∞ $=$ 300 K). Alternatively, (4.6) can also be obtained by balancing the buoyancy term ρ _∞ β $\Delta T$ g to the inertial term ρ _∞ $\boldsymbol{v}$ ·∇ $\boldsymbol{v}$ ∼ρ _∞ U ²/W in the Navier–Stokes equation for the velocity field $\boldsymbol{v}$ of the airflow around the droplet under the Boussinesq approximation (Tritton Reference Tritton1988). Since $\boldsymbol{v}$ is driven by the flame rather than the droplet, the flame width W serves as the more appropriate length scale in this context.

It is important to note that the droplet merely passively receives heat transfer from the up-flowing convective vapour stream resulting from this natural convection. In other words, from the droplet’s standpoint, since the buoyant convection here is not generated directly by the droplet but indirectly by the flame, the heat transfer to the droplet is of forced convection type. This is fundamentally different from the combustion of a solid particle, where heat transfer is driven by natural convection from a flame situated directly adjacent to the particle (Potter & Riley Reference Potter and Riley1980; Wang et al. Reference Wang, Cui, Bi, Liu, Dong, Xing, Li and Liu2020a ). In the latter case, the induced buoyant velocity scale is U∼ (gβ $\Delta \textit{TR}$ )^1/2, and heat transfer is typically characterised by the Grashof number Gr = gβΔTR ³/ν ² that measures the strength of buoyancy-induced flow relative to viscous dissipation, where ν is the kinematic viscosity of the gas phase. Since the flame size is comparable to the particle size 2R, Gr (which is proportional to R ³) may not be very large when the particle is small. As a result, buoyant convection may be relatively weak compared to that in droplet combustion where the flame is significantly greater than the droplet. A similar situation can also occur in droplet evaporation, where the droplet is typically small and buoyant convection is negligible, thereby allowing the classical D ²-law to remain valid.

Returning to our case, although Nu varies with Pe in a manner like (4.5), this relationship can be recast in terms of Gr. Specifically, by substituting the flame-induced buoyancy velocity scale from (4.6) into Pe = UR/α, we obtain Pe = (gβΔTR ² W/α ²)^1/2, which can be rewritten in terms of Gr and the Prandtl number Pr = ν/α as Pe = Gr ^1/2 Pr (W/R)^1/2. Here, the additional factor (W/R)^1/2 accounts for the influence of the flame structure on the droplet. This formulation highlights the indirect but important role of Gr in describing the buoyancy-enhanced heat transfer experienced by the droplet arising from the flame.

Since the flame is the place where the convection is generated, it acts like a permeable hot shell that entrains the surrounding air towards the droplet on the bottom portion of the flame where burning is the most vigorous. How far the flame is located with respect to the droplet, which is reflected by W/R, can also influence heat transfer to the droplet – the closer the flame, the greater the heat transfer rate to the droplet, similar to heat transfer by jet impingement (Liu, Lienhard & Lombara Reference Liu, Lienhard and Lombara1991). So we further add this flame size effect to (4.5) as

(4.7) \begin{equation} {\textit{Nu}}\sim {\textit{Pe}}^{\phi }\left(W/R\right)^{-\omega }, \end{equation}

with ω (>0) being the exponent to account for the degree of this flame size effect.

Since both ΔT and W in (4.6) are also part of the results waiting to be determined, how the heat flux term Nu ΔT in (4.3) depends on R for determining the shrinkage exponent n in (1.2) requires the knowledge of how ΔT and W vary with R. Hence additional conditions are needed to relate ΔT and W to each other.

First, we anticipate that W would tend to be suppressed by U. To find their relationship, we consider the lower envelope of the flame into which there is a strong air entrainment due to intense burning. When the vapourised fuel (with diffusivity $\mathscr{D}$ ) diffuses from the droplet (with the vapour phase concentration ρ _fuel at the surface) onto the flame, it is quickly consumed by the burning with fresh oxygen at a high supply rate from the air suction below. Hence on the flame, the fuel vapour’s diffusive flux $j_{fuel} \sim \mathscr{D}\rho_{\textit{fuel}}/W$ across the vapour zone has to match the consumed oxygen’s convective flux j _O2∼ U Δρ _O2 (times the stoichiometry ratio), where Δρ _O2 is the consumed oxygen concentration, which is typically very small because oxygen is quite abundant here, i.e. Δρ _O2<< ρ _fuel. As a result of j _fuel ∼ j _O2, the flame width is suppressed by the flow according to

(4.8) \begin{equation} W\sim \left({\mathscr{D}}/U\right)\left(\rho _{\textit{fuel}}/{\varDelta} \rho _{\mathrm{O2}}\right)\!. \end{equation}

As for the second condition, we consider the fuel mass balance over the entire flame. The fuel vapourised from the droplet has to be consumed by the oxygen supply from the upward airflow due to the burning. This requires the total fuel consumption rate −r _rxn V _f within the flame (of volume V _f and rate of reaction r _rxn < 0) to match the total convective flux of consumed oxygen Δρ _O2 UW ² across the flame width:

(4.9) \begin{equation} -r_{\textit{rxn}}V_{f}\sim \rho _{\mathrm{O2}}UW^{2}. \end{equation}

In other words, the rate of fuel consumption can be represented by that of oxygen. For the same reason, the heat generation from the flame can also be expressed in the same manner. This heat provides the heat flow needed to vapourise the droplet, albeit there is some heat loss through the upward airflow on the top of the flame. Therefore, we can relate the heat flow towards the droplet hΔTA = (k/R)Nu ΔTA in (4.3) to the convective heat flow carried by the airflow:

(4.10) \begin{equation} \left(k/R\right)\textit{Nu}{\Delta} T\,A\sim \left(-r_{\textit{rxn}}V_{f}\right)\,|{\Delta} H_{\textit{rxn}}|\sim {\Delta} \rho _{\mathrm{O2}}UW^{2}\,|{\Delta} H_{\textit{rxn}}|, \end{equation}

where $\Delta H$ _rxn (<0) is the heat of reaction.

Now having two additional conditions (4.8) and (4.10) to relate $\Delta T$ and W, we can determine how they vary with R. To do so, it is more convenient to scale both with respective inherent scales determined only by fluid properties. For $\Delta T$ , because it is generated from combustion through heat of reaction $\Delta H$ _rxn , it is natural to scale it with the temperature equivalent to the heat of burning reaction T _rxn = | $\Delta H$ _rxn | $ / $ C _p for measuring the fuel’s potential to release heat. As for W, it can be measured relative to some inherent length scale $\mathrm{\ell }$ . This length scale can be found by re-expressing the Péclet number Pe = UR/α after substitution of the velocity scale (4.6):

(4.11) \begin{equation} \textit{Pe}\sim \left({\Delta} T/T_{\textit{rxn}}\right)^{1/2}\left(W/\mathrm{\ell }\right)^{1/2}\left(R/\mathrm{\ell }\right), \end{equation}

wherein

(4.12) \begin{equation} \mathrm{\ell }=\big(g\beta T_{\textit{rxn}}/\alpha ^{2}\big)^{-1/3} \end{equation}

represents the inherent buoyancy length for a given fuel–gas system. It is approximately 100 μm at T _rxn of typical value 10⁴ K in air with α = 0.2 cm² s⁻¹. Since the droplets used in our experiments are on a mm scale, they are larger than this inherent length scale, within which the proposed theory is applicable.

Having scaled $\Delta T$ and lengths with T _rxn and $\mathrm{\ell }$ , respectively, we can rewrite the energy balance (4.10) in the dimensionless form

(4.13) \begin{equation} \textit{Nu}\left({\Delta} T/T_{\textit{rxn}}\right)\sim f\,\textit{Pe}\left(W/\mathrm{\ell }\right)^{2}\left(R/\mathrm{\ell }\right)^{-2}, \end{equation}

where f = Δρ _O2/ρ is the mass fraction of oxygen consumption. With (4.6), W in (4.8) can be related to $\Delta T$ through

(4.14) \begin{equation} W/\mathrm{\ell }\sim \left({\Delta} T/T_{\textit{rxn}}\right)^{-1/3}\big({\textit{Le}}^{-1}\chi \big)^{2/3}, \end{equation}

where χ = ρ _fuel /Δρ _O2 is the concentration ratio of the vapour fuel to the consumed oxygen, and $Le = {\alpha}/\mathscr{D}$ is the Lewis number. Substituting (4.7) and (4.11) into (4.13) and replacing W in terms of $\Delta T$ using (4.14), we can determine first $\Delta T$ and then W in terms of R as

(4.15a) \begin{align}&\qquad {\Delta} T/T_{\textit{rxn}}\sim {\Lambda} ^{\gamma }\left(R/\mathrm{\ell }\right)^{-\gamma \left(1+\omega +\phi \right)}, \\[-12pt] \nonumber \end{align}

(4.15b) \begin{align}& W/\mathrm{\ell }\sim {\Lambda} ^{-\gamma /3}\left(R/\mathrm{\ell }\right)^{\gamma \left(1+\omega +\phi \right)/3}\big({\textit{Le}}^{-1}\chi \big)^{2/3}, \\[9pt] \nonumber \end{align}

where Λ = f (Le ⁻¹ χ)^{(5 + 2ω− φ)/3} and γ = 3/(4+ω +φ). Equation (4.15a ) indicates that the temperature difference $\Delta T$ – approximately equal to the flame temperature – increases as the droplet shrinks, which is consistent with experimental observations (Shaw & Wei Reference Shaw and Wei2007; Won et al. Reference Verwey and Birouk2018; Meng et al. Reference Meng, Lai, Zhang, Willmott and Zhang2025). This behaviour can be understood by the fact that as the droplet shrinks, more fuel is vapourised and subsequently consumed on the flame, leading to an increase in the flame temperature. Equation (4.15b ) predicts that the flame width W decreases with the droplet size, as a shrinkage droplet releases less vapourised fuel, thereby reducing the amount of combustible material available to sustain the flame. This leads to a narrower flame envelope, consistent with observations during the later stages of combustion.

With (4.14), the Péclet number Pe = UR/α in (4.11) can be written in terms of $\Delta T$ :

(4.16) \begin{equation} \textit{Pe}\sim \left({\Delta} T/T_{\textit{rxn}}\right)^{1/3}\left(R/\mathrm{\ell }\right)\big({\textit{Le}}^{-1}\chi \big)^{1/3}. \end{equation}

With (4.14) again, Nu in (4.7) can also be written in terms of $\Delta T$ :

(4.17) \begin{equation} {\textit{Nu}}\sim \left({\Delta} T/T_{\textit{rxn}}\right)^{\left(\phi +\omega \right)/3}\left(R/\mathrm{\ell }\right)^{\phi +\omega }\big({\textit{Le}}^{-1}\chi \big)^{\big(\phi -2\omega \big)/3}. \end{equation}

Substituting (4.17) and (4.15a ) into the alternative form of (4.3) by writing $\Delta H$ $ _{v\textit{ap}} $ = C _p T $ _{v\textit{ap}} $ :

(4.18) \begin{equation} \frac{{\rm d}R^{2}}{{\rm d}t}=-2\alpha \left(\rho /\rho _{L}\right)\textit{Nu}\left(T_{\textit{rxn}}/T_{\textit{vap}}\right) \left({\Delta} T/T_{\textit{rxn}}\right), \end{equation}

we arrive at

(4.19a) \begin{equation} \frac{{\rm d}R^{2}}{{\rm d}t}=-2\alpha \left(\rho /\rho _{L}\right)\left(T_{\textit{rxn}}/T_{\textit{vap}}\right) \big({\textit{Le}}^{-1}\chi \big)^{\left(\phi -2\omega \right)/3}{\Lambda} ^{\left(\phi +\omega +3\right)\gamma /3}\left(R/\mathrm{\ell }\right)^{-m}, \end{equation}

where m is the deviation exponent from the D ²-law, found to be

(4.19b) \begin{equation} m=\frac{3}{4+\omega +\phi }. \end{equation}

After integration of (4.19a ) and writing the result in terms of D, we arrive at a non-square power law given by (1.2) with the shrinkage exponent

(4.20a) \begin{equation} n=2+m. \end{equation}

By integrating (4.19a ) to (1.2) in terms of D, the resulting burning rate constant is

(4.20b) \begin{equation} K=8\times2^{m}\left(1+m/2\right)c\mathrm{\ell }^{2+m}\left(\frac{\alpha }{\mathrm{\ell }^{2}}\right)\left(\frac{\rho }{\rho _{L}}\right)\left(\frac{T_{\textit{rxn}}}{T_{\textit{vap}}}\right)\big({\textit{Le}}^{-1}\chi \big)^{\left(2\phi +\omega +5\right)m/3}f^{\left(\phi +\omega +3\right)m/3}, \end{equation}

with c being a dimensionless numerical pre-factor.

As indicated by (4.20b ), K scales as $\mathrm{\ell }$ ^2+m (α/ $\mathrm{\ell }$ ²) and has units mm^2+m s⁻¹, consistent with the dimensional requirement of the non-square law D ^2+m = $D_{0}^{2+m}$ −Kt from (1.2). Unlike the classical K in the D ²-law, which has units mm² s⁻¹, this rate constant incorporates buoyant convection effects that lead to m>0, reflecting the non-quadratic nature of the law. This dimensional change is essential, as gravitational acceleration g must enter the shrinkage kinetics through the flame-induced buoyancy velocity scale (4.6), which in turn alters the shrinkage exponent n and thereby the units of K.

In other words, rather than introducing buoyant convection as a minor correction to the K in the D ²-law, as done in the study by Law & Williams (1972), the present work shows that buoyant convection fundamentally changes the burning behaviour. It leads to a new shrinkage kinetic law characterised by a different exponent n and a redefined rate constant K. This distinction is crucial, as it marks a shift from additive corrections to a structurally different formulation that captures the dominant role of buoyancy in shaping the burning dynamics of droplets under normal gravity.

It is worth noting that significant efforts have previously been made to develop theories of droplet vapourisation in both normal and burning regimes, most notably by Sirignano and Law. The theoretical framework presented in this study is fundamentally different from theirs.

Sirignano’s theory is mainly focused on spray combustion in which droplet burning occurs while the droplet is moving within a convective flow (Sirignano Reference Sirignano1983). In this case, the heat and mass transfers are of forced-convection type, with the transfer coefficients described by Ranz–Marshall-type correlations based on the classical Spalding model (Abramzon & Sirignano Reference Abramzon and Sirignano1989). This leads to a D ^3/2-law (Prakash & Sirignano Reference Prakash and Sirignano1980), which does not adequately capture the droplet vapourisation dynamics observed under buoyant convection. In contrast, our theory is fundamentally different: it explicitly considers buoyant convection by fully incorporating the coupled effects of momentum, heat and mass transfer. This results in a new Dⁿ -law, with n significantly greater than 2, enabling us to capture the enhanced vapourisation behaviour driven by buoyancy (see § 4.2).

Another aspect of Sirignano’s theory is that the shrinkage kinetics can be further influenced by internal conduction and convective effects within the droplet (Sirignano Reference Sirignano1983). However, our results suggest that such influence may not be as significant as previously assumed. In our analysis, the droplet interior is treated as rapid mixed (i.e. infinitely conductive), thereby neglecting internal thermal gradients. Despite this simplification, our model successfully captures the experimentally measured n values across a range of fuels (see § 4.2). This consistency indicates that under the conditions studied, the shrinkage kinetics is primarily governed by external transport mechanisms, and that internal conduction and convection play a secondary role.

In fact, internal vortical fluid motion within the droplet exhibits closed streamlines, which gives rise to very distinctive heat transfer characteristics. Acrivos and his co-workers have shown both theoretically and experimentally that for heat transfer within closed streamlines, Nu approaches a constant in the Pe → ∞ limit – a value several times higher than that in the pure conduction case (Pe = 0) (Pan & Acrivos Reference Pan and Acrivos1968; Robertson & Acrivos Reference Robertson and Acrivos1970; Acrivos Reference Acrivos1971). Unlike the typical power-law relationship Nu∼Pe ^φ described by (4.5), this behaviour arises because in the Pe → ∞ limit, the isotherms coincide with the streamlines across which heat transfer can be achieved only by conduction. This implies that internal vortical flow can enhance the droplet’s effective thermal conductivity (Abramzon & Sirignano Reference Abramzon and Sirignano1989). For a non-evaporative droplet, the overall heat transfer can be modelled as a combination of internal and external heat transfers in series (Abramzon & Borde Reference Abramzon and Borde1980). When external convection is strong – i.e. when the external Nusselt number exceeds the internal one – the internal vortex mechanism dominates the overall heat transfer. This internally dominated heat transfer behaviour, characterised by a constant Nusselt number in the high Péclet number regime, has been confirmed by Zhang, Yang & Mao (Reference Zhang, Yang and Mao2012) in their simulation study of a droplet in an external flow field.

However, in our case involving an evaporative droplet, the situation is different because internal and external heat transfers occur simultaneously at the interface, and act in parallel to determine the evaporation rate. As a result, when external convection is strong, the influence of internal vortex effects appears to be minor, as the predictions from the enhanced effective conductivity model arising from such effects do not differ significantly from those of the infinite conductivity model (Abramzon & Sirignano Reference Abramzon and Sirignano1989). Combined with the fact that our experimental results deviate substantially from the D ²-law and closely match theoretical predictions based on the assumption of rapid internal mixing, this strongly supports the validity of the infinite conductivity droplet approximation adopted in our theoretical model.

In fact, such infinite conductivity droplet approximation may be even more applicable in our case, as internal mixing can arise during droplet evaporation due to Marangoni thermal convection induced by non-uniform heating (Mandal & Bakshi Reference Mandal and Bakshi2012; Kazemi et al. Reference Kazemi, Saber, Elliott and Nobes2021), which further reinforces the D ²-law (Mandal & Bakshi Reference Mandal and Bakshi2012). During droplet combustion – where the D ²-law no longer holds – this internal convection can become further intensified and acts in the same direction as the upward buoyant convection, driven by higher temperatures at the droplet’s bottom where the flame burns more intensely. Similar to droplet evaporation under an imposed temperature gradient in a convective environment (Shih & Megaridis Reference Shih and Megaridis1996), such combined convection can further enhance the evaporation rate, thereby causing the droplet’s shrinkage behaviour to deviate even more from the classical D ²-law.

As for Law’s theory, it does account for natural convection effects (Law & Williams 1972), but these effects appear merely as a Grashof number (Gr) based correction to the rate constant K within the framework of the classical D² -law. Strictly speaking, such a correction is valid only in the limit of small Gr, which contrasts with the actual scenario where Gr is typically large. Our theory is fundamentally different, as it does not treat buoyant convection as a small perturbation to the classical framework or as a variant of the D ²-law. Instead, it fully incorporates buoyancy-driven flow into the governing momentum, heat and mass transfer processes, leading to a new Dⁿ -law with n significantly greater than 2. This new formulation captures the dominant role of buoyant convection in shaping the droplet burning behaviour under normal gravity.

4.2. Influence of thermal boundary layer structures on shrinkage kinetics

The above clearly shows that droplet burning under gravity can alter the D ²-law to a non-square law due to buoyant convection, as also displayed in figure 6 for the measured values of n from experiments. While the measured n values are found to span a narrow range, small but systematic variations can still be observed. Such variations are not unexpected as they can arise from different reaction pathways or thermochemistries, depending on chemical nature of fuels (Curran Reference Curran2019; Yang Reference Yang2021). Such fuel-chemistry-dependent nature can significantly affect flame structures and characteristics (Kohse-Höinghaus Reference Kohse-Höinghaus2021). Since the fuel droplet shrinkage kinetics at large are mainly determined by transport processes to display a nearly universal behaviour, the small variations in the measured n values for a variety of fuels may be explained by different flame and flow structures caused by different combustion chemistries of fuels. In this subsection, we account for such variations in terms of distinct thermal boundary layer structures according to different flame appearances seen in experiments.

We begin with the Nusselt number under the strong convection condition. It can be typically established from Nu ∼ (δ/R)⁻¹ (Deen Reference Deen1998) if the scale of the thermal boundary layer thickness δ can be identified for a given flow field. To do so, instead of using the standard boundary layer theory, we find it is much simpler to make use of δ ∼ (α t)^1/2 in transient heat conduction problems to establish the scaling for δ:

(4.21a) \begin{equation} \delta \sim \big({\alpha}\,t_{\textit{flow}}\big)^{1/2}, \end{equation}

with time t in δ ∼ (α t)^1/2 being replaced by the flow time t _flow based on the flow velocity scale u _δ within the boundary layer,

(4.21b) \begin{equation} t_{\textit{flow}}\sim R/u_{\delta }. \end{equation}

Using (4.21a ) together with (4.21b ) in Nu ∼ (δ/R)⁻¹ yields

(4.22) \begin{equation} {\textit{Nu}}\sim \left(\frac{u_{\delta }R}{k}\right)^{1/2}. \end{equation}

This simplified approach provides a more general recipe to determine Nu based on the local Péclet number Pe _δ = u _δ R/α within the boundary layer without having to involve a formal boundary layer analysis. The advantage of (4.22) is that it offers a unified way to determine Nu for three different flow types (i) slipping flow, (ii) shear flow, and (iii) straining flow commonly seen in practice. In fact, as illustrated in figure 9 as well as seen in the analysis below, these different flow types may happen in different burning situations, allowing us to explain why the measured values of n display small variations between approximately 2.6 and 2.7 shown in figure 6.

Flow type (i) corresponds to the situation where the droplet surface becomes slippery (figure 9 a) due to a possible existence of vapour film (Sirignano Reference Sirignano2014). In this case, u _δ ∼U, making (4.22) reduced to the well-known result like the one for bubbles (Deen Reference Deen1998):

(4.23) \begin{equation} {\textit{Nu}}\sim {\textit{Pe}}^{1/2}. \end{equation}

So with ν $=$ 1/2 and ω $=$ 0, (4.15a ) and (4.15b ) become

(4.24a) \begin{align}& {\Delta} T/T_{\textit{rxn}}\sim f^{2/3}\big({\textit{Le}}^{-1}\chi \big) \left(R/\mathrm{\ell }\right)^{-1}, \\[-12pt] \nonumber \end{align}

(4.24b) \begin{align}& W/\mathrm{\ell }\sim f^{-2/9}\big({\textit{Le}}^{-1}\chi \big)^{1/3}\left(R/\mathrm{\ell }\right)^{1/3}. \\[6pt] \nonumber \end{align}

The resulting deviation exponent (4.32b ) is m = 2/3, hence the shrinkage exponent from (4.20a ) is

(4.24c) \begin{equation} n=8/3\approx 2.67. \end{equation}

The burning rate constant from (4.20b ) is then

(4.24d) \begin{equation} K=8\times2^{2/3}(4c/3)\mathrm{\ell }^{8/3}\left(\frac{\alpha }{\mathrm{\ell }^{2}}\right)\left(\frac{\rho }{\rho _{L}}\right)\left(\frac{T_{\textit{rxn}}}{T_{\textit{vap}}}\right)\big({\textit{Le}}^{-1}{\chi} \big)^{4/3}f^{7/9}. \end{equation}

The predicted exponent 2.67 from (4.24c ) well captures the measured values n = 2.66 ± 0.21 (tetradecane) and 2.65 ± 0.17 (dodecanol) for less violent burning like figure 4(b).

Flow type (ii) may take place in the scenario where the droplet surface becomes rigidified due to contamination by soot particles (figure 9 b). It has been observed experimentally that such soot contamination on the surface of a droplet can become quite severe in practically used fuels such as kerosene and diesel (Rasid & Zhang Reference Rasid and Zhang2019). The local flow behaviour around the droplet surface in this case is therefore of shear flow type with u _δ ∼Uδ/R, rendering (4.22) to be like the equation for solid particles (Deen Reference Deen1998):

(4.25) \begin{equation} {\textit{Nu}}\sim {\textit{Pe}}^{1/3}, \end{equation}

giving φ = 1/3 and ω = 0 in (4.7), and m = 9/13 in (4.19b ). So (4.15a ), (4.15b ), (4.20a ) and (4.20b ) become

(4.26a) \begin{align}&\qquad\qquad\qquad\quad {\Delta} T/T_{\textit{rxn}}\sim f^{9/13}\big({\textit{Le}}^{-1}\chi \big)^{14/13}\left(R/\mathrm{\ell }\right)^{-12/13}, \\[-12pt] \nonumber \end{align}

(4.26b) \begin{align}&\qquad\qquad\qquad\qquad W/\mathrm{\ell }\sim f^{-3/13}\big({\textit{Le}}^{-1}\chi \big)^{4/13}\left(R/\mathrm{\ell }\right)^{4/13}, \\[-12pt] \nonumber \end{align}

(4.26c) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad n=35/13\approx 2.69, \\[-12pt] \nonumber \end{align}

(4.26d) \begin{align}& K=8\times2^{9/13}\left(35c/26\right)\mathrm{\ell }^{35/13}\left(\frac{\alpha }{\mathrm{\ell }^{2}}\right)\left(\frac{\rho }{\rho _{L}}\right)\left(\frac{T_{\textit{rxn}}}{T_{\textit{vap}}}\right)\big({\textit{Le}}^{-1}\chi \big)^{17/13}f^{10/13}. \\[6pt] \nonumber \end{align}

The predicted exponent 2.69 from (4.26c ) is slightly higher than 2.67 in (4.24c ), closely capturing the measured values n = 2.73 ± 0.22 (kerosene) and 2.72 ± 0.2 (diesel) for incomplete burning like figure 4(c).

Flow type (iii) may happen to the very violent burning situation for highly volatile fuels such as ethanol and hexane. This is because air entrainment beneath the droplet could become so strong that the uprising flow towards the droplet becomes of straining type impinging the droplet (figure 9 c). Such an impinging flow field around the bottom of the droplet is u _δ ∼ (U/W)R, with strain rate ∼U/W, which gives

(4.27) \begin{equation} {\textit{Nu}}\sim {\textit{Pe}}^{1/2}\left(W/R\right)^{-1/2}, \end{equation}

according to (4.22). With $\phi$ $=$ 1/2 and ω $=$ 1/2 in (4.7), and m = 3/5 in (4.19b ), the quantities of interest (4.15a ), (4.15b ), (4.20a ) and (4.20b ) are

(4.28a) \begin{align}&\qquad\qquad\qquad\quad {\Delta} T/T_{\textit{rxn}}\sim f^{3/5}\big({\textit{Le}}^{-1}\chi \big)^{11/10}\left(R/\mathrm{\ell }\right)^{-6/5}, \\[-12pt] \nonumber \end{align}

(4.28b) \begin{align}&\qquad\qquad\qquad\qquad W/\mathrm{\ell }\sim f^{-1/5}\big({\textit{Le}}^{-1}\chi \big)^{3/10}\left(R/\mathrm{\ell }\right)^{2/5}, \\[-12pt] \nonumber \end{align}

(4.28c) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad n=13/5=2.6, \\[-12pt] \nonumber \end{align}

(4.28d) \begin{align}& K=8\times2^{3/5}\left(13c/10\right)\mathrm{\ell }^{13/5}\left(\frac{\alpha }{\mathrm{\ell }^{2}}\right)\left(\frac{\rho }{\rho _{L}}\right)\left(\frac{T_{\textit{rxn}}}{T_{\textit{vap}}}\right)\big({\textit{Le}}^{-1}\chi \big)^{13/10}f^{4/5}. \\[6pt] \nonumber \end{align}

The exponent 2.6 from (4.28c ) predicted in this case well captures the measured values n = 2.56 ± 0.16 (ethanol), 2.58 ± 0.07 (hexane), 2.57 ± 0.18 (acetone) and 2.60 ± 0.12 (ethyl acetate) for violent burning like figure 4(a).

As shown above, different flow and boundary layer structures can make $\Delta T$ and W vary with R in different manners and hence yield different values of n. Table 1 summarises the various scaling relationships for $\Delta T$ , W, n and K given above.

Table 1.

Summary of distinct scaling relationships for the flame–droplet temperature difference ΔT, the flame width W, and the burning rate constant K at different values of the shrinkage exponent n resulting from different flow and thermal boundary layer structures illustrated in figure 8.

4.3. Burning rate constant

As we can establish the non-square power law (1.2) as above to describe the shrinkage behaviour of a burning droplet under the influence of buoyant convection, the burning rate constant K can be readily obtained from the slope of a plot of (D/D ₀) ⁿ versus t/D ₀ ⁿ if n is specified. Since the measured values of n only span a narrow range 2.6–2.7, we can use the D ^8/3-law to represent such a nearly universal power-law behaviour. This allows us to make use of the D ^8/3–t plots to determine the values of K from the slopes for the various fuels that we use in this study, as illustrated in figure 10 for representative fuels: ethanol, tetradecane and kerosene. Figure 11 plots the measured values of K against theory’s prediction

(4.29) \begin{equation} K\sim {\alpha}\mathrm{\ell }^{2/3}\left(\frac{\rho }{\rho _{L}}\right)\left(\frac{T_{\textit{rxn}}}{T_{\textit{vap}}}\right), \end{equation}

according to (4.24d ). We find that the measured K values for pure liquid fuels approximately follow (4.29), which supports our theory. The slope of the best-fitted line is approximately 9, which gives the value of the numerical pre-factor c ´ for the scaling expression (4.29). However, the measured K values for practically used fuels such as kerosene, diesel and petrol are significantly lower than those for pure fuels, implying that the corresponding values of the pre-factor c ´ are substantially smaller. If these practically used fuels still obey (4.29), then their much lower values of c ´ compared to those for pure fuels can only be explained by mixture effects in these fuels, rendering c ´ in this case to be composition-dependent. Perhaps because these fuels are originally distilled from crude oils and are far from pure, the degree of incomplete burning could become much more severe than that of pure liquid fuels, with a lot more soot particles being generated. This not only makes their flames become much extended compared to those of pure fuels, but also significantly hinder heat release from the flames to the droplets, explaining why the observed K values appear lower than for pure fuels.

Figure 10.

Plots of (D/D ₀)^8/3 versus t/D ₀ ^8/3 for determining the values of the burning rate constant K for fuel droplet combustion of (a) ethanol, (b) tetradecane and (c) kerosene, corresponding to the combustion sequences shown in figure 4.

Having a close look at how the measured K values distribute for pure fuels in figure 11, we find that they seem to approximately vary according to fuel types – the K values for alcohol, ketone and ester appear to be lower, whereas those for alkane are higher. For the same family of fuels, the K values do not seem to vary monotonically with carbon content. For alkanes, in particular, the lower carbon content fuel such as hexane can even have a comparable or greater value of K compared to values for higher carbon content, such as decane and tetradecane.

We also find that for a more volatile fuel like ethanol, while it appears to display a higher average flame temperature $\Delta T$ than a less volatile fuel like tetradecane (see figures 4 a,b), the former’s K value is smaller than the latter’s. This may be explained using the alternative form of (4.29) due to the inherent buoyancy length $\ell\propto T_{\textit{rxn}}^{-1/3}$ from (4.12):

(4.30) \begin{equation} K\sim \alpha ^{3}\mathrm{\ell }^{-7/3}g^{-1}\left(\frac{\rho }{\rho _{L}}\right)\left(\beta T_{\textit{vap}}\right)^{-1}\propto \frac{{T_{\textit{rxn}}}^{7/9}}{T_{\textit{vap}}}\propto \frac{\left| {\Delta} H_{\textit{rxn}}\right| ^{7/9}}{{\Delta} H_{\textit{vap}} }. \end{equation}

Equation (4.30) reveals that K grows with heat of reaction $\Delta H$ _rxn (<0) as $|$ $\Delta H$ _rxn $|$ ^7/9. This is because according to (4.3), vapourisation of a fuel droplet is driven by $\Delta T$ due to generation of this heat on the flame; the higher $|$ $\Delta H$ _rxn $|$ , the larger K. However, also according to (4.3), the rate of fuel mass vapourisation is diminished by the amount of latent heat $\Delta H$ $ _{v\textit{ap}} $ required to vapourise liquid into vapour, rendering a lower value of K when $\Delta H$ $ _{v\textit{ap}} $ is larger. Since $|$ $\Delta H$ _rxn $|$ typically increases with boiling temperature more rapidly than $\Delta H$ $ _{v\textit{ap}} $ , this explains why tetradecane has a greater value of K than ethanol.

For practically used fuels, however, their burning characteristics seem different compared to those of pure fuels. To illustrate this, we take kerosene as an example, and make comparison with tetradecane, whose volatility is the lowest among the pure fuels selected here. Although kerosene’s $|$ $\Delta H$ _rxn $|$ is comparable to tetradecane’s, its K value is significantly lower than tetradecane’s. This result may be attributed to formation of soot particles caused by incomplete combustion (Avedisan 2000; Wang Reference Wang2011). This may not only reduce the heat release from the flame, but also cause soot contamination on the droplet surface (Rasid & Zhang Reference Rasid and Zhang2019), thereby further lowering the vapourisation rate and hence the pre-factor c ´ in (4.29). Such soot particles are also responsible for the long and red flames observed in the experiments for this type of fuels (see figure 4 c). Perhaps also because the droplet becomes immobilised due to contamination of these particles on the droplet’s surface, the thermal boundary layer structure becomes of shear flow type with the Nusselt number scaling (4.25) to yield n ≈ 2.69 according to (4.26c ), reasonably capturing the measured n values at approximately 2.73 for kerosene, diesel and petrol.

It is worth mentioning that the burning rate constant (4.29) or (4.30) for droplet combustion under gravity is very different from that without gravity described by the D ²-law. In the latter case, the burning rate constant (4.4b) with $\Delta T$ ∼T _rxn is

(4.31) \begin{equation} K\left(g=0\right) \sim \alpha \left(\frac{\rho }{\rho _{L}}\right)\left(\frac{T_{\textit{rxn}}}{T_{\textit{vap}}}\right). \end{equation}

Compared to (4.29), one may simply say that (4.29) can be related to (4.31) as

(4.32) \begin{equation} K\sim K\left(g=0\right)\,\mathrm{\ell }^{2/3}. \end{equation}

However, as implied by $\mathrm{\ell }$ according to (4.12), T _rxn from combustion on the flame must work with g to set up buoyancy Δρg ∼ ρg β T _rxn to drive upward convection for oxygen supply. To illuminate this important dependence on T _rxn under the influence of g, (4.30) is a more appropriate form to represent the burning rate constant.

4.4. Possible sources of experimental variations

As demonstrated above, droplet shrinkage kinetics of various fuels approximately follows a universal law with n = 2.6–2.7, with the observed variations attributed to differences in flow and boundary layer structures. While the differences in the measured n values between fuels are found to be statistically significant, it is worth noting that additional factors may also contribute to these variations.

First, the most prominent factor is thermal radiation, which arises from high flame temperature during combustion. Whether this effect is important can be assessed by comparing the radiative heat flux q _rad to the convective heat flux q = h $\Delta T$ from the flame to the droplet in (4.1). The radiative heat flux is evaluated as

(4.33) \begin{equation} q_{rad}=\varepsilon \sigma _{B}\big({T_{f}}^{4}-{T_{s}}^{4}\big)\sim \varepsilon \sigma _{B}{T_{f}}^{4}, \end{equation}

where ε ∼ 1 is the emissivity of the droplet, σ_B = 5.67 × 10⁻⁸ W m⁻² K⁻⁴ is the Stefan–Boltzmann constant, and T _f is the flame temperature, typically much higher than the droplet’s surface temperature T _s. The convective heat flux q = h $\Delta T$ can be approximated as h T _f and evaluated in terms of the Nusselt number Nu = hR/k:

(4.34) \begin{equation} q\sim \left(k/R\right)\textit{Nu} T_{f}, \end{equation}

where Nu can be estimated as Nu ∼ (UR/α) ^φ (with $\phi\gt 0$ being the exponent), k ∼ 2.62 × 10⁻² W mK^-1 is the thermal conductivity of the gas phase, and α ∼ 0.2 cm² s^-1 is the corresponding thermal diffusivity.

With the above, the ratio between (4.33) and (4.34) is

(4.35) \begin{equation} q_{rad}/q\sim \left(\frac{\varepsilon \sigma _{B}T_{f}^{3}}{k/R}\right)\left(\frac{UR}{\alpha }\right)^{-\phi }. \end{equation}

Here, $ \epsilon \sigma_{B}T_{\kern-1pt f}^{3}/(k/R) $ represents the ratio of the radiative heat transfer coefficient ε σ _B T _f ³ due to thermal radiation to the conductive heat transfer coefficient k/R. Equation (4.35) indicates that the effect of thermal radiation is diminished by convection. With typical values R ∼ 1 mm and U ∼ 1 m/s, the ratio given by (4.35) falls in the range 0.3–0.6 for φ = 1/3–1/2. Therefore, thermal radiation is less significant compared to convection. Moreover, since the flux ratio is q _rad /q $ \propto $ R ^λ with ${\lambda} = 1-{\phi}\gt 0$ according to (4.35), and the above estimate is based on R taken as the initial radius R ₀ ∼ 1 mm, the ratio will become even smaller as R decreases during combustion. A similar conclusion can also be drawn for droplet evaporation based on the analogous dependence of q _rad / q on R (Abramzon & Sazhin Reference Abramzon and Sazhin2005). This is further supported by model calculations, showing that the difference in vapourisation rates with and without radiation is merely 10 % or less at elevated temperatures (Abramzon & Sazhin Reference Abramzon and Sazhin2006).

Variations can also arise from a range of additional factors. The most common ones include droplet asphericity, fibre heating effects, variations in fluid properties, and the formation of soot particles, all of which can introduce complexities into the combustion process.

Droplet asphericity results from deformation of the droplet by the support fibre. Since this issue can be appropriately addressed using the volume-equivalent approach to determine the effective droplet diameter (see § 2.2) with re-confirmation of the D ²-law, it should not pose a significant problem.

The primary concern regarding the support fibre is actually its potential contribution to additional heating of the droplet. However, since we have used a thinner fibre to conduct the experiment, and found that the measured n values consistently fall within a narrow range regardless of the measurement methods, this issue may not be of great concern. In fact, the fibre heating effect in droplet combustion can be shown to be negligible by evaluating the ratio of the heat flow from the fibre to the droplet, Q _fibre = h _fibre $\Delta T$ _fibre A _fibre , to the heat flow from the ambient gas to the droplet, Q = h $\Delta \textit{TA}$ :

(4.36) \begin{equation} \frac{Q_{\textit{fibre}}}{Q}=\frac{h_{\textit{fibre}}{\Delta} T_{\textit{fibre}}}{h{\Delta} T}\left(\frac{A_{\textit{fibre}}}{A}\right), \end{equation}

where h _fibre , $\Delta T$ _fibre and A _fibre denote the heat transfer coefficient, temperature difference and heat transfer area between the fibre and the droplet, respectively. In the fibre heating part, h _fibre $=$ k _fibre /a _fibre since heat transfer primarily occurs through conduction across the fibre, where a _fibre is the fibre radius and k _fibre is its thermal conductivity. Writing h in terms of the Nusselt number Nu as h $=$ (k/R) Nu ∼ (k/R) (UR/α) ^φ , with Nu ∼(UR/α) ^φ and using $\Delta T$ _fibre ∼ $\Delta T$ and A _fibre /A ∼ (a _fibre /R)², (4.36) can be approximated as

(4.37) \begin{equation} \frac{Q_{\textit{fibre}}}{Q}\sim \left(\frac{k_{\textit{fibre}}}{k}\right)\left(\frac{a_{\textit{fibre}}}{R}\right) \left(\frac{UR}{\alpha }\right)^{-\phi }. \end{equation}

With k _fibre /k∼5, a _fibre /R ∼ 10⁻¹ and UR/α ∼50 under U ∼ 1 m s^-1, the resulting ratio is 0.07–0.14 for φ = 1/3–1/2. This indicates that the fibre heating effect contributes approximately 10 % of the ambient gas heating, which supports the observation that the measured n values are largely unaffected by the presence of the fibre (see figures 6 a,b). Comparing results for fibre diameters 35 μm and 100 μm (see figure 7 a,b) shows that while the measured n values for 2a _fibre = 100 μm are slightly lower, they exhibit a slight increase with boiling point. The trend also persists for 2a _fibre = 35 μm, indicating that it can be attributed not to the fibre but to the intrinsic burning characteristics of the fuel.

It is worth noting that compared to droplet combustion, fibre effects in droplet evaporation may become more pronounced. It can be seen from figure 2 that the cross-fibre method (figure 2 b) yields an n value significantly greater than 2, in contrast to the suspended fibre case (figure 2 a), resulting in a more substantial deviation from the D ²-law. The greater sensitivity of droplet evaporation to fibre effects arises from the absence of convection to enhance the heat flow Q from the ambient gas to the droplet – convection would otherwise diminish the heat flow Q _fibre from the fibre. Since both Q _fibre and Q are mainly through conduction, their ratio becomes

(4.38) \begin{equation} \frac{Q_{\textit{fibre}}}{Q}\sim \left(\frac{k_{\textit{fibre}}}{k}\right)\left(\frac{a_{\textit{fibre}}}{R}\right), \end{equation}

which makes the fibre contribution more significant compared to that in (4.37) in the presence of strong convection. As also indicated by (4.38), the fibre heating effect is important only when the droplet size 2R approaches the fibre diameter 2a _fibre near the end of evaporation. This effect is further amplified in droplet evaporation under the cross-fibre method due to increased droplet–fibre contact, but becomes insignificant in droplet combustion due to suppression by strong buoyant convection. The latter also explains why the fibre only has a marginal influence on the measured n values in droplet combustion, regardless of the fibre supporting method, as shown in figure 6.

Regarding variations caused by non-constant fluid properties, it is evident that properties such as density and heat capacity not only change with temperature but also exhibit spatial variations during combustion, which are not accounted for in our theory. Nevertheless, given the fact that the measured n values closely align with theoretical predictions based on the assumption of constant fluid properties, these effects do not qualitatively alter the nature of the phenomenon. This also implies that variations in fluid properties are unlikely to exceed an order of magnitude sufficient to affect the core physics and key characteristics of the phenomenon. Indeed, numerical studies have shown that variations in droplet vapourisation rates due to temperature-dependent fluid properties are typically within 10–25 % (Kassoy & Williams Reference Kassoy and Williams1968; Megaridis Reference Megaridis1993), which falls within the range of experimental uncertainty. These findings suggest that although such effects are present, they do not compromise the validity of the present theoretical framework.

Formation of soot particles is a common occurrence in flame combustion, resulting from incomplete burning that is intrinsically linked to the fuel’s chemical composition and detailed combustion kinetics (Jin et al. Reference Jin, Yuan, Li, Yang, Zhou, Zhao, Li and Qi2023). The narrow range of measured n values across various fuels suggests that soot formation does not significantly affect the droplet shrinkage mechanism, which is mainly transport-controlled, as demonstrated here. Since soot formation is not entirely deterministic but influenced by stochastic events arising from flame instability (Chu et al. Reference Chu, Qi, Feng, Dong, Hong, Qiu and Han2023), the effects inevitably introduce variations in the droplet shrinkage process, resulting in some degree of scatter in the measured data, as observed in figure 5(a ). Since the extent of incomplete burning may be influenced by the initial fuel quantity, which is controlled by the starting droplet diameter D ₀, such soot formation effects could cause the coefficient in the expression (4.20b ) for the burning rate constant K to depend on D ₀. This might explain the slight variations in the observed shrinkage kinetics with D ₀.