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

HAL
open science

# Large-Eddy Simulation framework for two-phase flows with heat transfer

Damien Fouquet, Julien Carmona, Vincent Moureau

► To cite this version:

Damien Fouquet, Julien Carmona, Vincent Moureau. Large-Eddy Simulation framework for two-phase flows with heat transfer. European Conference for Aerospace Sciences, Jun 2025, Rome, Italy. hal-05272251

HAL Id: hal-05272251

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

Submitted on 22 Sep 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Large-Eddy Simulation framework for two-phase flows with heat transfer

Damien Fouquet*†, Julien Carmona*, and Vincent Moureau*

* CNRS, INSA Rouen Normandie, Univ Rouen Normandie, Normandie Univ, CORIA UMR 6614, F-76000, Rouen, France

†Corresponding author: damien.fouquet@coria.fr

## Abstract

Within the past decades, many methods have been developed to study two-phase flows with heat and mass transfer. Most of these methods rely on the Direct-Numerical-Simulation (DNS) approach. Although giving an accurate representation of the heat and mass transfer at the interface, a high resolution is needed to accurately compute the gradients at the interface, becoming unaffordable for industrial processes at high Prandtl or Schmidt number, where thermal or species boundary layers becomes very thin. In this work, a Large-Eddy Simulation (LES) framework is presented and applied to two-phase heat transfer. It relies on a multi-fluid approach where a temperature is defined in each phase, weighted by a sharp phase indicator, and then filtered. The resulting filtered phase indicator has the same profile as in the Accurate-Conservative-Level-Set¹² method, which enables to use Chiodi and Desjardins’ reinitialization equation⁷ to maintain the profile over time. The filtering of the phase-weighted temperature equations leads to a set of transport equations in which the interface temperature and heat flux appear and need to be modeled. In this work, a closure is proposed assuming the sub-grid temperature profile at the interface. Moreover, the filtered phase weighted temperatures should be consistent with the filtered phase indicator. To do so, the filtered phase-weighted temperatures are reinitialized at the interface with a scalar reinitialization equation specifically derived and aiming at reimposing the modeled interface temperature profile. The developed framework is applied to two configurations: a sloshing tank where the numerical pressure drops are compared to experimental pressure drops from Himeno et al.,¹⁶ then on a two-dimensional jet in cross-flow, showing the ability of the framework in a complex configuration.

## 1. Introduction

Interfacial two-phase flows with heat and mass transfer play an important role in various industrial processes such as fuel injection and evaporation in combustion chambers or sloshing in cryogenic tanks. Methods have been developed to study two phase-flows with heat and mass transfer, mainly within a Direct-Numerical-Simulation approach.¹⁸,²³⁻²⁵,²⁹ These methods rely on two key ingredients: i) capturing or tracking the interface position, and ii) an accurate calculation of heat/mass/momentum exchange terms. For the first ingredient, different methods have been developed and can be classified into two approaches. On the one hand, the Front-Tracking method,²⁷ an Eulerian-Lagrangian method, where Lagrangian markers are advected to follow the interface position whereas the other quantities are described on an Eulerian mesh. On the other hand, Front-Capturing methods, such as Volume-of-Fluid¹⁷ method, Level-Set²² method, or Accurate-Conservative-Level-Set¹² (ACLS) method, where the interface position is implicitly retrieved based on a color function. The different methods have been applied in various configurations, on both structured and unstructured grids. For the second ingredient, the calculation of heat/mass/momentum exchange terms requires a high resolution at the interface to estimate with sufficient accuracy the temperature, species mass fractions, and velocity gradients. While this resolution can be reached in a DNS framework but at the expanse of a high CPU cost, it is generally not possible in the simulation of industrial processes, where characteristic length scales can range from nanometers to meters, as encountered in cryogenic tanks.

To tackle this challenge, two approaches have been developed. On the one hand, as for single-phase flows, Large-Eddy Simulation (LES) equations have been derived for two-phase flows as introduced by Labourasse et al.²⁰ for isothermal two-phase flows. Compared to single-phase LES equations, new sub-grid terms arise when the one-fluid mass and momentum equations are filtered. Closure have been proposed for some terms, for example by Toutant et al.²⁶ based on 3D *a priori* tests, or by Hermann¹⁵ for the sub-grid contribution of surface tension. On the other hand, DNS equations can be employed directly but with sub-grid reconstruction to adjust convective and diffusive species fluxes or temperature fluxes within cells containing the liquid-gas interface. This approach has been applied in a Volume-of-Fluid framework,³,⁵,²⁸ and interface-tracking strategies.¹,²,⁹
In this work, an LES approach is proposed, named Filtered-Interface-Multi-Fluid (FIMF) approach, and relies on both approaches described before. The transport equations are obtained by explicitly filtering the continuity equations. The filtering operation introduces interface quantities in the equations: interface temperature and heat flux appear and need to be modeled. These terms are closed assuming the temperature profile at the interface at the sub-grid level based on an analytical solution. Section 2 introduces the governing equations, as well as the retained model for the interface quantities. Then, some details about the numerical implementation are given in Section 3. Finally, the proposed framework is applied to two cases in Section 4, a sloshing tank for which results are compared with experimental data regarding pressure drop, and a two-dimensional jet in cross-flow, showing the good behavior of the proposed framework when complex motions and deformations of the interface are considered.

## 2. Governing equations

### 2.1 DNS equations

The proposed framework relies on a multi-fluid approach to study heat transfer within the Low-Mach number hypothesis. First, the equations are presented assuming that the density does not depend on the temperature. For simplicity, the temperature is transported here. The extension to a gas with a variable density will be briefly described after the derivation of the LES equation (Sec. 2.3). In this section, the DNS equations are introduced and will be filtered in the next section to obtain the LES equations. Although the proposed framework can be applied for N phases, the equations are presented for two phases, a liquid and a gas, denoted by subscripts l and g.

#### 2.1.1 Phase indicator equation

To distinguish the two phases, a phase indicator is introduced. Denoted $\phi_i$ for phase i, the phase indicator is defined as:

$$ \phi_i(\mathbf{x}_i, t) = \begin{cases} 1 & \text{if position } \mathbf{x}_i \text{ is in phase } i \text{ at time } t \\ 0 & \text{otherwise,} \end{cases} \qquad (1) $$

with the additional constraint that $\phi_l + \phi_g = 1$ everywhere in space. Time evolution of a phase indicator $\phi_i$ is given by the following transport equation:

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

$\mathbf{u}_{\Gamma}$ being the velocity of the phase interface and t the time. Note that for an incompressible flow without discontinuous velocity between phases, Eq. (2) is equivalent to:

$$ \frac{\partial \phi_i}{\partial t} + \nabla \cdot (\phi_i \mathbf{u}) = 0. \qquad (3) $$

where $\mathbf{u}$ is the one-fluid velocity discussed hereafter.

#### 2.1.2 Mass and momentum equations

Considering an incompressible and isothermal flow in each phase, the mass and momentum equations within the one-fluid formalism read:

$$ \nabla \cdot \mathbf{u} = 0, \qquad (4) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot (\mathbf{u} \otimes \mathbf{u}) = -\frac{1}{\rho} \nabla p + \nabla \cdot (\nu (\nabla \mathbf{u} + \nabla^T \mathbf{u})) + \mathbf{g}, \qquad (5) $$

where $\rho$ the density, $p$ the pressure, $\nu$ the kinematic viscosity, and $\mathbf{g}$ the gravity vector. These equations (Eqs (4)-(5)) are solved for both the gas and liquid phases, assuming that the material properties remain constant. Density and viscosity are discontinuous across the interface $\Gamma$ and can therefore be expressed as a function of the phase indicators $\phi_i$:

$$ \rho = \phi_l \rho_l + \phi_g \rho_g, \qquad (6) $$

$$ \mu = \phi_l \mu_l + \phi_g \mu_g, \qquad (7) $$

where $\rho_i$ and $\mu_i$ are the density and the dynamic viscosity of the pure phase i = l, g, respectively. Moreover, considering an absence of phase change and viscous fluids, the velocity is continuous across the interface and thus the jump in velocity $[\mathbf{u}]_{\Gamma}$ writes:

$$ [\mathbf{u}]_{\Gamma} = \mathbf{u}_l - \mathbf{u}_g = 0. \qquad (8) $$
Due to the effects of surface tension and the normal component of viscous stress, a pressure jump $[p]_Γ$ occurs across the interface:¹¹

$$[p]_Γ = \sigma\kappa + 2[\mu]_Γ \mathbf{n}_l' \cdot \nabla \mathbf{u} \cdot \mathbf{n}_l, \quad (9)$$

where $\sigma$ refers to the surface tension coefficient and $\kappa$ the interface curvature. This pressure jump is explicitly handled using the Ghost-Fluid Method,¹⁴ a robust and accurate approach for addressing this condition within a one-fluid formulation.

### 2.1.3 Multi-fluid temperature equation

According to the multi-fluid formalism, the temperature of each phase $T_i$ with $i = l, g$ has to be transported. Each transport equation is obtained by multiplying the temperature transport equation considering a pure phase by its phase indicator $\phi_i$:

$$\phi_i \left( \frac{\partial T_i}{\partial t} + \nabla \cdot (T_i \mathbf{u}) \right) = \phi_i \nabla \cdot \left( \frac{\lambda_i}{\rho_i c_{p,i}} \nabla T_i \right), \quad (10)$$

where $\lambda_i$ and $c_{p,i}$ correspond to the thermal conductivity and the heat capacity at constant pressure of the phase $i = l, g$. Using the transport equation of the phase indicator $\phi_i$ (Eq. (3)), Eq. (10) can be simplified as:

$$\frac{\partial \phi_i T_i}{\partial t} + \nabla \cdot (\phi_i T_i \mathbf{u}) = \nabla \cdot \left( \frac{\lambda_i}{\rho_i c_{p,i}} \phi_i \nabla T_i \right) - \frac{\lambda_i}{\rho_i c_{p,i}} \nabla \phi_i \cdot \nabla T_i. \quad (11)$$

The second term of the right hand side of Eq. (11) corresponds to the interface heat flux. It must be pointed out that this term is only non-zero at the interface. Within the multi-fluid formalism, the temperature is defined over the entire domain as the sum of the contributions from each phase:

$$T = \sum_{i=l,v} \phi_i T_i. \quad (12)$$

## 2.2 LES equations

To obtain the LES equations, the previous equations have to be filtered. Two kinds of filtering can be considered: an implicit filtering associated to the mesh, or an explicit filtering of the equations, introducing a filtering kernel $G$. The explicit filtering of a flow quantity $\xi_i$ is written as:

$$\tilde{\xi}_i = G * \xi_i(\mathbf{x}) = \int G(\mathbf{x} - \mathbf{y}) \xi_i(\mathbf{y}) d\mathbf{y}. \quad (13)$$

For convenience, the Drew¹³ notation is introduced $\tilde{\xi}_i$, corresponding to the phase weighted filtered $\xi_i$:

$$\tilde{\xi}_i = \frac{\overline{\phi_i \xi_i}}{\overline{\phi_i}}, \quad (14)$$

where $\overline{\phi_i}$ is the filtered phase indicator.

The filtering kernel considered in this work is given for a one-dimensional case as:

$$G(x_i) = \frac{1}{4\epsilon \cosh^2(x_i/2\epsilon)}, \quad (15)$$

with $\epsilon$ a parameter controlling the filter width.

### 2.2.1 Filtered phase indicator equations

The filtering of Eq. (3) leads to the following equation:

$$\frac{\partial \bar{\phi}_i}{\partial t} + \nabla \cdot (\bar{\phi}_i \bar{\mathbf{u}}) + \nabla \cdot (\overline{\phi_i \mathbf{u}} - \bar{\phi}_i \overline{\mathbf{u}}) = 0, \quad (16)$$

where $(\overline{\phi_i \mathbf{u}} - \bar{\phi}_i \overline{\mathbf{u}})$ is a term representing the contribution of the sub-grid scale in the phase indicator transport. In this work, this contribution is neglected, the modeling of this sub-grid contribution being a complex topic in its own in the
two-phase community. For some details about the modeling, see for example the work of Toutant et al.²⁶ Neglecting this contribution, the transport equation of the phase indicator becomes:

$$ \frac{\partial \bar{\phi}_l}{\partial t} + \nabla \cdot (\bar{\phi}_l \bar{\mathbf{u}}) = 0. \qquad (17) $$

For the specific case of two-phase flow, the filtered gas phase indicator is simply given by: $\bar{\phi}_g = 1 - \bar{\phi}_l$. Given the filtering kernel defined by Eq. (15) and assuming a planar interface, the filtered phase indicator reads:

$$ \bar{\phi}_i = \frac{1}{2} \left( 1 + \tanh \left( \frac{x_i}{2\epsilon} \right) \right). \qquad (18) $$

For a multi-dimensional problem, this profile can be recast introducing the signed distance $\psi_i$ to the $\frac{1}{2}$ iso-contour of $\bar{\phi}_i$:

$$ \bar{\phi}_l = \frac{1}{2} \left( 1 + \tanh \left( \frac{\psi_l(\mathbf{x}, t)}{2\epsilon} \right) \right). \qquad (19) $$

where $\epsilon$ is linearly dependent to the local cell size $\Delta_x$: $\epsilon = \epsilon_{fac}\Delta_x$, and ensure that $\bar{\phi}_l$ is well discretized considering the mesh resolution. A common value of $\epsilon_{fac}$ is 0.5. Having $\bar{\phi}_l$ described by Eq. (19) transported with Eq. (17) corresponds to the transport equation of the ACLS method.¹² To avoid numerical diffusion of $\bar{\phi}_l$, its profile has to be reinitialized after each time-step. The robust and accurate reinitialization proposed by Chiodi and Desjardins⁷ and adapted to unstructured grids by Janodet et al.¹⁹ is used:

$$ \frac{\partial \bar{\phi}_l}{\partial \tau} = \nabla \cdot \left( \frac{1}{4 \cosh^2\left(\frac{\phi_{map,l}}{2\epsilon}\right)} \left( \nabla \phi_{map,l} \cdot \mathbf{n}_l - \mathbf{n}_l \cdot \mathbf{n}_l \right) \right) \mathbf{n}_l \qquad (20) $$

$$ \text{with } \phi_{map,l} = \epsilon \ln \left( \frac{\bar{\phi}_l}{1 - \bar{\phi}_l} \right). \qquad (21) $$

### 2.2.2 Filtered mass and momentum equations

As explained in Sec. 2.1.2, the mass and momentum equations are solved in the one-fluid formalism. The evaluation of density and viscosity is carried out as described in the work of Janodet et al.,¹⁹ i.e. the physical properties are averaged around the interface with $\rho_i$ and $v_i$ constant per phase. For the sake of clarity, the density field is denoted as $\bar{\rho}$ even if a sharp representation is used here:

$$ \bar{\rho} = \sum_i \mathcal{H}(\mathbf{x}_i) \rho_i. \qquad (22) $$

Contrarily to the density, the viscosity is a filtered field and is calculated as follows:

$$ \bar{v} = \sum_i \bar{\phi}_i v_i. \qquad (23) $$

The filtered mass equation is given by:

$$ \nabla \cdot \bar{\mathbf{u}} = 0, \qquad (24) $$

with $\bar{\mathbf{u}} = \bar{\phi}_l \tilde{\mathbf{u}}_l + \bar{\phi}_g \tilde{\mathbf{u}}_g$. To determine the filtered momentum equation, the following assumptions are considered:

$$ \left( \frac{\nabla p}{\rho} \right) = \frac{\nabla \bar{p}}{\bar{\rho}}, \qquad (25) $$

$$ \nu (\nabla \mathbf{u} + \nabla \mathbf{u}^T) = \bar{v} (\nabla \bar{\mathbf{u}} + \nabla \bar{\mathbf{u}}^T), \qquad (26) $$

$$ \bar{\mathbf{u}} \otimes \bar{\mathbf{u}} - \bar{\mathbf{u}} \otimes \bar{\mathbf{u}} = \nu_t (\nabla \bar{\mathbf{u}} + \nabla \bar{\mathbf{u}}^T), \qquad (27) $$

with $\nu_t$ the turbulent viscosity, based on Boussinesq⁴ assumption. It leads to the following one-fluid momentum equation:

$$ \frac{\partial \bar{\mathbf{u}}}{\partial t} + \nabla \cdot (\bar{\mathbf{u}} \otimes \bar{\mathbf{u}}) = - \frac{\nabla \bar{p}}{\bar{\rho}} + \nabla \cdot \left( (\bar{v} + \nu_t) (\nabla \bar{\mathbf{u}} + \nabla \bar{\mathbf{u}}^T) \right) + \mathbf{g}. \qquad (28) $$
### 2.2.3 Filtered phase temperature equations

Assuming a planar interface without transverse temperature gradient, the filtering of Eq. (11) leads to the following LES equation:

$$ \frac{\partial \overline{\phi_i \widetilde{T}_i}}{\partial t} + \nabla \cdot (\overline{\phi_i \widetilde{T}_i \bar{\mathbf{u}}}) = \nabla \cdot \left( \frac{\lambda_i}{\rho_i c_{p,i}} (\nabla \overline{\phi_i \widetilde{T}_i} - T^\Gamma \nabla \overline{\phi_i}) \right) - \frac{j}{\rho_i c_{p,i}} \nabla \overline{\phi_i} \cdot \mathbf{n} + \nabla \cdot (\overline{\phi_i \widetilde{T}_i \bar{\mathbf{u}} - \phi_i \widetilde{T}_i \bar{\mathbf{u}}}) , \quad (29) $$

where $j = \lambda_i \nabla T_i |_{\Gamma} \cdot \mathbf{n}$ is the interface heat flux and $T^\Gamma$ is the interface temperature. In a LES context, the thermal gradients are not resolved at the interface and the previous quantities have to be modeled. The next section will focus on this point. The sub-grid transport of temperature appears in this equation, and it is assumed that it can be closed using a turbulent diffusivity:

$$ \overline{\phi_i \widetilde{T}_i \bar{\mathbf{u}} - \phi_i \widetilde{T}_i \bar{\mathbf{u}}} = D_t \overline{\phi_i \nabla T_i} = D_t (\nabla \overline{\phi_i \widetilde{T}_i} - T^\Gamma \nabla \overline{\phi_i}) . \quad (30) $$

Far from the interface, a single phase flows is retrieved and the closure is classical. However, we highlight that in the interface vicinity, the closure might not be suited. As for the closure of the sub-grid interface transport, the closure of this term is a complex topic on its own and must be improved in future work. The phase weighted temperature equation finally reads:

$$ \frac{\partial \overline{\phi_i \widetilde{T}_i}}{\partial t} + \nabla \cdot (\overline{\phi_i \widetilde{T}_i \bar{\mathbf{u}}}) = \nabla \cdot \left( \left( \frac{\lambda_i}{\rho_i c_{p,i}} + D_t \right) (\nabla \overline{\phi_i \widetilde{T}_i} - T^\Gamma \nabla \overline{\phi_i}) \right) - \frac{j}{\rho_i c_{p,i}} \nabla \overline{\phi_i} \cdot \mathbf{n} . \quad (31) $$

## 2.3 Extension to variable density

The previous equations have been derived assuming a constant density in each phase. However, to retrieve some thermodynamic effects, as pressure increase or decrease due to heat transfer, it is necessary to consider that gaseous density depends on the temperature. The extension to variable density is performed assuming that the low-Mach number hypothesis is still verified. When gas phase has a variable density, the velocity field is no longer divergence free in the gaseous phase and pressure effects must be taken into account in the temperature equation. The filtered multi-fluid temperature equation becomes for the gaseous phase

$$ \frac{\partial \overline{\phi_i \widetilde{T}_i}}{\partial t} + \nabla \cdot (\overline{\phi_i \widetilde{T}_i \bar{\mathbf{u}}}) = \nabla \cdot \left( \left( \frac{\lambda_i}{\rho_i c_{p,i}} + D_t \right) (\nabla \overline{\phi_i \widetilde{T}_i} - T^\Gamma \nabla \overline{\phi_i}) \right) - \frac{j}{\rho_i c_{p,i}} \nabla \overline{\phi_i} \cdot \mathbf{n} + \frac{\overline{\phi_i}}{\rho_i c_{p,i}} \frac{dP}{dt} + \overline{\phi_i \widetilde{T}_i} (\nabla \cdot \bar{\mathbf{u}}) , \quad (32) $$

and for the liquid phase:

$$ \frac{\partial \overline{\phi_i \widetilde{T}_i}}{\partial t} + \nabla \cdot (\overline{\phi_i \widetilde{T}_i \bar{\mathbf{u}}}) = \nabla \cdot \left( \left( \frac{\lambda_i}{\rho_i c_{p,i}} + D_t \right) (\nabla \overline{\phi_i \widetilde{T}_i} - T^\Gamma \nabla \overline{\phi_i}) \right) - \frac{j}{\rho_i c_{p,i}} \nabla \overline{\phi_i} \cdot \mathbf{n} + \overline{\phi_i \widetilde{T}_i} (\nabla \cdot \bar{\mathbf{u}}) . \quad (33) $$

An additional term $dP/dt$ appears in the gas temperature equation, where $P$ is a thermodynamic pressure. This term is computed based on the density field with a prediction/correction method: first, the phase temperatures are transported with a predictor of $dP/dt$, then the new gaseous mass is computed using the density updated with the equation-of-state, finally the pressure associated to the given gaseous mass is computed, giving a new $dP/dt$, that is used to correct the gaseous temperature equation.

## 2.4 Interface heat transfer modeling

The filtered temperature equation, Eq. (31), has two under-resolved contributions that need to be modeled: the interface temperature $T^\Gamma$ and heat flux $j$. A closure for these two variables is proposed in this section. For two semi-infinite media into contact at $x = 0$, the phase $l$ occupying the space $\Omega_l = \{x | x > 0\}$ and the phase $g$ occupying the space $\Omega_g = \{x | x < 0\}$, the time-dependent analytical solution⁶ reads:

$$ T(x, t) = \begin{cases} T^\Gamma + \frac{e_g}{e_l+e_g} \Delta T^0 \text{erf} \left( \frac{x}{2\sqrt{D_l t}} \right) & \text{if } x > 0, \\ T^\Gamma - \frac{e_l}{e_l+e_g} \Delta T^0 \text{erf} \left( \frac{x}{2\sqrt{D_g t}} \right) & \text{if } x < 0, \end{cases} \quad \text{with } T^\Gamma = \frac{e_l T_l^0 + e_g T_g^0}{e_l + e_g}, \quad (34) $$

where $e_l$ is the effusivity of phase $i$, $\Delta T^0 = T_l^0 - T_g^0$ the initial temperature difference between the two phases. The time $t$ appearing in this expression corresponds to the time after the contact of the two media and it is an unknown of the problem here. This analytical solution is used to derive the model for the interface quantities, i.e. it is assumed that at
the sub-grid scale the temperature profile is given by Eq. (34). The filtered temperature profile is $\overline{T}(x, t) = \overline{\phi}_l \overline{T}_l + \overline{\phi}_g \overline{T}_g$ with:

$$ \left\{ \begin{aligned} \overline{\phi}_l \overline{T}_l &= \overline{\phi}_l T^\Gamma + \frac{e_g}{e_l+e_g} \Delta T^0 \overline{\phi}_l \operatorname{erf} \left( \frac{\psi_l}{2\sqrt{D_g t}} \right), \\ \overline{\phi}_g \overline{T}_g &= \overline{\phi}_g T^\Gamma - \frac{e_l}{e_l+e_g} \Delta T^0 \overline{\phi}_g \operatorname{erf} \left( \frac{\psi_g}{2\sqrt{D_g t}} \right). \end{aligned} \right. \qquad (35) $$

Now, the modeled interface temperature and heat flux can be determined using known information. In Eq. (35), three unknowns appear: $T^\Gamma$, $\Delta T^0$, and $t$. As only two relations are known (Eq. (35)), another relation should be added. In our case, we considered that $T^\Gamma$ can be determined as

$$ T^\Gamma = \frac{e_l \overline{T}_l |_\Gamma + e_g \overline{T}_g |_\Gamma}{e_l + e_g}. \qquad (36) $$

where $\overline{T}_l|_\Gamma$ is the interface value of $\overline{T}_l$ at the resolved scale. $\Delta T^0$, or equivalently $j$, and $t$ can now be determined. Combining the previous filtered Eq. (35)-(36) equations and assuming that in explicit time stepping a Fourier-like stability condition is satisfied

$$ t \ll \frac{\psi^2}{4D_i}, \qquad (37) $$

and,

$$ \overline{\phi}_i \operatorname{erf}(\psi) \approx \operatorname{erf}(\overline{\phi}_i \overline{\psi}), \qquad (38) $$

simple expressions are derived for the time $t$ and the heat flux $j$:

$$ \frac{1}{\sqrt{t}} = \frac{2 \sqrt{D_g}}{\overline{\phi}_l \overline{\psi}_l} \operatorname{erf}^{-1} \left( \left( 1 + \frac{e_l}{e_g} \right) \frac{\overline{\phi}_l \overline{T}_l - \overline{\phi}_l T^\Gamma}{\overline{T}_l - \overline{T}_g} \right) = \frac{2 \sqrt{D_g}}{\overline{\phi}_g \overline{\psi}_g} \operatorname{erf}^{-1} \left( \left( 1 + \frac{e_g}{e_l} \right) \frac{\overline{\phi}_g T^\Gamma - \overline{\phi}_g \overline{T}_g}{\overline{T}_l - \overline{T}_g} \right), \qquad (39) $$

$$ j = \frac{1}{\sqrt{\pi t}} \frac{\overline{T}_l - \overline{T}_g}{\frac{1}{\alpha_l} + \frac{1}{\alpha_g}} \quad \text{with} \quad \alpha_i = \frac{e_i}{\overline{\operatorname{erf}}\left(\frac{\psi_l}{2\sqrt{D_i t}}\right)}. \qquad (40) $$

To summarize, this model is able to compute the interface temperature and the heat flux at the interface from the liquid and gas filtered temperatures assuming an approximate heat conduction problem at the sub-grid scale. The different approximations will be assessed globally in the following.

## 2.5 Correction of numerical consistency errors

### 2.5.1 Scalar reinitialization equations

The reinitialization of the Level-Set is applied so that $\overline{\phi}_i$ keeps a hyperbolic tangent profile despite transport numerical errors. This reinitialization creates fluxes of the phase indicator that are not present in the scalar transport equation. A first step to regain consistency between the phase indicator and the scalar transport is obtained by applying a similar operation to $\overline{\phi}_i \overline{T}_i$. The method proposed in this work consists in a reinitialization of $\overline{\phi}_i \overline{T}_i$ in a narrow band around the interface. A profile consistent with the model for interface quantities should be reimposed at the interface. If we assume that the profile is given by the model presented in Sec. 2.4, the profile to reimpose in the narrow band is:

$$ \overline{\phi}_i \overline{T}_i = \overline{\phi}_i T^\Gamma \pm j \frac{\sqrt{\pi t}}{e_i} \overline{\phi}_i \operatorname{erf} \left( \frac{\psi_i}{2\sqrt{D_i t}} \right). \qquad (41) $$

To derive the reinitialization equation, it is first assumed that the temperature is constant at the interface. The generalization to a case with temperature variation will be presented after the derivation. The profile to reimpose thus becomes:

$$ \overline{\phi}_i \overline{T}_i = \overline{\phi}_i T^\Gamma. \qquad (42) $$

As for the Level-Set reinitialization, this scalar reinitialization has two objectives: 1) prevent $\overline{\phi}_i \overline{T}_i$ from diffusing into its complementary phase, 2) avoid scalar accumulation or overshoots at the interface. A difference between the reinitialization of $\overline{\phi}_i$ and of $\overline{\phi}_i \overline{T}_i$ is the fact that the one on the scalars must be done only at the interface and not elsewhere. To achieve this, a mask $\chi_i$ is introduced to apply reinitialization only where needed. Different tests have
been made, and the mask used in the following is $\chi_i = (4\bar{\phi}_i(1 - \bar{\phi}_i))^{1/2}$. Hence, the following equation has been considered to reinitialize the scalars:

$$ \frac{\partial \bar{\phi}_i \tilde{T}_i}{\partial \tau} + \nabla \cdot \left( \chi_i \left( \left(1 - \frac{\bar{\phi}_i \tilde{T}_i}{T^\Gamma}\right) \bar{\phi}_i \tilde{T}_i - \left(1 - \bar{\phi}_i\right) \bar{\phi}_i T^\Gamma \right) \mathbf{n}_i \right) = \nabla \cdot \left( \chi_i \epsilon \left( \nabla (\bar{\phi}_i \tilde{T}_i - \bar{\phi}_i T^\Gamma) \cdot \mathbf{n}_i \right) \mathbf{n}_i \right). \quad (43) $$

When heat flux at the interface is considered, the profile to reimpose becomes $\bar{\phi}_i \tilde{T}_i = \bar{\phi}_i T^\Gamma + \bar{\phi}_i \tilde{P}_i$, here with $\tilde{P}_i = \pm j \frac{\sqrt{\pi i}}{e_i} \operatorname{erf}\left(\frac{\psi_i}{2\sqrt{D_i T}}\right)$, and the following procedure is applied:

1. Compute $\bar{\phi}_i \tilde{T}'_i = \bar{\phi}_i \tilde{T}_i - \bar{\phi}_i \tilde{P}_i$,

2. Use $\bar{\phi}_i \tilde{T}'_i$ to converge Eq. (43) toward $\bar{\phi}_i \tilde{T}'_i = \bar{\phi}_i T^\Gamma$,

3. Compute $\bar{\phi}_i \tilde{T}_i$ as $\bar{\phi}_i \tilde{T}_i = \bar{\phi}_i \tilde{T}'_i + \bar{\phi}_i \tilde{P}_i$.

As the reinitialization procedure is not directly applied to $\bar{\phi}_i \tilde{T}_i$, it is necessary to ensure that it is conservative. Indeed, issues can arise when $\bar{\phi}_i \tilde{P}_i$ changes during reinitialization steps. This quantity depends on two parameters: $\bar{\phi}_i$, that is reinitialized before $\bar{\phi}_i \tilde{T}_i$ and not modified during scalar reinitialization; and the heat flux, computed at the beginning of the time step and not modified during the scalars reinitialization. So, $\bar{\phi}_i \tilde{P}_i$ is not modified during reinitialization steps and as Eq. (43) is conservative, $\bar{\phi}_i \tilde{T}_i$ is also conserved during the scalar reinitialization when a flux is considered.

### 2.5.2 Pseudo phase indicator

The previous section introduced Eq. (43) that has been added to maintain the consistency between $\bar{\phi}_i$ and $\bar{\phi}_i \tilde{T}_i$ over time. However, the full consistency of both fields is ensured only if the two reinitialization equations (Eq. (43) and Eq. (20)) are perfectly converged and with the same truncation errors. However, as the equations are iterated based on a fixed number of iterations and that the reinitialization equations are different, the full consistency is not exactly ensured. Hence, $\tilde{T}_i$, $i = l, g$, necessary to compute $T^\Gamma$ and j modeled in Sec. 2.4 can suffer from these consistency errors that can grow over time. They can lead to overestimation or underestimation of the heat flux at the interface, creating non-physical behavior of the solution. The accurate computation of the interface quantities required for the transport equation and the reinitialization equation is thus of great importance. To reduce again the consistency errors, an additional field is transported. Considering $I_i$, $i = l, g$, a constant field equal to one, transported with the following equation:

$$ \frac{\partial I_i}{\partial t} + \nabla \cdot (I_i \mathbf{u}) = 0. \quad (44) $$

The equation is multiplied by $\phi_i$ and filtered, leading to

$$ \frac{\partial \bar{\phi}_i \tilde{I}_i}{\partial t} + \nabla \cdot (\bar{\phi}_i \tilde{I}_i \tilde{\mathbf{u}}) = 0. \quad (45) $$

Theoretically, as $I_i$ is equal to one everywhere, $\bar{\phi}_i \tilde{I}_i = \bar{\phi}_i$. To verify this property numerically, $\bar{\phi}_i \tilde{I}_i$ and $\bar{\phi}_i$ must be transported and reinitialized with the same equation. To quantify consistency errors, $\bar{\phi}_i \tilde{I}_i$ is reinitialized with the same equation than $\bar{\phi}_i \tilde{T}_i$, which reads:

$$ \frac{\partial \bar{\phi}_i \tilde{I}_i}{\partial \tau} + \nabla \cdot \left( \chi_i \left( \left(1 - \bar{\phi}_i \tilde{I}_i\right) \bar{\phi}_i \tilde{I}_i - \left(1 - \bar{\phi}_i\right) \bar{\phi}_i \right) \mathbf{n}_i \right) = \nabla \cdot \left( \chi_i \epsilon \left( \nabla (\bar{\phi}_i \tilde{I}_i - \bar{\phi}_i) \cdot \mathbf{n}_i \right) \mathbf{n}_i \right). \quad (46) $$

In Eq. (46), the theoretical value of $I_\Gamma$ is used, which is $I_\Gamma = 1$. As for $\bar{\phi}_i \tilde{T}_i$, if the reinitialization equation is exactly converged, $\bar{\phi}_i \tilde{I}_i$ is consistent with $\bar{\phi}_i$ and $\bar{\phi}_i \tilde{I}_i = \bar{\phi}_i$. Consistency errors can now be corrected introducing a new quantity: for any field $\xi_i$, $\langle \xi_i \rangle$ is defined as:

$$ \langle \xi_i \rangle = \frac{\bar{\phi}_i \tilde{\xi}_i}{\bar{\phi}_i \tilde{I}_i} = \frac{1}{\tilde{I}_i} \tilde{\xi}_i. \quad (47) $$

$(\tilde{I}_i)^{-1}$ appears as a correction factor of $\tilde{\xi}_i$. When the reinitialization equations are perfectly converged, $\tilde{I}_i = 1$ and the correction is not necessary. For not perfectly converged reinitialization, $\tilde{I}_i$ might differ from 1, as well as $\tilde{T}_i$ might differ
from its physical value: the division by $\tilde{T}_i$ corrects this error. Using this new quantity, the interface model defined in Sec. 2.4 is rewritten as

$$j = \frac{1}{\sqrt{\pi t}} \frac{\langle T_l \rangle - \langle T_g \rangle}{\frac{1}{\alpha_l} + \frac{1}{\alpha_g}}, \quad \text{with } \alpha_i = \frac{e_i}{\langle \text{erf} \rangle \left( \frac{\psi_i}{2\sqrt{D_i t}} \right)} \qquad (48)$$

$$T^{\Gamma} = \frac{e_l \langle T_l \rangle |_{\Gamma} + e_g \langle T_g \rangle |_{\Gamma}}{e_l + e_g}. \qquad (49)$$

### 3. Numerical methods

The proposed framework is implemented in the multi-physics code YALES2,²¹ a pair-based finite-volume low-Mach number solver with node-centered control volumes, designed for unstructured meshes and to be massively parallel. The implementation of the developed framework is carried out as described in the work of R. Janodet et al.¹⁹ regarding the hydrodynamic part. Regarding the interface model, the heat flux and temperature are computed at the edges crossed by the interface, these pairs values are then extended at the nodes in a narrow band around the interface using the Geometrical Projection Markers Method.¹⁹ The transport and reinitialization are implemented in the same manner than in the work of R. Janodet et al.,¹⁹ with some specifics treatment for the reinitialization of $\bar{\phi}_i \tilde{T}_i$ and $\bar{\phi}_i \tilde{T}_i$. Equations (43) and (20) are solved using a fourth-order Runge-Kutta method in the following way: first, Eq. (20) is iterated for a number of steps $N_{\bar{\phi}_i}$ with a pseudo-time step $\Delta\tau_{\bar{\phi}_i}$. Then, Eq. (43) is iterated for a number of steps $N_{\bar{\phi}_i \tilde{T}_i}$ with a pseudo-time step $\Delta\tau_{\bar{\phi}_i \tilde{T}_i}$. $N_{\bar{\phi}_i \tilde{T}_i}$ is a user-defined parameter, which is usually chosen around five following the work of Janodet et al.¹⁹ The pseudo-time used to iterate Eq. (20) is determined as in Chiodi and Desjardins⁷

$$\Delta\tau_{\bar{\phi}_i} = \min\left(\text{CFL}\Delta x, \frac{F_0\Delta x^2}{4\epsilon}\right), \qquad (50)$$

where CFL and $F_0$ are the Courant-Friedrichs-Lewy and Fourier numbers, set to 0.7 and 0.5 respectively for all the simulations. Once Eq. (20) has been iterated $N_{\bar{\phi}_i}$ times, a total pseudo-time $\tau_{\bar{\phi}_i} = N_{\bar{\phi}_i} \Delta\tau_{\bar{\phi}_i}$ has been done. The pseudo-time used to iterate Eq. (43) is computed in a similar manner. However, to improve the consistency between $\bar{\phi}_i$ and $\bar{\phi}_i \tilde{T}_i$, the same total pseudo-time should be done during the reinitialization. This means that the pseudo-time steps and the number of steps are determined as

$$\Delta\tau_{\bar{\phi}_i \tilde{T}_i} = \min\left(\Delta\tau_{\bar{\phi}_i}, \frac{\text{CFL}\Delta x}{\left\|\left(1 - \frac{\bar{\phi}_i \tilde{T}_i}{T^{\Gamma}}\right)\mathbf{n}_i\right\|_{\text{max}}}, \frac{F_0\Delta x^2}{4\epsilon}\right), \qquad (51)$$

$$N_{\bar{\phi}_i \tilde{T}_i} = \frac{\Delta\tau_{\bar{\phi}_i}}{\Delta\tau_{\bar{\phi}_i \tilde{T}_i}} N_{\bar{\phi}_i}. \qquad (52)$$

The full solution procedure is summarized here:

1. Compute $T^{\Gamma}$, j from the model derived in Sec. 2.4,

2. Transport $\bar{\phi}_i$, $\bar{\phi}_i \tilde{T}_i$, $\bar{\phi}_i \tilde{T}_i$, for $i = l, g$,

3. Reinitialize $\bar{\phi}_i$,

4. Reinitialize $\bar{\phi}_i \tilde{T}_i$, $\bar{\phi}_i \tilde{T}_i$, for $i = l, g$,

5. Solve momentum equation using a projection method:⁸ first, momentum equation is solved without the pressure term. Then, a Poisson equation on pressure is solved to ensure the divergence-free condition. More details can be found in the work of Janodet et al.¹⁹.

### 4. Results

#### 4.1 Validation

The developed framework is first validated on a one-dimensional case. The test case is as follows: initially, the temperature is linear on each side of the interface, the thermal gradient verifies $\lambda_g \nabla T_g = \lambda_l \nabla T_l$. When it is transported,
it should not be modified but only shifted of **ut** at time *t* compared to the initial solution, where **u** is the constant velocity field imposed for the test case. However, due to numerical errors, the thermal gradients are not maintained at the interface. This test case highlights the ability of the scalar reinitialization equation to maintain the gradient at the interface. The domain is one meter long, defined by $\Omega = [-0.5, 0.5]$, and the interface is initialized at $x = -0.3$ m. A constant velocity is imposed: $\mathbf{u} = 1 \text{ m.s}^{-1}$. Temperature fields are transported during $\Delta t = 0.6$ s. The properties of the fluids (oil and water) are defined in Table 1, except for the thermal conductivity, where different ratios have been investigated: 4, 10, 100, and 1000. All the comparisons are made on meshes ranging from 101 to 1001 nodes with increments of 100 nodes. Five reinitialization steps are considered. The $L_1$ and $L_2$ norms of the error on $\bar{T}$ are shown in Fig. 1. A comparison between the analytical temperature profiles and the obtained ones is also presented in Fig. 2 for the finest mesh. A convergence order between one and two is obtained for all the thermal conductivity ratios. Also, it appears that the higher the thermal conductivity ratio, the smaller the errors are. The higher the thermal conductivity ratio, the more the temperature derivatives are discontinuous, making the reinitialization equations robust in such cases.

<table><caption>Table 1: Physical properties of the two phases.</caption><thead><tr><th></th><th>Oil</th><th>Water</th></tr></thead><tbody><tr><td>ρ[kg.m<sup>-3</sup>]</td><td>930</td><td>1000</td></tr><tr><td>μ[Pa s]</td><td>6.45 × 10<sup>-5</sup></td><td>1.4354 × 10<sup>-4</sup></td></tr><tr><td>c<sub>p</sub>[J.kg<sup>-1</sup>.K<sup>-1</sup>]</td><td>2000</td><td>4180</td></tr><tr><td>e[J.K<sup>-1</sup>.m<sup>-2</sup>.s<sup>-1/2</sup>]</td><td>490</td><td>1584</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 1: $L_1(\bar{T})$ and $L_2(\bar{T})$ for different thermal conductivity ratios. Black dashed line corresponds to first order convergence and black solid line corresponds to second order convergence.</figcaption></figure>

## 4.2 Pressure drop in a sloshing tank

The developed framework is now applied to a sloshing tank. The studied configuration comes from the work of Himeno et al.,¹⁶ where they studied experimentally the pressure drop in a sloshing tank for different lateral accelerations. The set-up is as follows: 30 seconds before the beginning of sloshing, oil at 284 K and air at 298 K are in the tank. They let the tank evolve, while keeping a valve open for the pressure to remain constant. At $t = 0$ s, the sloshing begins: during the first 0.25 s, a lateral acceleration $-a$ is imposed to the tank; between $t = 0.25$ s and $t = 0.5$ s, $+a$ lateral acceleration is imposed; finally, no lateral acceleration is imposed for the remaining time. They studied different lateral accelerations $+a$: 0.1g, 0.2g, 0.3g, 0.4g, 0.5g with $g = 9.81 \text{ m.s}^{-2}$. From the start of the sloshing to the end of the pressure measurement, slightly less than 2.5 seconds are studied. Only these 2.5 seconds will be simulated in this section. The air contained in the tank is considered as an ideal gas. The physical properties of the two fluids are given in Table 2. The tank is 0.10 m wide and 0.2 m tall.

During the numerical investigations, it is found out that the first order parameters regarding the pressure drop were the initial temperature distribution and the captured interface area. To be representative of the experiment, the
<figure><img src="image_4.png" /><figcaption>Figure 2: Comparison of the temperature profile obtained (red line) with the analytical profile (black dashed line) at $t = 0.3s$.</figcaption></figure>

<table><caption>Table 2: Physical properties of the considered oil and air</caption><thead><tr><th></th><th>Oil</th><th>Air</th></tr></thead><tbody><tr><td>ρ [kg.m<sup>-3</sup>]</td><td>818</td><td>Ideal gas</td></tr><tr><td>μ [Pa.s]</td><td>8.18 × 10<sup>-4</sup></td><td>1.716 × 10<sup>-5</sup></td></tr><tr><td>λ [J.s<sup>-1</sup>.m<sup>-1</sup>.K<sup>-1</sup>]</td><td>0.1</td><td>0.0242</td></tr><tr><td>c<sub>p</sub> [J.kg<sup>-1</sup>.K<sup>-1</sup>]</td><td>2000</td><td>1006</td></tr><tr><td>σ [N.m]</td><td>0.0169</td><td>/</td></tr></tbody></table>

initial temperature distribution is taken as the analytical solution for two semi-infinite media at $t = 30$ s. In the following, 2D simulations are carried out. The mesh is refined around the interface using Adaptive Mesh Refinement (AMR), the prescribed mesh size at the interface being $Δx_Γ = 0.4$ mm. AMR is performed using MMG Library¹⁰ following Janodet et al.¹⁹

The pressure drop for the different lateral accelerations are presented in Fig. 3, as well as the experimental results. To understand the differences between the lateral accelerations, snapshot of the temperature field for the different accelerations are given in Fig. 4. The evolution of the interface area (length here, as the domain is two-dimensional) is given in Fig. 5. Pressure variations are correlated with the interface area: as the interface area increases, the pressure decreases. Also, interface area peaks correspond to breaks in the slope of the pressure curves. It is consistent with the
temperature fields in Fig. 3: a higher interface area corresponds to more heat exchange. The increase of interface area for cases with high accelerations corresponds to the creation of liquid droplets and ligaments, that are dispersed into the gas phase, enhancing the temperature mixing and consequently, leading to higher pressure drops.

For the cases with the smallest lateral acceleration (0.1g, 0.2g, 0.3g), the pressure evolution is close to the experimental results. The results for the highest lateral acceleration, i.e. 0.4g and 0.5g, exhibit some differences. While capturing the same trend, numerical pressure drop is smaller than in the experiment. The discrepancies between the experimental and numerical results could be explained in different ways: first, the simulations are two-dimensional and some three-dimensional effects are not taken into account, that could lead to a higher pressure drop. Also, the conjugated heat transfer with the solid surrounding the domain is not taken into account. The heat exchange with the solid could explain two phenomena: the higher pressure drop at the beginning of the experiment as the gas will exchange heat with the solid where the liquid was, cooling down the gas phase, and the pressure increase seen at the end for the highest lateral accelerations, as the gas phase will exchange with the surrounding solid, which is at the ambient temperature. Finally, it is considered that the temperature is homogeneous at $t = -30$ s within the liquid and the gas, at 284 K and 298 K, following the description in Himeno et al.¹⁶ However, the experimental temperature profiles are not homogeneous and the minimal liquid temperature is around 282.5 K, which can lead to higher pressure drops.

<figure><img src="image_5.png" /><figcaption>Figure 3: Comparison of the pressure drops (left: 2D simulations, right: experimental results, extracted from Himeno et al.¹⁶)</figcaption></figure>

## 4.3 Two-dimensional jet in cross-flow

This third case highlights the abilities of the framework in a complex configuration, with high deformation of the interface. The studied case corresponds to a two-dimensional jet in cross-flow, where a liquid is injected in a cross-flow of gas. This configuration is representative of fuel injection in engines, and is of interest for industrial applications to understand heat and mass transfer in such systems. Here, only the heat transfer is studied and even if two-dimensional, it shows how the developed framework is able to deal with such configuration. The liquid and gas have the physical properties given in Table 3. The test case is purely fictive and does not aim at validating the developed framework. Hydrodynamic simulation of a jet in cross-flow has already been performed with YALES2 in Janodet et al.¹⁹ and
<figure><img src="image_6.png" /><figcaption>Figure 4: Temperature fields at $t = 1$ s for the different lateral accelerations.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Interface length for the different lateral accelerations</figcaption></figure>

the goal here is to see if the developed framework behaves correctly for this kind of configuration, with many under-resolved droplets and ligaments. For our simulation, liquid is injected at 293 K and $u = 17 \text{ m.s}^{-1}$ from the bottom of the domain, and gas is injected at 313 K and $u = 37 \text{ m.s}^{-1}$ from the left of the domain. Liquid inlet, on the bottom of the domain, is 1.4 mm wide; the gas inlet, on the left of the domain, is 25 mm wide. The domain length is 60 mm. As for the previous case, MMG library is used to refine the mesh around the interface, the imposed mesh size at the interface being $\Delta x_{\Gamma} = 20 \mu\text{m}$. The global temperature field and a zoom-in around droplets of the temperature field are given at $t = 1$ ms in Fig. 6. First, the developed framework behaves well for this complex configuration and results are in line with physics. Also, the temperature field at the scale of the under-resolved droplets behaves well: while being under-resolved, the droplets keep exchanging heat with the surrounding gas over long distances and the reinitialization maintains the droplet temperature over long distances, avoiding diffusion of the temperature field at the resolution limit. This property comes mainly from the conservative nature of the framework: unlike one-fluid methods, the gas and liquid temperatures are transported conservatively and separately, and the heat exchange is modeled explicitly with a transfer term.

## 5. Conclusion

A multi-fluid LES method is proposed for the transport of temperature in the simulation of interfacial flows. This framework differs from classical one-fluid methods as the gas and liquid temperatures are transported separately with full conservation. Heat exchange is explicitly modeled with a new model based on the solving of an unsteady diffusion problem. Equations are derived considering constant properties in each phase. The extension to variable gaseous density is also presented to account for thermodynamic effects. Consistency between the filtered phase indicator and
<table><caption>Table 3: Physical properties for the Jet in Cross-flow</caption><thead><tr><th></th><th>Liquid</th><th>Gas</th></tr></thead><tbody><tr><td>ρ [kg.m<sup>-3</sup>]</td><td>790</td><td>11.59</td></tr><tr><td>μ [Pa.s]</td><td>1.23 × 10<sup>-3</sup></td><td>1.845 × 10<sup>-5</sup></td></tr><tr><td>λ [J.s<sup>-1</sup>.m<sup>-1</sup>.K<sup>-1</sup>]</td><td>0.679</td><td>0.15</td></tr><tr><td>c<sub>p</sub> [J.kg<sup>-1</sup>.K<sup>-1</sup>]</td><td>1000</td><td>1000</td></tr><tr><td>σ [N.m]</td><td>22.77 × 10<sup>-3</sup></td><td>/</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 6: Results obtained for the Jet in Cross-Flow, the white line representing the interface.</figcaption></figure>

the multi-fluid temperature is maintained introducing two key treatments: scalar reinitialization around the interface and a pseudo-phase indicator enabling the quantification of consistency errors as well as the correction of these errors. The developed framework is then validated on a one-dimensional case, showing a convergence order between one and two. Once validated, the framework is used to study pressure drop in a sloshing tank at different sloshing amplitudes. The results for smaller lateral accelerations show good agreement with the experimental results, while the pressure drops for higher lateral acceleration do not completely agree with the experimental ones. As conjugated heat transfer with the surrounding solid is not taken into account, it could explain the discrepancies obtained for the highest lateral accelerations. In the future, three-dimensional simulations will also be carried out to provide more accurate results. Finally, the framework is applied on a two-dimensional jet in cross-flow, showing the ability of the framework to maintain the temperature gradient on a complex configuration, with many under-resolved droplets.

## 6. Acknowledgments

The CRIANN (Centre Regional Informatique et d'Applications Numeriques de Normandie) is gratefully acknowledged (project No. 2012006) for its CPU resources. The authors acknowledge the ANR - France (French National Research Agency) for its financial support of the MACROBOIL project No. ANR-23-CE50-0018.

## References

[1] B Aboulhasanzadeh, S Hosoda, A Tomiyama, and G Tryggvason. A validation of an embedded analytical description approach for the computations of high schmidt number mass transfer from bubbles in liquids. *Chemical Engineering Science*, 101:165–174, 2013.

[2] B Aboulhasanzadeh, S Thomas, M Taeibi-Rahni, and G Tryggvason. Multiscale computations of mass transfer from buoyant bubbles. *Chemical Engineering Science*, 75:456–467, 2012.

[3] Dieter Bothe and Stefan Fleckenstein. A volume-of-fluid-based method for mass transfer processes at fluid particles. *Chemical Engineering Science*, 101:283–302, 2013.
[4] Joseph Boussinesq. *Essai sur la theorie des eaux courantes*. Impr. nationale., 1877.

[5] Kangbei Cai, Guangyuan Huang, Yuchen Song, Junlian Yin, and Dezhong Wang. A sub-grid scale model developed for the hexahedral grid to simulate the mass transfer between gas and liquid. *International Journal of Heat and Mass Transfer*, 181:121864, 2021.

[6] Horatio Scott Carslaw, John Conrad Jaeger, and Herman Feshbach. Conduction of heat in solids. *Physics Today*, 15(11):74–76, 1962.

[7] Robert Chiodi and Olivier Desjardins. A reformulation of the conservative level set reinitialization equation for accurate and robust simulation of complex multiphase flows. *Journal of Computational Physics*, 343:186–200, 2017.

[8] Alexandre Joel Chorin. The numerical solution of the navier-stokes equations for an incompressible fluid. *Bulletin of the American Mathematical Society*, 73(6):928–931, 1967.

[9] Claire MY Claassen, Shafiul Islam, EAJF Peters, Niels G Deen, JAM Kuipers, and Maike W Baltussen. An improved subgrid scale model for front-tracking based simulations of mass transfer from bubbles. *AIChE Journal*, 66(4):e16889, 2020.

[10] C Dapogny, C Dobrzynski, P Frey, and A Froelhy. Mmg platform: Robust, open-source and multidisciplinary software for remeshing. *GitHub repository*, 2020.

[11] Jean-Marc Delhaye. Jump conditions and entropy sources in two-phase systems. local instant formulation. *International Journal of Multiphase Flow*, 1(3):395–409, 1974.

[12] Olivier Desjardins, Vincent Moureau, and Heinz Pitsch. An accurate conservative level set/ghost fluid method for simulating turbulent atomization. *Journal of computational physics*, 227(18):8395–8416, 2008.

[13] Donald A Drew. Mathematical modeling of two-phase flow. *Annual review of fluid mechanics*, 15(1):261–291, 1983.

[14] Ronald P Fedkiw, Tariq Aslam, Barry Merriman, and Stanley Osher. A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method). *Journal of computational physics*, 152(2):457–492, 1999.

[15] Marcus Herrmann. A sub-grid surface dynamics model for sub-filter surface tension induced interface dynamics. *Computers & Fluids*, 87:92–101, 2013.

[16] Takehiro Himeno, Daizo Sugimori, Katsutoshi Ishikawa, Yutaka Umemura, Seiji Uzawa, Chihiro Inoue, Toshi-nori Watanabe, Satoshi Nonaka, Yoshihiro Naruo, Yoshifumi Inatani, et al. Heat exchange and pressure drop enhanced by sloshing. In *47th AIAA/ASME/SAE/ASEE joint propulsion conference & exhibit*, page 5682, 2011.

[17] Cyril W Hirt and Billy D Nichols. Volume of fluid (vof) method for the dynamics of free boundaries. *Journal of computational physics*, 39(1):201–225, 1981.

[18] Muhammad Irfan and Metin Muradoglu. A front tracking method for direct numerical simulation of evaporation process in a multiphase system. *Journal of Computational Physics*, 337:132–153, 2017.

[19] Romain Janodet, Carlos Guillamón, Vincent Moureau, Renaud Mercier, Ghislain Lartigue, Pierre Bénard, Thibaut Ménard, and Alain Berlemont. A massively parallel accurate conservative level set algorithm for simulating turbulent atomization on adaptive unstructured grids. *Journal of Computational Physics*, 458:111075, 2022.

[20] Emmanuel Labourasse, Delphine Lacanette, Adrien Toutant, P Lubin, S Vincent, O Lebaigue, J-P Caltagirone, and P Sagaut. Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests. *International journal of multiphase flow*, 33(1):1–39, 2007.

[21] V. Moureau, P. Domingo, and L. Vervisch. Design of a massively parallel CFD code for complex geometries. *Comptes Rendus Mécanique*, 339(2-3):141–148, 2011.

[22] Stanley Osher and James A Sethian. Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations. *Journal of computational physics*, 79(1):12–49, 1988.
[23] Yohei Sato and Bojan Ničeno. A sharp-interface phase change model for a mass-conservative interface tracking method. *Journal of Computational Physics*, 249:127-161, 2013.

[24] Changxiao Shao, Kun Luo, Min Chai, Haiou Wang, and Jianren Fan. A computational framework for interface-resolved DNS of simultaneous atomization, evaporation and combustion. *Journal of Computational Physics*, 371:751-778, 2018.

[25] Sébastien Tanguy, Thibaut Ménard, and Alain Berlemont. A level set method for vaporizing two-phase flows. *Journal of Computational Physics*, 221(2):837-853, 2007.

[26] A Toutant, E Labourasse, O Lebaigue, and Olivier Simonin. DNS of the interaction between a deformable buoyant bubble and a spatially decaying turbulence: a priori tests for the two-phase flow modelling. *Computers & Fluids*, 37(7):877-886, 2008.

[27] Grétar Tryggvason, Bernard Bunner, Asghar Esmaeeli, Damir Juric, N Al-Rawahi, W Tauber, J Han, S Nas, and Y-J Jan. A front-tracking method for the computations of multiphase flow. *Journal of computational physics*, 169(2):708-759, 2001.

[28] Andre Weiner and Dieter Bothe. Advanced subgrid-scale modeling for convection-dominated species transport at fluid interfaces with application to mass transfer from rising bubbles. *Journal of Computational Physics*, 347:261-289, 2017.

[29] Samuel WJ Welch and John Wilson. A volume of fluid based method for fluid flows with phase change. *Journal of computational physics*, 160(2):662-682, 2000.