Analysis of two interactive burning droplets with different temperatures

Abstract

This paper presents a comprehensive theoretical analysis of the interaction between two quasi-steady burning droplets with differing temperatures, sizes and distances, building upon the mass-flux-potential model and flame-sheet assumption. In contrast to existing research, this study introduces a fresh perspective on droplet interactions by considering the different temperatures of the droplets. Utilizing the bispherical coordinate approach, theoretical solutions for the Stefan flow, scalar fields, droplet evaporation/burning rates, interaction coefficients and flame positions have been derived successfully. A comparison with extensive numerical simulations indicates a good agreement between the analytical and numerical results under a variety of conditions. It is revealed that proximity between the droplets causes non-uniform evaporation rates on their surfaces, and in some cases, leads to condensation on the cooler droplet. Notably, when the temperatures of the two droplets differ, this results in an uneven temperature distribution across the flame surface, and increasing the temperature of one droplet substantially elevates the temperature of the nearby flame. This study also establishes a criterion for the transition between different combustion modes, specifically between group and separated combustion. The findings of this study are crucial in deepening our understanding of evaporation and combustion processes, as well as the dynamics of flame spreading, local ignition, and extinction in systems involving multiple droplets.

1. Introduction

In combustion devices that utilize liquid fuels, the implementation of spray combustion offers significant advantages. This process involves atomizing the fuel into numerous fine droplets, thereby optimizing heat and mass transfer in the combustor. It enhances liquid evaporation and fuel–air mixing, contributing to improved combustion efficiency and increased power generation. In this context, the interactive combustion of droplets within sprays plays a critical role in determining the performance and emission characteristics of the combustor (Annamalai & Ryan 1992; Sirignano 2014).

Extensive research has been conducted to explore droplet interactions within systems involving two or more droplets (Annamalai & Ryan 1992; Umemura 1994; Sazhin 2022; Li, Zhang & Law 2023a). In static environments, Labowsky (1978, 1980) proposed a simplified model for studying the evaporation and combustion processes in multi-droplet systems. Using this model, the physical processes occurring in the gas phase are assumed to be quasi-steady (QS), and the mass flux in the gas phase is represented as a gradient of a scalar potential function. Consequently, the nonlinear momentum equation, which describes the Stefan flow resulting from droplet evaporation, can be simplified to a linear Laplace equation. Based on Labowsky's QS mass-flux-potential model, several subsequent studies have been conducted. For example, Imaoka & Sirignano (2005) applied the model to various liquid systems, including droplet arrays and liquid films, accounting for variable thermophysical properties. Sirignano (2007) investigated the combustion of dense sprays with non-unity Lewis numbers, incorporating transient variations in droplet temperature. Sirignano & Wu (2008) then examined the evaporation dynamics in droplet groups with multiple fuel components.

When forced and natural convection are negligible, and Stefan convection primarily drives gas movement around the droplets, the QS mass-flux-potential model proves to be effective. This model allows for the representation of scalar quantities such as gas-phase temperature and species concentration in three-dimensional space as one-dimensional functions of the potential function (Imaoka & Sirignano 2005; Sirignano 2014). Nevertheless, deriving theoretical solutions from these simplified equations remains challenging, mainly because of the intricate boundary conditions present at the surfaces of multiple droplets. Some studies have used direct numerical simulations to solve these equations (Imaoka & Sirignano 2005; Sirignano 2007; Sirignano & Wu 2008), while others have adopted the method of images (Labowsky 1980; Annamalai & Ryan 1992; Zheng et al. 2019), a technique commonly utilized in electrostatics studies. In contrast to directly solving partial differential equations, the method of images strategically places multiple point sources of varying strengths at different positions (Zheng et al. 2019). Through iterative calculations, this approach seeks to align the boundary values with the given conditions. For widely spaced droplets, the initial iteration often yields fairly accurate results, a process commonly referred to as the point-source method (Marberry, Ray & Leung 1984; Cossali & Tonini 2020). Nonetheless, the convergence rate diminishes with increasing inter-droplet spacing, requiring multiple iterations to achieve a reasonable result.

The method of images, primarily an iterative algorithm, is particularly effective for systems containing a limited number of droplets (Annamalai & Ryan 1992). It does not, however, generally provide a theoretical solution to the underlying problems. Fortunately, for the system with two droplets, which is the simplest case involving droplet interactions, analytical solutions are attainable using the bispherical coordinate approach (Moon & Spencer 2012). Under axisymmetric conditions, the bispherical coordinate system transforms the irregular space outside the two droplets into a regular rectangular region, facilitating the application of the separation of variables method for solving the Laplace equation. Several theoretical studies (Brzustowski et al. 1979; Umemura, Ogawa & Oshima 1981a,b) have employed this approach, along with the flame-sheet assumption, to analyse the QS burning of two droplets. This leads to theoretical solutions for the potential flow field, burning rates and flame positions. The Stefan flow is taken into account, and the model allows for the two droplets to be of arbitrary sizes. It is found that for closely spaced droplets, there is typically a single enveloping flame, while at greater distances, each droplet is encased in its separate flame. As the distance between two droplets decreases, the evaporation and burning rates of each droplet decrease due to inter-droplet interactions.

The aforementioned studies based on the QS mass-flux-potential model (Labowsky 1980; Umemura et al. 1981a,b; Imaoka & Sirignano 2005; Sirignano 2007; Sirignano & Wu 2008) typically focus on scenarios with equal droplet temperatures, and hence equal fuel vapour concentrations at the droplet surfaces under thermal equilibrium assumptions. However, It is crucial to recognize that in real droplet evaporation or combustion processes, the temperatures of different droplets are not always identical. For instance, droplets with different initial temperatures and sizes reach equilibrium at varying rates. Moreover, multi-droplet systems may experience additional heat sources, such as thermal radiation from the environment (Millán-Merino, Fernández-Tarrazo & Sánchez-Sanz 2021) or heat conduction from fibre lattices (Yoshida et al. 2019; Mikami et al. 2023), commonly employed in experiments to support and fix droplets. These factors can lead to temperature disparities across droplets. When the droplet temperatures differ, especially in closely spaced droplets, condensation may occur on parts of the surface of the cooler droplets, in contrast to evaporation. Additionally, the conclusion that flame temperatures are uniform across the flame surface, based on the analysis of droplets with identical temperatures (Labowsky 1980; Sirignano 2014), is no longer valid and requires further investigation. Considering the crucial role of flame temperature in flame limit phenomena (Law 2006), studying interactions between droplets with different temperatures can enhance understanding of flame spreading, as well as local ignition and extinction within droplet clouds (Brzustowski, Sobiesiak & Wojcicki 1981; Mikami et al. 2006, 2023; Yoshida et al. 2019; Li et al. 2023a).

Reflecting on the above discussions, this paper conducts a theoretical analysis of the QS burning of two interacting droplets with different temperatures based on the mass-flux-potential model and the flame-sheet assumption, using the bispherical coordinate approach. The paper is structured as follows. In § 2, the physical problem is defined and formulated. In § 3, the numerical approach and parameters are presented. In § 4, a detailed theoretical analysis of the problem is provided, with a comprehensive comparison and discussion of the theoretical findings against simulation results.

2. Formulation

As depicted in figure 1, this study investigates the QS interactions between two burning single-component fuel droplets, designated as Droplets 1 and 2, in a stationary oxygen-containing environment. The droplets have different radii $r_1$ and $r_2$, and temperatures $T_{s,1}$ and $T_{s,2}$, with centre-to-centre distance $H$. This study simplifies its analysis under some reasonable assumptions. The QS mass-flux-potential model (Labowsky 1980; Annamalai & Ryan 1992; Sirignano 2014) is applied to describe the flow field outside the droplets. A single-step gas-phase reaction $\nu _F F + \nu _O O \rightarrow \nu _P P$ is considered. Employing the flame-sheet assumption, commonly used in diffusion-controlled combustion studies, the reactions are constrained to an infinitely thin zone, a consequence of the rapid chemical kinetics. Additionally, this study assumes constant thermophysical properties across all species, with a unity Lewis number signifying equivalent mass and thermal diffusivities. The effects of gravity and thermal radiation, and the Soret and Dufour effects, are ignored. Additionally, given the uniform temperatures and liquid concentrations along the surfaces of each droplet, the Marangoni effect can be disregarded.

Figure 1. Schematic of the interaction between two burning droplets and the bispherical coordinate system.

The steady-state continuity equation for the gas phase reads

(2.1)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho_g \boldsymbol{v}) = 0, \end{equation}

where $\rho _g$ denotes the gas density, and $\boldsymbol {v}$ is the gas velocity. Three coupling functions $\beta _k$ (i.e. $\beta _{FO}$, $\beta _{TF}$ and $\beta _{TO}$) can be introduced using the Shvab–Zeldovich formulation (Law 2006) as

(2.2a)$$\begin{gather} \beta_{FO} \equiv \frac{Y_F }{ W_F \nu_F } - \frac{Y_O }{ W_O \nu_O }, \end{gather}$$

(2.2b)$$\begin{gather}\beta_{TF} \equiv \frac{c_{p,g} T_g }{ q_c } + \frac{Y_F }{ W_F \nu_F }, \end{gather}$$

(2.2c)$$\begin{gather}\beta_{TO} \equiv \frac{c_{p,g} T_g }{ q_c } + \frac{Y_O }{ W_O \nu_O }, \end{gather}$$

where $T_g$ is the gas temperature, $Y_j$, $W_j$ and $\nu _j$ are the mass fraction, molar mass and stoichiometric coefficient of species $j$ ($\kern 0.06em j = F, O$), respectively, $c_{p,g}$ is the specific heat capacity of the gas, and $q_c$ refers to the reaction heat. Under flame-sheet assumption and unity Lewis number conditions, the three coupling functions satisfy the convection–diffusion equation without the chemical source term as

(2.3)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho_g \boldsymbol{v} \beta_k - \rho_g D\, \boldsymbol{\nabla} \beta_k) = 0, \end{equation}

where $\rho _g D$ is the mass and thermal diffusivity, assumed to be a fixed value equal to $(\rho _g D)_\infty$. By solving (2.3) and considering the definition of $\beta _k$ in (2.2), the temperature and component fields inside and outside the flame surface can be determined straightforwardly.

Based on the mass-flux-potential assumption, the mass flux in the gas phase can be expressed as

(2.4)\begin{equation} \rho_g \boldsymbol{v} ={-} \boldsymbol{\nabla} \varphi, \end{equation}

where $\varphi$ represents the potential function. Substituting (2.4) into the continuity equation (2.1), we obtain the Laplace equation for $\varphi$ as

(2.5)\begin{equation} \nabla^2 \varphi = 0. \end{equation}

Similarly, the coupling function equation (2.3) can be rewritten as

(2.6)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} ( \beta_k\,\boldsymbol{\nabla} \varphi + \rho_g D\, \boldsymbol{\nabla} \beta_k) = 0. \end{equation}

Finally, to finalize the problem, the gas state equation $P = \rho _g R_g T_g$ should be considered, with $R_g$ the gas constant.

Selecting an appropriate coordinate system is critical for effective analysis of specific physical problems. In this context, the bispherical coordinate system – an orthogonal curvilinear coordinate system – is well-suited for problems involving two spheres (Moon & Spencer 2012). Figure 1 illustrates the distribution of the three coordinate lines $(\xi, \theta, \psi )$ in the bispherical coordinates. Surfaces of constant $\xi$ consist of non-intersecting spheres, each with a distinct radius. In contrast, surfaces of constant $\theta$ align with intersecting tori, again each with a differing radius. Assuming axisymmetry for the current problem, i.e. $\partial /\partial \psi = 0$, the transformation between the bispherical coordinate system $(\xi, \theta, \psi )$ and the cylindrical coordinate system $(\sigma, z, \psi )$ reads

(2.7a,b)\begin{equation} \sigma = \frac{a \sin \theta }{ \cosh \xi - \cos \theta }, \quad z = \frac{a \sinh \xi }{ \cosh \xi - \cos \theta }, \end{equation}

where $\sinh (\cdot )$ and $\cosh (\cdot )$ are the hyperbolic sine and cosine functions, respectively. The $z$ axis in the cylindrical coordinate system, serving as the axis of rotational symmetry, coincides with the line connecting the centres of the two droplets. The parameter $a$ in (2.7), as a function of the inter-droplet distance $H$ and droplet radii $r_1$, $r_2$, is defined as

(2.8)\begin{equation} a = \frac{\sqrt{ (H + r_1 + r_2) (H + r_1 - r_2) (H - r_1 + r_2) (H - r_1 - r_2)} }{ 2 H}. \end{equation}

In the current coordinate system, the surfaces of two droplets are located on two specific $\xi$-constant surfaces as

(2.9a,b)\begin{equation} \xi_1 ={-} \sinh^{{-}1} \frac{a }{ r_1 } < 0, \quad \xi_2 = \sinh^{{-}1} \frac{a }{ r_2 } > 0. \end{equation}

Therefore, the range of the $\xi$ coordinate under consideration is $[\xi _1, \xi _2]$, encompassing the region from the surface of Droplet 1 to that of Droplet 2. Comparatively, the range of the $\theta$ coordinate spans from $0$ to ${\rm \pi}$, with $\theta = 0$ and $\theta = {\rm \pi}$ corresponding to different parts of the symmetrical axis. Specifically, $\theta = 0$ refers to the part of the axis lying outside the two droplets. Conversely, $\theta = {\rm \pi}$ aligns with the axis segment located between the two droplets. Notably, in the bispherical coordinates, the region at infinity is mapped to a single point, represented as $(\xi = 0, \theta = 0)$.

A non-dimensional inter-drop distance can be introduced as $\mathcal {H} \equiv H / (r_1 + r_2)$. When $\mathcal {H} \to 1$, which indicates that the droplets are on the verge of contact, $a \propto \sqrt {\mathcal {H} - 1} \to 0$. In this case, $\xi _1 \to 0^{-}$ and $\xi _2 \to 0^{+}$. As $\mathcal {H}$ increases, both $a$ and $|\xi _i|$ increase accordingly. In the limit where the droplets are significantly distant from each other (i.e. $\mathcal {H} \to \infty$), $a \sim H / 2 \to \infty$, while $\xi _1$ and $\xi _2$ approach $-\infty$ and $+\infty$, respectively.

In the bispherical coordinate system, the Laplace equation for $\varphi$, as shown in (2.5), is expressed as

(2.10a,b)\begin{equation} \frac{\partial}{\partial \xi} \left( \mathcal{K}\,\frac{\partial \varphi}{\partial \xi} \right) + \frac{\partial}{\partial \theta} \left( \mathcal{K}\,\frac{\partial \varphi}{\partial \theta} \right) = 0, \quad \mathcal{K} (\xi, \theta) \equiv \frac{\sin \theta }{ \cosh \xi - \cos \theta }. \end{equation}

The coupling function equation (2.6) is represented as

(2.11)\begin{equation} \mathcal{K} \left( \frac{\partial \varphi}{\partial \xi}\,\frac{\partial \beta_k}{\partial \xi} + \frac{\partial \varphi}{\partial \theta}\,\frac{\partial \beta_k}{\partial \theta} \right) + \rho_g D \left[ \frac{\partial}{\partial \xi} \left( \mathcal{K}\,\frac{\partial \beta_k}{\partial \xi} \right) + \frac{\partial}{\partial \theta} \left( \mathcal{K}\,\frac{\partial \beta_k}{\partial \theta} \right) \right] = 0. \end{equation}

The relationship between the mass flux and $\varphi$ reads

(2.12)\begin{equation} (\rho_g v_\xi, \ \rho_g v_\theta) ={-} \left( \frac{\cosh \xi - \cos \theta }{ a } \right) \left( \frac{\partial \varphi}{\partial \xi}, \frac{\partial \varphi}{\partial \theta} \right), \end{equation}

where $v_\xi$ and $v_\theta$ represent the components of the velocity in the $\xi$ and $\theta$ directions, respectively.

The boundary conditions at infinity are $P = P_\infty$, $T_g = T_{\infty }$, $Y_F = Y_{F,\infty } = 0$, $Y_O = Y_{O,\infty }$. Here, $P$ is the ambient pressure. The boundary conditions on the surfaces of the two droplets are

(2.13)\begin{equation} \left.\begin{array}{@{}c@{}} T_g = T_{s,i}, \quad Y_F = Y_{F,s,i}, \\[6pt] \dot{m}_{ev, i} (1 - Y_{F,s,i}) ={-} \rho_g D \left( \dfrac{\cosh \xi - \cos \theta }{ a }\, \dfrac{\partial Y_F}{\partial \xi} \right)_{s,i},\\[16pt] \dot{m}_{ev, i} Y_{j,s,i} = \rho_g D \left( \dfrac{\cosh \xi - \cos \theta }{ a }\, \dfrac{\partial Y_j}{\partial \xi} \right)_{s,i}, \quad j \neq F, \end{array}\right\} \end{equation}

where $T_{s,i}$ and $Y_{F,s,i}$ are the temperature and fuel vapour mass fraction at the surface of droplet $i$, respectively, and $\dot {m}_{ev,i}$ represents the local evaporation rate, expressed as $\dot {m}_{ev,i} \equiv (\rho _g v_\xi ) |_{s,i}$. In the presence of an enclosed flame, the oxygen mass fraction at the droplet surface is zero, i.e. $Y_{O,s,i} = 0$. To investigate the influence of interactive burning droplets with different temperatures, the temperatures of both droplets are held constant at prescribed values. Given the absence of temperature and liquid concentration gradients along the respective surfaces of the two droplets, we can disregard the influence of the Marangoni effect. Assuming thermodynamic equilibrium at the liquid surface, the fuel vapour mass fraction at the surfaces of the two droplets can be calculated based on their temperatures using the Clausius–Clapeyron relation:

(2.14)\begin{equation} Y_{F,s,i} = \left( \frac{W_F }{ W_{mix} } \right) \exp \left[ \frac{q_L }{ R_u } \left( \frac{1 }{ T_{BP} } - \frac{1 }{ T_{s,i} } \right) \right], \end{equation}

where $W_{mix}$ is the average molar mass of the mixture, $R_u$ is the universal gas constant, and $q_L$ and $T_{BP}$ are the latent heat and boiling point of the liquid fuel, respectively.

3. Numerical approach and parameters

To verify the accuracy of the theoretical analysis presented in this paper, extensive simulations are conducted by using the second-order finite difference approach to iteratively solve the coupled gas-phase governing equations in the bispherical coordinate system as outlined in § 2. A variety of critical parameters for the dual-droplet system are explored, including different droplet temperatures, sizes and inter-drop distances. In the explicit iteration process, we set the relative convergence tolerance for both the potential function $\varphi$ and the coupling function $\beta _k$ to $10^{-6}$. A uniform grid division is generated for the rectangle region ($\xi \in [\xi _1, \xi _2]$ and $\theta \in [0, {\rm \pi}]$), with $200$ grid points in both coordinate directions. Considering that the bispherical coordinate system has denser grids near the droplets and sparser grids farther from the droplets (refer to figure 1), the generated mesh is sufficient to accurately resolve the interactions between droplets. Mesh-independence tests were conducted based on representative examples to assess and confirm the accuracy and reliability of the simulation.

Referring to the thermophysical parameters for liquid n-heptane/air systems (Lide 2004; Li, Zhang & Law 2023b), the subsequent results employ the following parameters unless stated otherwise: $P = 1$ atm, $(\rho _g D)_\infty = 5.8 \times 10^{-5}$ kg m$^{-1}$ s$^{-1}$, $c_{p,g} = 1100$ J kg$^{-1}$ K$^{-1}$, $q_L = 3.4 \times 10^4$ J mol$^{-1}$, $q_c = 6 \times 10^5$ J mol$^{-1}$, $T_{BP} = 371$ K, $T_{g, \infty } = 300$ K, $Y_{F,\infty } = 0$, $W_F = 100$ g mol$^{-1}$, $W_{i \neq F} = 30$ g mol$^{-1}$, $\nu _F = 1$, $\nu _O = 11$, $\nu _P = 15$ and $Y_{O,\infty } = 1$. Note that the theoretical analysis presented in the next section is applicable beyond the particular conditions outlined above.

4. Theoretical analysis

4.1. Solution for potential function $\varphi$

Considering the derivative relationship between gas-phase velocity and the scalar potential $\varphi$, as described in (2.4), we can, without compromising the generality of the solution, set the value of $\varphi$ to zero at infinity, i.e. $\varphi (\xi = 0, \theta = 0) = 0$. When $\varphi$ maintains constant values, $\varphi _{s,1}$ and $\varphi _{s,2}$, on the surfaces of two droplets, defined at $\xi = \xi _1$ and $\xi = \xi _2$, respectively, the Laplace equation of $\varphi$, given by (2.10), can be solved using the separation of variables method (Morse & Feshbach 1954). This approach yields the exact solution for $\varphi$ in the form of an infinite series as

(4.1)\begin{align} \varphi &= \sqrt{ 2(\cosh \xi - \cos \theta) } \sum_{n=0}^\infty \left[ \frac{( \varphi_{s,1} \,\mathrm{e}^{(2n+1)\xi_1} - \varphi_{s,2} )}{\mathrm{e}^{(2n+1)\xi_1} - \mathrm{e}^{(2n+1)\xi_2}} \,\mathrm{e}^{(n + 1/2) \xi}\right.\nonumber\\ &\quad\left. {}+ \frac{\mathrm{e}^{(2n+1)\xi_1}\,( \varphi_{s,2} - \varphi_{s,1} \,\mathrm{e}^{(2n+1)\xi_2} ) }{ \mathrm{e}^{(2n+1)\xi_1} - \mathrm{e}^{(2n+1)\xi_2} } \,\mathrm{e}^{-(n + 1/2) \xi} \right] P_n (\cos \theta), \end{align}

where $P_n (\cdot )$ represents the $n$th-order Legendre polynomial.

Because the mass flux at the droplet surface is affected by the local gradient of the fuel vapour mass fraction, the values of $\varphi _{s,i}$ have not yet been determined. In fact, even the uniformity of $\varphi _{s,i}$ across the surfaces of the droplets cannot be ensured. To address this, it is necessary to incorporate the convection–diffusion equation of the coupling function $\beta _k$, as specified in (2.11). The corresponding boundary conditions on both droplet surfaces should be applied to derive coupled solutions for $\varphi$ and $\beta _k$.

4.2. Solution for coupling function $\beta _{FO}$

Considering the boundary conditions for $T_g$, $Y_F$ and $Y_O$ at the surfaces of both droplets and at infinity, the boundary conditions for the three coupling functions, as defined in (2.2), read

(4.2)\begin{equation} \left.\begin{array}{@{}c@{}} \beta_{FO,s,1} = \dfrac{Y_{F,s,1} }{ W_F \nu_F }, \quad \beta_{TF,s,1} = \dfrac{c_{p,g} T_{s,1} }{ q_c } + \dfrac{Y_{F,s,1} }{ W_F \nu_F }, \quad \beta_{TO,s,1} = \dfrac{c_{p,g} T_{s,1} }{ q_c }, \\[10pt] \beta_{FO,s,2} = \dfrac{Y_{F,s,2} }{ W_F \nu_F }, \quad \beta_{TF,s,2} = \dfrac{c_{p,g} T_{s,2} }{ q_c } + \dfrac{Y_{F,s,2} }{ W_F \nu_F }, \quad \beta_{TO,s,2} = \dfrac{c_{p,g} T_{s,2} }{ q_c }, \\[10pt] \beta_{FO,\infty} ={-} \dfrac{Y_{O,\infty} }{ W_O \nu_O }, \quad \beta_{TF,\infty} = \dfrac{c_{p,g} T_\infty }{ q_c }, \quad \beta_{TO,\infty} = \dfrac{c_{p,g} T_\infty }{ q_c } + \dfrac{Y_{O,\infty} }{ W_O \nu_O }. \end{array}\right\} \end{equation}

One possible solution for $\beta _k$ to (2.6) can be represented in the form

(4.3a,b)\begin{equation} \beta_k = A_k + B_k \times \mathcal{L} (\varphi), \quad \mathcal{L} (\varphi) \equiv \exp \left[ - \frac{\varphi }{ (\rho_g D)_\infty } \right], \end{equation}

where $A_k$ and $B_k$ are undetermined coefficients. This form of solution indicates that $\beta _k$ is a one-dimensional function of $\varphi$. Since (4.3) contains only two undetermined coefficients, while there are three boundary conditions to satisfy (see (4.2)), the validity of this general solution hinges on the existence of a specific relationship among these three boundary conditions.

Given that only fuel vapour exhibits mass flux at the surfaces of the droplets, i.e.

(4.4)\begin{equation} \left.\begin{array}{@{}c@{}} \dot{m}_{F,s,i} = \boldsymbol{n} \boldsymbol{\cdot} ( \rho_g \boldsymbol{v} Y_F - (\rho_g D)_\infty\,\boldsymbol{\nabla} Y_F )_{s,i} = ( \rho_g \boldsymbol{v})_{s,i},\\[3pt] \dot{m}_{O,s,i} = \boldsymbol{n} \boldsymbol{\cdot} ( \rho_g \boldsymbol{v} Y_O - (\rho_g D)_\infty\,\boldsymbol{\nabla} Y_O )_{s,i} = 0, \end{array}\right\} \end{equation}

where $\boldsymbol {n}$ represents the unit normal vector on the droplet surface. Considering (2.12) and (4.4), we have

(4.5a,b)\begin{equation} \left. (\rho_g D)_\infty\,\frac{\mathrm{d} Y_F}{\mathrm{d} \varphi} \right|_{s,i} = 1 - Y_{F,s,i}, \quad \left. {-(\rho_g D)_\infty}\,\frac{\mathrm{d} Y_O}{\mathrm{d} \varphi} \right|_{s,i} = Y_{O,s,i} = 0. \end{equation}

Then by using the relation between $\beta _{FO}$ and $Y_j$, we obtain

(4.6)\begin{equation} \left. \frac{\mathrm{d} \beta_{FO}}{\mathrm{d} \varphi} \right|_{s,i} = \frac{1 }{ (\rho_g D)_\infty }\,\frac{1 - Y_{F,s,i} }{ W_F \nu_F }. \end{equation}

Alternatively, if $\beta _{FO}$ adheres to the solution form provided in (4.3), based on the boundary conditions of $\beta _{FO}$ as outlined in (4.2), then we deduce that

(4.7)\begin{equation} \left. \frac{\mathrm{d} \beta_{FO}}{\mathrm{d} \varphi} \right|_{s,i} ={-} \frac{1 }{ (\rho_g D)_\infty } \left( \frac{Y_{O,\infty} }{ W_O \nu_O } + \frac{Y_{F,s,1} }{ W_F \nu_F } \right) \frac{\mathcal{L} (\varphi_{s,i}) }{ \mathcal{L} (\varphi_{s,1}) - 1 }. \end{equation}

Combining (4.6) and (4.7), we obtain

(4.8)\begin{equation} \frac{1 - Y_{F,s,i} }{ W_F \nu_F } ={-} \left( \frac{Y_{O,\infty} }{ W_O \nu_O } + \frac{Y_{F,s,1} }{ W_F \nu_F } \right) \frac{\mathcal{L} (\varphi_{s,i}) }{ \mathcal{L} (\varphi_{s,1}) - 1 }, \end{equation}

which is the relation that must be satisfied on the surfaces of both droplets. Further combining the expression of $\mathcal {L} (\varphi )$ from (4.3), the values of $\varphi$ on the droplet surfaces can be derived as

(4.9)\begin{equation} \varphi_{s,i} = (\rho_g D)_\infty \ln ( 1 + B_{v,i} ), \end{equation}

where $\varphi _{s,i}$ and $B_{v,i}$ represent the evaporation potential and Spalding transfer number of droplet $i$, respectively, and $B_{v,i}$ is defined as

(4.10)\begin{equation} B_{v,i} = \frac{Y_{F,s,i} + \phi Y_{O,\infty} }{ 1 - Y_{F,s,i} }, \quad \phi \equiv \frac{W_F \nu_F }{ W_O \nu_O }, \end{equation}

where $\phi$ is the fuel-to-oxygen mass stoichiometric ratio. According to (4.10), $B_{v,i}$ is a function of $\phi Y_{O,\infty }$ and $Y_{F,s,i}$, which further depends on $T_{s,i}$. Ultimately, with the boundary conditions of $\beta _{FO}$, its solution can be determined as

(4.11)\begin{equation} \beta_{FO} (\varphi) ={-} \frac{Y_{O,\infty} }{ W_O \nu_O } + \left( \frac{Y_{O,\infty} }{ W_O \nu_O } + \frac{Y_{F,s,1} }{ W_F \nu_F } \right) \frac{\mathcal{L} (\varphi) - 1 }{ \mathcal{L} (\varphi_{s,1}) - 1 }. \end{equation}

In light of the preceding derivation, provided that the values of $\varphi$ on the droplet surfaces conform to (4.9), $\beta _{FO}$ has a solution given by (4.11). Concurrently, the infinite series solution for $\varphi$, as specified in (4.1), is applicable. Considering the uniqueness theorem for partial differential equations, the above provides the exact solution for the current problem. As a result, we have successfully derived the solutions for both $\varphi$ and $\beta _{FO}$, applicable to pairs of droplets with differing evaporation potentials. Our analysis confirms that in such scenarios, the coupling function $\beta _{FO}$ preserves a one-to-one correspondence with the potential function, which has been established previously in studies addressing scenarios with identical evaporation potentials (Imaoka & Sirignano 2005; Sirignano 2014). It is important to note, however, that the other two coupling functions, $\beta _{TF}$ and $\beta _{TO}$, no longer maintain a one-to-one correspondence with $\varphi$ in instances where $\varphi _{s,1} \neq \varphi _{s,2}$, which will be further discussed in this paper.

4.3. Asymptotic solution for potential function $\varphi$

The infinite series solution for the potential function $\varphi$, detailed in (4.1), exhibits a slow convergence rate, particularly notable at considerably small inter-droplet distances. Therefore, a substantial number of terms must be considered to achieve satisfactory accuracy. Under such conditions, the practicality of the series solution for analytical exploration becomes limited due to the complexity involved in assessing the impact of various parameters. In response to this issue, we have performed a perturbation analysis of (4.1) (detailed in Appendix A), leading to the derivation of the following asymptotic solution for $\varphi$:

(4.12)\begin{align} \varphi &\sim \left[ \frac{\mathcal{A} (\xi, \theta) }{ \mathcal{A} (\xi - 2 \xi_1, \theta) } - \frac{\mathcal{B} (\xi)\,\mathcal{A} (\xi, \theta) }{ \mathcal{A} (\xi_2 - 2 \xi_1, \theta) } \right] \varphi_{s,1}\nonumber\\ &\quad + \left[ \frac{\mathcal{A} (\xi, \theta) }{ \mathcal{A} (2\xi_2 - \xi, \theta) } - \frac{( 1 - \mathcal{B} (\xi) )\, \mathcal{A} (\xi, \theta) }{ \mathcal{A} (2\xi_2 - \xi_1, \theta) } \right] \varphi_{s,2}, \end{align}

where

(4.13a,b)\begin{equation} \mathcal{A} (\xi, \theta) \equiv \sqrt{ \cosh \xi - \cos \theta }, \quad \mathcal{B} (\xi) \equiv \frac{\xi - \xi_1 }{ \xi_2 - \xi_1 } \in [0, 1]. \end{equation}

The asymptotic solution, expressed in (4.12), satisfies the boundary conditions at both droplet surfaces ($\xi = \xi _i$) and at infinity ($\xi = \theta = 0$). It distinctively elucidates the influence of each droplet's evaporation potential on $\varphi$ in different external regions. This solution contrasts with the series solution in (4.1), which demands numerous terms to achieve a comparable level of accuracy. Notably simpler, the asymptotic solution in (4.12) employs only a single term and yet effectively approximates the exact solution. To the best of our knowledge, this particular solution has not been reported previously. In cases where the two droplets have identical radii ($r_1 = r_2 \Rightarrow \xi _1 = - \xi _2$), (4.12) can be simplified further as

(4.14)\begin{align} \varphi |_{r_1 = r_2}&\sim \mathcal{A} (\xi, \theta) \left[ \frac{\varphi_{s,1} }{ \mathcal{A} (2 \xi_2 + \xi, \theta) } + \frac{\varphi_{s,2} }{ \mathcal{A} (2\xi_2 - \xi, \theta) } \right.\nonumber\\ &\quad \left.{}- \frac{(\varphi_{s,1} + \varphi_{s,2}) \, \xi_2 + (\varphi_{s,1} - \varphi_{s,2}) \, \xi }{ 2 \xi_2 \, \mathcal{A} (3\xi_2, \theta) } \right]. \end{align}

Figures 2 and 3 present a comparison between the asymptotic solution (see (4.12)) and the numerical solution for $\varphi$, demonstrating a reasonably good agreement across various parameter conditions. Since the infinite series solution (see (4.1)) closely matches the numerical solution when up to $32$ terms in the series are considered, it is not illustrated separately in these contour plots. In these plots, the direction of the gradient descent in $\varphi$ indicates the direction of the local mass flux, with larger absolute values of the gradient denoting higher local mass flow rates. At relatively large distances between the droplets, as shown in figures 2(a,c) and 3(a,c), the distribution of $\varphi$ near each droplet's surface appears almost spherically symmetric, suggesting minimal droplet interaction and nearly uniform evaporation rates across each droplet surface. In contrast, at shorter inter-droplet distances, as seen in figures 2(b,d) and 3(b,d), the flow field around each droplet is significantly influenced by the presence of the other, particularly for the smaller droplet in the dual-droplet system.

Figure 2. Comparative distribution of the dimensionless potential function $\varphi / \varphi _{s,1}$ for different droplet sizes and distances, with $r_1 = 100$ ${\rm \mu}$m. The evaporation potentials of the two droplet surfaces are identical, i.e. $\varphi _{s,2} / \varphi _{s,1} = 1$. The asymptotic solutions (see (4.12)) are presented alongside corresponding numerical solutions for comparison. Plots are for (a) $r_2/r_1=1$, ${\mathcal {H}}=2$, (b) $r_2/r_1=1$, ${\mathcal {H}}=1.1$, (c) $r_2/r_1=2$, ${\mathcal {H}}=2$, and (d) $r_2/r_1=2$, ${\mathcal {H}}=1.1$.

Figure 3. Comparative distribution of the dimensionless potential function $\varphi / \varphi _{s,1}$ for different droplet sizes and distances, with $r_1 = 100$ ${\rm \mu}$m. The two droplet surfaces have differing evaporation potentials, specifically, $\varphi _{s,2} / \varphi _{s,1} = 1.5$. The asymptotic solutions (see (4.12)) are presented alongside corresponding numerical solutions for comparison. Plots are for (a) $r_2/r_1=1$, ${\mathcal {H}}=2$, (b) $r_2/r_1=1$, ${\mathcal {H}}=1.1$, (c) $r_2/r_1=2$, ${\mathcal {H}}=2$, and (d) $r_2/r_1=2$, ${\mathcal {H}}=1.1$.

As depicted in figures 3(b,d), in scenarios where $\varphi _{s,1} \neq \varphi _{s,2}$ and $\mathcal {H}$ is relatively small ($\mathcal {H} = 1.1$), it is observed that Droplet 1, which has a lower evaporation potential, exhibits a negative evaporation rate on its side facing the other droplet. This phenomenon essentially manifests as condensation on the surface of one droplet, a result of the high fuel vapour concentration on the surface of its adjacent counterpart. The occurrence of local condensation arises wherever the partial vapour pressure near a liquid surface segment surpasses the saturation pressure. To the best of our knowledge, this effect has not been addressed theoretically in the existing literature on interactive burning droplets, as it is absent in scenarios where both droplets possess identical surface temperatures and consequently the same evaporation potentials.

Our subsequent focus is the distribution of $\varphi$ in certain critical areas for the present problem. These areas of interest include primarily the region along the symmetrical axis, particularly the portion of the axis extending outside the two droplets (indicated by $\theta = 0$), as well as the portion of the axis situated between the two droplets (denoted by $\theta = {\rm \pi}$). The interaction is comparatively weak on the sides of the droplets that are furthest apart from each other. In contrast, the area directly between the droplets, where they face each other, is characterized by the most intense interaction. Based on (4.12), for the portion of the axis extending outside the two droplets ($\theta = 0$), we obtain

(4.15)\begin{align} \varphi (\theta = 0) &\sim |\mathrm{e}^{\xi} - 1 | \times \left[ \left( \frac{1 }{ \mathrm{e}^{\xi - 2\xi_1} - 1 } - \frac{\mathrm{e}^{(\xi_2 - \xi)/2} }{ \mathrm{e}^{\xi_2 - 2\xi_1} - 1 }\,\mathcal{B} (\xi) \right)\mathrm{e}^{- \xi_1}\,\varphi_{s,1} \right.\nonumber\\ &\quad \left.{}+ \left( \frac{1 }{ \mathrm{e}^{2\xi_2 - \xi} - 1 } - \frac{\mathrm{e}^{(- \xi_1 + \xi)/2} }{ \mathrm{e}^{2\xi_2 - \xi_1} - 1 }\,( 1 - \mathcal{B} (\xi) ) \right) \mathrm{e}^{\xi_2 - \xi}\,\varphi_{s,2} \right], \end{align}

while for the portion of the axis located between the two droplets ($\theta = {\rm \pi}$), we have

(4.16)\begin{align} \varphi (\theta = {\rm \pi}) &\sim ( \mathrm{e}^{\xi} + 1 ) \times \left[ \left( \frac{1 }{ \mathrm{e}^{\xi - 2\xi_1} + 1 } - \frac{\mathrm{e}^{(\xi_2 - \xi)/2} }{ \mathrm{e}^{\xi_2 - 2\xi_1} + 1 }\,\mathcal{B} (\xi) \right) \mathrm{e}^{- \xi_1}\,\varphi_{s,1} \right.\nonumber\\ &\quad \left.{}+ \left( \frac{1 }{ \mathrm{e}^{2\xi_2 - \xi} + 1 } - \frac{\mathrm{e}^{(- \xi_1 + \xi)/2} }{ \mathrm{e}^{2\xi_2 - \xi_1} + 1 }\,( 1 - \mathcal{B} (\xi) ) \right) \mathrm{e}^{\xi_2 - \xi}\,\varphi_{s,2} \right]. \end{align}

Figure 4 displays the asymptotic solutions for $\varphi (\theta = 0)$ and $\varphi (\theta = {\rm \pi})$, as detailed in (4.15) and (4.16). These solutions demonstrate good agreement with the numerical solutions, particularly for the portion of the axis extending outside the two droplets ($\theta = 0$), where the asymptotic solution in (4.15) aligns perfectly with the numerical result. As shown in figure 4(a), a smaller inter-droplet distance results in reduced $|\xi |_i$ values, leading to larger gradients of $\varphi$ on the droplet surfaces. This signifies a more intense mass flux on the surface of both droplets in the facing area. According to figure 4(b), when the droplets are closely spaced (i.e. $\mathcal {H} = 1.1$), the variation of $\varphi$ along the portion of the axis located between the two droplets ($\theta = {\rm \pi}$) is generally linear. In this case, the gradients of $\varphi$ on the surfaces of the two droplets exhibit opposite signs, suggesting evaporation on the droplet, with a higher evaporation potential and condensation on the other. As the distance between the droplets increases (e.g. $\mathcal {H} = 3$), this distribution evolves into a concave shape, with a minimum $\varphi$ value between the two droplets. Under these circumstances, the closest points on both droplet surfaces undergo evaporation, leading to the inference that condensation does not occur on either surface.

Figure 4. Distribution of the dimensionless potential function $\varphi / \varphi _{s,i}$ along the symmetrical axis: (a) the portion of the axis extending outside the two droplets (indicated by $\theta = 0$), and (b) the portion of the axis situated between the two droplets (denoted by $\theta = {\rm \pi}$). Here, $r_1 = 100$ ${\rm \mu}$m and $\varphi _{s,2} / \varphi _{s,1} = 1.5$. The scatter data represent numerical results, while the solid lines depict asymptotic solutions provided by (4.15) and (4.16).

Referring to the derivations presented in Appendix A, when the inter-droplet distance becomes substantially large ($\mathcal {H} \to \infty$), the asymptotic solution for $\varphi$ (4.12) can be simplified further. Specifically, for the areas near Droplets 1 and 2, as well as for the region not in close proximity to the surfaces of either droplet, we respectively have

(4.17)\begin{gather} \left.\begin{gathered} \varphi |_{\mathcal{H} \to \infty} ( \xi \to \xi_1 ) \sim \mathrm{e}^{- (\xi - \xi_1)}\, \varphi_{s,1} + ( 1 - \mathrm{e}^{- (\xi - \xi_1)/2} )\, \mathrm{e}^{- \xi_2}\,\varphi_{s,2}, \\ \varphi |_{\mathcal{H} \to \infty} ( \xi \to \xi_2 ) \sim ( 1 - \mathrm{e}^{(\xi - \xi_2)/2} ) \,\mathrm{e}^{\xi_1}\,\varphi_{s,1} + \mathrm{e}^{\xi - \xi_2}\,\varphi_{s,2}, \\ \varphi |_{\mathcal{H} \to \infty} ( \xi \not\to \xi_1, \xi_2 ) \sim ( 1 - \mathrm{e}^{-(\xi_2 - \xi)/2} + \mathrm{e}^{-\xi} ) \,\mathrm{e}^{\xi_1}\, \varphi_{s,1} + ( 1 - \mathrm{e}^{- (\xi - \xi_1)/2} + \mathrm{e}^{\xi}) \,\mathrm{e}^{- \xi_2}\,\varphi_{s,2}, \end{gathered}\right\} \end{gather}

which intuitively demonstrates the interaction between two distant droplets. As stated in (4.17), $\varphi$ is an exclusive function of $\xi$ and remains independent of $\theta$. This implies that when two droplets are significantly separated, the evaporation flux along their surfaces is nearly uniformly distributed. The uniform distribution arises because a large $\mathcal {H}$ value leads to the evaporation of one droplet generating a weak yet consistent concentration of fuel vapour around the other droplet. In such circumstances, the evaporation (or burning) behaviour of each droplet is predominantly governed by its own evaporation potential. The presence of the other droplet exerts a relatively minor effect, as reflected in the expression for $\varphi$ at the limits $\xi \to \xi _1$ and $\xi \to \xi _2$ in (4.17).

4.4. Droplet evaporation/burning rates and interaction coefficients

4.4.1. Droplet evaporation/burning rates

By employing the derived solution for $\varphi$, and leveraging the relationship between local mass flux and $\varphi$ (see (2.12)), the local evaporation rate on the two droplet surfaces can be obtained with

(4.18)\begin{equation} \dot{m}_{ev,i} (\theta) = ({-}1)^{i} \times \left( \frac{\cosh \xi - \cos \theta }{ a }\, \frac{\partial \varphi}{\partial \xi} \right)_{\xi = \xi_i}, \end{equation}

where $(-1)^i$ accounts for the direction of evaporation flux on droplet surfaces. Owing to the axial symmetry of the problem, the evaporation rate $\dot {m}_{ev,i}$ varies as a function of $\theta$, a relationship depicted in figure 5. Note that the horizontal axis in the figure represents the $\theta$ axis in the bispherical coordinate system $(\xi, \theta )$ (refer to figure 1). This axis does not exhibit a uniform variation across the droplet surface. The results shown in figure 5 correspond to a pair of droplets of identical size, where Droplet 2 possesses a slightly larger Spalding number, i.e. $B_{v,2} > B_{v,1}$. It is found that the asymptotic solutions from (4.12) and (4.18) effectively capture the local evaporation rates for both droplets.

Figure 5. Evaporation rate distribution on two equal-sized droplets at varying distances ($r_1 = r_2 = 100$ $\mathrm {\mu }$m; $B_{v,1} = 3$, corresponding to $T_{s,1} \approx 341$ K and $Y_{F,s,1} \approx 0.67$; $B_{v,2} = 4$, corresponding to $T_{s,2} \approx 347$ K and $Y_{F,s,2} \approx 0.74$). The scatter data represent numerical solutions, while the solid lines denote asymptotic solutions provided by (4.12) and (4.18).

Figure 5 reveals distinct $\dot {m}_{ev,i} (\theta )$ distributions due to different droplet temperatures. Droplet 1, at a lower temperature, exhibits a less uniform distribution of evaporation rate compared to the warmer Droplet 2, implying a more significantly impacted evaporation/burning process. When the droplets are distantly spaced (e.g. $\mathcal {H} > 1.5$, not shown in figure 5 for clarity), both possess positive evaporation rates across their surfaces, though the rate diminishes near the facing areas. As the droplets move closer, $\dot {m}_{ev,i} (\theta )$ becomes progressively non-uniform. According to figure 5, as $\mathcal {H}$ falls below a critical threshold, a portion of the Droplet 1 surface facing Droplet 2 experiences condensation instead of evaporation, a phenomenon that intensifies with decreasing $\mathcal {H}$. In contrast, Droplet 2 demonstrates a non-monotonic evaporation rate distribution, featuring an elevated rate on its side facing Droplet 1. Overall, the proximity of the droplets reduces the system's total evaporation rate, attributable primarily to the decreased evaporation rate of the cooler droplet.

By integrating $\dot {m}_{ev,i}$ from (4.18) over the surface of each droplet, their cumulative evaporation rates can be derived as

(4.19)\begin{equation} \dot{Q}_i = \oint_{s,i} \dot{m}_{ev,i} \, \mathrm{d} S = ({-}1)^{i} \times 2 {\rm \pi}a \int_0^{\rm \pi} \left( \mathcal{K} (\xi, \theta)\,\frac{\partial \varphi}{\partial \xi} \right)_{\xi = \xi_i} \mathrm{d} \theta. \end{equation}

Substituting the solution for $\varphi$ (see (4.1)) into (4.19) and integrating yields the cumulative evaporation rates for the two droplets:

(4.20a)$$\begin{gather} \dot{Q}_1 = 8 {\rm \pi}a \sum_{n=0}^\infty \frac{\mathrm{e}^{(2n+1)\xi_1}\,( \varphi_{s,2} - \varphi_{s,1} \,\mathrm{e}^{(2n+1)\xi_2} ) }{ \mathrm{e}^{(2n+1)\xi_1} - \mathrm{e}^{(2n+1)\xi_2} }, \end{gather}$$

(4.20b)$$\begin{gather}\dot{Q}_2 = 8 {\rm \pi}a \sum_{n=0}^\infty \frac{\varphi_{s,1} \,\mathrm{e}^{(2n+1)\xi_1} - \varphi_{s,2} }{ \mathrm{e}^{(2n+1)\xi_1} - \mathrm{e}^{(2n+1)\xi_2} }. \end{gather}$$

According to (4.20), the sign of $\dot {Q}_i$ is determined by the combined effect of the evaporation potentials of both droplet surfaces. Specifically, $\dot {Q}_i$ can be either positive, signifying evaporation, or negative, indicating condensation. This suggests that a droplet may experience either evaporation or condensation collectively, a variability that is absent when the surfaces of both droplets possess identical temperatures. Furthermore, it is crucial to note that the fuel vapour evaporating from each droplet might not fully react with oxygen; instead, it may condense on the surface of the adjacent droplet. Therefore, $\dot {Q}_i$ should be interpreted strictly as the evaporation rate of droplet $i$ and not as a measure of its burning rate.

It is noted that localized condensation on the surface of a droplet can occur when its temperature is lower than that of an adjacent droplet, as depicted in figure 5. However, for the cumulative evaporation rate of a droplet to become negative (i.e. $\dot {Q}_i < 0$), more stringent criteria, as indicated in (4.20), must be satisfied. These criteria necessitate the droplets being in close proximity, with one droplet being comparatively smaller and having a substantially lower temperature than the other.

Utilizing the individual evaporation rates of the two droplets as specified in (4.20), the total evaporation rate of the entire dual-droplet system can be calculated accordingly as

(4.21)\begin{equation} \dot{Q}_{tot} = 8 {\rm \pi}a \sum_{n=0}^\infty \frac{\mathrm{e}^{(2n+1)\xi_1}\,( 1 - \mathrm{e}^{(2n+1)\xi_2} ) \varphi_{s,1} + ( \mathrm{e}^{(2n+1)\xi_1} - 1 ) \varphi_{s,2} }{ \mathrm{e}^{(2n+1)\xi_1} - \mathrm{e}^{(2n+1)\xi_2} } > 0. \end{equation}

Taking into account the signs of $\xi _i$, we can deduce from (4.21) that $\dot {Q}_{tot}$ is invariably positive. This suggests that in scenarios where two adjacent droplets exhibit surface temperature disparities, the system predominantly remains in an evaporation state. This observation persists even when the cooler droplet experiences condensation collectively. Since $\dot {Q}_{tot}$ inherently accounts for the mass transfer between the two droplets, it can also be interpreted as the overall burning rate of the dual-droplet system.

For the special case where the evaporation potentials of the two droplet surfaces are equal ($\varphi _{s,1} = \varphi _{s,2} \equiv \varphi _s$), the individual cumulative evaporation rates of the two droplets, as expressed in (4.20), can be reduced to:

(4.22a)$$\begin{gather} \dot{Q}_1 = 8 {\rm \pi}a \varphi_s \sum_{n=0}^\infty \frac{\mathrm{e}^{(2n+1)\xi_1} ( 1 - \mathrm{e}^{(2n+1)\xi_2} ) }{ \mathrm{e}^{(2n+1)\xi_1} - \mathrm{e}^{(2n+1)\xi_2} }, \end{gather}$$

(4.22b)$$\begin{gather}\dot{Q}_2 = 8 {\rm \pi}a \varphi_s \sum_{n=0}^\infty \frac{\mathrm{e}^{(2n+1)\xi_1} - 1 }{ \mathrm{e}^{(2n+1)\xi_1} - \mathrm{e}^{(2n+1)\xi_2} }. \end{gather}$$

According to (4.22), when $\varphi _{s,1} = \varphi _{s,2}$, both $\dot {Q}_i$ values are greater than $0$, indicating the absence of condensation on each droplet.

Based on (4.20), figure 6 depicts the variation of the cumulative evaporation rates of the two droplets of equal size as their separation distance increases. The Spalding number for Droplet 2 either matches or surpasses that of Droplet 1, leading to the Droplet 2 total evaporation rate being at least as high as that of Droplet 1. In particular, when $B_{v,1} = B_{v,2} = 4$, the evaporation rates of both droplets equalize and converge to a finite value of approximately $\dot {Q}_i \approx 0.8$ kg s$^{-1}$, provided that the droplets are in close proximity ($\mathcal {H} = 1$). As $\mathcal {H}$ increases, their evaporation rates also rise, eventually approaching the limiting value for isolated droplet combustion, estimated at $\dot {Q}_i \approx 1.2$ kg s$^{-1}$.

Figure 6. Theoretical solution for the variation in evaporation rates for two equal-sized droplets at varying distances ($r_1 = r_2 = 100$ $\mathrm {\mu }$m). Here, $B_{v,1} = 4$ is held constant (corresponding to $T_{s,1} \approx 347$ K and $Y_{F,s,1} \approx 0.74$), while $B_{v,2}$ varies among $[4,5,6,7]$ (corresponding to $T_{s,2} \approx [347, 350, 353, 355]$ K and $Y_{F,s,2} \approx [0.74, 0.78, 0.81, 0.84]$).

On the other hand, in cases where $B_{v,2} > B_{v,1}$, the deviation of $\dot {Q}_1$ from $\dot {Q}_2$ increases as the gap between $B_{v,1}$ and $B_{v,2}$ widens. In close proximity, the evaporation rate of Droplet 2 tends towards positive infinity, while that of Droplet 1 approaches negative infinity, representing a singularity of infinite evaporation and condensation rates at their contact point. Despite this, the total evaporation/burning rate $\dot {Q}_{tot}$ of the system remains finite due to the mutual counteraction of these processes. Note that the singularity resulting from the interaction of two droplets with different temperatures is not sustainable in practical scenarios, and serves primarily as a theoretical construct.

4.4.2. Droplet interaction coefficients

In cases where each of the two burning droplets is independently situated in an infinite space, the respective cumulative evaporation rate for each droplet is determined as (Law 2006; Li et al. 2023b)

(4.23)\begin{equation} \dot{Q}_{iso, i} = 4 {\rm \pi}r_i \times (\rho_g D)_\infty \ln ( 1 + B_{v,i} ) = 4 {\rm \pi}r_i \, \varphi_{s,i}. \end{equation}

Employing the isolated droplet evaporation/burning rate $\dot {Q}_{iso, i}$ given in (4.23) as a reference, the individual interaction coefficients $\eta _i$ for the two droplets, as well as the total interaction coefficient $\eta _{tot}$ for the system, can be determined as

(4.24a,b)\begin{equation} \eta_i \equiv \frac{\dot{Q}_i }{ \dot{Q}_{iso,i} } = \frac{\dot{Q}_i }{4 {\rm \pi}r_i\, \varphi_{s,i}}, \quad \eta_{tot} \equiv \frac{\sum_i \dot{Q}_i }{ \sum_i \dot{Q}_{iso,i} }. \end{equation}

Based on the formulation of $\dot {Q}_i$ as presented in (4.20), the individual interaction coefficients for the two droplets are obtained as

(4.25a)$$\begin{gather} \eta_1 = \frac{2a }{ r_1 } \sum_{n=0}^\infty \frac{\mathrm{e}^{(2n+1)\xi_1}\,( \varphi_{s,2} / \varphi_{s,1} - \mathrm{e}^{(2n+1)\xi_2} ) }{ \mathrm{e}^{(2n+1)\xi_1} - \mathrm{e}^{(2n+1)\xi_2} }, \end{gather}$$

(4.25b)$$\begin{gather}\eta_2 = \frac{2a }{ r_2 } \sum_{n=0}^\infty \frac{(\varphi_{s,1} / \varphi_{s,2})\, \mathrm{e}^{(2n+1)\xi_1} - 1 }{ \mathrm{e}^{(2n+1)\xi_1} - \mathrm{e}^{(2n+1)\xi_2} }. \end{gather}$$

Then, utilizing the expression for $\dot {Q}_{tot}$ in (4.21), the total interaction coefficient for the system can be expressed as

(4.26)\begin{align} \eta_{tot} = \frac{2 a }{ r_1 \varphi_{s,1} + r_2 \varphi_{s,2} } \sum_{n=0}^\infty \frac{\mathrm{e}^{(2n+1)\xi_1}\,( 1 - \mathrm{e}^{(2n+1)\xi_2} ) \varphi_{s,1} + ( \mathrm{e}^{(2n+1)\xi_1} - 1 ) \varphi_{s,2} }{ \mathrm{e}^{(2n+1)\xi_1} - \mathrm{e}^{(2n+1)\xi_2} }> 0. \end{align}

According to (4.25) and (4.26), it is evident that the interaction coefficients $\eta _i$ and $\eta _{tot}$ are affected by the geometric parameters of the system, including the radii of the two droplets $r_i$ and their separation distance $H$. Additionally, these coefficients are influenced by the relative evaporation potential of the droplet surfaces ($\varphi _{s,2} / \varphi _{s,1}$). In particular, in instances where the evaporation potentials of the two droplets are equal, i.e. $\varphi _{s,2} / \varphi _{s,1} = 1$, the interaction coefficients are determined exclusively by the geometric parameters of the system.

Based on (4.25), figure 7 illustrates the variation of $\eta _{i}$ with respect to different droplet radii and distances. It is observed that as $\mathcal {H}$ approaches infinity, $\eta _i$ converges to $1$, signifying minimal interaction. Conversely, at typical droplet distances, $\eta _i$ is less than $1$ and generally decreases with a reduction in $\mathcal {H}$. Notably, when $B_{v,2} = B_{v,1}$, both $\eta _1$ and $\eta _2$ exhibit a gradual decrease as $\mathcal {H}$ reduces, approaching a finite value when $\mathcal {H} \to 1$. In contrast, for $B_{v,2} > B_{v,1}$, a distinct trend is observed: $\eta _1$ decreases rapidly, while $\eta _2$ correspondingly increases sharply as $\mathcal {H}$ descends below a certain threshold. As $\mathcal {H}$ approaches $1$, $\eta _1 \to -\infty$ and $\eta _2 \to +\infty$. This trend is attributed to localized extreme evaporation and condensation rates on the inner sides of the droplets when they are in close proximity, a phenomenon unique to situations where the two droplet surfaces have differing temperatures. A comparison between figures 7(a,b) reveals that as Droplet 1 becomes relatively smaller, its interaction coefficient $\eta _1$ diminishes, indicating a greater suppression of its evaporation rate. Conversely, the interaction coefficient for the larger Droplet 2, $\eta _2$, approaches closer to $1$. Therefore, droplets with higher temperatures and larger radii experience less impact from their neighbouring droplets.

Figure 7. Interaction coefficients for two droplets with varying radii and distances ($r_1 = 100$ $\mathrm {\mu }$m, $B_{v,1} = 4$): (a) $r_2/r_1=1$, (b) $r_2/r_1=2$. The theoretical solution given in (4.25), and the asymptotic solution for $H \to \infty$ given in (4.27) and (4.28), are presented. Here, $\eta = 1$ represents isolated droplet combustion.

4.4.3. Interaction coefficients for two distant droplets

When the two droplets are considerably far apart (i.e. $\mathcal {H} \to \infty$, resulting in $|\xi _i| \to \infty$), based on (4.25) and (4.26), and using the summation formula $\sum _{n=0}^{\infty } \mathrm {e}^{-(2n+1)|\xi |} = 1 / (2 \sinh |\xi |)$, the interaction coefficients for the two droplets can be derived respectively as

(4.27)\begin{equation} \eta_1 |_{\mathcal{H} \to \infty} \sim \frac{a }{ r_1 } \left[ \frac{\varphi_{s,2} / \varphi_{s,1} }{ \sinh (\xi_1 - \xi_2) } - \frac{1 }{ \sinh \xi_1 } \right] \sim 1 - \frac{r_2 }{ H }\,\frac{\varphi_{s,2} }{ \varphi_{s,1} }, \end{equation}

and

(4.28)\begin{equation} \eta_2 |_{\mathcal{H} \to \infty} \sim \frac{a }{ r_2 } \left[ \frac{\varphi_{s,1} / \varphi_{s,2} }{ \sinh (\xi_1 - \xi_2) } + \frac{1 }{ \sinh \xi_2 } \right] \sim 1 - \frac{r_1 }{ H }\,\frac{\varphi_{s,1} }{ \varphi_{s,2} }. \end{equation}

Additionally, the total interaction coefficient can be obtained as

(4.29)\begin{equation} \eta_{tot} |_{\mathcal{H} \to \infty} \sim 1 - \frac{r_1 r_2 }{ H }\, \frac{\varphi_{s,1} + \varphi_{s,2} }{ r_1 \varphi_{s,1} + r_2 \varphi_{s,2} }. \end{equation}

Based on (4.27) and (4.28), it is evident that the interaction coefficient for a droplet reduces as the size and relative evaporation potential of the other droplet increase. As the separation between the two droplets narrows (i.e. $\xi _2 - \xi _1$ decreases), the mutual interference between them intensifies. Notably, when the droplets are substantially separated, their interaction resembles that of point sources. Therefore, the results for $\eta _i$ and $\eta _{tot}$ in the limit of $\mathcal {H} \to \infty$, as determined from (4.27), (4.28) and (4.29), can also be deduced using the point-source method, a technique explored in various studies (Marberry et al. 1984; Annamalai & Ryan 1992; Cossali & Tonini 2020). However, it is worth mentioning that the existing literature has not explored the impact of differing $\varphi _{s,i}$. Figure 7 compares the asymptotic solutions of $\eta _i$ for two distant droplets with the exact solution from (4.25). It is observed that except in cases involving essentially small $\mathcal {H}$, the asymptotic solution aligns closely with the exact solution. This agreement is particularly apparent for droplets with relatively smaller sizes and Spalding numbers.

Upon comparing the expressions for $\eta _i$ of the two distant droplets, as specified in (4.27) and (4.28), the following relationship is derived:

(4.30)\begin{equation} r_i \varphi_{s,i}^2 (1 - \eta_i) = \frac{r_1 r_2 }{ H } \, \varphi_{s,1} \varphi_{s,2}. \end{equation}

The left-hand side of (4.30) delineates the characteristics of an individual droplet, including its radius, surface evaporation potential and interaction coefficient, while the right-hand side reflects the overall characteristics of the dual-droplet system. This relationship suggests that for a system with distant droplet pairs, the product $r_i \varphi _{s,i}^2 (1 - \eta _i)$ can be considered approximately constant. Consequently, a larger radius $r_i$ and a higher surface evaporation potential $\varphi _{s,i}$ lead to an increased $\eta _i$, simulating a scenario where droplet evaporation and combustion are minimally affected by neighbouring droplets. This finding, in line with our previous discussions, implies that in sparse sprays within a stationary environment, larger and warmer droplets experience less influence from nearby droplets compared to smaller and cooler droplets, which are more susceptible to such influences.

Based on the asymptotic solutions of $\eta _i$ for $\mathcal {H} \to \infty$, as outlined in (4.27) and (4.28), we can deduce that in a multi-droplet burning system with sparsely distributed droplets, their interactions can be linearly superimposed, similar to treating them as point sources. In this context, $\eta _i$ and $\eta _{tot}$ can be expressed as

(4.31a,b) \begin{equation} \eta_i \sim 1 - \sum_{\substack{j = 1 \\ j \neq i}}^{j=N} \frac{r_j }{ H_{ij} }\, \frac{\varphi_{s,j} }{ \varphi_{s,i} }, \quad \eta_{tot} \sim 1 - \frac{1}{ \displaystyle\sum_{i=1}^N ( r_i \varphi_{s,i}) } \times \sum_{i=1}^N \begin{pmatrix}\displaystyle r_i \sum_{\substack{j = 1 \\ j \neq i}}^{j=N} \frac{r_j \varphi_{s,j} }{ H_{ij} } \end{pmatrix}, \end{equation}

where $N$ represents the total number of droplets in the system, $r_j$ is the radius of the $j$th droplet, and $H_{ij}$ is the distance between the $i$th and $j$th droplets. The applicability of (4.31) is $( r_j / H_{ij} ) ( \varphi _{s,j} / \varphi _{s,i} ) \ll 1$ ($i, j = 1, 2, \ldots N$; $j \neq i$). When all droplets in the system maintain the same temperature, and thus the same evaporation potential $\varphi _{s,i}$, the condition simplifies to requiring that the inter-droplet distance be substantially larger than the size of any individual droplet. Conversely, in scenarios where $\varphi _{s,i}$ varies, the condition necessitates that droplets with a higher evaporation potential be positioned farther away from other droplets.

4.4.4. Interaction coefficients for two closely approaching droplets

As the inter-droplet distance reduces, the convergence rate for the infinite series solution of $\eta _i$, detailed in (4.25), becomes slower. In such instances, an integral approach can be employed as a practical approximation for the infinite series. Specifically, for Droplet 1, the value of $\eta _1$ can be approximated using the integral

(4.32a,b)\begin{equation} \eta_1 \approx \frac{2a }{ r_1 } \int_{0}^{\infty} F (n) \,\mathrm{d} n, \quad F (n) \equiv \frac{\mathrm{e}^{(2n+1)\xi_1}\,( \varphi_{s,2} / \varphi_{s,1} - \mathrm{e}^{(2n+1)\xi_2} ) }{ \mathrm{e}^{(2n+1)\xi_1} - \mathrm{e}^{(2n+1)\xi_2}}. \end{equation}

Considering the complexity of the integral in (4.32) under general circumstances, our discussion will focus on the specific case where the two droplets possess equal radii. By employing the defined expression for $F(n)$, the integral in (4.32) can be derived as

(4.33)\begin{equation} \int_{0}^{\infty} F (n) \,\mathrm{d} n = \left( \frac{\varphi_{s,2} }{ \varphi_{s,1} } - 1 \right) \left( \frac{\ln (\mathrm{e}^{2\xi_2} - 1) }{ 4 \xi_2 } - \frac{1 }{ 2 } \right) + \left( \frac{\ln (1 + \mathrm{e}^{\xi_2}) }{2 \xi_2} - \frac{1 }{ 2 } \right). \end{equation}

As two droplets with equal radii approach each other, we can deduce that $a \to r_1 \sqrt {2(\mathcal {H} - 1)}$ and $\mathrm {e}^{\xi _2} \to 1 + \sqrt {2(\mathcal {H} - 1)}$. Utilizing these relationships, the limit of $\eta _1$ as $\mathcal {H} \to 1$ can be determined as

(4.34)\begin{equation} \lim_{\mathcal{H} \to 1} \eta_1 |_{r_1 = r_2} = \ln 2 + \frac{1 }{ 2 } \left( \frac{\varphi_{s,2} }{ \varphi_{s,1} } - 1 \right) \times \lim_{\mathcal{H} \to 1} \ln ( \sqrt{2 (\mathcal{H} - 1)} ). \end{equation}

Similarly, for Droplet 2, the limit of $\eta _2$ as $\mathcal {H} \to 1$ reads

(4.35)\begin{equation} \lim_{\mathcal{H} \to 1} \eta_2|_{r_1 = r_2} = \ln 2 + \frac{1 }{ 2 } \left( \frac{\varphi_{s,1} }{ \varphi_{s,2} } - 1 \right) \times \lim_{\mathcal{H} \to 1} \ln ( \sqrt{2 (\mathcal{H} - 1)} ). \end{equation}

Additionally, the total interaction coefficient for the dual-droplet system can be calculated as

(4.36)\begin{equation} \lim_{\mathcal{H} \to 1} \eta_{tot} |_{r_1 = r_2} = \frac{\varphi_{s,1} \lim_{\mathcal{H} \to 1} \eta_1 |_{r_1 = r_2} + \varphi_{s,2} \lim_{\mathcal{H} \to 1} \eta_2 |_{r_1 = r_2} }{ \varphi_{s,1} + \varphi_{s,2} } = \ln 2. \end{equation}

In cases where $\varphi _{s,1}$ and $\varphi _{s,2}$ differ, (4.34) and (4.35) suggest that $\lim _{\mathcal {H} \to 1} \eta _i = O ( \ln (\mathcal {H} - 1) )$. As the droplets approach each other closely, the interaction coefficient for the droplet with the larger evaporation potential approaches infinity, while for the droplet with the smaller evaporation potential, it tends towards negative infinity. However, in the particular case where $\varphi _{s,1} = \varphi _{s,2}$, $\eta _i (\mathcal {H} \to 1)$ for two equal-sized droplets stabilizes at $\ln 2 \approx 0.693$. This finding aligns with the numerical solutions shown in figure 7(a), and is consistent with previous research (Annamalai & Ryan 1992; Sirignano 2014). Interestingly, as indicated in (4.36), the total interaction coefficient $\eta _{tot}$ for the system consistently equals $\ln 2$ in close proximity situations, regardless of the relative evaporation potentials. This outcome arises from the counterbalancing of infinite mass transfer at the contact point between the two droplets.

4.5. Flame position and combustion mode transition

The position and morphology of the flame surface are crucial in the analysis of interactive combustion within multi-droplet systems. These aspects play a key role in discerning whether the system exhibits a separated combustion mode or a group combustion mode (Sirignano 2014). Under the assumption of infinitely fast reactions, the mass fractions of fuel and oxygen on the flame sheet are both reduced to zero (i.e. $Y_F = Y_O = 0$). Therefore, in accordance with the definition of $\beta _{FO}$ in (2.2), it is deduced that $\beta _{FO,f} = \beta _{FO} (\varphi _f) = 0$. Utilizing the solution for $\beta _{FO}$ as presented in (4.11), the value of $\varphi _f$ at the flame surface can be established as

(4.37)\begin{equation} \varphi_f = (\rho_g D)_\infty \ln (1 + \phi Y_{O,\infty}). \end{equation}

Equation (4.37) implies that the flame surface aligns with an isosurface of $\varphi$, a finding consistent with earlier research on multi-droplet systems with identical temperatures (Labowsky 1980; Imaoka & Sirignano 2005). As per (4.37), $\varphi _f$ correlates with the oxygen content $Y_{O,\infty }$ in the far field, and remains unaffected by the surface temperature $T_{s,i}$ or the fuel vapour mass fraction $Y_{F,s,i}$ at the droplet surfaces. We introduce a dimensionless quantity $\mathcal {R}_{f,i}$ to quantify the relative distance between the flame surface and individual droplets:

(4.38)\begin{equation} \mathcal{R}_{f,i} \equiv \left( \frac{\varphi_f }{ \varphi_{s,i} } \right)^{{-}1} = \frac{\ln ( 1 + B_{v,i} ) }{ \ln ( 1 + \phi Y_{O,\infty} ) }. \end{equation}

Note that $\mathcal {R}_{f,i}$ represents the relative position of the flame surface within a one-dimensional $\varphi$-based coordinate system, as opposed to its physical location in three-dimensional Cartesian space. The formulation of $\mathcal {R}_{f,i}$ in (4.38) resembles the radius ratio between the flame and droplet surface during the burning of a single, isolated droplet (Sirignano 2014). According to (4.38), the distance between the flame surface and a droplet within the system increases in response to either a higher evaporation potential of the droplet or a reduced oxygen concentration in the surrounding environment. This observation is believed to apply not only to isolated droplets but also to systems encompassing two or more droplets.

In scenarios where two droplets are significantly distant from each other, each droplet is encircled by a nearly spherical flame, resembling combustion individually in an infinite space, characterizing the system as being in a separated combustion mode. As the droplets approach each other more closely, the shape of their flames starts to change. When the separation becomes sufficiently small, the two distinct flames coalesce into a single flame surface that encompasses both droplets, signifying a transition to group combustion mode. This transition of combustion mode alters the topological properties of the flame surface. Additionally, the operating combustion mode serves as an indicator of the interaction intensity between the burning droplets.

Given the axial symmetry of the problem, the initial point of contact between the two separate flames surrounding the droplets, occurring as $\mathcal {H}$ decreases, is situated along the axis between the surfaces of the two droplets. In cases where the droplets have identical radii and surface potentials, this contact point is precisely at their midpoint $(\xi = 0, \theta = {\rm \pi})$. For more general scenarios, as illustrated in figure 4(b), the criterion for the transition of combustion mode can be established as

(4.39)\begin{equation} \min_{\xi \in [\xi_1, \, \xi_2]} \varphi (\theta = {\rm \pi}) = \varphi_f. \end{equation}

By applying the asymptotic solution for $\varphi (\theta = {\rm \pi})$ from (4.16) into (4.39), the critical droplet distance $\mathcal {H}_{cri}$ that signifies the transition between combustion modes can be determined. Comparing the solution for $\varphi _{s,i}$ (see (4.9)) with the expression for $\varphi _f$ (see (4.37)), it is evident that $\varphi _{s,i}$ is consistently greater $\varphi _f$ across all valid values of $B_{v,i}$ and $\phi Y_{O,\infty }$. This indicates that the minimum value of $\varphi (\theta = {\rm \pi})$ is situated between the surfaces of the two droplets, rather than being directly on either surface. Considering the concave profile of the $\varphi (\theta = {\rm \pi})$ curve, as shown in figure 4, it is demonstrated that during the critical transition of combustion mode, evaporation occurs at the two closest points between the droplets, without any condensation. Accordingly, evaporation is the exclusive process taking place on all surface areas of both droplets. Consequently, the previously discussed phenomenon of partial surface condensation on a droplet, as illustrated in figure 5, is feasible exclusively within group combustion mode.

Based on (4.39), figure 8(a) displays the critical inter-droplet distance $\mathcal {H}_{cri}$ for the transition between combustion modes in a dual-droplet system with varying droplet radii. When $\mathcal {H}$ is greater than $\mathcal {H}_{cri}$, the dual-droplet system displays a separated combustion mode, as depicted in figure 8(b). Conversely, if $\mathcal {H}$ falls below $\mathcal {H}_{cri}$, then the system transitions to a group combustion mode, as shown in figure 8(d).

Figure 8. (a) Critical inter-droplet distance for combustion mode transition ($r_1 = 100$ $\mathrm {\mu }$m). The theoretical solution is calculated based on (4.12) and (4.39), while the asymptotic solution for $r_2 / r_1 \to \infty$ is obtained from (4.40). (b–d) The three typical combustion modes indicated in (a): separated combustion, critical combustion and group combustion modes.

Figure 8(a) shows that for $r_2 / r_1 = 1$, indicating droplets of equal size, the evaporation potentials of both droplets contribute equally to determining $\mathcal {H}_{cri}$. Increasing $B_{v,i}$ results in a greater $\mathcal {H}_{cri}$. Generally, $\mathcal {H}_{cri}$ decreases monotonically as the ratio $r_2 / r_1$ increases. However, in assessing the dimensional counterpart of $\mathcal {H}_{cri}$, represented by $H_{cri} = \mathcal {H}_{cri} (r_1 + r_2)$, it becomes evident that $H_{cri}$ actually increases with an increasing ratio $r_2 / r_1$. It suggests that all else being constant, enlarging the size of either droplet facilitates the interaction between their flames.

In situations with a marked size difference between the two droplets, the contact point of the two distinct flames, during the critical transition of combustion mode, is notably closer to the smaller droplet (see figure 12). In light of this, by applying the asymptotic solution for $\varphi$ from (4.17), the explicit expression for $\mathcal {H}_{cri}$ can be obtained as

(4.40a)$$\begin{gather} \mathcal{H}_{cri} |_{{r_2 }/{ r_1 } \to 0} \sim \left( 1 - \sqrt{1 - \frac{\varphi_f }{ \varphi_{s,2} } } \right)^{{-}1} \frac{\varphi_{s,1} }{ 2 \varphi_{s,2} }, \end{gather}$$

(4.40b)$$\begin{gather}\mathcal{H}_{cri} |_{{r_2 }/{ r_1 } \to \infty} \sim \left( 1 - \sqrt{1 - \frac{\varphi_f }{ \varphi_{s,1} } } \right)^{{-}1} \frac{\varphi_{s,2} }{ 2 \varphi_{s,1}}. \end{gather}$$

According to (4.40), when there is a substantial size difference between the two droplets, resembling a scenario of a droplet adjacent to a large liquid plane, the critical distance $\mathcal {H}_{cri}$ is not determined by the size of the droplets. Instead, it is determined exclusively by $\varphi _f$, which is influenced primarily by $\phi Y_{O,\infty }$, in addition to $\varphi _{s,i}$ of the two droplets. As shown in figure 8(a), when the ratio $r_2 / r_1$ approaches infinity, $\mathcal {H}_{cri}$ approximately converges to the limit specified in (4.40). It is found that this limit is more sensitive to variations in the surface potential of Droplet 2 compared to Droplet 1, indicating that the evaporation potential of the larger droplet exerts a greater influence on the transition between combustion modes.

It should be noted that the results presented in figure 8 are based on the condition where $Y_{O,\infty } = 1$, indicative of an oxygen-rich environment, which brings the flame closer to the droplet surface and results in a smaller $\mathcal {H}_{cri}$. For instance, with $B_{v,1} = B_{v,2} = 1$, $\mathcal {H}_{cri}$ is approximately $4.3$ when $Y_{O,\infty } = 1$, in contrast to approximately $21.8$ in an environment where $Y_{O,\infty } = 0.2$. To provide a clear illustration of the distribution of various physical quantities near the droplet and flame surfaces, subsequent results are obtained based on $Y_{O,\infty } = 1$.

4.6. Solutions for scalar fields and flame temperature distribution

In our above analyses, we have successfully derived both the exact solution, as presented in (4.1), and the asymptotic solution, as per (4.12), for the potential function $\varphi$. Furthermore, as indicated in (4.11), the coupling function $\beta _{FO}$ exhibits a one-to-one correspondence with $\varphi$. This correlation has led to the conclusion that the flame surface aligns with an isosurface of $\varphi$, as outlined in (4.37).

The distribution of fuel vapour and oxygen can be deduced by referencing the definition of $\beta _{FO}$, as presented in (2.2), in conjunction with the flame-sheet assumption. However, to determine the temperature field, it is imperative to find solutions for the additional two coupling functions, $\beta _{TF}$ and $\beta _{TO}$, also defined in (2.2). In cases where the temperatures of the two droplets differ, leading to $\varphi _{s,1} \neq \varphi _{s,2}$, these two coupling functions do not maintain a one-to-one correspondence with $\varphi$ as $\beta _{FO}$ does. Nonetheless, by employing a derivation approach similar to that used for the asymptotic solution of $\varphi$ in Appendix A, we have successfully derived asymptotic solutions for $\beta _{TF}$ and $\beta _{TO}$ as

(4.41)\begin{equation} \beta_k \sim \beta_{k,\infty} - \left[ \frac{\beta_{k,\infty} - \beta_{k,s,1} }{ \mathcal{L} (\varphi_{s,1}) - 1 }\,( 1 - \mathcal{B} (\xi) ) + \frac{\beta_{k,\infty} - \beta_{k,s,2} }{ \mathcal{L} (\varphi_{s,2}) - 1 }\,\mathcal{B} (\xi) \right] ( \mathcal{L} (\varphi) - 1 ), \end{equation}

where $\beta _{k,s,i}$ and $\beta _{k,\infty }$ are specified in the boundary conditions as given in (4.2). Equation (4.41) indicates that $\beta _{TF}$ and $\beta _{TO}$ can be approximated as functions reliant on both $\varphi$ and $\mathcal {B} (\xi )$, which is a bending function ranging from $0$ to $1$ (see (4.13)). Incorporating the boundary conditions for $\beta _{TF}$ and $\beta _{TO}$ from (4.2) into (4.41), the solutions for these two coupling functions can be derived as

(4.42)\begin{align} \beta_{TF} &= \frac{c_{p,g} T_\infty }{ q_c } - \left[ \frac{1 - \mathcal{B} (\xi) }{ \mathcal{L} (\varphi_{s,1}) - 1 } \left( \frac{c_{p,g} (T_\infty - T_{s,1}) }{ q_c } - \frac{Y_{F,s,1} }{ W_F \nu_F } \right) \right. \nonumber\\ &\quad \left.{}+ \frac{\mathcal{B} (\xi) }{ \mathcal{L} (\varphi_{s,2}) - 1 } \left( \frac{c_{p,g} (T_\infty - T_{s,2}) }{ q_c } - \frac{Y_{F,s,2} }{ W_F \nu_F } \right) \right]( \mathcal{L} (\varphi) - 1 ) \end{align}

and

(4.43)\begin{align} \beta_{TO} &= \left( \frac{c_{p,g} T_\infty }{ q_c } + \frac{Y_{O,\infty} }{ W_O \nu_O } \right) - \left[ \frac{1 - \mathcal{B} (\xi) }{ \mathcal{L} (\varphi_{s,1}) - 1 } \left( \frac{c_{p,g} (T_\infty - T_{s,1}) }{ q_c } + \frac{Y_{O,\infty} }{ W_O \nu_O } \right) \right. \nonumber\\ &\quad \left.{}+ \frac{\mathcal{B} (\xi) }{ \mathcal{L} (\varphi_{s,2}) - 1 } \left( \frac{c_{p,g} (T_\infty - T_{s,2}) }{ q_c } + \frac{Y_{O,\infty} }{ W_O \nu_O } \right) \right]( \mathcal{L} (\varphi) - 1 ). \end{align}

Utilizing the solutions for the three coupling functions, namely $\beta _{FO}$, $\beta _{TF}$ and $\beta _{TO}$, as presented respectively in (4.11), (4.42) and (4.43), the scalar fields including $T_g$, $Y_F$ and $Y_O$, for regions both inside and outside the flame surface, can be obtained. These computations are facilitated by integrating the definitions of these coupling functions provided in (2.2). In scenarios of infinitely fast chemical reactions, oxygen does not exist between the droplet and flame surfaces, and fuel vapour is absent beyond the flame surfaces. Consequently, in the region enclosed by the diffusion flame (where $\varphi > \varphi _f$), we have

(4.44a–c)\begin{equation} Y_F = (W_F \nu_F) \times \beta_{FO}, \quad Y_O = 0, \quad T_g = \frac{q_c }{ c_{p,g} } \times \beta_{TO}, \end{equation}

while for the region outside the flame surface (where $\varphi < \varphi _f$), we obtain

(4.45a–c)\begin{equation} Y_F = 0, \quad Y_O ={-} (W_O \nu_O) \times \beta_{FO}, \quad T_g = \frac{q_c }{ c_{p,g} } \times \beta_{TF}. \end{equation}

Upon substituting $\varphi = \varphi _f$ into the solutions for $\beta _{TF}$ and $\beta _{TO}$ given in (4.42) and (4.43), and considering the absence of fuel and oxygen on the flame surface, the temperature distribution along the flame surface can be determined as

(4.46)\begin{align} T_f (\xi)& = T_\infty + \phi Y_{O,\infty} \left[ \frac{1 - \mathcal{B} (\xi) }{ \phi Y_{O,\infty} + Y_{F,s,1} } \left( (T_{s,1} - T_\infty) + \frac{q_c }{ c_{p,g} }\,\frac{Y_{F,s,1} }{ W_F \nu_F } \right) \right. \nonumber\\ &\quad \left.{}+ \frac{\mathcal{B} (\xi) }{ \phi Y_{O,\infty} + Y_{F,s,2} } \left( (T_{s,2} - T_\infty) + \frac{q_c }{ c_{p,g} }\,\frac{Y_{F,s,2} }{ W_F \nu_F } \right)\right]. \end{align}

According to (4.46), contrasting with the results from existing studies based on equal surface temperatures for the droplets (Labowsky 1980; Umemura et al. 1981b; Imaoka & Sirignano 2005), it is evident that the flame temperature $T_f$ is not uniform when the droplet temperatures are unequal. Given that the chemical energy released during combustion substantially outweighs the internal energy of the evaporated fuel vapour, $T_f$ is influenced primarily by $Y_{F,s,i}$ rather than $T_{s,i}$. In accordance with the Clausius–Clapeyron relation (see (2.14)), which establishes a direct correlation between $Y_{F,s,i}$ and $T_{s,i}$ in thermodynamic equilibrium, elevating the surface temperature of a droplet notably increases its adjacent flame temperature, whereas lowering the droplet's temperature reduces the flame temperature on that side.

Based on (4.46), figure 9 presents the theoretical flame temperature distribution. For the case with $B_{v,1} = 1$ and $B_{v,2} = 2$, the respective surface temperatures of the two droplets are approximately $315\ {\rm K}$ and $333\ {\rm K}$. Despite the modest difference of $18\ {\rm K}$ in droplet surface temperatures, a substantial temperature variation becomes apparent on the flame surface, with the maximum discrepancy exceeding $200\ {\rm K}$. This substantial variation suggests that differences in surface fuel vapour concentrations, stemming from the distinct droplet surface temperatures, are the primary factors contributing to the observed uneven temperature distribution on the flame surface.

Figure 9. Theoretical distribution of flame temperature outside the two droplets at varying droplet distances and sizes ($r_1 = 100$ $\mathrm {\mu }$m; $B_{v,1} = 1$, corresponding to $T_{s,1} \approx 315\ {\rm K}$ and $Y_{F,s,1} \approx 0.35$; $B_{v,2} = 2$, corresponding to $T_{s,2} \approx 333\ {\rm K}$ and $Y_{F,s,2} \approx 0.57$). The theoretical solution is calculated using (4.12), (4.37) and (4.46).

In scenarios where $B_{v,2} > B_{v,1}$, as depicted in figure 9, the lowest temperature on the flame surface resides at the outermost side of Droplet 1, denoted as $T_{f,out,1}$, while the highest temperature is found at the outermost side of Droplet 2, identified as $T_{f,out,2}$. When $\mathcal {H}$ is smaller than $\mathcal {H}_{cri}$, a unified flame surface encompasses both droplets, resulting in a gradual transition of flame temperature $T_f$ from $T_{f,out,1}$ to $T_{f,out,2}$ along the flame surface. In contrast, when $\mathcal {H}$ surpasses $\mathcal {H}_{cri}$, the flame surface surrounding the droplets separates, leading to two distinct temperature discontinuities. These discontinuities are located at the nearest points of the two separate flames, with the temperatures at these points identified as $T_{f,in,1}$ and $T_{f,in,2}$, respectively. The $T_{f,in,i}$ value for each droplet rapidly approaches its corresponding $T_{f,out,i}$ value for increasing $\mathcal {H}$, indicating a more uniform temperature distribution across each flame surface. According to figure 9, as $\mathcal {H}$ approaches infinity, both $T_{f,in,i}$ and $T_{f,out,i}$ converge to $T_{f,iso,i}$, which represents the flame temperature of Droplet $i$ when it burns in isolation.

Figure 9 provides additional insight into the influence of the relative sizes of the two droplets on the range of flame temperatures. Specifically, within the group combustion mode, an increase in the size of Droplet 2 results in a higher maximum temperature $T_{f,out,2}$ on the flame surface. This increase in the maximum temperature also simultaneously raises the minimum temperature of the flame surface, $T_{f,out,1}$. However, when the system transitions into the separated combustion mode, the sensitivity of the variation in flame temperatures to the relative sizes of the droplets diminishes.

Figure 10 displays the contours of various physical quantities during the QS combustion of two interacting droplets, identical in both temperature and radius. A comparison with numerical solutions reveals that the asymptotic solutions, derived from (4.11), (4.12), (4.42) and (4.43), effectively capture variations in the flow, composition and temperature fields at different inter-droplet distances. It is seen that with droplets at equal temperatures, the flame temperature distribution remains uniform. When the droplets are in close proximity, as shown in figure 10(a) with $\mathcal {H} = 1.5$, they are enveloped by an ellipsoidal-shaped flame. As the inter-droplet distance increases, as depicted in figure 10(b) with $\mathcal {H} = 3$, the ellipsoidal flame stretches, and the Stefan flow intensifies around both droplets, particularly on their adjacent sides. Further increase of $\mathcal {H}$, as illustrated in figure 10(c) with $\mathcal {H} = 5$, leads to the division of the enveloping flame into two distinct flames, each surrounding one droplet. These two separate flames initially exhibit an inward-pointing egg shape around each droplet. As $\mathcal {H} \to \infty$, the two flames gradually transition into spherical shapes, with the evaporation and combustion dynamics of each droplet mirroring that of an isolated droplet.

Figure 10. Distributions of various physical quantities during the combustion of two droplets of equal size and temperature at varying distances ($r_1 = r_2 = 100$ $\mathrm {\mu }$m; $B_{v,1} = B_{v,2} = 1$, corresponding to $T_{s,i} \approx 315\ {\rm K}$ and $Y_{F,s,i} \approx 0.35$). Asymptotic solutions, calculated primarily using (4.11), (4.12), (4.42) and (4.43), are compared with numerical solutions.

In contrast to the special case of equal-temperature droplets shown in figure 10, figure 11 presents a more general scenario involving droplets with disparate temperatures, as indicated by their differing Spalding numbers. It is evident that our theoretical solutions maintain a close alignment with the numerical results. In figure 11, the Spalding number for Droplet 2, $B_{v,2}$, is increased from $1$ to $2$, in contrast to the case shown in figure 10. As a result, Droplet 2 exhibits a higher temperature and fuel vapour concentration on its surface, which is visible in figures 11(a iii), 11(b iii) and 11(c iii). This leads to a faster evaporation rate for Droplet 2, causing its adjacent flame surface to be located further from the droplet surface. A comparison between figures 10(c) and 11(c) reveals that a higher Spalding number for either of the droplets increases the probability of group combustion and reinforces the interaction between the droplets. In contrast to the uniform flame temperature distribution observed in figure 10, scenarios involving differing droplet temperatures, as shown in figure 11, display a non-uniform temperature distribution across the flame surface. This non-uniformity is characterized by notably lower flame temperatures near Droplet 1 in comparison to the flame temperatures observed on the side of Droplet 2.

Figure 11. Distributions of various physical quantities during the combustion of two droplets of equal size but unequal temperatures at varying distances ($r_1 = r_2 = 100$ $\mathrm {\mu }$m; $B_{v,1} = 1$, corresponding to $T_{s,1} \approx 315\ {\rm K}$ and $Y_{F,s,1} \approx 0.35$; $B_{v,2} = 2$, corresponding to $T_{s,2} \approx 333\ {\rm K}$ and $Y_{F,s,2} \approx 0.57$). Asymptotic solutions, calculated primarily using (4.11), (4.12), (4.42) and (4.43), are compared with numerical solutions.

While this study focuses primarily on the QS combustion of two interacting droplets with differing temperatures by assuming infinitely fast reactions, the insights gained are invaluable for understanding important phenomena associated with finite-rate reactions, such as flame spreading, as well as local flame extinction and ignition, in multi-droplet systems (Brzustowski et al. 1981; Mikami et al. 2006, 2023; Yoshida et al. 2019). In scenarios involving finite-rate reactions, the high sensitivity of combustion reactions to temperature variations is crucial. A reduction in temperature at any segment of the flame surface can potentially cause local flame extinction (Farouk & Dryer 2023). In the context of the present dual-droplet system, such local extinction is likely to occur near the cooler droplet, resulting in the leakage of unburned fuel vapour. Conversely, raising the temperature of either one or both droplets in the pair may lead to local ignition and flame spreading.

Indeed, in real-world scenarios, droplet temperatures are influenced by their surroundings, including heat transfer through thermal conduction and radiation from the flames. In this context, when a burning droplet approaches an unignited, cold droplet, two primary outcomes are possible. First, the cold droplet can be heated up substantially, eventually leading to its ignition, and potentially resulting in group combustion involving both droplets. Alternatively, the burning droplet could cool down considerably, with a significant reduction in its flame temperature, which may cause flame extinction. For simplicity, this paper focuses exclusively on scenarios with fixed droplet temperatures, enabling valuable discussions regarding the one-way influence of droplet temperatures on flame temperature distributions in dual-droplet systems. To comprehensively address the intricate aspects of transient ignition, extinction and flame spreading in such systems, further research and advanced modelling efforts are essential.

Note that figures 10 and 11 illustrate cases of droplets with equal radii; figure 12 broadens the scope to include a more general situation with droplets that have unequal radii. Specifically, Droplet 2 has a radius double that of Droplet 1. Figures 12(a–d) display the situation where $B_{v,2} / B_{v,1} = 0.5$, while figures 12(e–h) depict the scenario when the ratio is $2$. These figures reveal that a larger droplet radius enhances inter-droplet interaction, enlarges the flame, and delays flame separation as $\mathcal {H}$ increases. In figures 12(a–d), despite Droplet 2 being larger, its lower evaporation potential results in a flame size comparable to that of Droplet 1. Conversely, figures 12(e–h) show Droplet 2 as having both a larger radius and a greater evaporation potential, leading to a notably larger flame surface. In this instance, the evaporation and combustion process of the dual-droplet system is dominated primarily by the characteristics of Droplet 2. The flame shape and size around Droplet 2 remain relatively unchanged, irrespective of the inter-droplet distance or the combustion mode. Droplet 2 exhibits a comparatively uniform evaporation rate, characterized by spherically symmetric streamlines in the contours. In contrast, Droplet 1 shows a more variable evaporation rate with disturbed streamlines, particularly at small droplet distances (see figure 12e). Notably, regions exhibiting higher flame temperatures are consistently located closer to the droplet with a higher evaporation potential. This observation holds regardless of the relative sizes of the two droplets, in accordance with the theoretical solution presented in (4.46).

Figure 12. Theoretical solutions for the flow and temperature fields during the combustion of two droplets with dissimilar radii and temperatures at various inter-droplet distances ($r_1 = 100$ $\mathrm {\mu }$m, $r_2 = 200$ $\mathrm {\mu }$m). The theoretical solutions are calculated using (4.11), (4.12), (4.42) and (4.43).

It should be noted that more intricate evaporation/condensation and combustion phenomena are anticipated with multi-component droplets or those exhibiting non-uniform temperature distributions (Sirignano & Wu 2008; Cira, Benusiglio & Prakash 2015; Tonini & Cossali 2016; Diddens et al. 2017; Wang et al. 2022). However, these complexities are outside the scope of the current study.

5. Conclusions

This study presents a comprehensive theoretical analysis of the interaction between two quasi-steady burning droplets at different temperatures. The mass-flux-potential model and the flame-sheet assumption are employed in conjunction with bispherical coordinates. Theoretical solutions are derived for the gas-phase flow and scalar fields, droplet evaporation/burning rates and interaction coefficients, as well as the position and morphology of flames. The transition criterion for combustion modes is also obtained, specifically identifying the transition between group combustion and separated combustion modes. The comparison with extensive numerical simulations reveals a good alignment between our analytical findings and the numerical outcomes under various parameter conditions, such as the radii, distances and temperatures of the two droplets. It is shown that a reduction in the distance between droplets intensifies their interaction, resulting in an increasingly uneven distribution of evaporation rates on their surfaces. Notably, droplets with smaller evaporation potential and radius are more affected by their neighbouring droplets. In scenarios where a cooler droplet is in close proximity to a warmer one, part of the cooler droplet's surface may experience condensation instead of evaporation. This phenomenon occurs exclusively when two droplets with different temperatures engage in group combustion mode or pure evaporation state. For the current situation with unequal droplet temperatures, the distribution of temperature across the flame surface becomes non-uniform, contrary to the case with equal droplet temperatures. Increasing the temperature of one droplet significantly raises the temperature of the flame adjacent to it. This study effectively broadens existing research to encompass the general scenario of interactive burning droplets with differing temperatures. The insights provided are crucial for understanding local extinction, ignition and flame spreading in multi-droplet systems.