# Coupling Volume-Of-Fluid and Chemical Kinetics for Direct Numerical Simulations of Droplet Combustion

Edoardo Cipriano, Riccardo Caraccio, Alessio Frassoldati, Tiziano Faravelli, Alberto Cuoci

► To cite this version:

Edoardo Cipriano, Riccardo Caraccio, Alessio Frassoldati, Tiziano Faravelli, Alberto Cuoci. Coupling Volume-Of-Fluid and Chemical Kinetics for Direct Numerical Simulations of Droplet Combustion. International Journal of Heat and Mass Transfer, 2025, 253, pp.127586. 10.1016/j.ijheatmasstransfer.2025.127586. hal-04962085v2

HAL Id: hal-04962085

https://hal.science/hal-04962085v2

Submitted on 31 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Coupling Volume-Of-Fluid and Chemical Kinetics for
Direct Numerical Simulations of Droplet Combustion

Edoardo Cipriano<sup>a,*</sup>, Riccardo Caraccio<sup>a</sup>, Alessio Frassoldati<sup>a</sup>, Tiziano
Faravelli<sup>a</sup>, Alberto Cuoci<sup>a</sup>

<sup>a</sup>CRECK Modeling Lab, Department of Chemistry, Materials, and Chemical Engineering
"G. Natta", Politecnico di Milano, Piazza Leonardo da Vinci, 32, Milano, 20133, Italy

Abstract

Since the development of the d² law, in 1952, numerical modeling of droplet combustion has predominantly relied on the assumption of spherical symmetry. In this work, we relax this hypothesis by introducing a strategy to combine interface-resolved phase change models with the direct solution of the combustion chemistry. The multiphase system is described using the geometric Volume-Of-Fluid approach, previously extended to simulate the evaporation of mixtures with variable thermodynamic and transport properties. Gas-phase chemical reactions are incorporated using an operator splitting technique, addressing the stiffness and non-linearity of the reactive step. The model is validated by comparing its predictions with experimental data and with a spherically-symmetric model, proving convergence and capturing the burning rates accurately. The impact of natural convective fluxes at different ambient pressures is quantified by suspending the droplet on a solid fiber in normal gravity conditions. Lastly, the multicomponent formulation is employed to predict droplet extinction resulting from water absorption. The model's implementation and simulation setups are publicly available on the Basilisk code website.

Keywords: vof, droplet, evaporation, combustion, microgravity, buoyancy

<footnote>*Corresponding author.</footnote>

<footnote>Email address: edoardo.cipriano@polimi.it (Edoardo Cipriano)</footnote>
# 1. Introduction

Spray combustion of liquid fuels is a complex and multiphysics process involving a multitude of temporal and spatial scales. To simplify the analysis, several experimental studies have focused on the behavior of single, isolated droplets. This approach allows for a more refined understanding of the evaporation and combustion of fuel droplets, by eliminating the complexities arising from droplet interactions within the spray. Experimentally, droplets with diameters ranging between 0.1 to 2 mm have been studied, typically suspended on a solid fiber, both in microgravity and in normal gravity conditions. Nomura et al. [1] studied the evaporation of n-heptane droplets at high pressure and temperature conditions; Kumagai and Isoda [2] developed an accurate methodology to observe the behavior of burning droplets in microgravity conditions, and they applied this approach to n-heptane droplets at different operative conditions, in static and convective environments [3, 4, 5]. Nakaya et al. [6] built upon these studies by analyzing the behavior of isolated n-decane and ethanol droplets in a CO₂-enriched environment. Other studies focused on the effect of the solid suspender [7] and the oxygen concentration on the droplet burning rate [8]. Larger droplets, up to 5 mm, are investigated with experiments in microgravity achieved from orbiting spacecrafts. The notable work by Dietrich et al. [9] analyzed the more pronounced importance of radiative effects for larger droplets. These experiments, performed on the International Space Station, revealed radiative extinction effects and the formation of cool flames [10], which could not be observed for smaller droplets.

These studies are essential for determining the optimal operating conditions for combustion processes, which maximize energy efficiency while minimizing pollutant formation. In this context, mathematical models of droplet combustion provide valuable insights that guide experimental investigations and aid in designing cleaner combustion systems.

Detailed spherically-symmetric (1D) models have been extensively used to simulate isolated droplets in many literature works, including [11, 12, 13, 14, 15, 16, 17]. The models revealed particularly useful for investigating the effect of radiation [12], preferential evaporation [18], heat conduction from the solid suspender [14], soot formation from pure and multicomponent droplets [17], and the formation of cool flames [15]. The advantage of these models is their limited computational cost, which enables a detailed description of chemical kinetics, radiation physics, and species diffusion, while retaining
the solution of the interface equilibrium conditions without sub-grid-scale
correlations. However, the inherent assumption of spherical symmetry in
these models limits their ability to capture multidimensional fluid dynamics
phenomena.

With increasing availability of computational power, new droplet models
based on multiphase CFD have emerged, enabling the direct solution of the
two-phase system. These models can resolve complex multidimensional phe-
nomena such as droplet deformation, breakups, buoyancy-driven flows, liquid
internal recirculation, turbulent conditions, and interactions between multi-
ple droplets. In recent years, several incompressible phase change models
using multiphase CFD were proposed, based on the front tracking approach
[19, 20], the level set [21, 22, 23], the Volume-Of-Fluid model [24, 25, 26, 27],
or on alternative interface representation methods [28, 29]. Despite the abun-
dance of multiphase models for the phase change, just a limited number of
contributions proposed the combination of the evaporation model with the
combustion chemistry. The introduction of combustion reactions comes with
several challenges:

1. The reaction terms introduce stiffness, coupling, and non-linearity to the transport equations.

2. The solution of the reaction terms strongly increases the computational cost of the model.

3. Chemical reactions drastically change the temperature and composition of the system. Thus, realistic description of reactive systems require a phase change model with variable thermodynamic and transport properties.

4. Additional physics such as radiation heat losses in gas phase are not negligible for accurate simulations.

5. The phase change model must be able to handle strong vaporization rate values, which are due to the high flame temperature.

6. The generation of species in gas phase leads to possible condensation which requires handling multicomponent vapor-liquid equilibrium (VLE) at the interface.

To the authors’ knowledge, the first attempt to introduce combustion in a
multiphase CFD model was made by Shao et al. [30] who simulated the
combustion of static, falling, and colliding droplets in a level set based phase
change framework, considering n-heptane and a global kinetic scheme. They
included the chemical reaction terms as explicit sources in the transport
equations for the chemical species mass fractions and temperature. Irfan and Muradoglu [31] used a front tracking based phase change model, and they improved the discretization of the reaction terms by using an operator splitting approach, which was used for simulating n-heptane droplets with global and skeletal kinetics. Both works rely on two strong assumptions: constant thermodynamic and transport properties, and negligible radiative heat flux. These effects were included in the VOF model proposed by Saufi et al. [32], who successfully simulated the effect of the solid suspender on methanol droplets in normal gravity conditions. The limitation of the latter model is the use of numerical methods tailored for the single isolated droplet configuration. Additionally, Saufi et al. [32] employs an artificial volumetric force for the droplet suspension, neglecting the surface tension force.

In this work, we consider the general multicomponent VOF-based phase change model proposed by Cipriano et al. [27], which was recently extended to a low-Mach variable material properties formulation [33]. We build upon this framework by including the direct integration of the chemical reaction terms, and an optically-thin model to describe the radiation heat loss in the gas phase. This work combines cutting-edge discretization techniques for phase change with the direct integration of the combustion chemistry. The model is used for simulations of realistic droplet combustion in microgravity conditions and for the simulation of suspended droplets in normal gravity, including variable properties, gas phase radiation, and water absorption inside the droplet.

Unlike widespread spherically symmetric models, the present multidimensional formulation is not limited to isolated droplet combustion in microgravity conditions. It can be applied to a wide range of configurations, including: (i) droplet deformation due to surface tension; (ii) interactions between multiple droplets; (iii) coalescence and breakup phenomena; (iv) buoyancy-driven flows; (v) internal liquid recirculation; and (vi) droplets in forced convection and turbulent flows. All these effects require the detailed solution of the fluid dynamics. In a multiscale perspective, these simulations allow the development of sub-grid-scale correlations in order to correct results from 1D models, which can accommodate complex kinetics at a significantly lower computational cost.
## 2. Mathematical model

This section introduces the governing equations for the model, together with the underlying assumptions. The equations are integrated over a generic control volume which comprises two immiscible phases (gas and liquid), separated by a zero-thickness interface. We distinguish the different phases using a sharp indicator function $H$, which varies in space $\mathbf{x}$ and time $t$, and it is defined as:

$$H(\mathbf{x}, t) = \begin{cases} 1 & \text{if liquid phase} \\ 0 & \text{if gas phase} \end{cases} \qquad (1)$$

According to Drew [34], the transport of the gas-liquid interface by the interfacial velocity $\mathbf{u}_\Gamma$ is achieved from the solution of a transport equation on the indicator function:

$$\frac{DH}{Dt} = \frac{\partial H}{\partial t} + \mathbf{u}_\Gamma \cdot \nabla H = 0 \qquad (2)$$

which must be solved with appropriate numerical methods in order to preserve the sharpness of the interface and mass conservation. The sharp interface description is valid, in the context of droplet combustion problems, in sub-critical conditions [35].

### 2.1. Navier-Stokes equations

The velocity $\mathbf{u}$ and pressure $p$ are obtained from the solution of the Navier-Stokes equations. The two phases are assumed to be compressible according to the low-Mach formalism proposed by [36, 37] and extended to two-phase flows in [33, 38]. Following this approach, the thermodynamic pressure of the system is considered constant, and the density varies according to the temperature and chemical species mass fraction gradients, with the following continuity equation, written for each phase $k = l, g$:

$$\nabla \cdot \mathbf{u}_k = - \frac{1}{\rho} \frac{D\rho}{Dt} \Big|_k = \beta_k \frac{DT_k}{Dt} + M_{w,k} \sum_{i=1}^{NS} \frac{1}{M_{w,i}} \frac{D\omega_{i,k}}{Dt} \qquad (3)$$

where $\rho$ is the density, $T$ is the temperature, $\omega_i$ is the mass fraction of the chemical species $i$, $\beta$ is the thermal expansion coefficient, while $M_w$ and $M_{w,i}$ are the molecular weights of the mixture and of the single chemical species, respectively. The last term on the RHS (Right-Hand Side) of Equation 3 is derived under the hypothesis of ideal gases, and it is neglected for the liquid
phase, assuming that the thermal expansion dominates over composition effects. This is always valid for pure liquid droplets, and it is generally valid also for the systems investigated in this work. This approximation must be relaxed if multicomponent droplets with large density differences are simulated.

The momentum equation is written in non-conservative form:

$$ \rho_k \frac{D\mathbf{u}_k}{Dt} = -\nabla \cdot \boldsymbol{\tau}_k - \nabla p_{d,k} - \mathbf{g} \cdot \mathbf{x} \nabla \rho_k \quad (4) $$

where $p_d = p - \rho \mathbf{g} \cdot \mathbf{x}$ is the pressure without the hydrostatic contribution, while $\mathbf{g}$ is the gravitational acceleration. The gravitational term is re-written as a function of the density gradients and the dynamic pressure, in order to simplify the formulation of the boundary conditions by applying the effect of gravity just where the density changes. Note that, for combustion modeling, the Boussinesq approximation is not sufficiently accurate to account for the strong temperature and composition changes of the system [36]. The viscous stress tensor $\boldsymbol{\tau}$ is computed by neglecting the compressible part which, for the low-Mach systems under investigation, is typically small compared to the incompressible part. Considering the dynamic viscosity $\mu$, the viscous stress tensor is obtained as: $\boldsymbol{\tau} = -\mu(\nabla \mathbf{u} + (\nabla \mathbf{u})^T)$. Equation 4 is valid both for the gas and for the liquid phases.

## 2.2. Scalar transport equations

The transport equations for the chemical species mass fractions and temperature derive from a set of conservation laws on the mass of each chemical species and on the enthalpy of the system:

$$ \rho_k \frac{D\omega_{i,k}}{Dt} = -\nabla \cdot \mathbf{j}_{i,k} + \dot{\Omega}_{i,k} \quad (5) $$

The species equation is written in non-conservative form and solved for the chemical species mass fractions. The convective term is computed using the velocity of the phase being solved, while the diffusive flux $\mathbf{j}_i$ is obtained from a modified Fick's law, as a function of the mole fractions $x_i$:

$$ \mathbf{j}_{i,k} = -\rho_k \mathcal{D}_{i,k} \frac{M_{w,i}}{M_{w,k}} \nabla x_{i,k} + \omega_{i,k} \sum_{i=1}^{NS} \rho_k \mathcal{D}_{i,k} \frac{M_{w,i}}{M_{w,k}} \nabla x_{i,k} \quad (6) $$

where the first term on the RHS is the classic Fickian diffusion law, while the second term is the diffusion velocity correction proposed by Coffee and
Heimerl [39], which guarantees closure of the diffusive fluxes when the mixture components have different diffusivities $D_i$. The diffusivities are computed from the binary diffusion coefficients $D_{i,j}$ using the Hirschfelder and Curtiss approximation [40] for the gas phase, and the Leffler Cullinan approach [41] for the liquid phase:

$$D_{i,g} = \frac{\sum_{j \neq i}^{NGS} x_{j,g} M_{w,j}}{M_{w,g} \sum_{j \neq i} \frac{x_i}{D_{ij,g}}} \quad D_{i,l} = \frac{\Pi_{j \neq i}^{NLS} (D_{ij,g} \mu_{j,l})^{x_{j,g}}}{\mu_l} \quad (7)$$

Equation 6 is re-arranged in the following form in order to preserve the functional dependency from the mass fractions, thus keeping the implicit part of the numerical solution of the diffusion term in the species equation [33].

$$\mathbf{j}_{i,k} = -\rho_k D_{i,k} \nabla \omega_{i,k} + \left[ \left( \sum_{i=1}^{NS} \rho_k D_{i,k} \frac{M_{w,i}}{M_{w,k}} \nabla x_{i,k} \right) - \rho_k D_{i,k} \frac{\nabla M_{w,k}}{M_{w,k}} \right] \omega_{i,k} \quad (8)$$

The last term on the RHS of Equation 5 is the formation rate $\dot{\Omega}_i$, included just in gas phase and computed as explained in Section 2.5, thus neglecting liquid phase reactions.

The conservation of enthalpy is exploited to obtain an equation for temperature, assuming negligible viscous heating and pressure work terms, which are small for low-Mach number systems [42]:

$$\rho_k C p_k \frac{DT_k}{Dt} = -\nabla \cdot \mathbf{q}_k - \left( \sum_{i=1}^{NS} C p_{i,k} \mathbf{j}_{i,k} \right) \cdot \nabla T_k - \sum_{i=1}^{NS} h_{i,k} \dot{\Omega}_{i,k} \quad (9)$$

where $Cp$ and $Cp_i$ are the mixture and species heat capacities, $\mathbf{q}$ is the heat flux vector, $h_i$ is the enthalpy of each individual species, and $\sum_{i=1}^{NS} h_i \dot{\Omega}_i$ is the heat of reactions. The heat flux is computed as:

$$\mathbf{q} = -\lambda \nabla T + \mathbf{q}_{\text{rad}} + \mathbf{q}_{\text{spark}} \quad (10)$$

it comprises the heat conduction term, calculated using Fourier's law where $\lambda$ is the mixture thermal conductivity, the radiation heat flux, and the heat provided by the electric wire for the ignition of the mixture. The heat dissipated by the radiation in gas phase is described by an optically-thin model [43], assuming that the only radiative chemical species are $H_2O$, CO, $CO_2$, and $CH_4$. This model is valid for flames with small optical thickness, much smaller than
1, as it assumes that the medium emits radiative energy, but it does not absorb any of the emitted energy [44]. The radiation contribution can then be written as:

$$ \nabla \cdot \mathbf{q}_{rad} = -4\sigma_{SB}a_p(T^4 - 300^4) \quad (11) $$

where $\sigma_{SB}$ is the Stefan-Boltzmann constant, 300 K is the reference gas phase bulk temperature, while $a_p$ is the Plank mean absorption coefficient, computed as:

$$ a_p = \sum_{i=1}^{NRS} P_i a_{p,i} \quad (12) $$

where $P_i$ is the partial pressure of the chemical species, $a_{p,i}$ is the extinction coefficient, while $NRS$ is the number of radiative species, computed as in [43, 45]. The heat flux provided by the spark $\nabla \cdot \mathbf{q}_{spark}$ is set to a constant value in a small circle close to the droplet surface, for a small simulation time. The specific values used depend on the simulation and they are specified in Section 4.

It is important to note that Equations 5 and 9 are solved both for the gas and for the liquid phases. However, the chemical reaction terms, the radiation heat flux, and the heat provided by the spark are included only in gas phase. Although the method that we propose is general and it can be used also for the solution of the liquid phase kinetics, in this work we neglect chemical reactions in the condensed phase which do not occur in the systems under investigation. The introduction of liquid phase reactions can rely on the same methodology, but it requires handling the possible increase of temperature in liquid phase due to the heat of reactions. This can be challenging if the temperature overcomes the boiling temperature of the light components, triggering internal boiling effects. These phenomena require careful modeling, which will be developed in future works. Additionally, the current formulation neglects the Soret effects, as previous studies on similar configurations have shown that its impact on the numerical results is negligible [16]. However, it may become significant in the presence of light species, such as when helium is used instead of nitrogen as the inert diluent [46], and it will be investigated in the future. On the contrary, the Dufour effect is usually small and it can be neglected [see 47, pag. 768].

### 2.3. Interface jump conditions

The set of equations 3, 4, 5, 9 are valid for each individual phase, and they must be solved by including appropriate gas-liquid interface boundary
conditions. Such conditions are obtained by integrating the mass, momentum, chemical species, and enthalpy conservation laws along an infinitesimal portion of the interface. Using the jump notation $[\phi]_\Gamma = \phi_l - \phi_g$, we express the interface jump in a generic quantity $\phi$ between the gas and the liquid phase. By doing so, we obtain a set of jump conditions that characterize the interface:

$$[\mathbf{u}]_\Gamma \cdot \mathbf{n}_\Gamma = \dot{m} [1/\rho]_\Gamma \qquad (13)$$

$$[p]_\Gamma = \sigma\kappa - \dot{m} [\mathbf{u}]_\Gamma \cdot \mathbf{n}_\Gamma \qquad (14)$$

$$[\rho (\mathbf{u} - \mathbf{u}_\Gamma) \cdot \mathbf{n}_\Gamma]_\Gamma = 0 \qquad (15)$$

$$[\rho\omega_i (\mathbf{u} - \mathbf{u}_\Gamma) \cdot \mathbf{n}_\Gamma + \mathbf{j}_i \cdot \mathbf{n}_\Gamma]_\Gamma = 0 \qquad (16)$$

$$[\lambda \nabla T \cdot \mathbf{n}_\Gamma]_\Gamma = \sum_{i=1}^{NS} \dot{m}_i [h_i]_\Gamma + \dot{q}_{rad} \qquad (17)$$

where $\dot{m}$ is the total vaporization rate per unit of surface, i.e. the sum of the vaporization rate for each chemical species $\dot{m}_i$. The velocity jump (Equation 13) contains the Stefan flow contribution due to the phase change, the pressure jump (Equation 14) includes the surface tension (where $\sigma$ is the surface tension coefficient and $\kappa$ is the interface curvature) and the recoil pressure term [48, 38]. The tangential momentum jump is not reported because the Marangoni convection is neglected. This effect requires accurate numerical methodologies [49] and will be addressed in dedicated studies, as it can influence liquid recirculation under reduced gravity conditions. Equations 15, 16, and 17 express the conservation of mass, species, and energy across the interface. The energy balance includes the interface radiation term, which is computed as $\dot{q}_{rad} = \epsilon\sigma_{SB}(\hat{T}^4 - T_{bulk}^4)$, where $\epsilon$ is the emissivity of the liquid fuel, which is a material property, while $\hat{T}$ and $T_{bulk}$ are the interface and the gas phase initial environment temperatures, respectively.

The system of jump conditions is closed by including the thermodynamic equilibrium at the gas-liquid interface, which translates in the continuity of the temperature and fugacity for each chemical species: $[T]_\Gamma = 0$ and $[f_i]_\Gamma = 0$. Assuming an ideal gas phase, unitary activity coefficient, and negligible Poynting correction, the VLE condition is written:

$$P \hat{x}_{i,g} = P_{vap}(\hat{T}) \hat{x}_{i,l} \qquad (18)$$

where $\hat{x}$ are the interface mole fractions, $P$ is the thermodynamic pressure of the system, while $P_{vap}$ is the vapor pressure, which is computed from Antoine's law.
The combustion chemistry included in this work does not affect the interface jump conditions, because the chemical reactions occur in the bulk gas phase, and there are no (heterogeneous) chemical reactions at the interface. The presence of heterogeneous chemistry would require including the reaction terms in equations 16 and 17 (see for example [50]).

## 2.4. Thermodynamic and transport properties

The thermodynamic and transport properties of each phase are calculated as a function of the thermodynamic state of the mixture: $\phi = f(P, T, x_i)$, selecting the pressure, temperature, and mole fractions as the state variables. This operation is performed using the OpenSMOKE++ library [51], which computes the gas phase properties using the CHEMKIN formalism. The liquid properties cannot be obtained from a unified theory, similarly to the gas phase. Therefore, empirical correlations reported in engineering handbooks are used for the liquid phase properties, such as those reported in [52]. The only properties which are kept constant in space and time during the simulation are the emissivity of the liquid $\epsilon$, and the surface tension coefficient $\sigma$, which is constant because Marangoni effect is neglected in this study. Small variations in the droplet shape due to the droplet heating and consequent change in the surface tension coefficient do not affect the results in this work appreciably.

## 2.5. Chemical kinetics

The formation rates $\dot{\Omega}_i$ in Equations 5 and 9 are calculated from the reaction rates $\dot{r}_j$ and the stoichiometric coefficients $\nu_{j,i}$:

$$ \dot{\Omega}_i = M_{w,i} \sum_{j=1}^{NR} \dot{r}_j \nu_{j,i} \quad (19) $$

where $j$ is the index of a generic reaction, while $NR$ is the total number of reactions. The expression for the reaction rates depends on the choice of the kinetic scheme. In this work, we use the two-step schemes proposed by Westbrook and Dryer [53], consisting of the following reactions:

$$ C_n H_m + \left(\frac{n}{2} + \frac{m}{4}\right) O_2 \xrightarrow{k_1} nCO + \frac{m}{2} H_2O $$

$$ CO + \frac{1}{2} O_2 \xrightleftharpoons{k_2} CO_2 $$
The partial oxidation of the fuel is described by the following reaction rate:

$$ \dot{r}_1 = k_1 C_{C_n H_m}^a C_{O_2}^b \qquad (20) $$

where C indicates the concentration of the chemical species, $k_1$ is the kinetics constant of the first reaction, while the exponents a and b are constants that depend on the fuel. The kinetics constant is described by an Arrhenius type equation<fnref n="1" />:

$$ k_1 = A \exp(-E_a/RT) \qquad (21) $$

where A is the pre-exponential factor, $E_a$ is the activation energy, while R is the gas constant. The values of A and $E_a$ depend on the fuel, and they are reported in Section 4 for each test case. The rate of CO oxidation is described by the following forward and backward reactions [53]:

$$ \dot{r}_2 = 10^{14.6} \exp(-40/RT) C_{CO}^{1} C_{H_2O}^{0.5} C_{O_2}^{0.2} \qquad (22) $$

$$ \dot{r}_2 = 5 \times 10^8 \exp(-40/RT) C_{CO_2}^{1} \qquad (23) $$

The OpenSMOKE++ library [51] is employed to compute reaction rates and heat of reaction from a kinetic scheme defined in CHEMKIN format [54]. This framework is flexible and capable of handling arbitrarily complex kinetic mechanisms. However, significant limitation arises from the computational cost, which increases non-linearly with the number of species in the system. Therefore, the choice of the two-step kinetics allows the model to be validated by simplifying the chemistry in order to minimize the computational time. This approximation restricts the range of chemical phenomena that can be accurately captured; for instance, it is not suitable for droplet combustion involving significant soot formation or for simulations of auto-ignition and cool flames. More generally, two-step kinetic schemes are unable to describe phenomena that require complex radical chemistry to achieve accurate results. While future work will employ more complex kinetic mechanisms, the current approach prioritizes computational efficiency and focuses on the coupling between the combustion and the phase-change model.

<footnote n="1">¹Note that the OpenSMOKE++ [51] framework, used to manage the kinetics, thermodynamic and transport properties, can accommodate kinetics of any complexity, including the modified form of the Arrhenius equation: $k_j = A_j T^{n_j} \exp(-E_{a,j}/RT)$. For the kinetic schemes used in this work $n_j = 0$.</footnote>
## 2.6. Overview of model assumptions

For completeness, this section provides a list of the model’s main assumptions, which may limit its applicability under certain operating conditions.

1. **Sub-critical conditions:** the systems under investigation are characterized by a clearly defined gas-liquid interface dynamics.

2. **Thermodynamic vapor-liquid equilibrium:** the characteristic time of the system allows interface thermodynamic equilibrium to be achieved.

3. **Low-Mach-number systems, with constant thermodynamic pressure:** we filter-out acoustic waves but we consider density changes due to temperature and mass fraction gradients.

4. **Ideal gas and liquid mixtures:** needs to be relaxed to achieve more accurate results at high pressure or for non-ideal liquid mixtures.

5. **Negligible Soret and Dufour effects:** although the Dufour term is usually small [47], the Soret effect is important when small species (i.e. hydrogen or helium) are present, with soot formation, and when the flame interacts with hot solid walls.

6. **Negligible Marangoni effect:** it enhanced liquid mixing, and it can explain the recirculation observed in methanol-water systems in micro-gravity conditions.

7. **Negligible thermal influence of the solid suspender:** large quartz fiber may influence the burning rate significantly, triggering liquid re-circulation through thermal Marangoni convection, promoting natural convective fluxes in the gas phase, and also quenching the flame.

8. **Optically-thin radiation:** accurate for small droplets and low pressure values. More detailed models are required if the optical thickness is not small and the absorption of the emitted energy is non-negligible.

9. **Absence of liquid phase reactions:** this assumption is reasonable for relatively small hydrocarbons, but may become critical when modeling the combustion of polymeric liquids and heavy fuel oils.

10. **Absence of bubble formation within the liquid phase:** observed for liquid mixtures whose components have a wide range of volatility.

11. **Newtonian fluids:** this is a common assumption in droplet combustion models.

The use of complex kinetics models is not included in the list of assumptions, because the current formulation of the model is already able to simulate complex kinetic mechanisms. Assumptions 1-3 form the foundation of the
present model. Relaxing any of these assumptions would require a fundamentally different framework and is therefore beyond the scope of this work. In compressible and supercritical regimes, where thermodynamic equilibrium cannot be assumed, a different formulation must be employed. The other assumptions can be relaxed depending on the systems under investigations.

### 3. Numerical methods

The transport equations introduced in Section 2 are discretized using the Finite Volume Method [55], in a quad/octree grid. This section briefly describes the numerical solution of the equations, focusing on the strategy used to resolve the gas-liquid interface, and to impose the boundary conditions on the governing equations. The solution of the reactive step, which is the peculiarity of this work, is described in detail in Section 3.4. Further details on the discretization techniques can be found in [27] and [33].

#### 3.1. Geometric Volume-Of-Fluid

The objective of the geometric VOF approach is to transport the interface while maintaining the sharpness of the indicator function $H$, which is discretized by defining the volume fraction $c$ in each computational cell of the domain, with cell volume $V$:

$$c = \frac{1}{V} \int_V H dV \qquad (24)$$

obtaining pure gas cells ($c = 0$), pure liquid cells ($c = 1$), and interfacial cells with intermediate volume fraction values. This definition is applied to Equation 2, obtaining the following transport equation for the volume fraction:

$$\frac{\partial c}{\partial t} + \frac{1}{V} \oint_S H \mathbf{u}_\Gamma \cdot \mathbf{n} dS = \frac{1}{V} \int_V H \nabla \cdot \mathbf{u}_\Gamma dV \qquad (25)$$

which is solved using the directionally-split method described by Scardovelli and Zaleski [56] and improved by Weymouth and Yue [57] to guarantee discrete volume conservation in incompressible flows. When the same approach is applied to a low-Mach phase change model, the liquid volume changes according to the phase change and to the density variations. This approach results in the sharp and conservative transport of the gas-liquid interface.
### 3.2. Solution of the interface jump condition

The solution of the interface jump conditions is necessary to obtain the vaporization rates $\dot{m}_i$, the interface mass fractions $\hat{\omega}_{i,g}$, $\hat{\omega}_{i,l}$, and the interface temperature $\hat{T}$. These quantities are used to enforce the gas-liquid interface boundary conditions in the governing equations, as explained in Section 3.3 and 3.4. The system of jump conditions is solved by exploiting the decoupled procedure explained in [27]. The liquid phase composition is assumed to be equal to the interface mass fraction on the liquid phase side $\hat{\omega}_{i,l}$. The mass fractions are converted into mole fractions $\hat{x}_{i,l}$ and the thermodynamic equilibrium (Equation 18) is used to calculate $\hat{x}_{i,g}$, which is converted back to $\hat{\omega}_{i,g}$. The total vaporization rate is obtained as:

$$ \dot{m} = \frac{\sum_{i=1}^{NLS} \mathbf{j}_{i,g} \cdot \mathbf{n}_{\Gamma}}{1 - \sum_{i=1}^{NLS} \hat{\omega}_{i,g}} \qquad (26) $$

where the diffusive fluxes are calculated as a function of the interface gradients, obtained using the same approach reported in [58, 59, 27]. The vaporization rate for each chemical species is obtained from:

$$ \dot{m}_i = \dot{m}\hat{\omega}_{i,g} + \mathbf{j}_{i,g} \cdot \mathbf{n}_{\Gamma} \qquad (27) $$

which is equal to zero for non-condensable species. The interface temperature is finally updated by zeroing Equation 17, given the vaporization rates:

$$ f(\hat{T}) = \sum_{i=1}^{NLS} \dot{m}_i \Delta h_{ev,i} + \dot{q}_{rad} - \lambda_l \nabla T_l \cdot \mathbf{n}_{\Gamma} + \lambda_g \nabla T_g \cdot \mathbf{n}_{\Gamma} = 0 \qquad (28) $$

This procedure assumes that the evaporation process is driven by a gradient of chemical species mass fractions, and it does not include automatic transition to a boiling regime. This assumption reflects the conditions studied in this work, although it may require extension for simulations of very small droplets (i.e. $D \le 30$ µm for the methanol droplets studied by Awasthi et al. [60]). The solution of the interface jump condition for multicomponent liquids allows condensation to be directly included. In the condensation case, the interface gradients of the diffusive fluxes have opposite sign, leading to the absorption of gas in liquid phase.
### 3.3. One-field solution of the Navier-Stokes equations

The Navier-Stokes equations are solved using a one-field approach [61], meaning that a single set of equations is solved over the whole two-phase domain, by defining variables and properties which vary abruptly at the gas-liquid interface:

$$ \nabla \cdot \mathbf{u} = \dot{m} \left[ \frac{1}{\rho} \right]_{\Gamma} \delta_{\Gamma} - \frac{1}{\rho} \frac{D\rho}{Dt} \quad (29) $$

$$ \rho \frac{D\mathbf{u}}{Dt} = -\nabla \cdot \boldsymbol{\tau} - \nabla p_d - \mathbf{g} \cdot \mathbf{x} \nabla \rho + (\sigma\kappa - \dot{m}^2 [1/\rho]_{\Gamma}) \nabla c \quad (30) $$

where the properties and variables, like density $\rho$ and velocity $\mathbf{u}$ are weighted by the volume fraction using a simple arithmetic average:

$$ \phi = \phi_l c + \phi_g (1 - c) \quad (31) $$

while the phase properties are obtained from the thermodynamic state of the mixture: $\phi_k = f(P_k, T_k, \mathbf{x}_k)$, as explained in Section 2.4. Therefore, Equations 29 and 30 reduces to Equations 3 and 4 in the bulk of the gas and liquid phases, while the interface conditions are introduced, following the Whole Domain Formulation [56], as source terms localized in the interfacial cells.

Equations 29 and 30 are discretized using the collocated algorithm described by Popinet in [62, 63]. The advection step is solved in an explicit manner using the second order upwind Bell-Colella-Glaz approach [64]. The diffusion step is solved implicitly using a multigrid solver [65], while the projection step is modified by including the volume expansion due to the phase change and the density Lagrangian derivative [33].

When solving the Navier-Stokes equations using a one-field approach with phase change, a proper approach must be used to extract the liquid or gas phase velocities for the interface advection. In this work, we use the Double Pressure Velocity Coupling introduced in [27]. According to this approach, the solution of Equations 29 and 30 is repeated by neglecting the volume expansion due to the phase change $\dot{m}[1/\rho]_{\Gamma}$. Following this strategy, we obtain an extended velocity $\mathbf{u}_E$, representative for the liquid phase velocity, whose divergence is equal to the liquid density Lagrangian derivative. Finally, the interface velocity is obtained from Equation 15:

$$ \mathbf{u}_{\Gamma} = \mathbf{u}_{E} - \frac{\dot{m}}{\rho_{l}} \mathbf{n}_{\Gamma} \quad (32) $$
### 3.4. Two-field solution of temperature and chemical species with reactions

The transport equations for chemical species mass fractions and temperature are solved using a two-field approach: Equations 5 and 9 are multiplied by the indicator function and integrated over the whole domain, solving two different equations for the same fields in the liquid and in the gas phase, respectively [66]. The interface jump conditions appear into the governing equations as source terms localized in the interfacial cells. The resulting equations are solved by splitting the advection, diffusion, and chemical reactions terms, in order to approximate each contribution appropriately.

The advection term is discretized using the same geometric directional-split procedure used for the VOF field (Section 3.1), using the second order upwind Bell-Colella-Glaz approach [64] to approximate values of the scalar fields from the cell centers to the cell faces.

The diffusion term is re-arranged into a Poisson-Helmholtz equation, comprising the phase change source terms, and solved in an implicit manner using a multigrid solver [65].

The complications that arise when introducing the reaction step are mainly due to the high stiffness and the strong non-linearity of the reaction terms. These problems were tackled by different authors in the context of single phase flame simulations, by proposing a fully coupled solution of the governing equations, decoupled solutions with linearization of the source terms, and the operator-splitting method [67, 68]. The latter approach is adopted here because it allows the reactive terms for chemical species and temperature to be combined in an ODE system of equations:

$$ \left\{ \begin{aligned} \frac{dy}{dt} &= f(y, t) \\ y(t_0) &= y \end{aligned} \right. \qquad (33) $$

where $y = (\omega_i, T)$ is the vector of mass fractions and temperature values in the computational cell, while $f$ is the RHS of Equations 5 and 9, considering only the chemical reactions. This ODE system is integrated implicitly by a full time step $\Delta t$, using the native ODE solver implemented in the OpenSMOKE++ library [51] and specifically conceived for stiff reactive systems. The algorithm employs a variable coefficients backward differentiation formula (BDF) [51]. The advantages of this approach are the following:

1. The non-linearity of the reactive step is preserved by advancing the chemical species and temperature together.
2. This method does not require the integration time step of the fluid dynamics to be adjusted based on the stability conditions of the chemical reactions.

3. The ODE system integration depends only on the local composition and temperature of the cell, therefore the resulting algorithm can be easily implemented and efficiently parallelized.

Despite the many advantages, the solution of the ODE system represents by far the most time demanding task of each time step. To overcome this problem, many approaches to reduce the computational cost of the chemistry were proposed, although their implementation in a multiphase framework are beyond the aim of this work. Possible methods for the speed-up of the reactive step are the cell agglomeration technique [69], the tabulated chemistry [70], and the In-Situ Adaptive Tabulation (ISAT) [71].

# 4. Results and discussion

The model described at the previous sections is implemented in the open-source code Basilisk [72], a powerful tool for direct numerical simulations of gas-liquid flows. The code implements state-of-the-art methodologies for surface tension modelling in a quad/octree grid, which allows accurate and time efficient simulations of capillary driven flows, such as small droplets and bubbles. This framework was extended in [27, 33] including multicomponent phase change and a low-Mach formulation for variable thermodynamic and transport properties. In this study, we extend the capabilities of the code by incorporating chemical reactions in the gas phase, utilizing the coupling between Basilisk [72] and the OpenSMOKE++ library [51]. Validation is performed by comparing the model's results with experimental data in microgravity conditions [5, 6], and with a simplified spherically-symmetric (1D) model [73] (Sections 4.2 and 4.3). These test cases demonstrate the model's convergence and its validation under various operational conditions for different fuels. In Section 4.4 we discuss the different behavior of methanol droplets in microgravity and in normal gravity conditions, measuring the variation of the consumption rate when natural convective fluxes are present. Lastly, in Section 4.5, the multicomponent model is applied to calculate the droplet extinction resulting from water condensation effects during methanol combustion. All simulations use global kinetic schemes to focus on the methodology while minimizing computational time.
4.1. *Simulation Setups*

The simulations presented in this work are performed using two different configurations, depicted in Figure 1. The combustion of droplets in microgravity conditions is simulated using the initial setup in Figure 1 (a), where a droplet is initialized at the left-bottom corner of the square domain. This configuration exploits the symmetry axis to reduce the computational time and to avoid lateral drifts of the droplets that may break the symmetry of the problems. Despite the model used for surface tension is well-balanced, and it respects the equilibrium pressure jump (Laplace equation) [74], it does not guarantee conservation, meaning that the sum of the surface tension contributions in each interfacial cell of the droplet may not cancel out, leading to lateral droplet displacement. This effect is worsened by the disturbance given by the introduction of the phase change. To resolve this issue, we exploit the symmetry axis of the domain.

For droplets in microgravity conditions, we set symmetry boundary conditions on the bottom and left boundaries in contact with the droplet, while outflow boundary conditions are imposed on the top and right sides. This configuration runs in 2D with axial symmetry around the left side. Simulations in normal gravity exploit similar square domain with axial symmetry around the left axis, but introducing a model for the suspension by the action solid fiber. As shown in Figure 1 (b), half droplet is initialized along the left boundary of the domain, with gravity **g** pointing downward. The left side is shifted from the axis of symmetry by a length which corresponds to the fiber radius. The droplet suspension strategy developed by Cipriano et al. [33] is used, allowing the droplet to be suspended at the pinning point **x**<sub>p</sub> by the surface tension force. Outflow boundary conditions are used for the top, bottom, and right boundaries, while no-slip boundary conditions are used for the left boundary, mimicking the surface of the solid suspender. The solid fiber's presence is considered purely geometrically (i.e., for droplet suspension), while heat transfer between the solid and the gas-liquid interface is neglected; this will be addressed in future studies.

In both the configurations, we use an adaptive quadtree grid, with wavelet-
based error estimation [75]. The target fields for refinement and coarsening
are the fuel mass fraction, the temperature, and the velocity fields, while
the 3x3 stencil around the interface is always maintained at the maximum
level of refinement. This approach significantly reduces simulation time by
minimizing the number of cells required for the direct integration of combus-
tion chemistry. The grid is always refined around the droplet interface and
<figure><img src="image_2.png" /><figcaption>Figure 1: Sketch of the simulation setups in microgravity (a), and in normal gravity conditions (b). The blue line indicates the droplet interface, the red line is the flame front, **g** is the gravity, while **x**<sub>p</sub> is the pinning point at which the droplet is suspended.</figcaption></figure>

the flame front, while the domain length is maintained large to avoid the
influence of the boundaries on the droplet combustion process.

4.2. N-Heptane combustion in microgravity

As a first test case, we consider the combustion of a n-heptane droplet in microgravity conditions. This test aims at verifying the convergence and the accuracy of the numerical model. To this purpose, we compare the simulation results with experimental data from Hara and Kumagai [5], and with the 1D model by Cuoci et al. [73], which performs the direct solution of the system under the assumption of spherical symmetry. We consider this model to be reliable as it was extensively validated in several operative conditions [13, 16, 73, 76]. The simulation setup is reported in Figure 1 (a), considering a n-heptane droplet in air, and neglecting the effect of gravity. The initial droplet diameter is D₀ = 0.38 mm, the initial temperature of both the droplet and the environment is 300 K, while the ambient pressure is 0.1 MPa. The surface tension value is σ = 0.03 N/m, and we select a domain length which is 40 times larger than the initial droplet diameter (L₀ = 40D₀). The simulation runs at three different levels of refinement: 10, 11, and 12. Using this notation, the maximum number of cells along each dimension of the domain is 2<sup>LEVEL</sup>, which corresponds to approximately 26, 51, and 102
<figure><img src="image_3.png" /><figcaption>Figure 2: Maps of the temperature and mass fraction fields for the combustion of n-heptane in microgravity conditions at time 0.1 s.</figcaption></figure>

cells per diameter, respectively. The chemical reactions are activated by a spark, which is included by setting a constant heat flux ($q_{spark}$ in Equation 9) in a small circle with diameter $0.2D_0$, located at 1.2 mm from the droplet surface, from time 0 s to time 0.008 s and with intensity $\nabla \cdot q_{spark} = 100$ MW/m$^3$/K. Note that this value is obtained from optimal numerical tuning, and not from the experimental setup. The two-step kinetics described in Section 2.5 is considered, where the kinetic parameters of the fuel oxidation are: $A = 6.3 \times 10^{11}$, $E_a = 30000$, $a = 0.25$, $b = 1.5$ [53]. The same kinetic scheme and the same radiation model are used by the 1D simulations.

Figure 2 shows the temperature and mass fraction fields at simulation time 0.1 s. The black line represents the gas-liquid interface, while the white line is the flame front, defined in this section as the temperature peak isoline. The maps show the formation of the spherical flame, due to the absence of gravity, and its strong expansion with respect to the droplet diameter, which justifies the use of a large domain. The trends of the mass fractions in Figure 2 show that the n-heptane evaporates and it reaches its concentration peak around the droplet surface. It is consumed by the oxygen in the environment leading to the formation of $H_2O$, and CO in the rich part of the
flame, which is the one with larger amount of carbon atoms. Further oxidation of CO produces CO₂ which exhibits a trend similar to that of H₂O (not shown), while the inert species remains unaffected by the combustion chemistry.

The time evolutions of the square droplet diameter and of the flame diameter are plotted in Figure 3. We observe droplet ignition in the very first instants of the simulation, after which the square diameter decreases almost linearly throughout the entire droplet lifetime. The temperature of the environment is sufficiently hot to drastically increase the vaporization rate, hiding unsteadiness in the droplet consumption, which is typically observed for pure evaporation cases [33]. The trend of the square diameter matches well the experimental data already at level 10, and it converges toward the 1D model with increasing mesh resolution. The values of the burning rate, which is the slope of the square diameter line with opposite sign, are 0.784 mm<fnref n="2" />/s, 0.782 mm<fnref n="2" />/s, and 0.780 mm<fnref n="2" />/s, respectively for the three different levels of refinement. Increasing the grid resolution, the burning rates tend to the value predicted by the 1D model, whose burning rate is 0.769 mm<fnref n="2" />/s and to the experimental data (0.71 mm<fnref n="2" />/s)<fnref n="2" />. Minor discrepancies between the experimental and numerical diameters arise from the use of simplified kinetics.

Increasing the grid resolution, also the flame diameter tends to the 1D model, as shown in Figure 3. The flame diameter exhibit a sudden increase right after the ignition period, and it keeps expanding with respect to the droplet diameter throughout the entire simulation. Although the numerical results match very well, there is a discrepancy between the models and the experimental data. This behavior has already been observed in other numerical works [73], and it can be explained by considering that the numerical and the experimental definitions of the flame front are substantially different. Experimentally, the flame front is obtained from the droplet combustion pictures, by considering the region of the flame with highest luminosity [5]. This region may differ from the temperature peak, which instead is used by the numerical simulations to define the position of the flame front [73].

It may seem surprising that the model accurately captures the consump-

<footnote n="2">²The numerical burning rates are obtained from linear interpolation of the square diameter in time, excluding the initial region of the plot which contains unsteadiness due to the droplet ignition.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3: Plots of the square droplet diameter and the flame diameter in time. The points are experimental data from Hara and Kumagai [5], the triangles are results from the reference spherically-symmetric model [73], while the lines are numerical results from this work at different levels of refinement.</figcaption></figure>

tion dynamics despite significant simplifications in the chemical kinetics. This can be explained by the operative conditions used in this test case, which minimize soot formation. The droplet studied here is sufficiently small, such that the time required for its consumption is shorter than the characteristic time for soot formation. This trend is confirmed by Hara and Kumagai [5], which does not report any soot shell diameter, unlike the study by the same authors in [4], where non-negligible soot formation was observed for larger n-heptane droplets ($D_0 \approx 0.9$ mm). Simulations of larger n-heptane droplets with this setup tend to underestimate the burning rate, leading to a longer droplet lifetime. The use of more detailed kinetic schemes, including soot formation, will eventually resolve the problem with consequent increase of the computational time. For these reasons, in the following sections, we limit our investigation to fuels and operative conditions that minimize soot formation.

The simulation times are approximately 11 hours, 42 hours, and 144 hours for the cases with 25, 51, and 102 cells per diameter, respectively. These simulations were run in parallel, using at most 48 cores (AMD EPYC 7H12) in order to minimize communication overhead between different cores. In the
following sections, we run the simulations just at the intermediate level of
refinement, in order to balance accuracy with computational cost.

4.3. *Ethanol combustion in microgravity*

In this second validation case we challenge our model by simulating differ-
ent operative conditions and a different fuel. To this purpose, we consider the
combustion of ethanol droplets in microgravity conditions, with various ini-
tial diameters ranging from 0.36 to 0.71 mm, in CO2 enriched environments.
The simulation setup is reported in Figure 1 (a). The initial temperature
is 300 K for both the droplet and the environment, and all simulations are
performed at atmospheric pressure (0.1 MPa). The initial droplet diameter
varies depending on the case, with exact values provided in Figure 4. In
each case, the domain length is set to 40D₀ to maintain the same number of
cells per diameter, regardless of the initial droplet dimension. Simulations
are performed using an adaptive grid with maximum level of refinement 11,
which is a good compromise between accuracy and time efficiency as demon-
strated in Section 4.2. The kinetics parameters for the partial oxidation of
ethanol are: A = 1.8 × 10¹², Eₐ = 30000, a = 0.15, b = 1.6 [53]. The
surface tension is set to σ = 0.02 N/m, and the combustion starts because
of the heat provided by an electric spark, which is included as discussed in
Section 4.2. The simulations are performed with four different initial droplet
diameters and at various air compositions: the oxygen content is fixed at
the level of pure air, while the N₂ and CO₂ concentrations are adjusted to
achieve environments with 0%, 30%, and 60% CO₂ in volume, respectively.

Figure 4 shows the time evolution of the square droplet diameter for all simulations in this section. Similarly to the n-heptane combustion case, the droplet diameter proceeds linearly and it obeys the d² law. By looking at the square diameters we observe a remarkable agreement between our model and the experimental data [6]. The experimental data show an initial plateau, which corresponds to the transients before the spark ignites the droplet. Instead of simulating this initial period, which is not of interest for the discussion in this work, we impose the spark directly at time 0 s (for a duration of 4 ms) and we shift our results in order to compare the model and the experiments just in the combustion region, right after the ignition. In this region, the square diameter history in Figure 4 does not show any important thermal expansion effects, but we observe a decreased burning rate with increasing initial amount of CO₂ in the atmosphere, for each initial diameter value. To understand this effect we look at the flame properties (temperature
<figure><img src="image_5.png" /><figcaption>Figure 4: Plots of the square droplet diameter in time (first layer) the flame diameters (second layer), the flame temperature (third layer), and the dimensionless square diameters (last layer) for the ethanol droplet combustion at varying initial diameters and CO₂ enriched environments. In all the plots, the circles are experimental data by Nakaya et al. [6], while the triangles are numerical results from the simplified 1D model [73].</figcaption></figure>
and diameter), which are reported in Figure 4. Because these values are not provided by Nakaya et al. [6], we use the 1D model [73] introduced at the previous section in order to compare the prediction of our model. We observe that the results obtained using the multidimensional approach are able to capture the prediction from the 1D model accurately. The flame diameter decreases with the initial amount of CO₂, and the temperature of flame, which is considered as the peak temperature of the gaseous environment, decreases with increasing ambient CO₂ as well. Since a simplified kinetic scheme is employed, due to the non-sooting tendency of ethanol [6], the slower consumption of the droplets in rich CO₂ environment is due to the lower ambient (flame) temperature. This temperature value is caused by the variable material properties, in particular the factor ρCp which slightly increases with the CO₂ concentration, and by the effect of radiation, which increases with the CO₂, decreasing the ambient temperature, and therefore the flame temperature. On the contrary, the thermal diffusivity (λ/ρ/Cp) decreases with increasing amount of CO₂, limiting the dispersion of the flame.

If we fix the CO₂ concentration, and we compare the burning rates of the droplets at different initial diameters, we observe that, in any case, the initial diameter decreases the burning rate. This effect is confirmed by Kazakov et al. [77], who compared numerical and experimental burning rates, highlighting that simplified models, which do not include radiative heat losses, cannot capture the decreased burning rate with increasing initial droplet diameter. Additionally, in this work we notice that the decreased burning rate is more evident with increasing CO₂ concentration, since it enhances the importance of non-luminous radiative heat losses.

From the simulations in this section we confirm the correct implementation of the model, which captures the experimental data and the prediction from the 1D model for droplets of different fuels, and in different operative conditions. We demonstrate the importance of including gas phase radiation, and its utility in explaining experimental trends.

### 4.4. Methanol combustion in normal gravity

In this section, we leverage the model's capability to resolve multidimensional effects by simulating the combustion of methanol droplets in both microgravity and normal gravity conditions, emphasizing the differences between the two regimes. Microgravity simulations are performed using the same setup shown in Figure 1 (a), while normal gravity simulations exploit the setup depicted in Figure 1 (b), which accounts for non-spherical flames
including the effect of gravity and the droplet suspension (or pinning) due to surface tension. The initial droplet diameter is set to $D_0 = 0.5$ mm, which is an intermediate value between the small droplets produced during atomization ($D_0 < 0.1$ mm) and the larger droplets typically studied in isolated droplet combustion experiments ($D_0 > 1$ mm). In the microgravity case, the domain length is set to $L = 40D_0$, whereas in the normal gravity case it is set to $L = 80D_0$. The same spatial resolution is maintained for both cases by applying a maximum refinement level of 11 in microgravity and 12 in normal gravity, which lead to the same number of cells per diameter in the two configurations. The gravitational acceleration is $g = 9.81 \text{ m/s}^2$, the suspending fiber diameter is set to $0.1D_0$, and the surface tension is $\sigma = 0.022 \text{ N/m}$. The fiber diameter is taken as the minimum value which allows the droplet to remain suspended. In a real experimental setups, more complex fiber configuration (i.e. cross, or rod with a spherical tip), and the complexity of the solid material surface may limit the fiber diameter. These effects will be considered in future studies. Water absorption is neglected in this section in order to focus only on the effect of buoyancy on the burning rates. All simulations are performed with an initial droplet and ambient temperature of 300 K, at pressure values of 0.1, 0.5, 1, and 2.0 MPa (subcritical). The kinetic parameters for methanol combustion are: $A = 3.2 \times 10^{12}$, $E_a = 30000$, $a = 0.25$, and $b = 1.5$ [53].

The chemical reactions are initiated by a spark, similarly to the previous test cases. Once combustion starts, the microgravity simulations exhibit the formation of a spherical flame, as observed in Figure 2. In contrast, under normal gravity conditions, the flame significantly reduces the gas-phase density around the droplet, creating density gradients, which are responsible for upward buoyancy-driven flows. These flows transport the scalar fields, as illustrated in Figure 5, and they can be quantified by the Grashof number, which represents the ratio between buoyancy and viscous forces:

$$Gr = \frac{(\rho_s - \rho_b)}{\rho_b \nu_b^2} gD^3 \quad (34)$$

where the subscripts s and b denote surface and bulk properties, respectively, while $\nu$ is the kinematic viscosity. The functional dependence of the Grashof number reveals that the initial droplet diameter has a significant effect on buoyancy-driven flows. For instance, increasing the droplet diameter from 0.1 mm to 1 mm results in a 1000 times larger Grashof number, with substantial influence on the system dynamics. The dependence of buoyancy-driven
<figure><img src="image_6.png" /><figcaption>Figure 5: Maps of the temperature and chemical species mass fractions at time $t = 0.1$ s for the methanol combustion in normal gravity conditions at atmospheric pressure.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Maps of the temperature at time $t = 0.1$ s for the methanol combustion in normal gravity conditions at different ambient pressures.</figcaption></figure>

flows on pressure ($Gr \propto P^2$) is less evident but still noteworthy, as it en-
hances natural convective fluxes quadratically. The effect of pressure on the
flame geometry is shown in Figure 6, which displays the temperature field
and flame front at the same time for the various pressures considered. As
ambient pressure increases, the flame becomes thinner, the upward stretch-
ing intensifies, and the flame front moves closer to the droplet surface. The
reduced diffusivity of the temperature in Figure 6 is due to the different
properties, and it is observed also for the microgravity droplets (not shown).
In particular, the strong increase of the gas phase density with the ambient
pressure limits the thermal diffusivity, while the thermal conductivity and
the specific heat are unaffected by the ambient pressure. Moreover, radiative
heat losses increase with pressure according to Equation 12.

Quantitative results for this test cases are reported in Figure 7, which
shows the trend of the dimensionless square diameter in time. We observe
that, in microgravity conditions, the pressure introduces slight discrepancies
from the $d^2$ law, due to the increased thermal expansion with increasing
ambient pressure. Subsequently, the droplet consumption proceeds linearly
and with a higher burning rates for the high pressure cases. In normal gravity
conditions, the discrepancy from the $d^2$ is larger, especially after the ignition,
due to the initial stabilization of the flame around the droplet surface, until
a stable *envelope flame* is formed. For every pressure value, the burning
rate in normal gravity $K_{\text{ng}}$ is always larger compared to the same case in
microgravity conditions $K_{\text{mg}}$. This discrepancy is measured as the ratio
between the two cases, as reported in Figure 7, and it can not be attributed
to the higher flame temperatures, which are very similar for each simulation
in this section (in the range between 1900 and 2000 K). As predicted by the
functional dependency of the Grashof number (Equation 34), the effect of
the buoyancy-driven flows justifies the increased gap between the burning
rates at higher ambient pressure. The ratio between the burning rates is an
important parameter which can be used to correct results from 1D model in
order to account for the effect of buoyancy-driven flows. These discrepancies
can be observed in Figure 7 using the 1D model from the OpenSMOKE++
suite [73].

Figure 7 also reports the comparison between the flame front in microgravity and in normal gravity conditions. Conversely from sections 4.2 and 4.3, in these simulations the flame front is calculated from the mixture fraction $Z$, computed as in [32], and defined as the isoline with mixture fraction equal to the stoichiometric mixture fraction $Z_{st}$. This approach is more convenient than the temperature peak profile for multidimensional simulations, because it allows a level set function $Z$ to be obtained, with a rich ($Z > Z_{st}$) and a lean region ($Z < Z_{st}$). This definition of flame front is more practical for multidimensional simulations, but it corresponds to the temperature peak only if the species diffusivities are equal. This is not a big issue because, in this section, we want to compare the flame fronts in microgravity and in normal gravity conditions, based on the same definition. The flame front for the normal gravity simulations is considered as the horizontal coordinate of the flame front.

Figure 7 shows that the flame front in microgravity approaches the droplet surface as the pressure increases. In normal gravity, the flame front along the horizontal direction shows an initial peak due to the stabilization of the flame. After the formation of the envelope flame, the expansion of the front is always limited with respect to the microgravity case, due to the influence of buoyancy-driven flows on the flame geometry. As a consequence, the flame remains always closer to the droplet surface leading to higher droplet consumption rates.

This test case confirms the ability of the code to describe the combustion
of suspended droplets at different pressure values. It highlights that large
droplets and high pressure conditions cannot be simulated accurately with
a 1D model without introducing errors on the burning rates and on the
flame geometry. However, simplified models can be corrected accounting for
buoyancy-driven flows using information from the detailed fluid dynamics
model following a multi-scale approach.

4.5. *Methanol combustion with water condensation*

The dynamics of methanol combustion is known to be affected by the condensation of water produced by the chemical reactions. This phenomenon is due to low boiling temperature of methanol (337.85 K, at atmospheric pressure) with respect to the water. During methanol combustion, the gas-liquid interface will always be at a temperature which is smaller (or equal) to methanol's boiling point, promoting condensation of vapor produced by the combustion process.

In this section, we simulate this phenomenon by utilizing the model's
ability to handle both evaporation and condensation of liquid mixtures. To
achieve this, we simulate the combustion of methanol droplets in microgravity
conditions using the simulation setup shown in Figure 1. The droplet initial
diameter is equal to D₀ = 0.5 mm, simulation are performed at atmospheric
pressure and with different initial water content inside the droplet: 0%, 10%,
and 20% in mass. We consider the same simulation parameters listed in
Section 4.4 for the microgravity simulations. For simplicity, we neglect the
non-ideality of the methanol-water mixture, in order to focus only on the
influence of water absorption on the burning rate. The implementation of
more complex non-ideal thermodynamics and mixing rules will be the subject
of future studies.

Figure 8 reports the dimensionless square diameter in time for the methanol droplets at different initial water contents. The reference case is the pure methanol combustion neglecting water absorption (from the previous section), which displays the maximum burning rate. Adding water, we observe a decreased burning rate for every case under investigation, which is proportional to the initial amount of water inside the droplet. The reference case has burning rate 0.61 mm²/s, which decreases to 0.58, 0.53, and 0.48 mm²/s respectively, with increasing water content. The slower consumption rate is caused by a decreased flame temperature, and changed liquid properties such as density, heat capacity, and latent heat of vaporization [78]. Figure 8 also shows that, after the ignition, the square diameter decreases linearly for each test case, until eventually the water accumulation extinguishes the flame, leading to a sudden change in the slope of the square diameter. The
<figure><img src="image_8.png" /><figcaption>Figure 7: Plots of the dimensionless square diameter and the flame diameter for the methanol combustion in microgravity and in normal gravity conditions. The lines are results from this model, while the triangles are the results from the 1D benchmark model by Cuoci et al. [73]. The ratio between the burning rate in normal gravity $k_{ng}$ and in microgravity conditions $k_{mg}$ measures the displacement between the two simulations.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Plots of the dimensionless square diameter for methanol combustion with water absorption, with different initial amounts of water in the droplet.</figcaption></figure>

extinction diameter is anticipated for the cases with higher initial amount of water, as observed by Marchese and Dryer [12].

This test case confirms the model's ability in predicting the effect of water absorption during the combustion of methanol droplets. Detailed modeling of this phenomenon and comparison with experimental data requires the introduction of liquid internal recirculation, which may be due to either the droplet generation procedure, and by thermocapillary and diffusocapillary Marangoni effects [78].

**5. Conclusions and future steps**

We have proposed and validated the first multiphase CFD model that combines: (i) interface-resolved multicomponent phase change, (ii) variable thermodynamic and transport properties, (iii) chemical reactions and radiation, and (iv) droplet suspension via surface tension.

1. The model couples a low-Mach VOF based phase change model with the direct integration of the combustion chemistry.

2. The intrinsic coupling, stiffness, and non-linearity of the reaction terms are correctly managed by an operator splitting approach using an implicit adaptive ODE solver.
3. By using quadtree adaptivity, we can maintain a large simulation domain, thereby avoiding boundary effects on flame dynamics while simultaneously reducing computational costs.

4. Validation using microgravity conditions demonstrate the convergence of the model toward the experimental data and the accurate prediction of the burning rate.

5. The model accurately captures experimental burning rates for varying initial droplet diameters and CO₂-enriched environments, explaining experimental trends and confirming the impact of non-luminous radiative heat losses.

6. Results of intermediate size ($D_0$ = 0.5 mm) droplet combustion in normal gravity conditions exhibit strong deviation from the microgravity conditions, which increase with increasing ambient pressure. These results demonstrate inaccuracies of the spherically-symmetric models for simulations of intermediate sized droplets in normal gravity conditions.

7. The multicomponent formulation can directly account for condensation of water produced by the combustion chemistry, and the resulting droplet extinction.

Future research will be devoted to relaxing the assumptions 4-6 in Section 2.6, in order to simulate methanol-water systems in detail, including multidimensional recirculation effects, and using a complex kinetic mechanism. The recent implementation of the Marangoni effect in the Basilisk framework [49] represents an important starting point to model systems characterized by strong liquid internal recirculation. The recent advancement in the numerical methods for VOF with complex solid boundaries and surface tension [79] will be exploited to include the solution of the conjugate heat transfer between the solid and the gas-liquid system, thus resolving the thermal effect of the solid fiber suspender. Most importantly, we will introduce strategies to reduce the computational time of the chemistry, in order to apply the model to complex simulation scenarios, such as interactions among droplet arrays, breakup and coalescence phenomena, and droplet-solid interactions. These simulations can be performed using this multidimensional model, in order to develop correlations for simplified droplet combustion models, thereby merging the detailed multiphase fluid dynamics with the computational efficiency of 1D models.
References

[1] H. Nomura, Y. Ujiie, H. J. Rath, J. Sato, M. Kono, Experimental study on high-pressure droplet evaporation using microgravity conditions, in: Symposium (International) on Combustion, volume 26, Elsevier, 1996, pp. 1267–1273.

[2] S. Kumagai, H. Isoda, Combustion of fuel droplets in a falling chamber, in: Symposium (International) on Combustion, volume 6, Elsevier, 1957, pp. 726–731.

[3] S. Okajima, S. Kumagai, Experimental studies on combustion of fuel droplets in flowing air under zero-and high-gravity conditions, in: Symposium (International) on Combustion, volume 19, Elsevier, 1982, pp. 1021–1027.

[4] H. Hara, S. Kumagai, Experimental investigation of free droplet combustion under microgravity, in: Symposium (International) on Combustion, volume 23, Elsevier, 1991, pp. 1605–1610.

[5] H. Hara, S. Kumagai, The effect of initial diameter on free droplet combustion with spherical flame, in: Symposium (International) on combustion, volume 25, Elsevier, 1994, pp. 423–430.

[6] S. Nakaya, D. Segawa, T. Kadota, Y. Nagashima, T. Furuta, Combustion behaviors of isolated n-decane and ethanol droplets in carbon dioxide-rich ambience under microgravity, Proceedings of the Combustion Institute 33 (2011) 2031–2038.

[7] Y. Liu, Y. Xu, C. Avedisian, M. Hicks, The effect of support fibers on micro-convection in droplet combustion experiments, Proceedings of the combustion institute 35 (2015) 1709–1716.

[8] A. K. Yadav, A. Chowdhury, A. Srivastava, Interferometric investigation of methanol droplet combustion in varying oxygen environments under normal gravity, International Journal of Heat and Mass Transfer 111 (2017) 871–883.

[9] D. L. Dietrich, J. B. Haggard Jr, F. L. Dryer, V. Nayagam, B. D. Shaw, F. A. Williams, Droplet combustion experiments in spacelab, in: Symposium (International) on Combustion, volume 26, Elsevier, 1996, pp. 1201–1207.
[10] D. L. Dietrich, V. Nayagam, M. C. Hicks, P. V. Ferkul, F. L. Dryer, T. Farouk, B. D. Shaw, H. K. Suh, M. Y. Choi, Y. C. Liu, et al., Droplet combustion experiments aboard the international space station, Microgravity Science and Technology 26 (2014) 65-76.

[11] S. Cho, R. A. Yetter, F. Dryer, A computer model for one-dimensional mass and energy transport in and around chemically reacting particles, including complex gas-phase chemistry, multicomponent molecular diffusion, surface evaporation, and heterogeneous reaction, Journal of Computational Physics 102 (1992) 160-179.

[12] A. Marchese, F. Dryer, The effect of non-luminous thermal radiation in microgravity droplet combustion, Combustion science and technology 124 (1997) 371-402.

[13] A. Cuoci, M. Mehl, G. Buzzi-Ferraris, T. Faravelli, D. Manca, E. Ranzi, Autoignition and burning rates of fuel droplets under microgravity, Combustion and Flame 143 (2005) 211-226.

[14] T. Farouk, F. Dryer, Microgravity droplet combustion: effect of tethering fiber on burning rate and flame structure, Combustion Theory and Modelling 15 (2011) 487-515.

[15] T. Farouk, M. Hicks, F. Dryer, Multistage oscillatory “cool flame” behavior for isolated alkane droplet combustion in elevated pressure microgravity condition, Proceedings of the Combustion Institute 35 (2015) 1701-1708.

[16] A. Cuoci, A. Frassoldati, T. Faravelli, E. Ranzi, Numerical modeling of auto-ignition of isolated fuel droplets in microgravity, Proceedings of the Combustion Institute 35 (2015) 1621-1627.

[17] A. Cuoci, C. T. Avedisian, J. D. Brunson, S. Guo, A. Dalili, Y. Wang, M. Mehl, A. Frassoldati, K. Seshadri, J. E. Dec, et al., Simulating combustion of a seven-component surrogate for a gasoline/ethanol blend including soot formation and comparison with experiments, Fuel 288 (2021) 119451.

[18] A. Stagni, L. Esclapez, P. Govindaraju, A. Cuoci, T. Faravelli, M. Ihme, The role of preferential evaporation on the ignition of multicomponent
fuels in a homogeneous spray/air mixture, Proceedings of the Combustion Institute 36 (2017) 2483-2491.

[19] D. Juric, G. Tryggvason, Computations of boiling flows, International journal of multiphase flow 24 (1998) 387-410.

[20] M. Irfan, M. Muradoglu, A front tracking method for direct numerical simulation of evaporation process in a multiphase system, Journal of Computational Physics 337 (2017) 132-153.

[21] S. Tanguy, T. Ménard, A. Berlemont, A level set method for vaporizing two-phase flows, Journal of Computational Physics 221 (2007) 837-853.

[22] L. R. Villegas, R. Alis, M. Lepilliez, S. Tanguy, A ghost fluid/level set method for boiling flows and liquid evaporation: application to the leidenfrost effect, Journal of Computational Physics 316 (2016) 789-813.

[23] K. Luo, C. Shao, M. Chai, J. Fan, Level set method for atomization and evaporation simulations, Progress in Energy and Combustion Science 73 (2019) 65-94.

[24] S. W. Welch, J. Wilson, A volume of fluid based method for fluid flows with phase change, Journal of computational physics 160 (2000) 662-682.

[25] N. Scapin, P. Costa, L. Brandt, A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows, Journal of Computational Physics 407 (2020) 109251.

[26] B. Boyd, Y. Ling, A consistent volume-of-fluid approach for direct numerical simulation of the aerodynamic breakup of a vaporizing drop, Computers & Fluids 254 (2023) 105807.

[27] E. Cipriano, A. E. Saufi, A. Frassoldati, T. Faravelli, S. Popinet, A. Cuoci, Multicomponent droplet evaporation in a geometric volume-of-fluid framework, Journal of Computational Physics 507 (2024) 112955.

[28] T. Long, J. Pan, S. Zaleski, An edge-based interface tracking (ebit) method for multiphase flows with phase change, Journal of Computational Physics 513 (2024) 113159.
[29] S. T. Kumaran, R. Kumar, B. Premachandran, A multi-directional advection based moment of fluid method for phase change problems, Journal of Computational Physics 508 (2024) 113011.

[30] C. Shao, K. Luo, M. Chai, H. Wang, J. Fan, A computational framework for interface-resolved dns of simultaneous atomization, evaporation and combustion, Journal of Computational Physics 371 (2018) 751–778.

[31] M. Irfan, M. Muradoglu, A front tracking method for particle-resolved simulation of evaporation and combustion of a fuel droplet, Computers & Fluids 174 (2018) 283–299.

[32] A. E. Saufi, A. Frassoldati, T. Faravelli, A. Cuoci, Interface-resolved simulation of the evaporation and combustion of a fuel droplet suspended in normal gravity, Fuel 287 (2021) 119413.

[33] E. Cipriano, A. Frassoldati, T. Faravelli, S. Popinet, A. Cuoci, A low-mach volume-of-fluid model for the evaporation of suspended droplets in buoyancy-driven flows, International Journal of Heat and Mass Transfer 234 (2024) 126115.

[34] D. A. Drew, Mathematical modeling of two-phase flow, Annual Review of Fluid Mechanics 15 (1983) 261–291.

[35] N. Ly, M. Ihme, A regularized-interface method as a unified formulation for simulations of high-pressure multiphase flows, Journal of Computational Physics 518 (2024) 113310.

[36] R. G. Rehm, H. R. Baum, The equations of motion for thermally driven, buoyant flows, Journal of research of the National Bureau of Standards 83 (1978) 297.

[37] R. B. Pember, L. H. Howell, J. B. Bell, P. Colella, W. Y. Crutchfield, W. Fiveland, J. Jessee, An adaptive projection method for unsteady, low-mach number combustion, Combustion Science and Technology 140 (1998) 123–168.

[38] G. Mialhe, S. Tanguy, L. Tranier, E.-R. Popescu, D. Legendre, An extended model for the direct numerical simulation of droplet evaporation. influence of the marangoni convection on leidenfrost droplet, Journal of Computational Physics 491 (2023) 112366.
[39] T. Coffee, J. Heimerl, Transport algorithms for premixed, laminar steady-state flames, Combustion and Flame 43 (1981) 273-289.

[40] J. O. Hirschfelder, C. F. Curtiss, R. B. Bird, The molecular theory of gases and liquids, John Wiley & Sons, 1964.

[41] J. Leffler, H. Cullinan Jr, Variation of liquid diffusion coefficients with composition. dilute ternary systems, Industrial & Engineering Chemistry Fundamentals 9 (1970) 88-93.

[42] T. Poinsot, D. Veynante, Theoretical and numerical combustion, RT Edwards, Inc., 2005.

[43] A. Cuoci, A. Frassoldati, T. Faravelli, E. Ranzi, Numerical modeling of laminar flames with detailed kinetics based on the operator-splitting method, Energy & fuels 27 (2013) 7730-7753.

[44] M. F. Modest, S. Mazumder, Radiative heat transfer, Academic press, 2021.

[45] A. Cuoci, A. Frassoldati, T. Faravelli, E. Ranzi, A computational tool for the detailed kinetic modeling of laminar flames: Application to c2h4/ch4 coflow flames, Combustion and Flame 160 (2013) 870-886.

[46] B. Shaw, V. Dee, Influence of the gas phase inert on reduced-gravity combustion of decane/hexadecane droplets, Microgravity-Science and Technology 17 (2005) 26-34.

[47] R. B. Bird, Transport phenomena, Appl. Mech. Rev. 55 (2002) R1-R4.

[48] M. S. Dodd, D. Mohaddes, A. Ferrante, M. Ihme, Analysis of droplet evaporation in isotropic turbulence through droplet-resolved dns, International Journal of Heat and Mass Transfer 172 (2021) 121157.

[49] P. K. Farsoiya, S. Popinet, H. A. Stone, L. Deike, Coupled volume of fluid and phase field method for direct numerical simulation of insoluble surfactant-laden interfacial flows and application to rising bubbles, Physical Review Fluids 9 (2024) 094004.

[50] G. Standart, The mass, momentum and energy equations for heterogeneous flow systems, Chemical Engineering Science 19 (1964) 227-236.
[51] A. Cuoci, A. Frassoldati, T. Faravelli, E. Ranzi, Opensmoke++: An object-oriented framework for the numerical modeling of reactive systems with detailed kinetic mechanisms, Computer Physics Communications 192 (2015) 237-264.

[52] R. C. Reid, J. M. Prausnitz, B. E. Poling, The properties of gases and liquids, McGraw Hill Book Co., New York, NY, 1987.

[53] C. K. Westbrook, F. L. Dryer, Simplified reaction mechanisms for the oxidation of hydrocarbon fuels in flames, Combustion science and technology 27 (1981) 31-43.

[54] R. Kee, F. Rupley, J. Miller, Chemkin: The chemkin thermodynamic database. rep. sand87-8215, Sandia National Laboratories, Livermore (CA), USA (1987).

[55] J. H. Ferziger, M. Perić, R. L. Street, Computational methods for fluid dynamics, 4 ed., springer, 2019.

[56] R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Annual review of fluid mechanics 31 (1999) 567-603.

[57] G. D. Weymouth, D. K.-P. Yue, Conservative volume-of-fluid method for free-surface simulations on cartesian-grids, Journal of Computational Physics 229 (2010) 2853-2865.

[58] S. Fleckenstein, D. Bothe, A volume-of-fluid-based numerical method for multi-component mass transfer with local volume changes, Journal of Computational Physics 301 (2015) 35-58.

[59] G. Gennari, R. Jefferson-Loveday, S. J. Pickering, A phase-change model for diffusion-driven mass transfer problems in incompressible two-phase flows, Chemical Engineering Science 259 (2022) 117791.

[60] I. Awasthi, G. Gogos, T. Sundararajan, Effects of size on combustion of isolated methanol droplets, Combustion and flame 160 (2013) 1789-1802.

[61] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct Numerical Simulations of Gas-Liquid Multiphase Flows, Cambridge university press, 2011.
[62] S. Popinet, Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries, J. Comput. Phys. 190 (2003) 572-600.

[63] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, J. Comput. Phys. 228 (2009) 5838-5866.

[64] J. B. Bell, P. Colella, H. M. Glaz, A second-order projection method for the incompressible Navier-Stokes equations, J. Comput. Phys. 85 (1989) 257-283.

[65] S. Popinet, A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations, J. Comput. Phys. 302 (2015) 336-358.

[66] G. H. Yeoh, J. Tu, Computational techniques for multiphase flows, Butterworth-Heinemann, 2019.

[67] E. S. Oran, J. P. Boris, J. P. Boris, Numerical simulation of reactive flow, volume 2, Citeseer, 2001.

[68] A. Cuoci, Numerical modeling of reacting systems with detailed kinetic mechanisms, in: Computer Aided Chemical Engineering, volume 45, Elsevier, 2019, pp. 675-721.

[69] G. M. Goldin, Z. Ren, S. Zahirovic, A cell agglomeration algorithm for accelerating detailed chemistry in cfd, Combustion Theory and Modelling 13 (2009) 721-739.

[70] U. Maas, S. B. Pope, Implementation of simplified chemical kinetics based on intrinsic low-dimensional manifolds, in: Symposium (international) on combustion, volume 24, Elsevier, 1992, pp. 103-112.

[71] S. B. Pope, Computationally efficient implementation of combustion chemistry using in situ adaptive tabulation, Combustion Theory and Modelling (1997).

[72] S. Popinet, Basilisk C reference manual., http://basilisk.fr/Basilisk%20C, 2013.
[73] A. Cuoci, E. Cipriano, A. E. Saufi, A. Frassoldati, A numerical framework for modeling evaporation and combustion of isolated, spherically-symmetric, multi-component fuel droplets, Journal of Computational Science 100 (2024) 10–100.

[74] S. Popinet, Numerical models of surface tension, Annu. Rev. Fluid Mech. 50 (2018) 49–75.

[75] J. A. Van Hooft, S. Popinet, C. C. Van Heerwaarden, S. J. Van der Linden, S. R. De Roode, B. J. Van de Wiel, Towards adaptive grids for atmospheric boundary-layer simulations, Boundary-layer meteorology 167 (2018) 421–443.

[76] A. Cuoci, A. E. Saufi, A. Frassoldati, D. L. Dietrich, F. A. Williams, T. Faravelli, Flame extinction and low-temperature combustion of isolated fuel droplets of n-alkanes, Proceedings of the Combustion Institute 36 (2017) 2531–2539.

[77] A. Kazakov, J. Conley, F. L. Dryer, Detailed modeling of an isolated, ethanol droplet combustion under microgravity conditions, Combustion and Flame 134 (2003) 301–314.

[78] A. J. Marchese, F. L. Dryer, The effect of liquid mass transport on the combustion and extinction of bicomponent droplets of methanol and water, Combustion and Flame 105 (1996) 104–122.

[79] M. Tavares, C. Josserand, A. Limare, J. M. Lopez-Herrera, S. Popinet, A coupled vof/embedded boundary method to model two-phase flows on arbitrary solid surfaces, Computers & Fluids 278 (2024) 106317.