# Numerical study of droplet/flame interaction using a diffuse interface method

Benoît Péden, Pierre Boivin, Nicolas Odier

► To cite this version:

Benoît Péden, Pierre Boivin, Nicolas Odier. Numerical study of droplet/flame interaction using a diffuse interface method. 5th International Seminar on Non-Ideal Compressible Fluid Dynamics for Propulsion & Power, Ecole Centrale de Lyon, Oct 2024, Lyon, France. 10.1007/978-3-031-92695-2_12. hal-05123420

HAL Id: hal-05123420

https://hal.science/hal-05123420v1

Submitted on 4 Dec 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
Numerical study of droplet/flame interaction using a
diffuse interface method

Benoît Péden¹, Pierre Boivin², and Nicolas Odier³

¹CECI, Université de Toulouse, CERFACS/CNRS, 31100, Toulouse, France

²Aix Marseille Univ, CNRS, Centrale Med, M2P2, Marseille, France

³CERFACS, 42, Avenue Gaspard Coriolis, 31057 Toulouse, France

Abstract

The interaction between a fuel droplet and a premixed laminar flame is investigated in this study using a multi-fluid Diffuse-Interface approach. The influence of two parameters on the reacting zone is examined: the equivalence ratio of the premixed laminar flame and the relative velocity between the droplet and the stabilized flame front. After validating the equation of state's extension to temperature-dependent heat capacities, direct numerical simulations of droplet-flame interaction are performed. The results are qualitatively analyzed, revealing a significant influence of the operating conditions.

Keywords: Two-phase flows, Combustion, Evaporation, Diffuse Interface Method

# 1 Introduction

The combustion of sprays is characteristic of aeronautical gas turbines, leading to intricate spray/combustion interactions. Liquid injection and atomization processes are crucial as they determine the distribution of fuel vapor within a combustion chamber [1, 2, 3], significantly impacting the flame regime, burner efficiency, and pollutant emissions [4]. In such combustion chambers, liquid droplets of high inertia often cross the flame front [5]. In order to understand the interactions between an evaporating fuel droplet and a premixed flame, numerical simulations are mandatory.
Several approaches exist such as Direct Numerical Simulations (DNS) coupled with Discrete Particle Simulation (DPS) [1, 5, 6]. The later method does not explicitly resolve the liquid phase. Several approaches have been proposed in the literature to conduct numerical simulation of fully eulerian two-phase flows. The method used in this study falls within this category and in particular in Multi-Fluid methods [7, 8].

The paper is structured as follows, (i) the methodology employed is described in Section 2 with the system of conservation equation, the thermodynamic closure and the phase-transition solver. (ii) The results are given in Section 3 where the method is validated on a canonical one dimensional flame and the study of droplet/flame interaction follows.

# 2 Methodology

This section presents the models for the simulation of compressible, reactive multi-component flows in the AVBP solver [9, 10] from CERFACS. Initially, the system of conservation equations is introduced. Subsequently, the non-ideal thermodynamic closure and its extension for reactive flows are explained. Finally, the employed phase-transition solver is described.

## 2.1 The 4-equation Homogeneous Relaxation Method

The 4-equation, or Homogeneous Relaxation Method, model [11] is employed to simulate compressible reactive and multicomponent flows. This model assumes mechanical and thermal equilibrium between phases in a computational cell, i.e., all components share the same pressure, velocity, and temperature. Consequently, a single set of conservation laws for mass, momentum, and energy is solved throughout the numerical domain, regardless of the fluid state (liquid or gas). This model, which belongs to the multi-fluid method [7, 8], treats the interfacial area as a numerically diffused zone. Under these assumptions, the 4-equation model considers a mono-component liquid ($k = 1$), its associated vapor ($k = 2$), and a multi-component, non-condensable gas mixture ($k = [3, ..., N]$), with $N$ being the total number of species. The equilibrium conditions
lead to the following variable definitions,

$$ \left\{ \begin{aligned} T &= T_k, & \forall k, \\ p &= p_k, & \forall k, \\ v &= 1/\rho = \sum_{k=1}^{N} Y_k v_k, \\ e &= \sum_{k=1}^{N} Y_k e_k, \end{aligned} \right. \qquad (1) $$

where $Y_k$ represents the mass fraction of species $k$, $T$ is the temperature, $p$ is the pressure, $\rho$ is the density, and $e$ is the internal energy. The resulting four equations, with non-Euler right-hand side terms, are

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0, \qquad (2) $$

$$ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I}) = \mathbf{F}_{\sigma} + \nabla \cdot (\mu \boldsymbol{\tau}), \qquad (3) $$

$$ \frac{\partial \rho E}{\partial t} + \nabla \cdot ((\rho E + p)\mathbf{u}) = \mathbf{u} \cdot \mathbf{F}_{\sigma} + \nabla \cdot (\mu \boldsymbol{\tau} \cdot \mathbf{u}) - \nabla \cdot \mathbf{q}, \qquad (4) $$

$$ \frac{\partial \rho Y_k}{\partial t} + \nabla \cdot (\rho Y_k (\mathbf{u} + \mathbf{v}_k)) = \dot{\omega}_{c,k} + \dot{\omega}_{p,k}, \qquad (5) $$

where $\mathbf{u}$ is the velocity vector and $E$ is the total mixture energy ($E = e + \mathbf{u}^2/2$). The right-hand side non-Euler term expressions are described in [12, 13]. The $\mathbf{F}_{\sigma}$ term represents the surface tension force, modelled using the Continuum Surface Force of Brackbill et al. [14], and reads $\mathbf{F}_{\sigma} = \sigma \kappa \nabla z_{\ell}$. The $\sigma$ parameter is the surface tension coefficient, $\kappa$ is the curvature, computed using a color function [14], the liquid volume fraction here $z_{\ell}$. The average mixture viscosity and heat conductivity are based on algebraic averages between the phases [12]. The source terms of Eq. (5) are denoted $\dot{\omega}_{c,k}$ for chemical reactions and $\dot{\omega}_{p,k}$ for phase transitions, detailed in Section 2.3. The diffusion effect for each component $k$ is accounted for by introducing the diffusion velocity $\mathbf{v}_k$.

## 2.2 The extended Noble-Abel Stiffened-Gas Equation of State

Each component follows the Noble-Abel Stiffened-Gas (NASG) equation of state (EoS) [15]. This EoS is an improvement over the classic Stiffened-Gas EoS, accounting for repulsive molecular effects.

For a component $k$, the NASG EoS with temperature-variable heat capacities [16] is
expressed as:

$$ v_k(p, T) = \frac{(C_{p,k}(T) - C_{v,k}(T))T}{p + p_{\infty,k}} + b_k, \quad (6a) $$

$$ h_k(p, T) = \int C_{p,k}(T) dT + b_k p + q_k, \quad (6b) $$

$$ e_k(p, T) = \int C_{v,k}(T) dT + \frac{(C_{p,k}(T) - C_{v,k}(T))Tp_{\infty,k}}{p + p_{\infty,k}} + q_k, \quad (6c) $$

where $C_{p,k}$ and $C_{v,k}$ are the constant pressure and constant volume heat capacities, respectively; $b_k$ is the covolume; and $p_{\infty,k}$ is a correction term to account for attraction (both terms are present only in the liquid phase, i.e., $k = 1$). The reference energy is denoted $q_k$. All these EoS parameters are determined by the thermodynamic properties of the fluid [16], not using the saturation curves [15, 17].

When using temperature-variable heat capacities, the mixture density $\rho$ and internal energy $e$ can be derived from the intensive variables $p$ and $T$. However, the $\int C_{v,k}(T) dT$ term in the energy equation of System (6) precludes explicit expressions for $p$ and $T$ from $\rho$ and $e$. To address this issue, Boivin et al. [16] proposed a hybrid closure, stating that the constant heat capacity assumption only fails at high temperatures, where the liquid phase is no longer present. Heat capacities are assumed constant below a reference temperature ($T_{ref} = 700$ K), a temperature sufficiently high to ensure a mono-phase gaseous flow (far from the liquid saturation temperature). The hybrid closure posits that the EoS is NASG with constant heat capacities below $T_{ref}$ and the classic perfect gas (PG) EoS above. The NASG EoS coefficients must be rederived to ensure $C_1$ continuity at $T_{ref}$, in order to use the Janaf tables implemented in the AVBP solver. Originally, the NASG coefficients are derived from the NIST database [18], but the continuity of the enthalpy at $T_{ref}$ is not guaranteed due to the choice of reference temperature, which is set to $T = 0$ K in the AVBP solver. The NASG parameters used in this paper are given in Table 1. The density $\rho$ and the enthalpy $h$ obtained with the NASG EoS coupled with the previously described methodology is plotted against NIST database curves [18] in Fig. 1, showing an accurate prediction in both liquid and gaseous states. The temperature $T_{ref} = 700$ K is thus justified as it is far from the saturation point (visible as a jump in thermodynamic quantities).
<table><caption>**Table 1.** NASG EoS parameters employed in this study, with $\gamma = C_p/C_v$.</caption><thead><tr><th>Species</th><th>C<sub>p</sub> [J.kg<sup>-1</sup>.K<sup>-1</sup>]</th><th>γ</th><th>q [J.kg<sup>-1</sup>]</th><th>b [m<sup>3</sup>.kg<sup>-1</sup>]</th><th>p<sub>∞</sub> [Pa]</th></tr></thead><tbody><tr><td>O<sub>2</sub></td><td>1005</td><td>1.349</td><td>-41 783</td><td>0</td><td>0</td></tr><tr><td>H<sub>2</sub>O</td><td>2023</td><td>1.295</td><td>-78 766</td><td>0</td><td>0</td></tr><tr><td>N<sub>2</sub></td><td>1078</td><td>1.380</td><td>-18 777</td><td>0</td><td>0</td></tr><tr><td>KERO_LUCHE</td><td>2771</td><td>1.022</td><td>-754 009</td><td>0</td><td>0</td></tr><tr><td>CO</td><td>1090</td><td>1.374</td><td>-24 166</td><td>0</td><td>0</td></tr><tr><td>CO<sub>2</sub></td><td>1079</td><td>1.212</td><td>-139 174</td><td>0</td><td>0</td></tr><tr><td>C<sub>10</sub>H<sub>22,ℓ</sub></td><td>2386</td><td>1.216</td><td>-0.855 × 10<sup>6</sup></td><td>9.4 × 10<sup>-4</sup></td><td>2.875 × 10<sup>8</sup></td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Validation of the NASG parameters against the NIST database [18]. Profiles of density and enthalpy against temperature at $p = 1$ bar.</figcaption></figure>
## 2.3 Phase transition model

The phase transition occurs between the liquid phase ($k = 1$) and its vapor phase ($k = 2$). The source terms apply exclusively to species 1 and 2 and are defined as:

$$ \begin{cases} \dot{\omega}_{p,1} = \rho\nu(g_2 - g_1), \\ \dot{\omega}_{p,2} = -\rho\nu(g_2 - g_1), \end{cases} \qquad (7) $$

where $g_k$ denotes the Gibbs free energy of phase $k$, and $\nu$ is a function dependent on the specific interfacial area, pressure, and temperature. This parameter represents the relaxation rate of the mixture toward thermodynamic equilibrium. Within this framework, an instantaneous relaxation of Gibbs free energies is assumed [19]. Finite-rate models are also documented in the literature [20, 21]. Instead of directly solving Eq. (7), the fast phase transition solver proposed by Chiapolino et al. [19] is utilized.

# 3 Results

## 3.1 Validation of the implemented hybrid thermodynamic closure – Kerosene surrogate laminar pre-mixed flame

A one-dimensional, freely propagating pre-mixed flame is examined using fuel closely related to the targeted aeronautical engines applications. The chemical mechanism employed is the well-known two-step kerosene BFER model [22], which includes 6 + 1 species, with the liquid $C_{10}H_{22}$ added. The initial conditions are specified in Table 2, at atmospheric pressure. A small amount of liquid kerosene $C_{10}H_{22,\ell}$ is uniformly distributed throughout the domain. The species KERO\_LUCHE represents the vapor phase associated with the liquid $C_{10}H_{22}$.

<table><caption>**Table 2.** Initial conditions for the one-dimensional free-propagating $C_{10}H_{22}$–Air pre-mixed flame.</caption><thead><tr><th>Variable</th><th><i>T</i> [K]</th><th><i>u</i> [m.s<sup>-1</sup>]</th><th><i>Y</i><sub>C<sub>10</sub>H<sub>22,ℓ</sub></sub></th><th><i>Y</i><sub>KERO\_LUCHE</sub></th><th><i>Y</i><sub>CO<sub>2</sub></sub></th><th><i>Y</i><sub>CO</sub></th><th><i>Y</i><sub>O<sub>2</sub></sub></th><th><i>Y</i><sub>H<sub>2</sub>O</sub></th><th><i>Y</i><sub>N<sub>2</sub></sub></th></tr></thead><tbody><tr><td>Fresh gases</td><td>298</td><td>0.319</td><td>1 &times; 10<sup>-8</sup></td><td>5.775 &times; 10<sup>-2</sup></td><td>0</td><td>0</td><td>2.20 &times; 10<sup>-1</sup></td><td>0</td><td>1 &minus; &sum;<sup>N&minus;1</sup><sub>k=1</sub> <i>Y</i><sub>k</sub></td></tr><tr><td>Burnt gases</td><td>2217</td><td>2.5</td><td>1 &times; 10<sup>-8</sup></td><td>0</td><td>1.733 &times; 10<sup>-1</sup></td><td>5.015 &times; 10<sup>-3</sup></td><td>2.482 &times; 10<sup>-2</sup></td><td>7.416 &times; 10<sup>-2</sup></td><td>1 &minus; &sum;<sup>N&minus;1</sup><sub>k=1</sub> <i>Y</i><sub>k</sub></td></tr></tbody></table>

The resulting solution is depicted in Fig. 2 and is compared with results obtained from the AVBP solver using the perfect gas EoS. The two sets of simulation results are
superimposed, and no jump on the temperature/heat realease curve appears at $T_{\text{ref}}$, thereby validating the implementation.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Free-propagating kerosene flame. Temperature, species and heat release profiles for both equations of state: NASG without and PG with markers.</figcaption></figure>

## 3.2 Droplet/laminar flame front interaction

In this section, the interaction of an evaporating fuel droplet with a laminar premixed flame is investigated. The numerical framework, validated in the previous subsection is employed. The results are analyzed qualitatively. The laminar premixed flame equivalence ratio, as well as the relative velocity between the fuel droplet and the flame front, ($u^* = u_d/S_l$ where $u_d$ is the droplet initial velocity and $S_l$ is the laminar flame speed), are the two parameters of this study.

### 3.2.1 Configuration

The computational domain, as depicted in Fig. 3, comprises a [4.5 mm × 2.5 mm] two-dimensional region constructed with 537 × 10³ triangular elements. The boundary conditions, as detailed in [23], consist of a non-reflecting inlet on the left, a non-reflecting outlet on the right, both extended to NASG EoS [13], and two symmetries on the top and bottom boundaries. The chemistry is described using the 2-step mechanism, validated above, and reported in [22], encompassing 6 + 1 species and 2 reactions. A third-order Galerkin scheme in space and time, TTGC [24], is applied. Additionally, the Localized Artificial Diffusion (LAD) method, as described in [25], is employed. The surface tension coefficient for n-decane is set at $\sigma = 24 \times 10^{-3} \text{ N.m}^{-1}$. The gaseous mixture is injected at the flame front velocity. A liquid n-decane droplet of diameter $d = 80 \text{ µm}$ is placed in front of the reacting zone, after the establishment of the purely gaseous premixed flame. The composition of the gaseous premixed flame, for each
equivalence ratio, is determined considering the thermodynamic equilibrium condition with Cantera [26] for this chemical mechanism.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Schematic of the configuration.</figcaption></figure>

### 3.2.2 Influence of the premixed flame equivalence ratio

A droplet with a velocity equal to 10 times the premixed flame velocity is considered for two equivalence ratios, $\phi = 0.7$ and $\phi = 1.3$. For each initial equivalence ratio, the temperature and heat release solutions are presented in Fig. 4.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Influence of the equivalence ratio of the premixed flame. Left: $\phi = 0.7$, right: $\phi = 1.3$. The top half of each figure represents the temperature field and the bottom half represents the heat release field. For the case initially at $\phi = 0.7$, iso-contours of $\phi = 0.7$ and $\phi = 1$ are shown. The red iso-contour represents the liquid volume fraction ($z_{\ell} = 0.3$).</figcaption></figure>

The behavior differs between the two cases: for $\phi = 0.7$, a non-zero mass fraction of $O_2$ ($Y_{O_2}$) is observed in the burnt gases zone. Consequently, the interaction of the
evaporating fuel with this zone generates a second flame front, characterized by a lower
heat release, located at the $\phi = 1$ iso-contour (Fig. 4 bottom left). This topology is
reminiscent of diffusion flames. Conversely, for $\phi = 1.3$, no oxidizer is present in the
burnt gases zone, causing the reaction to cease. It is also noteworthy that the passage
of the evaporating droplet significantly affects the flame front for the $\phi = 0.7$ case;
the flame front is deformed and exhibits a higher heat release compared to the purely
gaseous premixed case ($\dot{Q} = 2 \times 10^9$ W.m$^{-3}.s^{-1}$). The reaction in the burnt gases is
evident in Fig. 5 (left): the O$_2$ mass fraction in the burnt gases, between the $\phi = 1$
iso-contour and the droplet, is zero. Additionally, the Takeno flame index (blue $\rightarrow$
non-premixed) evidences this diffusion mode in the burnt gases. Finally, the iso-contours
of liquid volume fraction $z_l = 0.3$ (in red) in Fig. 5 exhibit a higher evaporation rate for
the $\phi = 0.7$ case as the contour zone is smaller.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Influence of the equivalence ratio of the premixed flame. Left: $\phi = 0.7$, right: $\phi = 1.3$. The top half of each figure represents the vapor fuel mass fraction and the bottom half represents the oxidizer mass fraction. For the case initially at $\phi = 0.7$, iso-contours of $\phi = 0.7$ and $\phi = 1$ are shown. The red iso-contour represents the liquid volume fraction ($z_l = 0.3$).</figcaption></figure>

### 3.2.3 Influence of the relative velocity of the droplet to the flame front

For the case $\phi = 0.7$, a droplet with the same velocity as the flame front ($u^* = 1$) is
placed in front of the reacting zone. The global topologies of the temperature and the
heat release field are retrieved, the temperature and heat release solutions are presented in Fig. 6. The same observation regarding a second heat release in the burnt gases zone can be made: a high temperature zone is exhibited on the $\phi = 1$ iso-contour. However, as pointed out by the liquid volume fraction iso-contour ($z_{\ell} = 0.3$), the droplet does not cross the flame front due to its initial velocity.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Influence of the relative velocity of the droplet to the premixed flame for $\phi = 0.7$. (a) The top half represents the temperature field and the bottom half represents the heat release field. (b) The top half represents the Takeno flame index and the bottom half represents the $O_2$ species mass fraction field. Iso-contours of $\phi = 0.7$ and $\phi = 1$ are shown. The red iso-contour represents the liquid volume fraction ($z_{\ell} = 0.3$).</figcaption></figure>

Two instantaneous fields of liquid volume fraction ($z_{\ell}$) are shown in Fig. 7. An iso-contour of $\dot{Q} = 2 \times 10^9$ W.m$^{-3}$ is also shown, which corresponds to the initial premixed flame heat release. It means that for the case at $u^* = 1$, the droplet is almost completely vaporized before crossing the flame front. However, for the $u^* = 10$ case, the droplet crosses the flame front before being completely vaporized. For the latter case, the flame front retrieves the initial premixed value while the droplet finishes vaporizing in the burnt gases.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Influence of the relative velocity of the droplet to the premixed flame for $\phi = 0.7$. Liquid volume fraction $z_\ell$ fields with $\dot{Q} = 2 \times 10^9 \text{ W.m}^{-3}$ iso-contours (white lines).</figcaption></figure>

## 4 Conclusions and future works

The extension of the NASG equation of state (EoS) [15] to temperature-dependent heat capacities [16] has been implemented and validated within the AVBP solver. The interactions between droplets and flames have demonstrated significantly different behaviors when varying the following parameters: the premixed flame equivalence ratio $\phi$ and the relative velocity $u^*$ between the fuel droplet and the initial premixed flame front. Physical insights were obtained: when the premixed flame equivalence ratio $\phi$ is less than 1, a diffusion flame is present in the burnt gases zone. This phenomenon is naturally absent when $\phi$ is greater than 1. For droplets with an initial velocity greater than the flame front, the droplet traverses the reacting zone and continues to evaporate in the burnt gases zone; the reacting zone then relaxes towards the initial purely gaseous premixed flame configuration. Conversely, when the droplet's velocity matches that of the flame front, the droplet nearly completely vaporizes before crossing the reaction zone. These encouraging results indicate that the method captures physical phenomena accurately and can be applied to future, more complex cases and more detailed DNS studies.

## 5 Acknowledgments

This work was supported by the French National Research Agency through the project MADIRAN (grant ANR-21-CE50-0043) and benefited from GENCI-CINES HPC resources (grant 2023 - A0152B10157).
References

[1] J. Réveillon and L. Vervisch. "Analysis of weakly turbulent dilute-spray flames and spray combustion regimes". In: *Journal of Fluid Mechanics* 537.1 (2005), pp. 317-347.

[2] K. Luo, H. Pitsch, M. Pai, and O. Desjardins. "Direct numerical simulations and analysis of three-dimensional n-heptane spray flames in a model swirl combustor". In: *Proceedings of the Combustion Institute* 33.2 (2011), pp. 2143-2152.

[3] A. L. Sánchez, J. Urzay, and A. Liñán. "The role of separation of scales in the description of spray combustion". In: *Proceedings of the Combustion Institute* 35.2 (2015), pp. 1549-1577.

[4] C. Beck, R. Koch, and H.-J. Bauer. "Identification of droplet burning modes in lean, partially prevaporized swirl-stabilized spray flames". In: *Proceedings of the Combustion Institute* 32.2 (2009), pp. 2195-2203.

[5] F. Shum-Kivan, J. M. Santiago, A. Verdier, E. Riber, B. Renou, G. Cabot, and B. Cuenot. "Experimental and numerical analysis of a turbulent spray flame structure". In: *Proceedings of the Combustion Institute* 36.2 (2017), pp. 2567-2575.

[6] B. Rochette, E. Riber, and B. Cuenot. "Effect of non-zero relative velocity on the flame speed of two-phase laminar flames". In: *Proceedings of the Combustion Institute* 37.3 (2019), pp. 3393-3400.

[7] R. Saurel and R. Abgrall. "A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows". In: *Journal of Computational Physics* 150.2 (1999), pp. 425-467.

[8] A. Kapila, R. Menikoff, J. Bdzil, S. Son, and D. Stewart. "Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations". In: *Physics of Fluids* 13.10 (2001), pp. 3002-3024.

[9] T. Schönfeld and M. Rudgyard. "Steady and Unsteady Flow Simulations Using the Hybrid Flow Solver AVBP". In: *AIAA Journal* 37.11 (1999), pp. 1378-1385.

[10] L. Y. Gicquel, N. Gourdain, J.-F. Boussuge, H. Deniau, G. Staffelbach, P. Wolf, and T. Poinsot. "High performance parallel computing of flows in complex geometries". In: *Comptes Rendus Mécanique* 339.2-3 (2011), pp. 104-124.

[11] R. Saurel, P. Boivin, and O. Le Métayer. "A general formulation for cavitating, boiling and evaporating flows". In: *Computers & Fluids* 128 (2016), pp. 53-64.
[12] X. Deng and P. Boivin. "Diffuse interface modelling of reactive multi-phase flows applied to a sub-critical cryogenic jet". In: *Applied Mathematical Modelling* 84 (2020), pp. 405-424.

[13] B. Péden, J. Carmona, P. Boivin, T. Schmitt, B. Cuenot, and N. Odier. "Numerical assessment of Diffuse-Interface method for air-assisted liquid sheet simulation". In: *Computers & Fluids* 266 (2023), p. 106022.

[14] J. Brackbill, D. Kothe, and C. Zemach. "A continuum method for modeling surface tension". In: *Journal of Computational Physics* 100.2 (1992), pp. 335-354.

[15] O. Le Métayer and R. Saurel. "The Noble-Abel Stiffened-Gas equation of state". In: *Physics of Fluids* 28.4 (2016), p. 046102.

[16] P. Boivin, M. A. Cannac, and O. Le Métayer. "A thermodynamic closure for the simulation of multiphase reactive flows". In: *International Journal of Thermal Sciences* 137 (2019), pp. 640-649.

[17] O. Le Métayer, J. Massoni, and R. Saurel. "Elaboration des lois d'état d'un liquide et de sa vapeur pour les modèles d'écoulements diphasiques". In: *International Journal of Thermal Sciences* 43.3 (2004), pp. 265-276.

[18] P. J. Linstrom and W. G. Mallard. "The NIST Chemistry WebBook: A Chemical Data Resource on the Internet". In: *Journal of Chemical & Engineering Data* 46.5 (2001), pp. 1059-1063.

[19] A. Chiapolino, P. Boivin, and R. Saurel. "A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows". In: *Computers & Fluids* 150 (2017), pp. 31-45.

[20] D. Furfaro and R. Saurel. "Modeling droplet phase change in the presence of a multi-component gas mixture". In: *Applied Mathematics and Computation* 272 (2016), pp. 518-541.

[21] M. Pelanti. "Arbitrary-rate relaxation techniques for the numerical modeling of compressible two-phase flows with heat and mass transfer". In: *International Journal of Multiphase Flow* 153 (2022), p. 104097.

[22] B. Franzelli, E. Riber, M. Sanjosé, and T. Poinsot. "A two-step chemical scheme for kerosene-air premixed flames". In: *Combustion and Flame* 157.7 (2010), pp. 1364-1373.

[23] T. J. Poinsot and S. K. Lele. "Boundary conditions for direct simulations of compressible viscous flows". In: *Journal of Computational Physics* 101.1 (1992), pp. 104-129.
[24] O. Colin and M. Rudgyard. "Development of High-Order Taylor-Galerkin Schemes for LES". In: *Journal of Computational Physics* 162.2 (2000), pp. 338-371.

[25] T. Schmitt. "Large-Eddy Simulations of the Mascotte Test Cases Operating at Supercritical Pressure". In: *Flow, Turbulence and Combustion* 105 (2020), pp. 159-189.

[26] D. G. Goodwin, H. K. Moffat, and R. L. Speth. *Cantera - An Object-oriented Software Toolkit for Chemical Kinetics, Thermodynamics, and Transport Processes*. http://www.cantera.org. 2017.