# Thermal boundary layer modelling for bubbles at saturation: A posteriori analysis

Mathis Grosso, Guillaume Bois, Adrien Toutant

► To cite this version:

Mathis Grosso, Guillaume Bois, Adrien Toutant. Thermal boundary layer modelling for bubbles at saturation: A posteriori analysis. International Journal of Heat and Mass Transfer, 2025, 242 (1), pp.126744. 10.1016/j.ijheatmasstransfer.2025.126744. hal-04956101

HAL Id: hal-04956101

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

Submitted on 19 Feb 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
<figure><img src="image_2.png" /></figure>

# Thermal boundary layer modelling for bubbles at saturation: A posteriori analysis

Mathis Grosso<sup>a,b,*</sup>, Guillaume Bois<sup>a</sup>, Adrien Toutant<sup>b</sup>

## ARTICLE INFO

**Keywords:**

Two-phase flows

Heat transfer

Front-Tracking

Ghost Fluid Method

Boundary-layer modelling

Sub-grid modelling

## ABSTRACT

This study investigates different temperature and flux coupling strategies in Direct Numerical Simulations (DNS) of bubbles at saturation, employing local one-dimensional thermal boundary layer sub-resolutions. Specifically, a laminar radial sub-resolution (LRS) near the interface is employed to address challenges in capturing sharp temperature variations, which is crucial for liquid-vapour heat transfer correlations. State-of-the-art techniques use analytical profiles to capture very thin boundary layers around single-rising objects for very high Prandtl or Schmidt numbers. The original approach proposed in Grosso et al. (2024) relies on a more general embedded sub-resolution still applicable at low Prandtl numbers and coarse grids. To accurately integrate the sub-layer variations into the CFD grid, the literature recommends using the sub-grid profiles to evaluate the Eulerian face fluxes instead of correcting cell temperature. From experience, it avoids excessive flux leakage from the sub-layer region at high Prandtl numbers. The present article investigates these coupling methods while proposing adaptations for thick boundary layers and very coarse grids in the context of LRS. Two test cases, pure diffusion acting around a sphere and a single rising bubble configuration, are explored, measuring heat flux at the interface and its transmission to the fluid domain serving as figures of merit for each coupling method. In low Prandtl bubbly flows (Pr<sub>f</sub> ≤ 5), and on coarse and affordable grids (< 20 cells per bubble diameter), temperature coupling is found to be more stable though not conservative compared to flux coupling approaches. On the other hand, classical flux coupling strategies can exhibit artefacts and introduce potential instabilities with LRS. To overcome such problems, an improved local flux balance approach is proposed, demonstrating both robustness and efficiency in predicting and transmitting interfacial flux across the tested thermal layers' thickness ranges.

## 1. Introduction

The development and application of numerical tools for investigating nucleation, droplet physics, and particle-laden flows are highly active and rapidly evolving fields of research [1–9]. In the numerical solving procedures of these problems, one common challenge is the accurate treatment of the interfaces where heat, mass and momentum transfers occur between two fluids. While experimental measurements in these areas have primarily been conducted in simple configurations, a significant gap remains in the precise quantification of complex phenomena, particularly interfacial heat and mass transfers, which are crucial for understanding evaporation and condensation rates in boiling flows. Recent research efforts have focused on quantifying the effects of condensation [10] and evaporation [7] on the thermal boundary layer and investigating how various dimensionless numbers influence bubble heat transfer in vertical pipes [11]. However, collective phenomena

such as thermal wake interactions between inclusions remain inadequately understood. Some progress has been made in exploring similar interactions in particle-laden flows [12], yet further investigations are needed.

Interface-resolved simulations of two-phase flows are facing difficulties in capturing precise variations of quantities. The physics of boiling or chemical transport phenomena imply a mix of large and micro-scale regions which cause numerical bottlenecks with conventional methods. In the case of bubbly flows, we can identify a micro-scale region that corresponds to the boundary layer developing around the moving bubble interface and separating the two phases. On the contrary, large-scale regions correspond to the carrying phase where both convective and diffusive transfers are less intense.

To conduct multiple simulations, there is a critical need to enhance numerical approaches, particularly in regions near interfaces. Such advancements will improve the accuracy of simulations and facilitate

<footnote><sup>a</sup> Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides, Gif-sur-Yvette, 91191, France</footnote>

<footnote><sup>b</sup> PROMES - CNRS (UPR 8521), Université de Perpignan Via Domitia, Perpignan, 66100, France</footnote>

<footnote>* Corresponding author.
E-mail addresses: mathis.grosso@cea.fr (M. Grosso), guillaume.bois@cea.fr (G. Bois).</footnote>
the derivation of interfacial models for Euler–Euler simulations [13,14].

The resolution of these flows needs very high refinements in the vicinity of the interfaces where the heat and mass transfers take place. Unfortunately, the conventional methods do not allow such local refinements. The interface capturing or tracking methods (such as Volume of Fluid or Front-Tracking, later abbreviated FT) usually used to study complex bubbles and droplets' physics are very robust but have certain limitations. By privileging the freely moving aspect of the interface, the interfacial region is often spread over one or multiple cells on a uniformly refined grid. If compressive algorithm [15] or geometric advection [16] of extensive quantities can handle mass and energy conservation, it does not guarantee a good convergence of interfacial quantities, especially fluxes.

To address these challenges, solvers like Basilisk [17] have integrated Adaptive Mesh Refinement (AMR) techniques to enhance computational efficiency and accuracy. For instance, AMR has been successfully applied to study droplet evaporation [18]. However, the time step in such simulations remains constrained by the finest grid size, limiting computational efficiency. In a different approach, Vreman [19] explored the use of an overlapping spherical grid to model flow around rigid particles. This technique provided accurate results, particularly in capturing the velocity field near interfaces. More recently, in the context of phase change, Bures et al. [20] proposed an innovative method to correct mass flow rate predictions near interfaces. They employed a Robin boundary condition to sharpen the temperature field, ensuring that the temperature gradient is more accurately represented. Finally, the sub-grid scale approach introduced in [21] and recently used by Maarek et al. [22] employs analytical solutions to describe thin scalar boundary layers and has marked a significant advancement in the precise modelling of heat and mass transfer near interfaces at high Schmidt or Prandtl numbers.

For convective flows characterised by moderate Prandtl numbers (Pr ≥ 1), the sub-grid scalar region remains under-resolved compared to the hydrodynamics, as described by the Batchelor scale λB [23], which relates the diffusion of scalar quantities to molecular diffusion, characterised by the Kolmogorov length scale η = λB√PrI. At the discrete level, the local heat flux measured at the interfaces is more likely to be underestimated. On the other hand, the fluxes, computed in the interface's vicinity and transmitted to the carrying fluid domain can be underestimated or overestimated depending on the boundary layer thickness variations within interfacial cells. Our previous work has established a sub-grid refinement methodology [24]. This method alleviates issues encountered when the sub-grid model relies on analytical solutions, including challenges with boundary conditions and the limited representation of sub-grid physics used to derive the solution. It has been partially inspired by the work on velocity and thermal boundary layers around solid objects and referred to as Turbulent Boundary Layer Equations (TBLE) implemented in LES solvers [25–27].

For an isolated bubble configuration or homogeneous bubble swarms, the mixed cells are a minority in the computational domain. As a consequence, the numerical approaches have been focusing on correcting or substituting the conventional Finite-Volume numerical schemes in interfacial regions. At the discrete level, two regions need specific refinements. On the one hand, the interfacial quantities are usually ignored in the discretisation but rather computed a posteriori thanks to neighbouring phase values or ghost extensions [28,29]. Specific methods are often introduced in the conventional Computational Fluid Dynamics (CFD) resolution to enhance the measurement of these interfacial quantities (e.g. mass transfer evaluation $\dot{m}_v$ ∝ $\|\nabla T \cdot n_F\|_F$ at the exact interface position). On the other hand, other numerical methods focus on accurately predicting the liquid-vapour fluxes over the fixed grid in the interface vicinity by off-centring operators [30]. Some complex frameworks inherently answer these two stakes by replacing the CFD solver in the boundary layer. It thus ensures an accurate interfacial flux prediction to be transferred to the pure phase domain [21,31,32].

In more detail, the numerical enhancements brought for these two distinct purposes (interfacial and staggered grid fluxes) are the following:

(a) **Interfacial flux and temperature computation:** Locally, it is necessary to align with the normal interface direction and compute the normal variations of a given quantity (e.g. concentration or temperature). Then, the jump between two-sided temperature gradients can be accurately computed to predict the local heat transfer. In phase-change simulations, the rate of phase change, $\dot{m}_v$, is closely related to the interfacial gradient, as shown in Eq. (1) where $L^{vap}$ corresponds to the latent heat of vaporisation.

$$ \dot{m}_v L^{vap} = \lambda_I \nabla T_I \cdot n_F |_I - \lambda_v \nabla T_v \cdot n_F |_I \quad (1) $$

It can be done using various methods:

* The Ghost-Fluid extension Method (GFM) [29]. Authors use Ghost-Fluid extensions to establish two distinct continuous fields across the interface [28,33]. During the process one-sided interfacial gradients ($\lambda_k \nabla T_k \cdot n_F |_I$) are obtained and the interface temperature $T_I$ is prescribed.

* A simple Finite-Difference scheme aligned with the interface normal [6], as shown in Fig. 1.a.

* A sharpening method for the temperature field is used to accurately predict the mass flow rate and bubble growth [20].

* An advanced sub-grid fitting of an analytical profile, as illustrated in Fig. 1.b [21,31].

* A dedicated local sub-resolution technique, such as the Laminar Radial Sub-resolution (LRS) approach from our previous work [24], depicted in Fig. 1.c.

(b) **Fixed grid fluxes computation in mixed cells:** In interface resolved simulations, one of the difficulties is to capture the fluxes within the mixed cells. Faces of mixed cells can be purely liquid (l) or vapour (v) and are characterised by their surfaces $S_k$ ($k \in \{v; l\}$). Faces that are cut by the interfaces are referred to as mixed faces (m) and are separated into wetted surface areas $S_{k,m}$ as shown in Fig. 2. The interface $S_I$ and the union ($S_k \cup S_{k,m}$) of both surfaces delimitates the sub-volume of phase k. Vectors pointing outward the volume $V_I$ are denoted n

Ideally, the objective is to discretise and solve the energy balance in the liquid sub-volume, which is written in a conservative form as given by Eq. (2):

$$ \begin{aligned} & \rho_l C_{pl} \frac{\partial}{\partial t} \left( \iiint_{V_l} T_l dV \right) + \rho_l C_{pl} \iint_{S_l \cup S_{l,m}} T_l u \cdot n dS \\ &= \iint_{S_l \cup S_{l,m}} \lambda_l \nabla T_l \cdot n dS - \iint_{S_I} \lambda_l \nabla T_l \cdot n_I dS \end{aligned} \quad (2) $$

However, discretising the fluxes across faces using pure phase quantities is difficult. Considering the One-fluid approach [34], the mixed cell temperature consists of a mix of both phase temperatures which is not coherent with the face fluxes computation that requires one-sided temperature values and the level of accuracy is often unsatisfying [14]. As stated in (a), the GFM allows a separate treatment of each one-sided temperature field but without consideration for mixed cell surfaces $S_{l,m}$. In boiling flows, we assume a constant vapour and interface temperature $T^{sat}$. Face fluxes are then approximated straightforwardly with classical single-phase operators across pure liquid faces $S_l$ delimiting the liquid sub-volume as shown in Fig. 2. During each iteration, this volume is treated as a buffer region where discrete energy increments are assumed to be zero ($\partial_t \rho_l C_{pl} T_l = 0$), and values are obtained solely by extrapolating instantaneous quantities. On the other hand, the Two fluid approaches aim to resolve the velocity or scalars separately in the two phases in
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** (a) Interfacial gradient computation in a probe's manner (Finite-Difference) by [6] (b) Sub-grid fitting of an analytical function (erfc) reproduced and adapted from [32]. Comparison between the CFD grid fluxes calculation and the one derived from the local sub-grid solution. (c) Laminar Radial Sub-resolution approach [24]: the Eulerian field (E) is coupled to the thermal boundary layer (BL) through $T_p = T_{CFD}$ (●). $\nabla T \cdot n_f|_\Gamma$ (◆) is the interfacial heat transfer. $T_c$ is the cell centre temperature. $T_f$ and $\nabla T \cdot n_f|_f$ are computed to obtain the convective and diffusive faces fluxes.</figcaption></figure>

a cut-cell manner i.e. two unknowns are resolved in the mixed cells [30,35,36]. They are located at the respective sub-volume barycentres. Eq. (2) is thus discretised with fidelity but the numerical schemes employed in the immediate vicinity must be offset to use pure fluid phase quantities (see Fig. 3). Discrete points are no longer aligned. Moreover, a given scalar quantity is not discretised at the interface but rather computed a posteriori using the gradient value on both sides of the interface along with the specific transmission condition, such as a continuous heat flux, as discussed in [14,37], or a discontinuous heat flux, as shown in Eq. (1). In that case, the interfacial condition is still weakly imposed. In the case where two-sided Dirichlet Boundary conditions must be imposed (such as saturation temperature), this approach still requires highly refined grids or additional sub-grid models [32].

(c) **Combination of (a) and (b):** A last approach consists of combining the two previous steps by substituting the fixed grid solver locally with a uni-directional sub-grid solution thus achieving a complete two-way coupling. The sub-grid model evaluation is fed with fixed grid quantities (extensive quantities or local interpolated quantities). The sub-grid model then transmits temperature at the mixed cell centre [38] (or fluid sub-volume

barycentres) or fluxes to the neighbouring cell faces [32]. It should be noticed that the sub-grid methodology introduced in [21] and enhanced by Weiner et al. [32] is applied on the concentration field. Without the presence of chemical reactions, it is considered a passive scalar. It is then comparable to solving the temperature equation without phase change.

The previously developed quasi-static methodology also referred to as LRS underwent initial testing on post-processed temperature fields from Direct Numerical Simulations (DNS) at steady-state. It successfully demonstrated its ability to rectify interfacial gradients and accurately predict the surrounding temperature field [24]. This initial validation marked the completion of the first step in the process (a).

In the scope of achieving a complete two-way coupling between the newly sub-resolved region and the CFD grid solution (b), the LRS approach has been integrated into TRUST/TrioCFD a DNS-FT framework developed at the French Atomic Energy Commission (CEA) to study boiling flows. To achieve such coupling, the temperature field or the thermal fluxes (convective and diffusive) must be transmitted to the underlying CFD grid.

Several grid coupling strategies have been exposed in the literature and should be adapted to answer our needs. The most straightforward
<figure><img src="image_4.png" /><figcaption>Fig. 2. Surfaces delimiting a liquid control volume ($V_l$) in mixed cells for exact energy balance.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Off-centred schemes applied to scalar Y for the prediction of face fluxes $\beta$ (from [30]).</figcaption></figure>

coupling strategy consists of correcting the mixed cell centre temperature using the sub-resolved uni-directional temperature profile. In other words, the cell centre temperature becomes the local temperature value of the fluid as opposed to integrated quantities. The fluxes on the mixed cells' faces are then computed using this corrected temperature [31]. This approach can be employed on both sides of the interface i.e. in both the vapour and liquid phases. It can be directly coupled with the LRS approach. If the boundary layer is thick enough, it is expected to work well but when the latter becomes thinner, the temperature and fluxes are known to vary strongly across the cell and its faces which could cause a misprediction of the fluxes.

As a consequence, authors such as Weiner et al. [32] have first recommended using the sub-grid analytical concentration profile to compute face flux values along with the cell concentration in a Finite-Volume manner. They went further by integrating the uni-directional concentration profile and derivative over the faces to correct their fluxes. Their profile, described by an analytical function, representative of the boundary layer physics, is fitted and integrated analytically along with geometric variables that depend on the interface to cell intersections (Plane portion of interface).

Additionally, Weiner et al. discussed and investigated exhaustively the different aspects of the coupling approach to take maximum advantage of the method while mitigating numerical artefacts: e.g. treating inconsistencies in nearly pure cells, ensuring bounded concentration values and effects caused by micro to macro timescale difference. The phenomenon of unboundedness arises when an excessive amount of energy (or mass) is extracted from pure liquid cells, leading to a temperature decrease beyond a critical threshold $\Delta T$ (in sub-cooling configuration). All these considerations allow for reducing the dispersion in the prediction of heat flux and ensure the numerical stability of this coupling.

The mapping between two discretisations i.e. from 1D to 3D and vice versa is complex to solve. The flux-based methods rely on the evaluation of a heat flux density ($W m^{-2}$) from the one-dimensional

<figure><img src="image_6.png" /><figcaption>Fig. 4. Schematic of the two-way coupling.</figcaption></figure>

profiles and weighted by a given cell face surface. Any one-dimensional resolution considers uni-directional variations within the mixed cell. It thus provides a unique heat flux value by interface portion. By weighting this value by the surface of the interface portion, it leads to the prediction of the heat leaving the interface (in W). Due to the quasi-static hypothesis, and considering the solution of a one-dimensional problem, the heat injected at the interface must be balanced by the heat transmitted to the domain. A wrong prediction of the local heat flux at the liquid face will inevitably lead to a wrong energy balance and the introduction of an error that could destabilise the sub-resolution itself.

To align with state-of-the-art techniques, a rich investigation of the face fluxes coupling approaches has been undertaken to answer our needs and constraints (One-fluid solver for momentum, Ghost-Fluid Method for temperature and interfaces tracked by Front-Tracking) and understand details about each leverage that could enhance the method. After a brief description and reminder of the DNS-LRS approach in Section 2.1 and Section 2.2, the two main classes of approaches (temperature and fluxes) are described more precisely in Section 2.3. Some additional ideas are proposed, starting from the observation that some interface portions do not contribute to the overall flux coupling strategy. Finally, a novel approach, based on the evaluation of a local error and corrective term to ensure a complete and local flux balance is proposed.

To test our methodology, two configurations are simulated: a static spherical bubble diffusing heat and a single bubble rising in an initially quiescent fluid (Section 3.2 and Section 3.3). It provides valuable insights into the effect of each coupling method. Their results are discussed in Section 4.

Finally, the overall conclusions emphasise the best coupling approach that handles an improved interfacial heat transfer prediction while granting conservative transmission of the fluxes to the rest of the domain. More complex simulation setups and the potentiality of the method are briefly discussed as perspectives in boiling and engineering applications.

## 2. Laminar Radial Subresolution: sub-grid temperature and flux transmission strategies for a complete two-way coupling

The term LRS stands for Laminar Radial Sub-resolution. This term covers the entire two-way coupling procedure as described in Fig. 4. The Eulerian ghost-fluid temperature field $T^*$ from the CFD resolution is transferred to the boundary layer solver through boundary conditions (B.Cs) and adequate source terms. They are obtained using tri-linear interpolations. Then the refined temperature field $T$ and its derivative are solved and sent back to the CFD field to correct the face fluxes on the fixed grid.

### 2.1. Direct Numerical Simulations-Front Tracking framework

The Direct Numerical Simulations and Front-Tracking framework (abbreviated DNS-FT) have been described in [24]. A fixed Eulerian grid is used to solve for scalar and vector quantities along with Lagrangian markers to represent the interfaces and compute geometrical
<figure><img src="image_7.png" /><figcaption>Fig. 5. Coarse ($D_b/\Delta = 16$) to fine ($D_b/\Delta = 64$) mesh comparison. The Lagrangian mesh locating the interfaces must be refined along the Eulerian grid resolution.</figcaption></figure>

properties. Two uniform grid refinements are shown in Fig. 5. The intersection of the Lagrangian mesh and the Eulerian grid is utilised to compute an exact liquid indicator function ($I_l$) on the fixed grid. The velocity is solved in a one-fluid manner (arithmetic averaged of both viscosities and densities thanks to the indicator function $I_l$). This approach is robust but fails when the viscosity ratio becomes significant. The investigated viscosity ratio remains within reasonable limits so that hydrodynamics are always resolved.

Incompressibility is considered for both fluids and mass conservation and is ensured through a pressure projection step solved by an in-house multi-grid solver. The mean curvature $\kappa$ is computed accurately on the Lagrangian mesh using a technique based on differentials of volume and surface [39]. It is provided for surface tension forces calculation (with constant surface tension coefficient $\sigma$) and sub-grid resolution in order to scale the osculating sphere frame of reference. Convective (Quick, 3rd order) and diffusive (Centred, 2nd order) Cartesian grid operators for the temperature are fully explicit. The time scheme is first-order explicit (Euler).

The original temperature resolution in the interface vicinity is ensured by an in-house Ghost-Fluid Method (GFM) embedding a sub-grid pure diffusion model. It has proved to be robust, but it has a slow grid convergence on interfacial heat transfer. The solution embeds hypotheses which are not representative of the physics in the near interface region. The temperature is resolved only on the liquid side because the vapour is considered to be at saturation temperature. The advection-diffusion equation Eq. (3) in mixed and liquid cells is computed from the Ghost Fluid extended temperature field denoted $T^*$ while $T_{CFD}$ designates the Eulerian temperature field before the extrapolation procedure. In pure liquid cells, we have $T^* = T_{CFD}$.

$$ \frac{\partial T^*}{\partial t} + \nabla \cdot (\mathbf{u}T^*) = \alpha_l \nabla \cdot (\nabla T^*) \quad (3) $$

We solve for the temperature difference between the fluids and the saturation temperature i.e. $T^* \leftarrow (T^* - T^{sat})$. Thus, at the interfaces, the saturation condition writes $T^{sat} = 0$. If the liquid is subcooled, $T^* < 0$, whereas $T^* > 0$ for an overheated fluid.

In Fig. 5, the comparison between coarse and fine meshes is illustrated. The convergence of the temperature field is obtained on the finer grid while LRS is applied to the coarse grid.

## 2.2. LRS overall methodology

**Description of the method.** The first step of the Laminar Radial Sub-resolution technique has been extensively introduced in [24]. For the sake of clarity, a few key elements are reminded in the following paragraphs. Some additional details relative to transient simulations and geometrical aspects that have not been discussed yet are given.

**A simplified governing equation for the temperature in the boundary layer.** In the boundary layer region, a pseudo-three-dimensional temperature equation, evaluated in the radial direction, and written in a spherical frame of reference is used (see local and instantaneous Eq. (4a) in axisymmetric configuration and its numerical matrix formulation Eq. (4b)). The temperature at the interface $T_r$ is constant and equals the saturation temperature $T^{sat}$ at each time step (n). The probe tip temperature $T_p$ is obtained by interpolation of the CFD grid temperature. The radial coordinate $r$ corresponds to the radial distance to the osculating sphere centre. The osculating radius at the interface $R_o$ is obtained through the mean curvature $\kappa$ computed from the discontinuous interface $R_o = 2/\kappa$ ($\kappa > 0$ for a convex interface). A modified radial coordinate $\tilde{r} = r - R_o$ such that $\tilde{r} = 0$ at the interface completes the definitions of the local sub-problem variables. Considering only the mean curvature is sufficient for the bubble shapes aimed. The error introduced by this hypothesis is confined to the term $2/r = 2/(R_o + \tilde{r})$ in front of the diffusion term $\partial_r T$ in Eq. (4a). We prove that this error is proportional to the squared difference between the two radii of curvature, $|\tilde{r}(R_1 - R_2)|^2$ (see [24], App. B). The error is zero at the interface where $\tilde{r} = 0$.

$$ \frac{\partial^2 T}{\partial r^2} + \frac{2}{r} \frac{\partial T}{\partial r} - \frac{u_r}{\alpha_l} \frac{\partial T}{\partial r} = \frac{u_\theta}{\alpha_l} \frac{1}{r} \frac{\partial T^*}{\partial \theta} - \frac{1}{r^2} \frac{\partial^2 T^*}{\partial \theta^2} \quad (4a) \\ \underbrace{\boxed{D^r}}_{\text{Source terms } S^\theta} \quad \underbrace{\boxed{C^r}}_{\text{Source terms } S^\theta} \quad \underbrace{\boxed{C^\theta}}_{\text{Source terms } S^\theta} \quad \underbrace{\boxed{D^\theta=0}}_{\text{Source terms } S^\theta} $$

$$ [A^{D'} + A^{C'}]^{(n)} T_{\text{LRS}}^{(n)} = [S^{\theta}]^{(n)} \quad (4b) $$

$$ T_{\Gamma}^{(n)} = T^{sat} = 0 \quad (4c) $$

$$ T_p^{(n)} = T_{CFD}^{(n)} \quad (4d) $$

The temporal term Eq. (5), accounting for bubble rising velocity $u_{rise}$ and radial deformation $u_r^d|_r$ of the interface, is neglected. This consideration makes the sub-resolution quasi-static. The term is exactly zero at the interface, and we expect it to cancel out rapidly with the rapid settlement of the layer in the regions of the bubble top and equator.

$$ \frac{DT}{Dt} = \frac{\partial T}{\partial t} + (u_{rise} + u_r^d|_r \cdot \tilde{e}_r) \cdot \nabla T \quad (5) $$

A uniform probe length and discretisation are set for the overall simulation. For each sub-problem, the matrices for convective and diffusive operators are scaled by the velocity components and local curvature. A length of two cell diagonals ($\tilde{r} = 2\sqrt{3}\Delta$ with uniform grid spacing) is used to ensure that the interpolation of the temperature $T_p$ as a B.C is computed with single-phase values (i.e. the CFD temperature $T^*$ is not affected by the extrapolation procedure during the iteration (n)). This B.C remains constant during the sub-resolution process. LRS solution is then used to enhance fluxes calculation on the CFD grid. Several methods are later described in Section 2.3. The Ghost-Fluid temperature field $T^*$ is incremented to iteration (n+1) using explicit convection and diffusion operators. During this step, the Eulerian grid temperature, $T^{*(n+1)}$, in liquid cells overlapped by the probes can evolve freely.

Velocity components $u_r$, $u_\theta$ are obtained through interpolations of the one-fluid velocity field $\mathbf{u}$. The velocity components are then corrected to align with the bubble frame of reference (see Appendix A.1 for more details). The source terms consist of the tangential convective and diffusive effects. Tangential convection relies on an extended temperature inherited from the GFM whereas the tangential diffusion term $D^\theta = 1/r^2\partial_{\theta\theta}^2 T$ is neglected throughout this work in agreement with previous studies in [24]. In fact, second-order derivatives are poorly reconstructed on coarse grids using this technique. Both source terms could be represented by a surrogate model in further work. By
correcting the fluxes transmitted to the CFD grid, the temperature near the interface is altered, which in turn affects the boundary condition and the source terms interpolated at the probe. A complete two-way coupling between CFD and LRS is achieved. Some dispersion of the interfacial values has been observed in the first a priori tests due to the reconstruction of the tangential derivatives. This issue arises from the computation of derivatives from the Ghost-Fluid temperature field. GFM does not strongly enforce the condition $\partial_\theta T|_R = 0$ at the interface; its accuracy is sensitive to the interface position within the mixed cell. It is particularly visible in the equator region because of the magnitude of the convective term. However, direct assessment of tangential derivatives from the GFM temperature has provided promising results without any specific treatment. As a consequence, we have decided to put effort into the coupling step, to prove the feasibility of the LRS a posteriori i.e. directly within a DNS solver. Further work to bring the method to a higher level of accuracy and smoothness is still of interest.

**Strength and challenges of the method.** The LRS does not ensure perfect conservation of quantities in the first step of the two-way coupling procedure because the energy contained within the boundary layer is not constrained by the underlying Eulerian grid quantities ($E \to BL$ in Fig. 1.c and Fig. 4). This first step completely replaces the solver in the boundary layer region in contrast to analytical methods developed for thin boundary layers, which are fitted according to the underlying grid energy. The method is quasi-static by neglecting the temporal term while updating the boundary conditions at the probe tip. The validity of this assumption has been confirmed, as demonstrated in the time-dependent diffusion problem outlined in Section 4.1.

The strength of LRS lies in its ability to capture the interfacial gradient and better approximate the heat flux transmitted to the carrying fluid. The conservation of the heat flux transmitted from the boundary layer to the CFD field is a major concern and is aimed at studying several coupling strategies.

## 2.3. Coupling strategies within mixed cells

Some coupling strategies related to one-dimensional sub-grid solutions are available in the literature. Their adaptation to suit the LRS implementation is then discussed. Some details about the implementation are given in Appendix A.3.

The fluxes transmitted to the CFD field are corrected by taking advantage of the highly refined LRS profile (see Fig. 1.c). In each approach, the uni-directional sub-grid profiles ($T, \partial_\theta T...$) are used in interpolation/extrapolation i.e. using the normal distance $d_\perp$ between the interface and a given neighbouring point (cell and face centres) as shown in Fig. 6.

From the literature, two major coupling strategies have been found:

(A) A coupling strategy based on a local cell temperature correction [21]. From this newly computed temperature, fluxes crossing the mixed cells are better approximated.

(B) A coupling strategy based on fluxes evaluation on pure liquid faces [32]. Each face flux is computed using the sub-grid profile (in our case the quantities predicted by the LRS) (see Fig. 8).

### 2.3.1. Temperature coupling strategy (A)

The most obvious coupling strategy consists of correcting the temperature in mixed cells using the LRS while preserving both convective and diffusive Eulerian schemes straddling the interface. The cell centre position requires to be projected onto the probe so that the radial distance $d_\perp$ can be measured. $d_\perp$ is then used to evaluate the temperature $T_c$ from the LRS profile (Fig. 6). This is equivalent to considering the CFD solution is uni-axial within the mixed cell.

This approach is not inherently conservative but exhibits convergence properties. This family of coupling strategies will be referred to as case (A) in the rest of the paper.

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Correction of the cell temperature using the LRS profile at coordinate $\bar{r} = d_\perp$. Probes have been offset tangentially for visibility. The red line marks the reconstructed equivalent interface in each cell.</figcaption></figure>

LRS necessitates quantities from the CFD field that are computed by the GFM. Another possibility is to switch from LRS to GFM in case of poor interfacial predictions by LRS. In presence of convection and at steady-state, the coarse grid simulation employing GFM always underestimates the interfacial gradient. LRS is expected to enhance this interfacial gradient in comparison to the under-resolved GFM but the opposite is observed. Then, as the GFM has proved to be very robust both the interfacial heat flux and cell temperature are locally overridden when $|VT \cdot n_r|_r^{GFM} > |VT \cdot n_r|_r^{LRS}$. Further motivations for this functionality are given in Appendix A.2. This particular coupling case is introduced as ($A_{corr,gfm}$).

### 2.3.2. Fluxes coupling strategy (B)

**Challenges in flux coupling.** Weiner et al. [32] discussed with precision the importance of correcting fluxes for an accurate transmission between thin layers and coarser CFD grids. We observe both insufficient flux transmission and significant overprediction of the transmitted fluxes, depending on the Prandtl number $Pr_l$, in approaches based on temperature correction. Specifically, for low Prandtl numbers, the transmitted flux is underestimated, whereas, as the boundary layer thins, we significantly overpredict the flux on the pure liquid side as later exposed in Section 4.2.1 and Appendix E.1.

In our simulations, moderate Prandtl numbers, lead to a thin boundary layer at the bubble's top which reaches its minimum thickness just above the equator; then, it significantly expands below the equator. Based on the most advanced methodology exposed in the literature, the fluxes coupling strategy appears to be the best candidate to transfer heat to the interface vicinity with minor loss. It has been utilised in both interface-capturing and boundary-fitted methods [32,40].

The heat contributions ($F_f$ in W) crossing cell faces $f$ are obtained by computing heat fluxes ($F''_f$ in W m$^{-2}$) using sub-grid profiles and weighting them with the surface of the face (see Eq. (6)).

$$F_f = F''_f S_f \qquad (6)$$

The error in evaluating $F_f$ arises from applying the flux $F''_f$ to the whole surface $S_f$.

The 1D sub-resolution is only discretised on the probe without consideration for the liquid control volume contained in the mixed cell (Finite-Difference). Consequently, computing the face fluxes without enforcing their balance makes the LRS approach non-conservative during the coupling step, as will be explained by Eq. (15).
<figure><img src="image_9.png" /><figcaption>Fig. 7. Flux coupling problem illustrated using a continuous interface (a) and a discrete interface (b). The remaining discrete face portions ( ) need to be associated with a given sub-problem. (c) Flux coupling simplifications to overcome the problem of bijection.</figcaption></figure>

There are challenges in accurately associating each local sub-problem with the correct portions of the pure liquid faces. This complexity must be addressed to establish a coherent global solution that maintains flux balance.

Within the framework of a continuous and seamlessly integrated interface set within a Cartesian grid (Fig. 7.a), one can partition the interface into segments. Surrounding the interface, a convex envelope can be delineated by the liquid faces closest to the interface. A bijective relationship can be established between each segment of the interface and its corresponding liquid faces to delimit coherent control volumes. At the discrete level, correcting pure liquid faces may raise interrogation as there is no clear bijection between a given liquid face and the sub-problem candidates to correct the latter (see Fig. 7.b).

Restraining the mixed cell fluxes computation to the immediate local sub-problem is a way to solve this issue (Fig. 7.c).

These properties are analysed in depth in Section 4. In addition, in the case of GFM-based solvers, the part of the energy that flows tangentially to the interface is lost through mixed faces. These mixed cells are problematic to ensure a perfect flux balance and global liquid energy conservation because their energy conservation equation is not resolved in time. Yet, they receive fluxes from pure liquid cells. Therefore, only the liquid faces are available to transmit fluxes instantaneously to the liquid field.

Such concerns about correcting the liquid face fluxes to ensure very precise conservation are raised in the present article. Weaknesses of each coupling method are identified to isolate the best candidate. Non-conservative flux coupling has a huge impact on the thermal wake. A proposal to ensure exact conservative transmission of the interfacial heat flux towards the liquid domain is later exposed in Section 2.3.3.

The fluxes coupling necessitates two steps that are reminded.

* As the fluxes crossing liquid portions of mixed faces $f_{l,m}$ are not computed, it is necessary to identify pure liquid cell neighbours and faces $f_l$.

* Both convective $F_{f_l}^C$ and diffusive $F_{f_l}^D$ fluxes whom expressions are given by Eq. (7a) and Eq. (7b) are corrected.

<figure><img src="image_10.png" /><figcaption>Fig. 8. Fluxes coupling strategies on pure liquid faces. The interface portion (×) is ignored in (B). The red flux (→) could potentially receive a contribution from the unused portion ($B_{avg}$).</figcaption></figure>

$$F_{f_l}^C = u_l T \cdot n_{f_l} \Big|_{f_l} S_{f_l} \qquad (7a)$$

$$F_{f_l}^D = \alpha_l \nabla T \cdot n_{f_l} \Big|_{f_l} S_{f_l} \qquad (7b)$$

It is necessary to understand that in available flux coupling approaches, a liquid face is only associated with one sub-grid problem. Some interface-attached methods have been developed in the case of Front-Tracking [41,42] but do not apply fluxes on faces as they rather introduce source terms on the Eulerian grid instead of overriding its fluxes.

This second coupling approach and its variations are referred to as case (B). A discrete numerical integration of the face fluxes ($B_{int}$) is also considered. It theoretically improves the coupling when there are strong variations of the flux over the face due to the variation of the distance. In Section 4.1.1 and Section 4.2.1, it is demonstrated that the improvement is minimal. The diffusive face fluxes are evaluated using a combination of the radial gradient refined with LRS ($\partial_r T$) and the tangential one obtained using the extended temperature field ($\partial_\theta T^*$). Further information about its implementation and details about fluxes calculation are given in Appendix A.3.

Weighting the face fluxes contributions ($B_{avg}$). The computation of the face fluxes uses conventionally a single candidate (B), i.e. the single sub-problem associated with a mixed cell [32,43]. From Fig. 8, the inclusion of the unconsidered interface portion (×) contributing to the liquid face ( ) is of interest to achieve a coherent flux transmission.

To address this challenge, flux averaging to account for several neighbouring probes has been tested and is referred to as ($B_{avg}$) whereas the last coupling method employing a colinearity-based weighting procedure is referred to as ($B_{wavg}$). Details about the geometrical computations as well as the neighbours averaging procedure are given in Appendix A.4.

### 2.3.3. A strategy to overcome fluxes unbalance

Previous approaches ($B_{int}$, $B_{avg}$ and $B_{wavg}$) apply a heat flux on a given liquid face $f_l$. Applications of these methods were not entirely successful on coarse grids. Locally, the hypothesis stating that the one-dimensional field is valid over the entire cell introduces error.

Initially, the LRS temperature field is resolved only at the probe and lacks an explicit control volume. During the coupling step, constraining the one-dimensional solution to the liquid control volume shown in Fig. 2 and Fig. 7.c results in a failure to maintain strict flux balance over mixed cells in agreement with the temporal evolution of energy
<figure><img src="image_11.png" /><figcaption>Fig. 9. Surfaces of flux transmission (without a two-scalar approach).</figcaption></figure>

in that cell. The flux coupling stage must address this issue.

Since temporal terms are not considered in mixed cells (as values are continuously overridden without temporal consideration), energy balance is not naturally satisfied. As a result, heat transfer from the interface may be partially lost through mixed faces and within the grey buffer region represented in Fig. 7.

Despite the significant enhancement of the one-dimensional temperature profile and the temperature gradient at the interface brought by the first stage of LRS (finer discrete resolution), we initially found challenging to utilise this improvement during the flux coupling stage. A numerical adaptation based on the quasi-static hypothesis and energy conservation has been explored to ensure a conservative transmission of the flux leaving the interface towards the liquid domain.

**Ensuring a conservative coupling.** On the basis that the LRS approach is quasi-static, the flux leaving the boundary layer region at the continuous level should exactly match the flux leaving the interface (over a closed portion). This property is difficult to achieve at the discrete level because of the reasons stated in Section 2.3.2.

In Fig. 9, it can be seen that LRS raises a challenge to transmit instantaneously the global interfacial heat $Q_Γ$ to the liquid domain through the convex envelop formed by the red pure liquid faces (heat $Q_{f_l}$). The coupling strategy must bypass efficiently the buffer region in grey where the energy is not resolved in time. This issue does not occur in existing sub-grid analytical strategies in two-scalar simulations because they compute energy increment within mixed cell [32,43]. The two-scalar approach introduces a degree of freedom and relaxes the constraint on the control volume to ensure naturally the local flux balance (see Eq. (2)).

In LRS, by constraining the local control volumes around the probes to the grey liquid sub-volumes in Fig. 9, the one-dimensional solutions are not well recombined to form a coherent three-dimensional quasi-static solution verifying a global flux balance. This statement is expressed through Eq. (8), where $\circ_c$ denotes a mixed cell.

$$ \underbrace{\sum_c \mathcal{R}_{\Gamma,c}^{\prime\prime} \,^{(d)} S_{\Gamma,c}}_{\text{Interfacial flux}} \neq \underbrace{\sum_c \sum_{f_l} \mathcal{R}_{f_l}^{\prime\prime} \,^{(d)} S_{f_l}^{c}}_{\text{Liquid face fluxes}} \quad (8) $$

To elaborate the last strategy, the instantaneous flux balance across the buffer region has been exploited to formulate a compensation term.

**Exhibiting the flux balance on a control volume.** The starting point is the local quasi-static temperature equation in Finite Volume form, assuming a constant liquid thermal diffusivity $\alpha_l$ (see Eq. (9)). To determine fluxes at the liquid faces of the mixed cell, the equation is integrated over the liquid sub-volume $V_l$ (see Eq. (10) and Fig. 2).

$$ \nabla \cdot (\mathbf{u}T) = \nabla \cdot (\alpha_l \nabla T) \quad (9) $$

$$ \iiint_{V_l} \nabla \cdot (\mathbf{u}T) dV = \iiint_{V_l} \nabla \cdot (\alpha_l \nabla T) dV \quad (10) $$

Using the divergence theorem, the fluxes can be split into several faces contributions (see Eq. (11) and Fig. 10.a where $n$ is the normal pointing outward the liquid control volume $V_l$ while $n_Γ$ points towards $V_l$). The flux balance is similar to that found in two-scalar solvers, except that temporal variations of temperature within the liquid sub-volume are neglected [32]. A diffusive contribution comes from the interface surface ($S_Γ$) and both convective and diffusive contributions come from the union ($S^*$) of the pure liquid faces ($S_{f_l}$) and liquid wetted portions of mixed faces ($S_{f_l,m}$). Knowing that $n = -n_Γ$ and $T_Γ = 0$ on the surface $S_Γ$, we get:

$$ \iint_{S^*} \mathbf{T} \mathbf{u} \cdot \mathbf{n} dS = \iint_{S_Γ} \alpha_l \nabla T \cdot \mathbf{n} dS + \iint_{S^*} \alpha_l \nabla T \cdot \mathbf{n} dS \quad (11) $$

Finally, terms can be gathered as in Eq. (12) and further developed in Eq. (13).

$$ \iint_{S^*} (\alpha_l \nabla T - \mathbf{T} \mathbf{u}) \cdot \mathbf{n} dS = \iint_{S_Γ} \alpha_l \nabla T \cdot \mathbf{n}_Γ dS \quad (12) $$

$$ \underbrace{\iint_{S_{f_l}} (\alpha_l \nabla T - \mathbf{T} \mathbf{u}) \cdot \mathbf{n} dS}_{\text{C: Pure liquid faces}} = \underbrace{\iint_{S_Γ} \alpha_l \nabla T \cdot \mathbf{n}_Γ dS}_{\text{D: Interface portions}} - \underbrace{\iint_{S_{f_l,m}} (\alpha_l \nabla T - \mathbf{T} \mathbf{u}) \cdot \mathbf{n} dS}_{\text{E: Liquid wetted surfaces}} \quad (13) $$

At the discrete level and within a mixed cell, the energy conservation Eq. (13) writes:

$$ \sum_{f_l} (\mathcal{R}_{f_l}^{D} \,^{(d)} - \mathcal{R}_{f_l}^{C} \,^{(d)}) = \alpha_l \left. \frac{\partial T}{\partial r} \right|_{\Gamma} S_{\Gamma,c} - \underbrace{\sum_{f_l,m_l} (\mathcal{R}_{f_l,m_l}^{D} \,^{(d)} - \mathcal{R}_{f_l,m_l}^{C} \,^{(d)})}_{\text{Hyp. } E^{(d)}=0} \quad (14) $$

The surface of the local interface portion is denoted $S_{\Gamma,c}$. Convective and diffusive fluxes are denoted by $\mathcal{R}_{f_l}^C$ and $\mathcal{R}_{f_l}^D$ on a given pure liquid face $f_l$. Finally, convective and diffusive fluxes across liquid-wetted portions of mixed faces are expressed with the subscript $f_l,m_l$.

Eq. (14) only ensures local energy conservation over the liquid layer trapped in mixed cells. To facilitate the implementation, it is commonly assumed that $E^{(d)} = 0$, as in coupling case (B).

The main strategy to enforce global and local conservation is to extend coupling case (B) to correct the evaluation of the discrete fluxes (superscript (d)) by a compensation term $e^{(d)}$ redistributed on the pure liquid cell faces as follows:

$$ \sum_{f_l} \left[ (\mathcal{R}_{f_l}^{C} \,^{(d)} - \mathcal{R}_{f_l}^{D} \,^{(d)}) + e_{f_l}^{(d)} \right] = \alpha_l S_{\Gamma,c} \left. \frac{\partial T}{\partial r} \right|_{\Gamma} \quad (15) $$

This method requires handling the cases when there are less than three pure liquid faces within the mixed cell i.e. $\sum_{f_l} 1 < 3$. We need as many faces with different Cartesian directions as there are non-zero components in the interface normal vector. This effectively distributes the interfacial heat flux. In that case, the error contributions should be transmitted or relocated to neighbouring mixed cells and applied on their respective pure liquid faces to transmit the heat to the carrying phase properly (see Fig. 10.b). All the complexity lies in the establishment of a directional weighting coefficient $w_{f_l}$ taking into account the geometry and the principal direction of the transfers. It is required to reformulate the single equation for the error Eq. (15) into a solvable system:

$$ e_{f_l}^{(d)} = w_{f_l} \left[ \alpha_l S_{\Gamma,c} \left. \frac{\partial T}{\partial r} \right|_{\Gamma} - \sum_{f_l} \left( \mathcal{R}_{f_l}^{C} \,^{(d)} - \mathcal{R}_{f_l}^{D} \,^{(d)} \right) \right] \quad (16) $$

A strategy to determine these weighting coefficients ($w_{f_l}$) is described in the following paragraph.
<figure><img src="image_12.png" /><figcaption>Fig. 10. Reallocating fluxes and applying directional error terms: (a) Flux balance over the liquid sub-volume. (b) Transfer of the compensation term to the neighbouring pure liquid faces. (c) Directional splitting, relocation and application of the compensation term.</figcaption></figure>

**Directional splitting and application of the compensation term.** The relocation and correction of the fluxes should be guided by the physics of transfers as much as possible. The transfers primarily occur perpendicular to the interface. The contributions to fluxes are either relocated through mixed faces or applied to pure liquid faces whose normals $n_{f_i}$ align with the direction of the interface normal (see Eq. (17)).

$$n_{f_i} \cdot n_{f_i} \ge 0 \qquad (17)$$

Considering a normal face vector $n_{f_i}$ aligned with the Cartesian grid direction $\hat{e}_i$, the directional weight straightforwardly writes:

$$w_{f_i} = \frac{n_{f_i} \cdot n_{f_i}}{\sum_j n_{f_i} \cdot n_{f_j}} \qquad (18)$$

Here, index $i$ designates the face across which the flux is applied or relocated, which can be either a pure liquid face or a mixed one. Index $j$ denotes the total number of flux contributions within a cell (Eq. (17)) leads to a maximum of three possible contributions).

The calculation of the weighting coefficients is realised in two steps. For instance in Fig. 10.c, the compensation term $E^{(d)}$ splits for cell $c_1$ into $e_{f_1}^{(d)}$ and $e_{f_2}^{(d)}$ through the calculation of two weighting coefficients denoted $w_{f_1}$ and $w_{f_2}$. $e_{f_1}^{(d)}$ is directly applied on face $f_1$ whereas $e_{f_2}^{(d)}$ is relocated to cell $c_2$ across face $f_2$. The flux $e_{f_2}^{(d)}$ must then be applied to the available liquid faces of the neighbouring cell $c_2$. Using a criterion based on the alignment of the face normal with the normal to the interface portion $n_{f_i}$, the weighting values are determined as $w_{f_{2,1}} = 0$ and $w_{f_{2,2}} = 1$ for this example. It restrains the flux transmission to close pure liquid faces. All these rules are generalisable for any mixed cell configurations.

This coupling approach is referred to as (B$_{cons}$) throughout the rest of the paper.

### 2.3.4. Summary and nomenclature of the coupling approaches

An exhaustive study of coupling strategies is performed throughout the paper. Table 1 summarises the tested strategies and underlines

their pros and cons identified from the literature and observed by the authors. Appendix A.5 provides indications on the algorithm for each coupling method.

## 3. Numerical setups for the assessment of various two-way coupling strategies

### 3.1. Test cases objectives

A priori assessment of the LRS in [24] has been realised on post-processed quantities coming from steady rising bubble DNS at the final time-step. Transmission of information from the boundary layer towards the Eulerian fields (BL → E in Fig. 1.c) presents some obstacles.

Some numerical instabilities can occur due to the coupling of two discretisations. The discretised probe's domain is not necessarily aligned with the Cartesian direction and the quasi-static hypothesis may introduce inconsistencies in the field. By feeding back the probe with a boundary condition, it is important to keep the Eulerian field free from these artefacts that may be cascading and that may ultimately overflow the simulation with unphysical values. There is a clear retroaction and any error made during the coupling step will affect the interface regions.

To separate effects, a first thermal simulation performed with a stationary interface is investigated. Its ability to mimic a time-dependent problem is appraised. Then, the LRS is tested on a dynamic simulation where it exhibits its capacity to predict the interfacial thermal quantities properly. Simulation cost is then discussed. A relative error on the Nusselt number of less than 5% has been achieved between the reference simulation and the coarse grid simulations embedding the LRS.

Various quantities at several locations are post-processed. In particular, some interfacial quantities, and radial quantities attached to the bubble are collected to facilitate the comparison with the refined solution. Integral quantities have also been tracked and are compared to correlations at steady-state. The interfacial temperature gradient has been a major metric in assessing the ability of LRS in the previous one-way coupling step [24]. It is still important to track this value depending on the coupling approach because the prediction of this quantity depends on the temperature in the interface vicinity. It will then provide important insights into the effect of each coupling method. As heat is constantly injected at the interface, an exact conservation between the liquid outlet flux and the interfacial heat flux must be measured at steady-state. To check this statement for each coupling strategy, two-dimensional slices are post-processed behind the bubbles. Slices must follow the bubble over time to perform the flux balance in the bubble frame of reference.

### 3.2. Thermal diffusion acting around a sphere

The first case consists in a static 3D bubble of diameter $D_b$ = 2mm diffusing heat in a periodic domain of dimensions (5 × 5 × 5) $D_b$ and composed of quiescent fluid with a thermal diffusivity $\alpha_l$ = $1.15 \times 10^{-7} \text{ m s}^{-2}$.

The exact analytical solution for the temperature and its derivative is well-known and reminded by Thiam et al. [12] (see Eq. (19) where $R_b$ stands for the bubble radius).

$$T = T_{\infty} + (T_{\Gamma} - T_{\infty}) \frac{R_b}{r} \left[ 1 - \erf \left( \frac{r - R_b}{2\sqrt{\alpha_l t}} \right) \right] \qquad (19a)$$

$$\frac{\partial T}{\partial r} = (T_{\infty} - T_{\Gamma}) \left[ \frac{R_b}{r} \frac{1}{\sqrt{\alpha_l \pi t}} e^{-\left(\frac{r - R_b}{2\sqrt{\alpha_l t}}\right)^2} \right] \qquad (19b)$$

$$+ \frac{R_b}{r^2} \left[ 1 - \erf \left( \frac{r - R_b}{2\sqrt{\alpha_l t}} \right) \right] \qquad (19c)$$
<caption>**Table 1**</caption>

Summary of the tested coupling approaches, their theoretical and numerical pros, and cons. All except (Bcons) are theoretically non-conservative when employed in the context of LRS.

<table><thead><tr><th></th><th>Cases</th><th>Description</th><th>Pros</th><th>Cons</th><th>Section</th><th>Ref.</th></tr></thead><tbody><tr><td rowspan="6">Temperature</td><td>GFM</td><td>Sub-grid pure diffusion model</td><td>Robust</td><td>Sub-grid physics validity is mesh dependant</td><td>Section 2.1</td><td>[24]</td></tr><tr><td>A</td><td>LRS temperature profile</td><td>Accuracy in sub-grid temperature variations</td><td rowspan="2">Directional aspects are lost in the coupling phase</td><td>Section 2.3.1</td><td rowspan="2">[38]</td></tr><tr><td>Acorrgfm</td><td>Substitute A by the GFM in regions where predictions are inaccurate</td><td>Better accuracy in the bottom and equator regions</td><td>Section 2.3.1 Appendix A.2</td></tr><tr><td>B</td><td>Single flux value at the face centre</td><td>Ease of implementation</td><td>Altered the boundedness of the solution</td><td>Section 2.3.2</td><td rowspan="4">[32,43]</td></tr><tr><td>Bint</td><td>Discrete flux integration</td><td>Accuracy in sub-grid face flux variations</td><td>Geometrical computations and implementation</td><td>Section 2.3.2 Appendix A.3</td></tr><tr><td>Bbis</td><td>Averaging flux contributions</td><td>Smoothing the flux prediction</td><td>Geometrical and parallel computations and implementation</td><td>Section 2.3.2</td></tr><tr><td>Flux</td><td>Bbiscol</td><td>Weighting average of flux contributions</td><td>Prioritise the closest probe to the face</td><td>Geometrical and parallel computations and implementation</td><td>Section 2.3.2 Appendix A.4</td><td>New</td></tr><tr><td></td><td>Bcons</td><td>Correct flux evaluation to ensure the interface to liquid face flux balance</td><td>Conservative</td><td>Arbitrary correction, relocation rules and implementation</td><td>Section 2.3.3</td><td></td></tr></tbody></table>

The analytical solution is applicable for diffusion acting around a sphere in an infinite domain, while the simulation involves a periodic array of bubbles due to the boundary conditions. Despite this difference, the analytical solution effectively predicts the temperature field for a certain duration, as long as the spheres remain sufficiently separated. Although more suitable boundary conditions exist [44], our solver is optimised for periodic swarms and has proven sufficient for evaluating our methodology.

Eq. (19) is used to initialise a smooth temperature field in the 3D cases (CFD grid nodes) at an initial time $t_{ini}$ ensuring the temperature field is sufficiently developed. Particular attention has been given to temperature initialisation in this purely diffusive case. If the thermal field is not sufficiently developed, the quasi-static hypothesis introduces a bias in the fluid domain and is not compensated over time.

Fig. 11 shows the initial interpolated profile (blue) in comparison to the one resolved with the LRS approach (orange). The temperature along the probe does not match exactly the analytical solution because the surface mesh is extremely coarse. The normal vectors of the interface segments do not perfectly align with those of an ideal, perfectly spherical bubble, for which the analytical solution is computed. The quasi-static temperature field computed by the LRS is a good approximation of the time-dependent field from the simulation start.

The global heat transfer is measured to assess the coupling strategies. The liquid-to-vapour heat transfer for a perfect sphere corresponds to the Nusselt number $Nu^{th}$ whose expression is the following:

$$Nu^{th} = - \frac{D_b \int_T \lambda_l \nabla T \cdot n_T dS}{\lambda_l \Delta T} = \frac{D_b (T_T - T_\infty)}{\Delta T} \left[ \frac{1}{R_b} + \frac{1}{\sqrt{\pi \alpha_l t}} \right] \quad (20) \\ = 2 + \frac{2R_b}{\sqrt{\pi \alpha_l t}}$$

### 3.3. Single rising bubble in an initially quiescent fluid

Once the bubble is in movement, the LRS ability to correct the gradient may vary. The orientation, position and rate of displacement

<figure><img src="image_13.png" /><figcaption>Fig. 11. Radial temperature profiles interpolated and resolved at physical time $t = 3.68s$. The underlying field has been initialised using Eq. (19a).</figcaption></figure>

of the interface portions within the mixed cells could cause dispersion of the predicted values. This is especially true on very coarse grids. First, the major stake of this configuration is to measure the global heat transfer over the bubble surface and check that it is not perturbed by these aspects. Secondly, the local temperature prediction over the probe length can be observed and compared to fine grid. Predictions of the LRS temperature field reveals the influence of the tangential source term ($\propto \partial_\theta T$).

At steady-state, the Nusselt number correlation developed by Feng et al. (see Eq.23a in [45]) is used to assess the ability of the LRS to predict the correct interfacial heat flux. At moderate Reynolds number, there is no analytical solution so the local temperature fields post-processed from coarse grids can only be compared to the reference simulation performed at $\Delta/D_b = 64$.
<caption>**Table 3**</caption>

Geometrical and physical parameters for the set of simulations are defined by $Ar^* = 50$ and $Pr_l \in \{1; 2.5; 5\}$. These conditions correspond to a liquid-vapour system at saturation under 155 bar, with the liquid's thermal conductivity adjusted to modify the Prandtl number from 1 to 5.

<table><caption>**Table 4**</caption><thead><tr><th>Parameters</th><th>Notation</th><th>Unit</th><th>Values - Range</th></tr></thead><tbody><tr><td rowspan="3">Geometrical</td><td>Bubble diameter</td><td>D<sub>b</sub></td><td>mm</td><td>2</td></tr><tr><td>Cells per bubble diameter</td><td>D<sub>b</sub>/&Delta;<sub>x,y,z</sub></td><td>-</td><td>[16, 32, 64] for Pr<sub>l</sub> &in; {1; 2.5; 5.}</td></tr><tr><td>Domain size</td><td>&Omega;<sub>x,y,z</sub></td><td>D<sub>b</sub></td><td>[4, 4, 6]</td></tr><tr><td rowspan="4">&Omega;<sub>l</sub></td><td>Viscosity</td><td>&mu;<sub>l</sub></td><td>Pas</td><td>6.82&times;10<sup>-5</sup></td></tr><tr><td>Density</td><td>&rho;<sub>l</sub></td><td>kg m<sup>-3</sup></td><td>594.4&times;10<sup>0</sup></td></tr><tr><td>Conductivity</td><td>&lambda;<sub>l</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>[0.69&times;10<sup>0</sup>, 0.276&times;10<sup>0</sup>, 0.138&times;10<sup>0</sup>]</td></tr><tr><td>Specific heat capacity</td><td>C<sub>&rho;l</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>10 110&times;10<sup>0</sup></td></tr><tr><td rowspan="2">&Omega;<sup>g</sup></td><td>Viscosity</td><td>&mu;<sub>v</sub></td><td>Pas</td><td>2.3&times;10<sup>-5</sup></td></tr><tr><td>Density</td><td>&rho;<sub>v</sub></td><td>kg m<sup>-3</sup></td><td>101.9&times;10<sup>0</sup></td></tr><tr><td>L<sub>l</sub></td><td>Surface tension</td><td>&sigma;</td><td>N m<sup>-1</sup></td><td>9.79&times;10<sup>-5</sup></td></tr><tr><td>&Omega;</td><td>Gravity</td><td>g</td><td>m s<sup>-2</sup></td><td>4.97&times;10<sup>-3</sup></td></tr></tbody></table>

Comparison in terms of computational costs between coarse and fine grids single rising bubble simulations at several Prandtl numbers. The simulation cost is approximated for a duration of 10 s.

<table><thead><tr><th rowspan="2">Parameters</th><th colspan="3">Coarse</th><th colspan="2">Fine</th></tr><tr><th>1</th><th>2.5</th><th>5</th><th>1</th><th>2.5</th><th>5</th></tr></thead><tbody><tr><td>Resolution (D<sub>b</sub>/&Delta;)</td><td colspan="3">16</td><td colspan="2">64</td></tr><tr><td>Number of cells (M)</td><td colspan="3">0.4</td><td colspan="2">25.2</td></tr><tr><td>Number of CPUs (-)</td><td colspan="3">8 (2<sup>3</sup>)</td><td colspan="2">512 (8<sup>3</sup>)</td></tr><tr><td>CPU time per step (s)</td><td colspan="3">&asymp; 1</td><td colspan="2">&asymp; 2.3</td></tr><tr><td>Time-step (s)</td><td colspan="3">1&times;10<sup>-3</sup></td><td colspan="2">1&times;10<sup>-4</sup></td></tr><tr><td>Simulation cost (h.CPU)</td><td colspan="3">&asymp; 22.2</td><td colspan="2">&asymp; 33,000</td></tr></tbody></table>

proves advantageous, as excessively fine simulations encounter challenges with the scalability of geometric routines for the interface.<fnref n="3" /> The degrees of freedom (d.o.f.) of the interfaces are proportional to the CFD mesh refinement, which significantly slows down smoothing routines.

Finally, the interface-related stability time-step [24] is eight times larger in the coarse simulation. On the largest boundary layer, the reduction in computational cost is around a factor of 1300 which justifies the interest of the method. The method becomes more interesting as the thickness of the boundary layer decreases. The comparison of the coarse to fine simulations' costs for the single rising bubble case is summarised in Table 4.

## 4. A posteriori evaluation of the boundary-layer coupling strategies

($A_{corr,gfm}$●) and (GFM<sub>ref</sub>■) are not exposed in the spherical diffusion case. Case (GFM<sub>ref</sub>) is replaced by an analytical solution while ($A_{corr,gfm}$) is not interesting as it could potentially replace the overall LRS prediction due to the spherical symmetry of the problem (i.e. $\nabla T \cdot n_r|_{T^{GFM}} \ge \nabla T \cdot n_r|_{T^{A}}, \forall \theta \in [-\pi/2; \pi/2]$).

### 4.1. Thermal diffusion acting from a sphere

#### 4.1.1. Integral heat transfer measurements

Heat transfer at the interface ($Nu_r$). The total vapour to liquid heat transfer at the interface has been plotted over time in Fig. 13. Eq. (19a) is used as a reference (black line ■). The Nusselt number relative error is plotted in Fig. 14 to emphasise the difference between each method.

At the beginning of the simulation, a small offset between the theoretical Nusselt number value remains for every case. It is attributed to the quasi-static approach whose outputs underestimate the interfacial

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Nusselt number values over time for the different coupling strategies. Initial physical time is non-zero.</figcaption></figure>

gradient. Then, the Nusselt number decreases rapidly until reaching a rather constant decreasing rate towards the asymptotic value of 2 (see Eq. (20)). Every prediction is globally localised around the theoretical solution which confirms our trust in each implementation. Some of them overestimate the interfacial transfer while others underestimate it. The Nusselt number error reaches a plateau around $t \approx 20$s and then increases because of the periodic boundary conditions.

At this level of discretisation $D_b/\Delta = 9$, the Ghost Fluid Method performs well (●). The Sub-grid pure diffusion analytical model replicates the underlying physics. The relative error stagnates at 1%.

The most straightforward temperature coupling method (A-▲) is very close to the reference solution. It performs well against the original GFM reaching a relative error of -1%. On the contrary, the flux-based coupling strategies (B-★) and (B<sub>int</sub>+<sup>+</sup>) lead to an overestimation of the Nusselt number value above 2%.

To explain why the Nusselt number can be overestimated, it should be remembered that heat is constantly injected at the interface. The Nusselt number overestimation arises from an underestimation of the fluxes to be transmitted to the CFD grid. It then results in the misprediction of the temperature field increment $\Delta T^{(n+1)}$ and the boundary condition at the probe's tip in the next iteration. $T_p^{(n+1)}$ is then lower than the theoretical temperature $T_p^{th}$, and so forth.

<footnote n="3">³ Geometric routines maintain sufficient quality, such as regularisation, mesh balancing and volume preservation.</footnote>
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Nusselt number relative error $\frac{Nu_f - Nu_{qh}}{Nu_{qh}} \times 100$ (%). The Nusselt number is evaluated at the interface $\Gamma$.</figcaption></figure>

By averaging the diffusive fluxes (see case $B_{\text{avg}}$-●), there is a clear effect on the transmitted fluxes which is now overestimated from the beginning of the simulation. Once in the pseudo-constant regime of transfer, the bias introduced in the liquid field stops increasing. The relative error is slightly greater than -2%. Using a colinearity weighting to smooth the fluxes values looks beneficial as ($B_{\text{wavg}}$) nearly reaches the results of (A-▲) and shows a smaller error during the initial transient. In summary, all methods perform similarly to predict interfacial heat transfer in pure diffusion.

Finally, the coupling method ($B_{\text{cons}}$♦) gives the best results while ensuring conservation. It reaches an error of -0.8%. The relocation algorithm is coherent with the physics of transfer and corrects the fluxes in the regions where the interface normal is not aligned with the Cartesian directions.

Heat transfer on pure liquid faces ($Nu_{f_i}$). The fluxes leaving the convex envelope and the relative error are plotted for each coupling case in Fig. 15 and Fig. 16. Cases (GFM) and (A) are almost superimposed. They rely on temperature correction within the mixed cell. A small modification of the cell centre temperature has thus a small effect on the transmitted flux.

Cases (B) and ($B_{\text{int}}$) underestimate the flux transmitted to the domain by nearly 5% whereas ($B_{\text{avg}}$) is close to the reference. Adding a particular colinearity weighting ($B_{\text{wavg}}$) is beneficial as it constitutes the best prediction of the fluxes leaving the sub-layer region (less than 0.5%). Finally, enforcing the flux conservation in case ($B_{\text{cons}}$) gives similar results to temperature coupling cases. From the overall heat transfer point of view (Fig. 14 and Fig. 16), case ($B_{\text{cons}}$) constitutes the best compromise to measure the interfacial heat flux while transmitting it accordingly. Additional integral measurements are provided in Appendix B. To confirm these observations at the local scale, measures of transfers have been obtained and are exposed in Section 4.1.2 and Section 4.1.3.

### 4.1.2. Local interfacial heat transfer measurements

As it has been done in [24], it is interesting to study the dispersion of the values predicted by the LRS. In Fig. 17, the interfacial gradient $\nabla T \cdot n_{\Gamma}|_{\Gamma}$ has been plotted against the spherical coordinates $\theta$ at simulation time $t = 20$ s. At a given $\theta$ range, there are several points corresponding to various azimuthal positions. In pure diffusion, dispersion is observable and is mainly dictated by the local interface position

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Nusselt number computed from the fluxes leaving the boundary layer region i.e. the convex envelope delimited by the pure liquid faces.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Nusselt number relative error $\frac{Nu_{f_i} - Nu_{qh}}{Nu_{qh}} \times 100$ (%). The Nusselt is evaluated at the liquid faces $f_i$.</figcaption></figure>

within a cell. As expected, the LRS applies the interfacial boundary conditions more precisely within the mixed cell, thus reducing this dispersion. The sharp application of the saturation B.C at the interface thus gives a uniform gradient prediction around the bubble surface in case of spherical symmetry. It avoids artefacts appearing in Cartesian directions. On the other side, the flux feedback locally influences the probe's end B.C. and is thus also responsible for a certain uniformity in the gradient prediction.

The GFM involves multiple smoothing iterations when computing the field's extension [24]. Still, it appears that the temperature gradient is affected by the interface position. Spherical symmetry is lost and the gradient is minimum (in absolute value) when the interface is almost aligned with the Cartesian direction i.e. when $\theta = \{-\pi/2, 0, \pi/2\}$. The same observations can be drawn for coupling cases (B-★), ($B_{\text{int}}$-+). The flux prediction is greater when the colinearity between the interface normal and the face normal decreases (diagonal directions). Although
<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Dispersion of the normal temperature gradient values predicted at the interface. The LRS is compared to the gradient computed from the GFM.</figcaption></figure>

case ($B_{\text{cons}}^{\text{◆}}$) is conservative, it presents dispersed values. It is not surprising as the flux relocation is arbitrary and inevitably decreases the smoothness of the CFD field solution. In the presence of convection like in [31], dispersion is primarily observed at the bubble top and is linked to the interface location within the cell. This effect becomes particularly significant near the bubble's top because it corresponds to a region of intense heat transfer.

On the contrary, this figure emphasises a major advantage of the flux averaging methods ($B_{\text{avg}}^{\text{◆}}$) and ($B_{\text{wavg}}^{\text{▶}}$). They can capture the spherical symmetry of the problem. Similar observations are made with coupling case (A-▲). Imposing a zero<sup>th</sup> order quantity ensures the flux smoothness and the easiest transition between the two discretisations especially in regions where they are not aligned.

Finally, these figures demonstrate that the temporal behaviour is accurately represented.

### 4.1.3. Temperature and gradient profiles in the interface vicinity

The Nusselt number gives access to an integral measure of the error i.e. $Nu \propto \int_S \lambda_l \nabla T \cdot n_r |_{n_r} dS$. The simulations can give decent interfacial quantities outputs although the Eulerian field is badly impacted by the LRS. As a consequence, it is fundamental to study the local temperature profile and its derivatives by taking advantage of the test case's analytical solution (see Fig. 18.a and Fig. 18.b).

The worst probe's position has been chosen i.e. the probe is pointing outward at the bubble in the diagonal position. In that case $(\theta, \phi) \approx (45, -45)^\circ$.

The temperature profile resolved on the probe matches well with the analytical solution whatever the coupling cases. In cases (B-★) and ($B_{\text{int}}^{+}$), the temperature is slightly lower which translates into an underestimation of the flux transmitted to the CFD grid. These results confirm the same counter-intuitive observations discussed in Section 4.1.1. Finally, the temperature is very well predicted for the rest of the cases (A-▲, $B_{\text{cons}}^{\text{◆}}$, $B_{\text{avg}}^{\text{◆}}$ and $B_{\text{wavg}}^{\text{▶}}$). The difference in temperature between the one interpolated from the GFM and the one solved by the ($B_{\text{cons}}$) case is almost indistinguishable.

Concerning the radial temperature gradient profile Fig. 18.b at the interface, observations are slightly different. Diagonal cells are cut so that they have less pure liquid faces. Flux relocation acts thus principally in the diagonal cells which benefits the solution obtained through

$(B_{\text{cons}})$. Fluxes predicted by (A), ($B_{\text{avg}}$) and ($B_{\text{wavg}}$) are very close to ($B_{\text{cons}}$). Finally, flux coupling methods (B and $B_{\text{int}}$) overestimate the flux both locally and globally. It is interesting to see that the radial profile reconstructed from the GFM extended temperature field (●) gives a linear piece-wise profile. It means that higher-order extrapolation techniques [33] are not likely to enhance sufficiently the temperature field without refining the grid. This profile is post-processed but does not serve directly in the GFM solver.

## 4.2. Rising bubble in an initially quiescent fluid

### 4.2.1. Integral interfacial heat transfer measurements

Heat transfer at the interface ($Nu_r$). The interfacial Nusselt number measured for each coupling case at $Pr_l = 1$ is shown in Fig. 19. The reference calculation is above the value predicted by Feng et al. [45]. The same conclusions were already drawn in [24]. It is due to our restricted domain size which tends to slightly compress the boundary layers. The LRS is a function of the velocity, the prediction of LRS should match the reference solution and cannot compensate for the effects of the domain's size.

Although the Peclet number is moderate, the original Ghost Fluid approach (GFM) applied on a coarse grid fails to predict the Nusselt number. The relative error between the coarse and fine grids reaches -17% at the steady state. Compared to other coupling approaches it is one of the worst predictions. It is getting worse at higher Prandtl numbers (see Appendix E for more results).

Local and integral flux coupling cases (B-★) and ($B_{\text{int}}^{+}$) are superimposed which means that the integration of the face fluxes is superfluous. The variations of the solution over the faces to correct is thus low or linear (i.e., in that case, the evaluation of the flux at the face centre equals the integral value). The wrong prediction of the interfacial heat transfer is once again due to a wrong transmission of the fluxes to the carrying fluid. The relative error stays between -17 to -16%. It acts back on LRS and explains poor predictions. At this level of coarseness, the error made in the liquid face flux prediction is rather due to the flux balance problem exposed in Section 2.3.2.

Averaging the fluxes contributions is beneficial for the evaluation of the interfacial heat transfer (see ◀ and ▶). The type of averaging has less effect on the prediction of fluxes compared to the pure diffusion case. The relative error drops to -11% and -12%, respectively.

Finally, three methods provide results below 10% error. The temperature coupling (A-▲) enhances the prediction of the interfacial flux passing above the correlation of Feng et al. [45]. Resorting to the maximum gradient between GFM and LRS (see $A_{\text{corr,gfm}}^{\text{◆}}$), the prediction nearly reaches the reference solution (GFM<sub>ref</sub><sup>■</sup>) up to a 2.5% relative error. It is also interesting to notice that the transient is well predicted. However, these two coupling cases do not inherently ensure conservation (Fig. 20).

Case ($B_{\text{cons}}$) is the only flux coupling strategy that enhances notably the Nusselt number prediction achieving an error below 5% while ensuring conservation. To validate this coupling strategy among others, it is necessary to study in detail the fluxes' conservation as well as the local dispersion of interfacial quantities.

Heat transfer on pure liquid faces ($Nu_f$). To emphasise the error made on the transmitted fluxes, the interfacial to liquid face Nusselt numbers relative error ($Nu_r - Nu_f$)/$Nu_r$ is plotted in Fig. 20. The fluxes computed according to the extended temperature field of (GFM-●) lead to a relative error of -10% at steady state. The relative error stabilises on a plateau. Increasing the refinement mitigates the energy loss during the transmission step (see GFM<sub>ref</sub><sup>■</sup>).

Flux coupling strategies (B-★ and $B_{\text{int}}^{+}$) present a relative error of less than +5% in conserving the interfacial heat flux transmitted to the liquid domain (see Fig. 19). While this is acceptable, these cases fail to accurately predict the interfacial heat flux shown in the previous figure Fig. 19, resulting in an incorrect amount of heat being transmitted to
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Post-processing of radial quantities at $t \approx 20$ s, $(\theta, \phi) \approx (45, -45)^\circ$. (a) Radial temperature profile over the probe length. (b) Radial temperature gradient $\frac{\partial T}{\partial r}$. In the case of GFM (●), only the interfacial gradient $\left.\frac{\partial T}{\partial r}\right|_r$ is directly corrected and so not visible in this figure.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Overall Nusselt number $Nu$ at $Pr_f = 1$. and $Re_b \in \{67; 72\}$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Relative error between the interfacial transfers $Nu_f$ and the liquid face fluxes $Nu_{f_l}$ in the convective case at $Pr_f = 1$.</figcaption></figure>

the liquid domain.

Averaging the fluxes contributions in cases (◀) and (▶) does not lead to greater accuracy. The error increases up to a range of +12 to +14 %. These methods gave interesting results in diffusion but have revealed deceiving in convective flow. The error increases mainly during the transient which confirms this weakness. Since the temperature and flux coupling strategies (A), ($A_{corr,gfm}$), and ($B_{cons}$) yield more promising results, immediate improvements for (◀) and (▶) are not pursued.

If ($A_{corr,gfm}$) works extremely well in predicting the interfacial transfer $Nu_f$, but it fails to transmit the correct amount of energy to the CFD grid.

Surprisingly, (A-▲) gives a very low conservation error at steady state (< 0.5%), while the heat flux value measured at the interface is accurate. The combination of these two metrics makes this simple coupling strategy usable. However, given the drawbacks observed in the

pure diffusion case, we expect that this seemingly good performance results from a balance between over and under-estimated regions.

Finally, ($B_{cons}$) has been formulated such that $Nu_f$ and $Nu_{f_l}$ match. The error is not exactly zero because a single relocation is performed. Some rare contributions cannot be relocated in a single step as there are no pure liquid faces within their immediate neighbours. The relative error has proved to never exceed ±0.2%. It is more conservative than the refined GFM solution (GFM$_{ref}$).

Once again, ($B_{cons}$) answers both criteria. The interfacial heat transfer is notably improved on a coarse grid while the extra benefit of conserving the transmitted fluxes is achieved.

To keep results clear, and based on supplementary observations (available in Appendix E.1), only (A), ($A_{corr,gfm}$) and ($B_{cons}$) are commented on in the next part. These three methods provide results below 10% error while the others are less accurate and proved to be unstable.
<figure><img src="image_22.png" /><figcaption>**Fig. 21.** Local Nusselt number $D_b \nabla T \cdot n_r|_F / \Delta T$ at $Pr_l = 1$, $Re_b \in \{67; 72\}$ at steady-state.</figcaption></figure>

The discussion regarding (A) and ($A_{corr,gfm}$) at higher Prandtl numbers reveals slight differences in accuracy; however, the error remains acceptable within the studied range of boundary layer thicknesses. Visible artefacts and stability of other coupling methods are discussed in Section 4.3.

### 4.2.2. Interfacial heat transfer measurements

The local Nusselt number at $Pr_l = 1$ is plotted in Fig. 21 according to the elevation parameter $\theta$. The original (GFM-●) applied on a coarse grid gives a poor prediction compared to the reference solution (GFM$_{ref}$-■). The error is larger as $\theta$ increases i.e. from the equator ($\theta = 0^\circ$) to the top-half region seeing the upcoming fluid ($\theta = 90^\circ$). However, the (GFM) predicts a smooth gradient profile as the variations along the azimuthal axis are low.

On the contrary, the LRS approach gives more dispersed values although the overall heat transfer coefficient is well predicted. It confirms our observations in [24]. The coupling case (A-▲) gives an overall good prediction in every region (bottom, top and equator) but it fails to keep the right prediction in regions that were already sufficiently resolved.

To overcome this issue, ($A_{corr,gfm}$-●) switches between LRS and GFM prediction to retain the best of the two. Doing so, it is visible that the prediction is enhanced and smoothed in the bottom and equator regions. It reduces the dispersion. Difficulties in reconstructing the source term in the equator region are still visible as the predicted values (●) are below the reference (■).

Finally, the case ($B_{cons}$-♦), presented as the best compromise for predicting integral quantities, exhibits a higher dispersion at the local level. It is due to the relocation of fluxes which diminishes the smoothness of the temperature field at a given instant (n) acting back to the LRS at time (n+1). Globally, the predicted values are not excessively dispersed, but the variance of the predicted value is increased.

Similar observations are made at higher Prandtl numbers. Figures are available in Appendix E.2. Increasing the Prandtl conducts to increase the dispersion. The integral of the heat transfer remains smooth, indicating that these irregularities do not significantly affect the overall heat transfer.

### 4.2.3. Temperature and gradient profiles in the interface vicinity

The radial temperature as well as its gradient have been post-processed in three distinct bubble regions. The location of the probe is given by the legends presented separately in Fig. 22. The profiles are

<figure><img src="image_23.png" /><figcaption>**Fig. 22.** Positions of the radial probes. Legends are common to Fig. 23 and Fig. 24 (as well as Fig. E.41 and Fig. E.42 in Appendix E.3).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 23.** LRS temperature profile at the bottom, equator and bubble's top regions at $Pr_l = 1$. and $Re_b \in \{67; 72\}$.</figcaption></figure>

plotted in Fig. 23 and Fig. 24, respectively.

The temperature varies rapidly in the top region (see ■ in Fig. 23). The variations are still intense in the equator region and then decrease drastically in the bottom regions. The probes in the reference simulation are shorter and have been emphasised (see ■). The (GFM) case gives a piecewise linear profile over the probe length. The coarse GFM predictions (●) are far from the reference simulation results. On the contrary, the use of coupling cases (A-▲), ($A_{corr,gfm}$-●) and ($B_{cons}$-♦) conducts accurate predictions of the temperature field in the interface vicinity. For these cases, the curves are almost superimposed. Very small discrepancies between each coupling method are visible at the equator (■). It is expected that the source term has an important contribution to that region. The coupling approach influences the temperature in the interface vicinity and the prediction of this source term and so forth. Some variations are also visible in the bottom region and at the probe's end (■).

The same observations can be drawn for the radial temperature gradient in Fig. 24. The LRS can capture the reference field appropriately on the entire range of the reference probe (GFM$_{ref}$-■). The (GFM) critically fails to predict the temperature gradient. In the Equator region, discrepancies are visible between temperature coupling strategies (A-▲), ($A_{corr,gfm}$-●) and flux coupling strategies ($B_{cons}$-♦). It is attributed to the flux relocation for two geometrical and physical reasons:

* The flux relocation acts more strongly in the bubble's diagonal direction i.e. in the least colinear direction to the Cartesian grid. It corresponds to the immediate upstream and downstream regions around the equator.

* The tangential temperature gradients are high and receive upstream fluxes. The relocation creates roughness in the temperature field, especially in the equator region which alters the approximation of $C^\theta \propto \frac{\partial T}{\partial \theta}$.

According to each figure of merit presented in previous sections (see Section 4.2.1 and Section 4.2.2), ($B_{cons}$) remains the best coupling method to achieve the right integral prediction of the heat transfer.
<figure><img src="image_25.png" /><figcaption>**Fig. 24.** LRS temperature gradient profile at the bottom, equator and bubble's top regions at $Pr_l = 1$. and $Re_b \in \{67; 72\}$.</figcaption></figure>

It has proven to lead to dispersion in the prediction of quantities but results in Fig. 23 and Fig. 24 demonstrate that locally, the temperature field is well captured in comparison to other solutions. Temperature coupling strategies (A and $A_{corr,gfm}$), although non-conservative, tend to smooth the temperature field. It is beneficial to reconstruct the source terms and feed the LRS. However, the conservation issues are too prohibitive to use these methodologies.

Supplementary material is provided in Appendix E.3, highlighting the effectiveness of the flux coupling method ($B_{cons}$). In summary, at higher Prandtl numbers, dispersion increases across all cases, with ($B_{cons}$) showing notable variability. However, these irregularities do not affect the integral heat transfer value. As the boundary layer thickness decreases, the advantages of LRS over a coarse grid GFM become more apparent, aligning with expectations outlined in [24].

### 4.3. Coupling effects on the Eulerian grid

If the predictions of the interfacial temperature gradient lead to the right order of magnitude for every coupling case, the thermal wake establishing downstream holds the cumulative footprints of the errors introduced by each method (see Fig. 25). These observations complete the analysis of the fluxes transmission through pure liquid faces.

#### 4.3.1. Qualitative observations

For the lowest and the largest Prandtl number simulations, slices of the temperature field are plotted in Fig. 25. Coupling case ($B_{cons}$) is shown at two grid refinements (16 and 32 $D_b/\Delta$) and compared to the reference solution (GFM$_{ref}$). At a Prandtl number of 1, two regions can be distinguished in the wake. The first one corresponds to the thermal iso-contours close to the interface, following the bubble shape, and stretched behind the bubble. The second region is the region detaching immediately after the bubble's equator region ($\theta = 0$). These two regions are captured from the coarser grid ($D_b/\Delta = 16$) for the two extreme values of the Prandtl numbers tested range. Qualitatively, a similar intensity can be observed in the wake between the coarse, medium and fine solutions.

Some irregularities are visible in the top and equator regions at $Pr_l = 5$. The temperature is not strictly symmetrical. This is due to the discrete interface which is not inherently symmetrical and aligned with the fixed grid. At a constant spatial resolution, the top region becomes more irregular as the boundary layer thickness decreases. These irregularities are responsible for the dispersion observed on local

quantities (see Fig. 21). However, they do not affect significantly the integral quantities as diffusion tends to smooth the solution over time.

Given the coarseness of the grids, it can be concluded that the method ($B_{cons}$) captures properly the thermal wake even though some discrepancies are visible in the top and equator regions. It is very promising because it means the bubble's wake interactions could be reproduced with fidelity using this approach ($B_{cons}$).

#### 4.3.2. Local flux distribution within the thermal wake

Once transmitted to the liquid, the energy coming from the interface can be transported downstream. Once injected within the domain, the Finite-Volume operators conserve the energy. As a consequence a flux balance must be found between the interface, the convex envelope and any downstream plane at steady-state. Local post-processing of the wake fluxes allows a complete comprehension of the transfer and how to mitigate the bias introduced in the temperature field.

In the downstream region, the directional convective and diffusive flux densities denoted $F''_{C,w}$ and $F''_{D,w}$ are assessed in the direction of gravity. For case (GFM$_{ref}$), these fluxes can be observed in Fig. 26. The convective fluxes (●) dominate the thermal wake transfers compared to the downstream diffusion (■). Immediately after the bubble, the convective fluxes decrease within the wake core (see □). By translating the plane away from the bubble, the fluxes contribution coming from the upstream transfer reaches the core and the fluxes distribution tends to a bell-shaped curve (■). It then diffuses and reduces its amplitude (see □).

This reference results for convective flux are compared with case ($B_{cons}$) in Fig. 27. If qualitatively the temperature intensity in the wake is well captured (Fig. 25), the measurements of Fig. 27 show that the amplitude of the flux profiles is less important in the coarse grid simulation along with ($B_{cons}^-$). The curve shape matches well with the reference solution (■). Profiles are slightly lower with ($B_{cons}$) as the overall Nusselt number $Nu_f$ is approximately 5% lower than for case (GFM$_{ref}$) (see Fig. 19). Besides, It is not surprising that profiles are slightly more diffused on the coarse grid as the downstream evolution is directly related to the precision of the standard numerical scheme used on the regular grid. The measurements are very encouraging. Local flux maxima (see □ and ■) are precisely located and it proves that the physics is globally captured from the coarse grid.

Supplementary integral measurements across the wake are available in Appendix D and they confirm every previous flux balance figure. To finish with, some coupling strategies introduced instabilities which are discussed qualitatively.

#### 4.3.3. Instabilities

Each coupling category is reported in Table 5 depending on the stability they exhibited in the three studied configurations.

The temperature coupling strategies are very robust and work in many cases despite being non-conservative. The sub-grid diffusion model of GFM ensures a very smooth transition between the two grids and avoids the apparition of dispersed values. (B) is stable on coarse and medium grids, and moderate and thin thermal layers. Methods using local (B, $B_{int}$) and average fluxes ($B_{avg}$, $B_{wavg}$) calculation remain stable at the lowest Prandtl number. Increasing the Prandtl for a given grid refinement creates instabilities in the equator region which renders these methods unusable (see Fig. 28). In the equator region, the temperature in the sub-cooled liquid drops below the minimum possible value; this indicates either excessive energy loss from these specific liquid cells or insufficient heat flux transferred to them to maintain an energy increment that upholds solution boundedness. ($B_{cons}$) addresses this issue with a compensation term that supplements the missing flux contribution.

This exhaustive description and assessment of each method lead to the identification of the best coupling method i.e. the conservative flux-based approach ($B_{cons}$). Its accuracy and its stability observed qualitatively, compensates for its lower performance in terms of dispersion. There is still an opportunity to reduce the irregularities introduced by the relocation algorithm by ways of smoothing and modelling the tangential source terms.
<figure><img src="image_26.png" /><figcaption>Fig. 25. Slices of the temperature field at steady state for (B<sub>cons</sub>). Two grid refinements are represented are compared to the reference solution computed with the GFM. We solve for the temperature field in the sub-cooled liquid relative to the saturation temperature.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 26. Convective and diffusive flux density at several downstream positions for the reference solution (GFM<sub>ref</sub>) - Ar<sup>s</sup> = 50, Pr<sub>l</sub> = 1.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 27. Convective flux density at several downstream positions and at steady-state. Results for the reference calculation (GFM<sub>ref</sub>) and (B<sub>cons</sub>) are compared at Ar<sup>s</sup> = 50, Pr<sub>l</sub> = 1.</figcaption></figure>

## 5. Conclusions

This paper introduces the development, implementation, and testing of a two-way coupling strategy named the Laminar Radial Sub-resolution (LRS). The LRS framework integrates a thermal sub-resolution model of the boundary layer developing around bubbles within a Cartesian CFD grid resolution. This strategy is designed to facilitate the investigation of boiling flows.

The initial case, which involved studying thermal diffusion around a spheroid with LRS, allowed us to demonstrate that the LRS method successfully captures a time-dependent solution by updating the boundary conditions of the sub-grid quasi-static problems whilst preserving the underlying physics.

In the second stage, a huge gain in the heat transfer prediction around a rising bubble has been observed on very coarse grids, answering the objective of formulating a versatile sub-grid approach. While the gain is huge at the interface, it has been difficult to transmit these fluxes to the surrounding domain in a conservative way. Establishing a flux-based coupling approach without struggling with

<figure><img src="image_29.png" /><figcaption>Fig. 28. Temperature artefacts in the equator region for coupling case (B) at Pr<sub>l</sub> = 5 and for a medium mesh (D<sub>b</sub>/Δ = 32). Mixed cells are coloured in red. The liquid is sub-cooled with ΔT = -1, and numerically, the temperature is zero in the vapour phase.</figcaption></figure>

instabilities and visible artefacts and avoiding any grid refinements has been challenging.

The temperature coupling proves to be more robust and stable to the detriment of the flux conservation. Although the conservation is not inherently guaranteed, the flux leakage remains reasonable in the
<table><caption>**Table 5**
Temperature and flux coupling stability observed at the three tested liquid Prandtl numbers. The symbol (~) means that the simulation remains stable while presenting strong unphysical artefacts. (X) means that the simulation is unstable. Clipping the temperature on the latter simulations has proved to ensure boundedness and avoid instabilities despite altering energy conservation. Assessment has been made on two grids: $D_h/\Delta = 16$ and $D_h/\Delta = 32$.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="3">Coarse</th><th colspan="3">Medium</th></tr><tr><th>1.</th><th>2.5</th><th>5.</th><th>1.</th><th>2.5</th><th>5.</th></tr></thead><tbody><tr><td rowspan="3">Temperature</td><td>GFM</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>A</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>A<sub>corr, gfm</sub></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td rowspan="5">Flux</td><td>B</td><td>✓</td><td>~</td><td>~</td><td>✓</td><td>✓</td><td>~</td></tr><tr><td>Bint</td><td>✓</td><td>~</td><td>~</td><td>✓</td><td>✓</td><td>~</td></tr><tr><td>Bcons</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Bbis</td><td>✓</td><td>~</td><td>~</td><td>✓</td><td>✓</td><td>~</td></tr><tr><td>Bbiscol</td><td>✓</td><td>~</td><td>~</td><td>✓</td><td>✓</td><td>~</td></tr></tbody></table>

range of study. Its ease of implementation makes it a good candidate to be applied to more complex configurations for rapid estimation of heat transfer coefficients.

Flux coupling strategies are not conservative as they involve only approximating the heat flux density (not weighted by surface) without performing any local flux balance over control volumes partly bounded by interface portions. Once used with LRS, these coupling approaches exhibit strong artefacts and can potentially introduce instabilities within the solver. They are very sensitive to the heat flux density computation; the diverse variations of the method considered, comfort us in the introduction of a more generic flux balance approach.

The latter has proven its efficiency in both predicting the interfacial flux and transmitting it to the domain on the whole thermal thickness range. It finally constitutes the best tested two-way coupling strategy. Particular attention to the artificial dispersion of numerical predictions for interfacial fluxes may be necessary when applying the method in the presence of phase change, to ensure that the velocity field remains stable and free from potential instabilities. Employing a cut-cell/two-scalar approach alongside a two-fluid solver for temperature may address this issue by incorporating temporal terms in mixed cells.

Overall, LRS has demonstrated that the coarse grid solution inherently contains sufficient information to perform interfacial refinements, while being three orders of magnitude less computationally expensive than simulations using the original Ghost-Fluid Method.

### CRediT authorship contribution statement

**Mathis Grosso:** Writing – review & editing, Writing – original draft, Investigation, Formal analysis, Data curation, Conceptualization.
**Guillaume Bois:** Writing – review & editing, Validation, Supervision, Project administration, Funding acquisition, Formal analysis.
**Adrien Toutant:** Writing – review & editing, Validation, Supervision, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgement

The authors would like to thank CEA/DES and TGCC for providing reliable computational resources for performing the DNS calculations. This work was granted access to the IRENE supercomputer's HPC resources under the allocations A1602B07712.

## Appendix A. Additional LRS and flux coupling functionalities

### A.1. Rising bubble frame of reference

Previous work has been performed on a stationary bubble i.e. there were no theoretical corrections to be performed on the velocity field. However, the accurate evaluation of the radial velocity component $u_r$ has proved to be essential and a correction such that $u_r = 0$ at $\dot{r} = 0$ was systematically undertaken. In the case of a bubble's displacement, the radial velocity $u_r$ is constructed by subtracting both the radial component of the rising velocity $u_{rise}$ and the radial deformation of the interface $u_r^d|_r$:

$$u_r \leftarrow u_r - u_{rise} \cdot \tilde{e}_r - u_r^d|_r = u_r - u_r|_r \quad (A.1)$$

This comes down to taking away the radial velocity component $u_r|_r$ evaluated at the interface to the rest of the probe (as shown in Eq. (A.1)). It is important to notice that the frame of reference (in which the quasi-static assumption is used) is attached to the geometrical interface i.e. it is not advected tangentially. As a consequence, the local tangential velocity must also be corrected by subtracting the tangential component of the rising bubble velocity $u_{rise}$ (see Eq. (A.2)) from the initial interpolated component $u_\theta$.

$$u_\theta \leftarrow u_\theta - u_{rise} \cdot \tilde{e}_\theta \quad (A.2)$$

### A.2. Sub-resolution deactivation

As mentioned earlier, in the bubble wake, the temporality is cut off by the quasi-static hypothesis that has been employed. Moreover, the tangential source terms are sometimes not reconstructed properly, especially in regions where they vary strongly like at the equator. Because of these two effects, it has been observed that the normal temperature gradient is sometimes poorly predicted by the LRS in comparison to the original GFM prediction. From this assessment and because the interfacial flux is often underestimated in coarse simulations, an additional option has been implemented to deactivate the transmission of information from the sub-layer to the CFD grid. In that case, the GFM takes over the LRS. The case is referred to as ($A_{corr,gfm}$).

Finally, the LRS can be deactivated in the case where the local curvature is negative or the probe length exceeds the osculating radius. In that case, a singularity could appear in the scale factor $1/r = (\dot{r} - |R_o|)^{-1}$. The quality of the curvature computation is ensured by in-house smoothing routines of Trust/TrioCFD.

### A.3. Face fluxes numerical integration

To enhance the fluxes transmitted to the carrying fluid, a discrete flux evaluation has been developed. To achieve the discrete numerical integration of the LRS quantities into face fluxes, a quadtree algorithm has been implemented in Trust/TrioCFD. Starting from the reconstructed interface $\bar{\Gamma}$, the normal distance $d_f^\perp$ to the pure fluid face is computed (see Fig. A.29.a). Then, knowing the normal face direction $n_f$, the displacement vector towards the face portions (denoted $f_p$) of surface $S_{f_p}$ are straightforwardly retrieved and projected onto the probe (see Fig. A.29.b). It results in a positive or negative normal distance increment $\Delta d_f^{(f_p)}$. The fact that the probe is longer than a mesh diagonal avoids the computation of a distance larger than the probe length i.e. $d_f^\perp + \Delta d_f^{(m)} < l_p$. Finally, the integral contributions $\phi_f^{(f_p)}$ are summed over the face portion to give $\phi_f^{(d)}$ (see Eq. (A.3)).

$$\phi_f^{(d)} = \frac{1}{S_{f_p}} \sum_{f_p} \phi_f^{(d)} (\dot{r} = d_f^\perp + \Delta d_f^{(f_p)}) S_{f_p} \quad (A.3)$$
<figure><img src="image_30.png" /><figcaption>Fig. A.29. (a) Computation of the interface to face centre normal distance $d_f^\perp$. A mean planar interface $\bar{\Gamma}$ is reconstructed from a surface average of the interface portions $\Gamma$. (b) Increment of normal distance $\Delta d_f^{(i)}$ between the cell face centre (■) and the points where the fluxes are evaluated (♦). This coupling is referred to as (B$_{int}$).</figcaption></figure>

### A.4. Averaging face fluxes contributions

Instead of computing solely the first pure fluid faces that are reachable by the probe within the mixed cell, pure face neighbours at a larger distance are considered in methods ($B_{avg}$) and ($B_{wavg}$). In detail, starting from a sub-problem, the maximum directional probe length is computed by projecting the probe length vector $l_f$ onto the Cartesian directions (Eq. (A.4)).

$$l_f = l_p \mathbf{n}_\Gamma = l_{x,f} \mathbf{e}_x + l_{y,f} \mathbf{e}_y + l_{z,f} \mathbf{e}_z \quad (A.4)$$

Then, the neighbouring face centres lying in the quadrangle of dimensions $[0, l_{x,f}] \times [0, l_{y,f}] \times [0, l_{z,f}]$ are found (see Fig. A.30). The choice of the particular quadrangle is defined by the signs of the components of $\mathbf{n}_\Gamma$. Then the flux can be evaluated using the LRS temperature or flux profile. A mean of every flux contributing to the face flux is referred to as case $B_{avg}$.

More complex weighting coefficients can be looked for to satisfy geometrical and physical considerations as in Eq. (A.5).

$$\phi_f^{\text{avg}} = \left( \sum_i w_{f_i} \right)^{-1} \sum_i w_{f_i} \phi_{f_i} \quad (A.5)$$

The colinearity between the normal probe direction and the probe centre-to-face centre direction $\mathbf{O}_i \mathbf{F}$ has been used as a weight (see Eq. (A.6)).

$$\phi_f^{\text{col}} = \left( \sum_i \frac{\mathbf{O}_i \mathbf{F} \cdot \mathbf{n}_\Gamma}{\|\mathbf{O}_i \mathbf{F} \cdot \mathbf{n}_\Gamma\|} \right)^{-1} \sum_i \phi_{f_i} \left( \frac{\mathbf{O}_i \mathbf{F} \cdot \mathbf{n}_\Gamma}{\|\mathbf{O}_i \mathbf{F} \cdot \mathbf{n}_\Gamma\|} \right) \quad (A.6)$$

Eq. (A.6) embeds the relative position of the interface portion compared to the face to correct. A similar colinearity-based weighting has already been used in interfacial mass flow rate calculation in [5].

### A.5. Pseudo-code for the coupling strategies

This section presents the flux coupling algorithms written in pseudo-code. Faces index are shown in Fig. A.31.

## Appendix B. Flux coupling variations and balance in pure diffusion

Details about sub-categories of flux coupling strategies applied in the pure diffusion case are given in this section (B$_{int}$ and $B_{wavg}$).

The use of the discrete integral to evaluate the liquid face fluxes has proved to be inefficient in the spherical diffusion case as shown in Fig. B.32. It may be relevant in convective flow at high Prandtl numbers when the boundary layer is entirely lying within a mixed cell. The selection of the interface portion to correct a given face has a greater influence on the flux prediction. Enforcing the flux conservation is beneficial for the prediction of the interfacial heat flux.

<figure><img src="image_31.png" /><figcaption>Fig. A.30. Finding neighbouring faces barycentres on which fluxes can be evaluated. Fluxes can be evaluated if $0 \le d_\perp \le l_p$. However, the region is restrained to the quadrangle. Shadowed regions are ignored.</figcaption></figure>

### Algorithm 1 Temperature coupling (A, $A_{corr,gfm}$)

**Require:** Interfacial gradient $\nabla T \cdot \mathbf{n}_\Gamma|_\Gamma$, Temperature profile $T(r)$ from LRS

**for each** Mixed cell **do**

if $|\nabla T \cdot \mathbf{n}_\Gamma|_\Gamma^{\text{GFM}}| > |\nabla T \cdot \mathbf{n}_\Gamma|_\Gamma^{\text{LRS}}|$ and ($A_{corr,gfm}$) **then**

else

Compute distance $d_\perp^\perp$ between cell-centre and interface

Evaluate $T(d_\perp^\perp)$ from LRS

Correct cell centre temperature $T^*$

end if

end for

### Algorithm 2 Flux coupling (B, $B_{int}$)

**Require:** Temperature and gradient profiles $T(r)$, $\partial_r T(r)$ from LRS

**for each** Mixed cell **do**

Locate pure liquid faces $f_i$

Reconstruct the gradient profile $\nabla T(r)$ on the probe

**for each** Face $f_i$ in $\Omega_i$ **do**

Compute face centre distance $d_\perp^\perp$ to interface

if Discrete integral (B$_{int}$) **then**:

Divide the face (see Fig. A.29.b)

**for each** Face portions $f_p$ **do**

Compute distance from face centre $\Delta d_\perp^\perp$

Evaluate $T(d_\perp^\perp + \Delta d_\perp^\perp)$ and $\nabla T(d_\perp^\perp + \Delta d_\perp^\perp)$

Evaluate fluxes $\mathcal{F}_{f_p}^C$ and $\mathcal{F}_{f_p}^D$

end for

Get fluxes $\mathcal{F}_{f_i}^C$ and $\mathcal{F}_{f_i}^D$ by averaging fluxes contributions

$\mathcal{F}_{f_p}^C$ and $\mathcal{F}_{f_p}^D$

else:

Evaluate $T(d_\perp^\perp)$, $\nabla T(d_\perp^\perp)$

Evaluate fluxes $\mathcal{F}_{f_i}^C$ and $\mathcal{F}_{f_i}^D$

end if

Apply LRS fluxes $\mathcal{F}_{f_i}^C$ and $\mathcal{F}_{f_i}^D$ at faces $f_i$

end for

end for
**Algorithm 3 Flux coupling (B<sub>avg</sub>, B<sub>wavg</sub>)**

**Require:** Temperature and gradient profiles T(r), ∂<sub>r</sub>T(r)

**for each** Mixed cell **do**

Locate pure liquid faces *f*<sub>*k*</sub>

Reconstruct gradient profile ∇T(r) on the probe

**for each** Liquid face *f*<sub>*k*</sub> **do**

Compute face centre distance *d*<sub>*f*</sub><sup>⊥</sup> to interface

Evaluate *T*(*d*<sub>*f*</sub><sup>⊥</sup>) and ∇*T*(*d*<sub>*f*</sub><sup>⊥</sup>)

Evaluate fluxes *F*<sub>*f*</sub><sup>*C*</sup> and *F*<sub>*f*</sub><sup>*D*</sup>

Compute weighting coefficient *w*<sub>*f*</sub><sub>*k*</sub> by Eq. (A.5)

Map sub-grid fluxes and weighting coefficients to the target

face *f*<sub>*k*</sub>

**end for**

**end for**

**for all** Faces *f*<sub>*l*</sub> **do**

Average fluxes contributions

Apply the sub-grid fluxes *F*<sub>*f*</sub><sup>*C*</sup> and *F*<sub>*f*</sub><sup>*D*</sup>

**end for**

**Algorithm 4 Conservative flux transmission (B<sub>cons</sub>)**

**Require:** Temperature and gradient profiles T(r), ∂<sub>r</sub>T(r) from LRS

**for each** Mixed cell **do**

Locate pure liquid faces *f*<sub>*l*</sub>

Reconstruct gradient profile ∇T(r) on the probe

**for each** Face *f*<sub>*l*</sub> in Ω<sub>*l*</sub> **do**

Compute distance *d*<sub>*f*</sub><sup>⊥</sup> between the face-centre and interface

Evaluate *T*(*d*<sub>*f*</sub><sup>⊥</sup>) and ∇*T*(*d*<sub>*f*</sub><sup>⊥</sup>)

Evaluate fluxes *F*<sub>*f*</sub><sup>*C*</sup> and *F*<sub>*f*</sub><sup>*D*</sup>

Map sub-grid fluxes and weighting coefficients to the target

face *f*<sub>*l*</sub>

**end for**

Balance ∇T · *n*<sub>*Γ*</sub> |<sub>*Γ*</sub> = ∂<sub>*r*</sub>T(0) with sub-grid fluxes

Compute compensation term *E*<sup>(d)</sup> in Eq. (14)

Compute weighting coefficient *w*<sub>*f*</sub><sub>*l*</sub>

Split *E*<sup>(d)</sup> directionally into *e*<sub>*f*</sub><sub>*l*</sub> as in Eq. (16)

**for each** Corrected face *f*<sub>*l*</sub> **do**

**if** Mixed face **then**

Relocate *e*<sub>*f*</sub><sub>*l*</sub> to neighbour cell across *f*<sub>*l*</sub>

Second splitting of *e*<sub>*f*</sub><sub>*l*</sub> into *e*<sub>*f*</sub><sub>*l,j*</sub> on faces *f*'<sub>*l*</sub>

**else**

Apply correction *e*<sub>*f*</sub><sub>*l*</sub> to *F*<sub>*f*</sub><sup>*C*</sup> and *F*<sub>*f*</sub><sup>*D*</sup>

**end if**

**end for**

**end for**

# Loop on faces where relocation has occurred #

**for each** Faces *f*'<sub>*l*</sub> **do**

Apply second correction *e*<sub>*f*</sub><sub>*l,j*</sub> to fluxes *F*<sub>*f*</sub><sup>*C*</sup> and *F*<sub>*f*</sub><sup>*D*</sup>

**end for**

The colinearity weighting (B<sub>wavg</sub>) is beneficial to reduce the error between the flux measure at the interface and the one transmitted to the domain in comparison to the standard weighting (see Fig. B.33). The colinearity weighting plays a positive role in smoothing and imposing the spherical symmetry of the problem to the Cartesian grid.

**Appendix C. Reynolds number and drag coefficients**

*Reynolds number measurements.* To ensure that the dynamics were accurately captured on both coarse and fine grids, the Reynolds number has been plotted in Fig. C.34. It can be seen that the dynamics for

<figure><img src="image_32.png" /><figcaption>Fig. A.31. Nomenclature for pseudocode.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Fig. B.32. Nusselt number comparison between the local flux interpolation (B) and discrete integral flux evaluation for case (B<sub>int</sub>). There are no significant differences.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Fig. B.33. Nusselt number comparison between the standard fluxes weighting and the one based on colinearity for case (B<sub>avg</sub>).</figcaption></figure>

both grid refinements are nearly the same. As a consequence, the LRS approach is not influenced by an eventual under-resolved velocity field.
<figure><img src="image_35.png" /><figcaption>Fig. C.34. Reynolds number measured from the coarse grid GFM ($D_h/\Delta = 16$) and reference GFM calculations ($D_h/\Delta = 64$).</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Fig. C.35. Steady drag coefficient $C_d$ measured from the coarse grid GFM ($D_h/\Delta = 16$) and reference GFM calculations ($D_h/\Delta = 64$).</figcaption></figure>

It is interesting to notice that the Reynolds number does not vary significantly according to the characteristic velocity used to establish it. The liquid to vapour velocity $\langle U_v \rangle - \langle U_l \rangle$ tends to underestimate the Reynolds number. Instead, the bubble terminal velocity $\langle U_v \rangle - U_\infty$ is preferable for the study of single rising objects as drag correlations are usually established in an infinite and quiescent domain.

**Drag coefficient.** Using the measure of the bubble terminal velocity $\langle U_v \rangle - U_\infty$ and the liquid-vapour viscosity ratio. The steady-state drag coefficient can be measured for both coarse and fine spatial resolutions (as in [24]). The transient of the drag coefficient is not relevant as it has not been computed from local interfacial forces. It is shown in figure C.35 and compared to the correlation of Feng and Michaelides [49]. Once again, the dynamic is sufficiently converged, although the simulation is very coarse.

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

<figure><img src="image_38.png" /><figcaption>Fig. D.36. Nusselt number $Nu_w$ derived from the heat transported across the wake at steady-state compared to $Nu_r$ and $Nu_{f_l}$. The main coupling cases are presented at $Ar^* = 50$. (a) $Pr_l = 1$, (b) $Pr_l = 5$.</figcaption></figure>

## Appendix D. Nusselt number across the wake

*Deriving a Nusselt number measurement from the thermal wake.* Using the measure of the convective and diffusive heat fluxes over a plane of surface $S_w$ crossing the wake ($w$) (see Eq. (D.1)), an equivalent Nusselt number denoted $Nu_w$ can be established (Green downstream plane in Fig. 12).

$$F_{C,w} = F''_{C,w} S_w = \int_{S_w} \rho_l C_{pl} T u \cdot n_w dS \qquad (D.1a)$$

$$F_{D,w} = F''_{D,w} S_w = \int_{S_w} \lambda_l \nabla T \cdot n_w dS \qquad (D.1b)$$

At steady-state, the energy increment within a sub-volume delimited by the interface and inlet-outlet planes must be zero. The numerical
<figure><img src="image_39.png" /><figcaption>Fig. E.37. Common legends to Fig. E.38 and Fig. E.39.</figcaption></figure>

scheme of the Cartesian grid being conservative, and the heat being injected only at the interface, a balance between the interfacial heat flux and the one crossing the surface of control should be found. The incoming fluid field has no energy surplus so it does not count in the energy balance. From these fluxes, an indirect measure of the interfacial Nusselt number can be achieved. This Nusselt number is denoted $Nu_w$ and is derived in Eqs. (D.2).

$$ \overbrace{\mathcal{F}''_{w}S_{w}=\mathcal{F}''_{v \to l}S_{b}}^{\approx 0} \qquad (D.2a) $$

$$ \overbrace{\mathcal{F}''_{C,w} S_w + \mathcal{F}''_{D,w} S_w}^{\approx 0} = h_{v \to l} S_b \Delta T \qquad (D.2b) $$

$$ h_{v \to l} = \frac{S_w \mathcal{F}''_{C,w} + \mathcal{F}''_{D,w}}{S_b \Delta T} \approx \frac{L_x/D_b L_y/D_b}{\pi} \mathcal{F}''_{C,w} \qquad (D.2b) $$

$$ Nu_w = \frac{h_{v \to l} D_b}{\lambda_l} \approx \frac{L_x/D_b L_y/D_b D_b}{\lambda_l \pi} \mathcal{F}''_{C,w} \qquad (D.2c) $$

Where $L_x/D_b$ and $L_y/D_b$ are the domain length to diameter ratio i.e. $(4 \times 4)$.

Strengthening our trust of conservation with the wake Nusselt numbers. The settling time of the thermal wake can be long compared to that of the overall interfacial transfer. The latter stabilises more rapidly around its final value because upstream transfers are the most intense. As a consequence, the wake transfers have been measured and integrated at $-1 \times D_b$ downstream to ensure the establishment of the wake in that region.

The flux conservation across the wake at two Prandtl numbers (1 and 5) is illustrated in Fig. D.36.a and Fig. D.36.b. Results are given on coarse grids $16D_b/\Delta$. ($B_{cons}$) proves his ability to ensure an accurate flux transmission. At the higher Prandtl number ($Pr_l = 5$), the need to refine the reference solution computed with GFM becomes apparent. A further mesh convergence is not in our interest as the aim was mainly to identify the best method for specific configurations. Mesh convergence will be mandatory in more complex configurations.

# Appendix E. Supplementary integral, interfacial and local post-processed quantities

## E.1. Integral interfacial heat transfer measurements

The legends common to all figures presented in Appendix E.1 are shown in Fig. E.37.

Interfacial Nusselt number $Nu_f$ at $Pr_l = 2.5$ and $Pr_l = 5$ are given in Fig. E.38.a and Fig. E.38.b, respectively. ($B_{avg}$) and ($B_{wavg}$) are not represented because they become unstable from the transient phase. ($B_{cons}$) outperforms the reference GFM solution at steady-state while the performance of (A) and ($A_{corr,gfm}$) decreases with the boundary layer thickness.

The error in flux transmission can be assessed for higher Prandtl numbers in Fig. E.39. The GFM exhibits excessive flux leakage ($Nu_f - Nu_{fl} < 0$), especially on coarse grids. Simple flux coupling methods (B) and ($B_{int}$) also show a high level of error. Temperature coupling methods (A) and ($A_{corr,gfm}$) present a reasonable level of error; however, this is somewhat coincidental. As the Prandtl numbers increase, there is a clear tendency for the level of flux leakage towards the liquid domain to rise in case (A). The value of $Nu_f - Nu_{fl}$ becomes negative between $Pr_l = 1$. and $Pr_l = 2.5$, indicating a transition from underestimated to overestimated fluxes.

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

<figure><img src="image_41.png" /><figcaption>Fig. E.38. Overall Nusselt number ($Nu$) over physical time $t$ at $Re_b \in \{67; 72\}$. (a) $Pr_l = 2.5$, (b) $Pr_l = 5$.</figcaption></figure>

## E.2. Local interfacial heat transfer measurements

At higher Prandtl numbers, the LRS predict more dispersed values. The boundary layer thickness is getting smaller and the dispersion caused by the relative position of the interface within the cell becomes greater (see Fig. E.40.a and Fig. E.40.b). The same observations have been made on the tangential source term. Finally, the relocation algorithm adds extra dispersion to the solution of ($B_{cons}$). Increasing the mesh resolution to 32 cells per diameter proves to be sufficient to reduce the dispersion. Any effort should focus on improving the capture of tangential quantities first and complement the relocation algorithm with conservative smoothing routines.

## E.3. Temperature and gradient profiles in the interface vicinity

At higher Prandtl numbers, LRS is still capturing the temperature and gradient field properly (see Fig. E.41, Fig. E.42 and legends in Fig. 22). Some discrepancies appear for the coupling case ($B_{cons}$). At $Pr_l = 5$, even if the temperature in the near-interface region is well predicted ($\dot{r} < 1 \times 10^{-4}$ m), some overshoots of the temperature such as $T \le T_\infty$ are observed in the top region (■). The prediction is getting coarser and more irregular in the equator region (○). The cause once again involves relocating fluxes and calculating the tangential
<figure><img src="image_42.png" /><figcaption>**Fig. E.39.** Relative error between the interfacial transfers $Nu_r$ and the liquid face fluxes $Nu_{f1}$ in the convective case at $Re_b \in \{67; 72\}$. (a) $Pr_l = 2.5$, (b) $Pr_l = 5$.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>**Fig. E.40.** Interfacial temperature gradient $\nabla T \cdot n_f|_r$ at $Re_b \in \{67; 72\}$ and steady-state. (a) $Pr_l = 2.5$, (b) $Pr_l = 5$.</figcaption></figure>
<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /><figcaption>Fig. E.41. LRS temperature (a) and gradient (b) profiles at the bottom, equator and bubble's top regions at $Pr_l = 2.5$ and $Re_b \approx 72$.</figcaption></figure>

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

<figure><img src="image_47.png" /><figcaption>Fig. E.42. LRS temperature (a) and gradient (b) profiles at the bottom, equator and bubble's top regions at $Pr_l = 5$. and $Re_b \approx 72$.</figcaption></figure>
derivative to formulate the LRS source term.

## Data availability

Data will be made available on request.

## References

[1] A. Urbano, S. Tanguy, C. Colin, Direct numerical simulation of nucleate boiling in zero gravity conditions, Int. J. Heat Mass Transfer 143 (2019) 118521, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2019.118521.

[2] L. Bureš, Y. Sato, Direct numerical simulation of phase change in the presence of non-condensable gases, Int. J. Heat Mass Transfer 151 (2020) 119400, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2020.119400.

[3] L. Bureš, Y. Sato, Direct numerical simulation of evaporation and condensation with the geometric VOF method and a sharp-interface phase-change model, Int. J. Heat Mass Transfer 173 (2021) 121233, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2021.121233.

[4] M.S. Dodd, D. Mohaddes, A. Ferrante, M. Ihme, Analysis of droplet evaporation in isotropic turbulence through droplet-resolved DNS, Int. J. Heat Mass Transfer 172 (2021) 121157, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2021.121157.

[5] L.C. Malan, A.G. Malan, S. Zaleski, P.G. Rousseau, A geometric VOF method for interface resolved phase change and conservative thermal energy advection, J. Comput. Phys. 426 (2021) 109920, http://dx.doi.org/10.1016/j.jcp.2020.109920.

[6] S. Zhao, J. Zhang, M. Ni, Boiling and evaporation model for liquid-gas flows: A sharp and conservative method based on the geometrical VOF approach, J. Comput. Phys. (2021) 110908, http://dx.doi.org/10.1016/j.jcp.2021.110908.

[7] 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, J. Comput. Phys. 491 (2023) 112366, http://dx.doi.org/10.1016/j.jcp.2023.112366.

[8] L. Torres, A. Urbano, C. Colin, S. Tanguy, On the coupling between direct numerical simulation of nucleate boiling and a micro-region model at the contact line, J. Comput. Phys. 497 (2024) 112602, http://dx.doi.org/10.1016/j.jcp.2023.112602.

[9] F. Mangani, A. Roccon, F. Zonta, A. Soldati, Heat transfer in drop-laden turbulence, J. Fluid Mech. 978 (2024) http://dx.doi.org/10.1017/jfm.2023.1002.

[10] A. Orazzo, S. Tanguy, Direct numerical simulations of droplet condensation, Int. J. Heat Mass Transfer 129 (2019) 432-448, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.07.094.

[11] N. Balcázar-Arciniega, J. Rigola, C.D. Pérez-Segarra, A. Oliva, Numerical study of the drag force, interfacial area and mass transfer in bubbles in a vertical pipe, Chem. Eng. J. (2024) 153124, http://dx.doi.org/10.1016/j.cej.2024.153124, URL https://www.sciencedirect.com/science/article/pii/S1385894724046126.

[12] E.I. Thiam, E. Masi, E. Climent, O. Simonin, S. Vincent, Particle-resolved numerical simulations of the gas-solid heat transfer in arrays of random motionless particles, Acta Mech. 230 (2) (2019) 541-567, http://dx.doi.org/10.1007/s00707-018-2346-5.

[13] N.S. Panicker, A. Passalacqua, R.O. Fox, Computational study of the effect of homogeneous and heterogeneous bubbling flows on bulk gas-liquid heat transfer, J. Fluids Eng. 142 (10) (2020) http://dx.doi.org/10.1115/1.4047806.

[14] A. Sonolet, Simulation Et Modélisation Des Transferts Thermiques Au Sein Des Écoulements Turbulents Diphasiques À Bulles (Ph.D. thesis), (2023PERP0005) Université de Perpignan, 2023, URL https://theses.hal.science/tel-04072390.

[15] E. Olsson, G. Kreiss, A conservative level set method for two phase flow, J. Comput. Phys. 210 (1) (2005) 225-246, http://dx.doi.org/10.1016/j.jcp.2005.04.007.

[16] R. Scardovelli, S. Zaleski, Interface reconstruction with least-square fit and split Eulerian-Lagrangian advection, Int. J. Numer. Methods FFluids 41 (3) (2002) 251-274, http://dx.doi.org/10.1002/fld.431.

[17] S. Popinet, Basilisk, 2014, URL http://basilisk.fr.

[18] E. Cipriano, A. Frassoldati, T. Faravelli, S. Popinet, A. Ciuci, A low-mach volume-of-fluid model for the evaporation of suspended droplets in buoyancy-driven flows, Int. J. Heat Mass Transfer 234 (2024) 126115, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2024.126115.

[19] A.W. Vreman, Particle-resolved direct numerical simulation of homogeneous isotropic turbulence modified by small fixed spheres, J. Fluid Mech. 796 (2016) 40-85, http://dx.doi.org/10.1017/jfm.2016.228.

[20] L. Bures, M. Bucci, Y. Sato, M. Bucci, A coarse grid approach for single bubble boiling simulations with the volume of fluid method, Comput. & Fluids 271 (2024) 106182, http://dx.doi.org/10.1016/j.compfluid.2024.106182.

[21] D. Bothe, M. Kröger, H.-J. Warnecke, A VOF-based conservative method for the simulation of reactive mass transfer from rising bubbles, Fluid Dyn. Mater. Process. 7 (3) (2011) 303-316, http://dx.doi.org/10.3970/fdmp.2011.007.303, URL http://www.techscience.com/fdmp/v7n3/24537.

[22] J. Maarek, S. Popinet, S. Zaleski, Subgrid scale modeling of the scalar transport equation in multiphase flows, Bull. Am. Phys. Soc. (2023).

[23] G.K. Batchelor, Small-scale variation of convected quantities like temperature in turbulent fluid Part 1. General discussion and the case of small conductivity, J. Fluid Mech. 5 (01) (1959) 113, http://dx.doi.org/10.1017/s002211205900009x.

[24] M. Grosso, G. Bois, A. Toutant, Thermal boundary layer modelling for heat flux prediction of bubbles at saturation: A priori analysis, Int. J. Heat Mass Transfer (2024) http://dx.doi.org/10.1016/j.ijheatmasstransfer.2023.124980.

[25] A. Chatelain, Simulation Des Grandes Echelles D'écoulements Turbulents Avec Transferts De Chaleur (Ph.D. thesis), Institut National Polytechnique de Grenoble - INPG, 2004, URL https://theses.hal.science/tel-00006845, Président : Pierre Sagaut ; Rapporteurs : Michel Deville et Eric Lamballais ; Examinateurs : Frédéric Ducros et Franck Nicoud.

[26] Y. Benarafa, O. Cioni, F. Ducros, P. Sagaut, RANS/LES coupling for unsteady turbulent flow simulation at high Reynolds number on coarse meshes, Comput. Methods Appl. Mech. Engrg. 195 (23-24) (2006) 2939-2960, http://dx.doi.org/10.1016/j.cma.2005.06.007.

[27] W. Bizid, Développement De Méthodes De Pénalisation Pour La Simulation De L'écoulement Turbulent Autour D'obstacles (Ph.D. thesis), Université de Bordeaux, 2017.

[28] S. Tanguy, M. Sagan, B. Lalanne, F. Couderc, C. Colin, Benchmarks and numerical methods for the simulation of boiling flows, J. Comput. Phys. 264 (2014) 1-22, http://dx.doi.org/10.1016/j.jcp.2014.01.014.

[29] R.P. Fedkiw, T. Aslam, B. Merriman, S. Osher, A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method), J. Comput. Phys. 152 (2) (1999) 457-492, http://dx.doi.org/10.1006/jcph.1999.6236.

[30] J. Palmore, O. Desjardins, A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows, J. Comput. Phys. 399 (2019) 108954, http://dx.doi.org/10.1016/j.jcp.2019.108954.

[31] D. Bothe, S. Fleckenstein, A volume-of-fluid-based method for mass transfer processes at fluid particles, Chem. Eng. Sci. 101 (2013) 283-302, http://dx.doi.org/10.1016/j.ces.2013.05.029.

[32] A. Weiner, D. Bothe, Advanced subgrid-scale modeling for convection-dominated species transport at fluid interfaces with application to mass transfer from rising bubbles, J. Comput. Phys. 347 (2017) 261-289, http://dx.doi.org/10.1016/j.jcp.2017.06.040.

[33] T.D. Aslam, A partial differential equation approach to multidimensional extrapolation, J. Comput. Phys. 193 (1) (2004) 349-355, http://dx.doi.org/10.1016/j.jcp.2003.08.001.

[34] I. Kataoka, Local instant formulation of two-phase flow, Int. J. Multiph. Flow 12 (5) (1986) 745-758, http://dx.doi.org/10.1016/0301-9322(86)90049-2, URL https://www.sciencedirect.com/science/article/pii/0301932286900492.

[35] D. Ingram, D. Causeux, C. Mingham, Developments in Cartesian cut cell methods, Math. Comput. Simulation 61 (3-6) (2003) 561-572, http://dx.doi.org/10.1016/s0378-4754(02)00107-6.

[36] Z. Xie, T. Stoesser, A three-dimensional cartesian cut-cell/volume-of-fluid method for two-phase flows with moving bodies, J. Comput. Phys. 416 (2020) 109536, http://dx.doi.org/10.1016/j.jcp.2020.109536.

[37] A. Sonolet, M. Grosso, G. Bois, A. Toutant, Conservative formulations for heat transfer simulation with sharp interfaces and discontinuous properties in two-phase flows, Comput. & Fluids 254 (2023) 105798, http://dx.doi.org/10.1016/j.compfluid.2023.105798.

[38] C. Ma, D. Bothe, Numerical modeling of thermocapillary two-phase flows with evaporation using a two-scalar approach for heat transfer, J. Comput. Phys. 233 (2013) 552-573, http://dx.doi.org/10.1016/j.jcp.2012.09.011.

[39] B. Mathieu, Études Physique, Expérimentale Et Numérique Des Mécanismes De Base Intervenant Dans Les Écoulements Diphasiques En Micro-Fluidique (Ph.D. thesis), Université de Provence, 2003.

[40] C. Pesci, A. Weiner, H. Marschall, D. Bothe, Computational analysis of single rising bubbles influenced by soluble surfactant, J. Fluid Mech. 856 (2018) 709-763, http://dx.doi.org/10.1017/jfm.2018.723.

[41] B. Abouhasanzadeh, S. Thomas, M. Taeibi-Rahni, G. Tryggvason, Multiscale computations of mass transfer from buoyant bubbles, Chem. Eng. Sci. 75 (2012) 456-467, http://dx.doi.org/10.1016/j.ces.2012.04.005.

[42] C.M.Y. Claassen, S. Islam, E.A.J.F. Peters, N.G. Deen, J.A.M. Kuipers, M.W. Baltussen, An improved subgrid scale model for front-tracking based simulations of mass transfer from bubbles, AIChE J. 66 (4) (2019) http://dx.doi.org/10.1002/aic.16889.

[43] K. Cai, G. Huang, Y. Song, J. Yin, D. Wang, A sub-grid scale model developed for the hexahedral grid to simulate the mass transfer between gas and liquid, Int. J. Heat Mass Transfer 181 (2021) 121864, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2021.121864.

[44] D. Bothe, M. Koebe, K. Wielage, J. Prüss, H.-J. Warnecke, Direct numerical simulation of mass transfer between rising gas bubbles and water, in: Bubbly Flows, Springer Berlin Heidelberg, 2004, pp. 159-174, http://dx.doi.org/10.1007/978-3-642-18540-3_13.

[45] Z. Feng, E.E. Michaelides, Drag coefficients of viscous spheres at intermediate and high Reynolds numbers, J. Fluids Eng. 123 (4) (2001) 841-849, http://dx.doi.org/10.1115/1.1412458.
[46] J. Hadamard, Mouvement permanent lent d'une sphere liquide et visqueuse dans un liquide visqueux, C. R. Hebd. Seances Acad. Sci. 152 (1911) 1735-1738.

[47] R. Satapathy, W. Smith, The motion of single immiscible drops through a liquid, J. Fluid Mech. 10 (04) (1961) 561, http://dx.doi.org/10.1017/s0022112061000366.

[48] I. Gledhill, K. Forsberg, P. Eliasson, J. Baloyi, J. Nordström, Investigation of acceleration effects on missile aerodynamics using computational fluid dynamics, Aerosp. Sci. Technol. 13 (4-5) (2009) 197-203, http://dx.doi.org/10.1016/j.ast.2009.04.008.

[49] Z. Feng, E.E. Michaelides, Heat and mass transfer coefficients of viscous spheres, Int. J. Heat Mass Transfer 44 (23) (2001) 4445-4454, http://dx.doi.org/10.1016/s0017-9310(01)00090-4.