# On the coupling between direct numerical simulation of nucleate boiling and a micro-region model at the contact line

Loric Torres, Annafederica Urbano, Catherine Colin, Sébastien Tanguy

► To cite this version:

Loric Torres, Annafederica Urbano, Catherine Colin, Sébastien Tanguy. On the coupling between direct numerical simulation of nucleate boiling and a micro-region model at the contact line. Journal of Computational Physics, 2024, 497, pp.112602. 10.1016/j.jcp.2023.112602. hal-04832248

HAL Id: hal-04832248

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

Submitted on 11 Dec 2024

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>

# On the coupling between Direct Numerical Simulation of Nucleate boiling and a Micro-Region model at the contact line

Loric Torresª, Annafederica Urbanoᵇ,¹, Catherine Colinª, Sébastien Tanguyª,∗

## ARTICLE INFO

Article history:

Keywords: Two-phase flows, Micro-region, Level-Set, Conjugate Heat Transfer

## ABSTRACT

During the past years, numerical simulations of Nucleate Boiling have attracted increasing attention. If significant advances have already been performed on the conception of accurate numerical schemes accounting for jump conditions at a non-isothermal interface, further studies are required on the modelling of the heat flux in the vicinity of the contact line. Especially, at very low scale (< 500 nm) the formation of a micro-region can be assumed, since the observed apparent contact angle can be higher than the microscopic contact angle and depend on the wall superheat. This variation of the apparent contact angle can be related to a temperature discontinuity between the saturation temperature and the local wall temperature at the contact line inducing a micro-region heat flux and an evaporation mass flow rate. Whereas this micro-region acts at very low scale, existing theoretical models suggest a strong influence of the micro-region heat flux on the global heat flux, and so on the overall bubble growth phenomenon. However, the numerical description of this micro-region in the framework of Direct Numerical Simulation (DNS) of nucleate boiling is a major issue due to the large range of spatial scales involved in the whole process. In the present paper, we propose a coupling methodology between an existing micro-region theoretical model and a larger scale DNS solver for Nucleate Boiling. Several benchmarks are proposed to show the consistency and the accuracy of the proposed method and to highlight the quantitative effect of the micro-region at the bubble scale.

© 2023 Elsevier Inc. All rights reserved.

## 1. Introduction

Nucleate Boiling is a critical phenomenon in many applications for which accurate numerical predictions are required to improve the efficiency and the reliability of the process. For example, the global heat transfer coefficient

<footnote>ªIMFT, Université de Toulouse, CNRS, INP, UPS, Toulouse, France</footnote>

<footnote>ᵇISAE SUPAERO, Université de Toulouse, Toulouse, France</footnote>

<footnote>*Corresponding author:
e-mail: tanguy@imft.fr (Sébastien Tanguy)</footnote>
of a heated wall immersed in a liquid can be significantly enhanced by increasing the number of bubbles nucleated on the wall. This can benefit to heat exchangers compactness which is crucial in many applications, as micro-electronics cooling devices. On the other hand, by increasing the wall temperature the number and the size of the nucleated bubbles also increase, which can lead to the formation of an insulated vapor layer between the liquid and the superheated wall. This phenomenon, well-known as the boiling crisis, is particularly feared in nuclear power plants, since it can lead to strong damages if the performances of the cooling system are altered. Nucleate Boiling is also of utmost importance in space applications for which the management of liquid propellant in micro-gravity is critical for the design of space launchers able to restart in flight conditions. However, in spite of these major applications, many issues must still be addressed to perform effective numerical simulations of Nucleate Boiling at the industrial scale. Indeed, the numerical simulations of nucleate boiling is made difficult by the multiscale coupling between different phenomena. The range of spatial scales involved in the whole process extends from ~ 10-100 nm, which is the typical size of the so-called micro-region, to several meters or even more for large industrial devices. Between these two extreme scales, several other typical scales can be introduced as the thickness of a micro-layer, which is approximately 1-10 µm or the bubble diameter which usually ranges from 1 mm to 10 mm in most cases. As a result, despite the continuous and huge progress in scientific computing during the last decades from the pioneering works in [1, 2] to more recent works, as in [3, 4, 5, 6, 7, 8, 9], a complete description based on first principle equations of all these phenomena remains unpractical for large scale applications. Therefore, different subgrid modelling strategies have to be developed in order to integrate the effect of the smallest scales in achievable large scale computations. The purpose of this paper is to propose a relevant coupling strategy between a Direct Numerical Simulation (DNS) solver with a micro-region model acting at a scale below 100 nm. Many experimental studies have highlighted high heat fluxes in the vicinity of the contact line. However, many issues are involved and a careful analysis is required to clearly understand the underlying physical phenomena. In particular, as pointed out in [10], two distinct phenomena associated with two different scales must be discriminated, i.e. the micro-region which is dominated by nanoscale effects as the interfacial resistance in the close vicinity of the contact line (order of 100 nm)and the micro-layer which is a much longer (about 100 micrometers) and thicker liquid film which development is driven by a coupling between hydrodynamics and thermal effects. Both phenomena have a significant effect on the global bubble growth rate and they require distinct numerical modelling strategies to be taken into account. The interested reader may find more details on experimental evidences of micro-layer in the following references [11, 12, 10, 13], and also recent numerical investigations in [14, 15, 16, 17, 18, 19, 20]. It is worth mentioning that the micro-layer does not appear in each nucleate boiling configuration, but rather develops in some specific conditions, corresponding to the so-called micro-layer regime highlighted in [10]. A correlation based on numerical simulations has been proposed in [16, 18] to determine if the nucleated bubbles will develop in the micro-layer or contact line regime.

The contact line regime can be observed essentially for lower wall superheat and higher contact angle. It is referred as the contact line regime since, during its growth, the bubble maintains a quasi-steady shape essentially driven by gravity, capillary and wall adhesion effects, as shown in [6]. Within the framework of nucleate boiling, either in micro-layer or in contact line regimes, a thermal singularity appears at the contact line since the imposed wall temperature is different from the interface temperature which is assumed to be equal to the saturation temperature. Although this singularity can be partially relaxed by local heat conduction in the wall, it is generally recognized that a thermal resistance results from this local breakage of the thermodynamic equilibrium. This means that a temperature jump between the wall and the interface must appear at the contact line, from what results a contact line heat flux. This local heat flux will have a significant impact on the interface shape in a tiny region in the vicinity of the contact line. In this so-called micro-region, the strong evaporation leads to a liquid flow toward the contact line and a significant pressure gradient bending the interface. One should notice that in the case where the thermal resistance tends to zero, the local thermodynamic equilibrium conditions specify that the interface and the wall reach the same temperature at the contact line, resulting in a local heat flux which tends to zero due to the local cooling at the contact line, and next rapidly increases, as shown in the numerical simulations presented in [21]. From this consideration, one can conclude that the contact line thermal resistance related to the local thermodynamic conditions is a critical factor to our understanding of the heat flux and phase change process in the vicinity of the contact line. Whereas an expression of this thermal resistance can be deduced to some out of equilibrium thermodynamic arguments based on kinetic equations as demonstrated in [22], an undetermined accommodation coefficient, appears in this model. This accommodation coefficient is an unknown parameter in most cases since direct measurements of the temperature jump are challenging at such low scales. Consequently, only indirect highlights of the micro-region are possible, unlike micro-layers for which infra-red temperature measurements in the solid enable the visualization of intense local heat
flux related to a long liquid film, as shown in [10, 13]. One measurable manifestation occuring from a micro-region onset is a variation of the contact angle at a macroscopic scale with the local wall superheat, as it has been shown in [23] for instance. Indeed, all existing micro-region models, as the ones developed in [24, 25, 26, 27, 28], predict a modification of the interface slope from the contact line to a given distance, of a few hundred nanometers, which is approximately the size of this micro-region. As a matter of fact, these models require to define a microscopic contact angle as an input and next provide outputs on the apparent contact angle and the heat flux distribution related to this micro-region. Unlike micro-layer (see [16]), a full resolution of this micro-region is not practicable for grid resolution considerations. Consequently, defining a suitable coupling strategy based on a subgrid resolution of this micro-region is mandatory to perform DNS of nucleate boiling in conditions where a micro-region is involved at the contact line. Several anterior studies have already proposed to couple a DNS solver with such a micro-region model, as in the following references [2, 29, 4, 30]. However, none of these works are fully focused on how to carry out such a non-trivial coupling. Therefore, some pending issues still remain and require additional investigations to understand clearly how this can be achieved in a consistent way. In particular, we will propose in this paper detailed investigations on the following points, (i) how to distribute a strong local mass flux at the contact line without altering stability computations ? (ii) how to couple the micro-region heat flux to the conjugate heat transfer in the solid wall ? (iii) how to ensure that the temperature discontinuity at the contact line, and thus the associated variation of the apparent contact angle, will be maintained when conjugate heat transfer is taken into account ? (iv) what is the measured impact of this micro-region heat flux on the global bubble growth ? As our numerical investigations have identified a strong impact of the contact line thermal resistance and the thermal conductivity ratio between the solid and the liquid, quantitave data will be also provided to quantify the respective effects of both parameters on the overall solution.

## 2. Fundamental equations and micro-region model

In this section, we remind briefly the conservation equations for the simulation of two-phase flows and their related jump conditions. Next, a subsection will be dedicated to micro-region models.

### 2.1. Conservation laws and related jump conditions

Although compressible two-phase flows solvers can be efficient at low Mach number as recently shown in [31, 32, 33, 34], our formalism is based on the incompressible Navier-Stokes equations which is an usual assumption in most Nucleate Boiling applications, since moderate temperature variations and low Mach number flows are involved in both phases. We consider a finite volume $Ω$, divided in two parts, $Ω^+$ the liquid phase and $Ω^-$ the vapor phase, and $Γ$ is the interface between both. For incompressible flows, the mass conservation equation writes,

$$ \nabla \cdot \mathbf{V} = 0, \qquad (1) $$

with $V$ the velocity field vector, which is valid both in the liquid phase and in the vapor phase. At the interface, to maintain the mass conservation in a configuration with liquid-vapor phase change, a velocity jump condition must be taken into account due to the jump of density,

$$ [\mathbf{V}]_Γ = \dot{m} \left[ \frac{1}{\rho} \right]_Γ \mathbf{n}_Γ, \qquad (2) $$

with $\dot{m}$ the local phase change mass flow rate, $\rho$ the density, $\mathbf{n}$ the normal vector pointing from the liquid phase to the vapor phase, and $[\cdot]_Γ$ the operator jump defined as $[f]_Γ = f_{liq} - f_{vap}$ for a given field $f$. This implies a discontinuous normal velocity because of the mass transfer at the interface. For the momentum conservation, we solve the incompressible Navier-Stokes equations,

$$ \rho \left( \frac{\partial \mathbf{V}}{\partial t} + (\mathbf{V} \cdot \nabla) \mathbf{V} \right) = -\nabla p + \nabla \cdot (2\mu \mathbb{D}) + \rho \mathbf{g}, \qquad (3) $$

with $p$ the pressure, $\mu$ the dynamic viscosity, $\mathbb{D}$ the strain-rate tensor and $\mathbf{g}$ the gravity vector. The following jump condition is related to the momentum conservation at the interface:

$$ [(\rho \mathbb{I} - 2\mu \mathbb{D}) \cdot \mathbf{n}_Γ]_Γ = \sigma \kappa \mathbf{n}_Γ + \nabla_Γ \sigma - \dot{m}^2 \left[ \frac{1}{\rho} \right]_Γ \mathbf{n}_Γ, \qquad (4) $$
with $\mathbb{I}$ the identity matrix, $\sigma$ the surface tension coefficient, $\nabla_{\Gamma}$ the gradient operator along the tangential direction at the interface. The term $\nabla_{\Gamma}\sigma$ is related to surface tension variations due to elastic membranes as in [35], surfactant concentrations at the interface or interface temperature variations, see [36] and [37], for more details on the latter point. Within the framework of Nucleate Boiling, a constant interface temperature is considered in most cases since both the vapor and the liquid are single species. From this, it results that $T_{\Gamma} = T_{\text{sat}}$, unlike configurations of evaporating droplet in air where tangential variations of the interface temperature occur due to the variations on the local mass fraction field, as shown in [38, 39, 37]. Finally, as we do consider neither a contaminated surface by surfactants, nor elasticity at the interface, the term $\nabla_{\Gamma}\sigma$ will be taken at zero. By projecting Eq. (4) in the normal and tangential directions, one obtains the following two relations,

$$[p]_{\Gamma} = \sigma\kappa + [2\mu\mathbf{n}_{\Gamma} \cdot \mathbb{D} \cdot \mathbf{n}_{\Gamma}]_{\Gamma} - \dot{m}^2 \left[\frac{1}{\rho}\right]_{\Gamma}, \quad (5)$$

$$[-2\mu\mathbf{t}_{\Gamma} \cdot \mathbb{D} \cdot \mathbf{n}_{\Gamma}]_{\Gamma} = 0, \quad (6)$$

with $\mathbf{t}_{\Gamma}$ the tangential vectors at the interface defined as $\mathbf{t}_{\Gamma} = \mathbb{I} - \mathbf{n}_{\Gamma} \otimes \mathbf{n}_{\Gamma}$, and $\mathbf{n}_{\Gamma} \cdot \mathbb{D} \cdot \mathbf{n}_{\Gamma}$ can be written more straightforwardly as $\frac{\partial v_{\mu}}{\partial n}$, giving,

$$[p]_{\Gamma} = \sigma\kappa + 2 \left[ \mu \frac{\partial V_{n}}{\partial n} \right]_{\Gamma} - \dot{m}^2 \left[ \frac{1}{\rho} \right]_{\Gamma}. \quad (7)$$

For an isobaric phenomenon and by neglecting the heat release due to viscous friction, the energy conservation equation can be formulated from the mass enthalpy $h$, as follows,

$$\rho \left( \frac{\partial h}{\partial t} + \mathbf{V} \cdot \nabla h \right) = \nabla \cdot (\lambda \nabla T) \quad (8)$$

where $\lambda$ is the thermal conductivity. It is worth mentioning that all the thermo-physical variables can be considered constant since nucleate boiling involves moderate temperature variations ($\Delta T < 20$ K). As a result $h$ can be defined as $h(T) - h_0 = C_p(T - T_0)$, which is valid both for a thermally perfect gas and a liquid, with $C_p$ the heat capacity at constant pressure. This gives,

$$\rho C_p \left( \frac{\partial T}{\partial t} + \mathbf{V} \cdot \nabla T \right) = \nabla \cdot (\lambda \nabla T). \quad (9)$$

The following jump condition is related to the latter equation,

$$[-\lambda \nabla T \cdot \mathbf{n}]_{\Gamma} = -\dot{m}(L_{\text{vap}} + (C_{p,\text{liq}} - C_{p,\text{vap}})(T_{\text{sat}} - T_{\Gamma})). \quad (10)$$

with $L_{\text{vap}}$ the latent heat of vaporization. It enables to maintain the conservation energy across the interface by accounting for latent heat effect. Finally, as in nucleate boiling applications the interface temperature is equal to the saturation temperature, Eq. (11) can be simplified:

$$[\lambda \nabla T \cdot \mathbf{n}]_{\Gamma} = \dot{m} L_{\text{vap}}. \quad (11)$$

## 2.2. Micro-region models

Since the pioneering work of Wayner et al. in [40], an important literature has grown around the topic of micro-region models. In [24], Stephan and Busse proposed an analytical solution for the fluid flow in the micro-region based on the lubrication theory. They also introduce the Kelvin effect, related to the variation of the surface tension involved by capillary effects at very low scale, and the interfacial resistance of Schrage [22]. This model has been developed within the framework of nucleate boiling of perfectly wetting fluids. It assumes an adsorbed film thick of a few nanometers, maintained by a disjoining pressure, in which no heat transfer can occur. In [25], Mathieu extended this model to partially wetting fluids for which a non-zero value of the microscopic contact angle is imposed. As a result there is no adsorbed film and the model describes a nanometer scale region in which the interface slope varies from a microscopic contact angle at the contact line to an apparent contact angle outside the micro-region. Schematic views of the micro-region for perfectly wetting (as described by Stephan and Busse model) and partial wetting fluids (as described by Mathieu model) are shown in Fig. 1. In Mathieu’s model for partial wetting fluids, the disjoining pressure is absent since there is no adsorbed film. Moreover, since the considered contact angles are larger than those in Stephan
<figure><img src="image_3.png" /><figcaption>Fig. 1. Schematic of a micro-region for a perfectly wetting (left) and partially wetting (right) configuration.</figcaption></figure>

& Busse's model, the lubrication theory does not hold to describe the local flow in the vicinity of the contact line. Instead, this flow is described with an analytical solution of the Stokes equations in a wedge. Janeček & Nikolayev present in [41] another model based on similar elements but which takes into account the Marangoni effect and the vapor recoil in the vicinity of the contact line. Even though Janeček & Nikolayev and Mathieu's models are similar, they have different temperature boundary conditions at the contact line. In fact, in Mathieu's model, a temperature jump at the contact line related to the local interface resistance is introduced, whereas Janeček & Nikolayev assume the temperature continuity at the contact line. Other micro-region models have been proposed in the following references [28, 42]. Although most of these models rely on a common baseline, they can be discriminated on the basis of subtle arguments of nano-scale physics. However, as experimental measurements at such small scales are hardly practicable, there is little hope to obtain in the near future conclusive measurements able to determine which one of these models is the most appropriate to describe the heat transfer and the interface shape modification due to intense liquid vapor phase change in the vicinity of the contact line. In light of this, simulations can play a significant role in the future to assess the respective qualities of these different models by performing comparisons with experimental results at large scales, for instance in terms of bubble shape, bubble diameter temporal evolution and apparent contact angle, which can provide indirect information on phenomena occurring at much lower scales. However, an important issue must be overcome before undertaking such studies: a proper implementation of the coupling between the DNS solver and such a micro-region model is required. The present paper will address this important step. The work is based on the micro-region model proposed in [25] since it has several attractive properties that will be detailed in what follows. However, one should keep in mind that the general procedure which is proposed can be easily generalized to any other micro-region model, provided it gives as an output an apparent contact angle and a local heat flux related to the micro-region.

Boiling in the micro-region induces a liquid flow so as to bring fluid in the vaporizing area. As previously stated, in the models proposed by Stephan and Busse in [24], this liquid flow is estimated from a lubrication equation which is consistent with a zero microscopic contact angle (perfectly wetting fluids), while Mathieu solves the velocity field from Stokes equations in a wedge which is suitable for a non-zero microscopic contact angle and thus partially wetting fluids. For a pinned contact line, the following velocity components $u_r$ and $u_\theta$ in the radial and azimuthal directions can be deduced from Stokes equations in a wedge with an angle $\alpha$, as shown in [43]:

$$ u_r = \frac{1}{r} \frac{Q}{\rho_{\text{liq}} L_{\text{vap}}} \frac{2 \cos(2(\alpha - \theta)) - 2 \cos(2\alpha)}{2\alpha \cos(2\alpha) - \sin(2\alpha)}, \quad (12a) $$

$$ u_{\theta} = 0, \quad (12b) $$

where $Q$ is the total heat flux evaporated in the wedge (in Watt per meter) and $\frac{Q}{L_{\text{vap}}}$ the associated vapor mass flux induced by vaporization in the wedge. It is supposed that the whole evaporation occurs in the wedge, and not along the interface. For the sake of illustration, a typical velocity field based on Eq. (12) is plotted in Fig. 2. From these equations, the pressure field can be deduced:

$$ P(r, \theta) = \frac{\mu_{\text{liq}} Q}{r^2 \rho_{\text{liq}} L_{\text{vap}}} \frac{4 \cos(2(\alpha - \theta))}{2\alpha \cos(2\alpha) - \sin(2\alpha)}. \quad (13) $$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Solution of Stokes equation in a wedge</figcaption></figure>

As the vapor pressure is considered constant, we can obtain the following relation for the pressure jump across the interface, by deriving the previous equation for $\theta = \alpha$,

$$ \frac{\partial \Delta P}{\partial s} = \frac{\mu_{liq}}{r^3} \frac{Q}{\rho_{liq} L_{vap}} \frac{8}{2\theta \cos(2\theta) - \sin(2\theta)}, \quad (14) $$

s being the curvilinear abscissa along the interface. As seen previously, the velocity expression, Eq. (12), depends on the vaporization mass flux incoming at the contact line $\frac{Q}{L_{vap}}$. It can be computed with the following equation for the local heat flux:

$$ q = \lambda_{liq} \left( \frac{T_w - T_{\Gamma}}{\delta + \lambda_{liq} R_{\Gamma}} \right), \quad (15) $$

with $R_{\Gamma}$ the interfacial resistance defined by Schrage [22], from arguments based on the out-of-equilibrium kinetic theory, as,

$$ R_{\Gamma} = \frac{2 - f T_{\text{sat}} \sqrt{2\pi R T_{\text{sat}}}}{2f \rho_{\text{vap}} L_{\text{vap}}^2}, \quad (16) $$

where $f$ is an accommodation coefficient and $R = R_g/M$ with $R_g=8.314 \text{ J.mol}^{-1}\text{K}^{-1}$ the ideal gas constant and $M$ the molar mass of the fluid. It is worth mentioning that the value of $f$ can be subject to many debates, as discussed by Batzdorf in [44]. Indeed, as its measurement is a major issue, there is no general agreement on the correct values of this accommodation coefficient. However, this parameter remains of foremost importance in calculating the micro-region superheat and flux. Its influence on the overall model is investigated further. To compute the local heat flux in Eq. (15), the local interface temperature $T_{\Gamma}$ is required. However, in the vicinity of the interface, the local saturation temperature can be modified at very low scale by strong pressure variations along the interface due to capillary forces, which is well-known as the Kelvin effect. From the Clausius-Clapeyron law, this local variation of the saturation temperature can be defined as,

$$ T_{\Gamma} = T_{\text{sat}} \left( 1 + \frac{\Delta P}{L_{\text{vap}} \rho_{\text{liq}}} \right). \quad (17) $$

Then, the local heat flux between the interface and the wall in the micro-region writes,

$$ q = \frac{T_w - T_{\text{sat}} \left( 1 + \frac{\Delta P}{L_{\text{vap}} \rho_{\text{liq}}} \right)}{\frac{\delta}{\lambda_{\text{liq}}} + \frac{2 - f T_{\text{sat}} \sqrt{2\pi R T_{\text{sat}}}}{2f \rho_{\text{vap}} L_{\text{vap}}^2}}. \quad (18) $$

From this, results the set of equations proposed by Mathieu in [25] for a partially-wetting fluids and a static contact
line,

$$ \left\{ \begin{array}{l} \frac{\partial x}{\partial s} = \cos(\theta), \hfill (19a) \\ \frac{\partial y}{\partial s} = \sin(\theta), \hfill (19b) \\ \frac{\partial \theta}{\partial s} = \frac{1}{\sigma} \Delta P, \hfill (19c) \\ \frac{\partial \Delta P}{\partial s} = \mu_l \left( \frac{Q}{\rho_{\text{liq}} L_v} \frac{1}{r^3} \times \frac{8}{2\theta \cos(2\theta) - \sin(2\theta)} \right), \hfill (19d) \\ \frac{\partial Q}{\partial s} = \frac{T_w - T_{\text{sat}}\left(1 + \frac{\Delta P}{\rho_{\text{liq}} L_{\text{vap}}}\right)}{\frac{r\theta}{\lambda_{\text{liq}}} + R_{\Gamma}}, \hfill (19e) \end{array} \right. $$

with $\frac{\partial Q}{\partial s} = q$ and $r$ the linear abscissa modified with the slip length $l_s$ defined as $r = (y + l_s)/\sin\theta$. A typical value of $l_s$ is taken to 0.1 nm and found to have little impact on the overall results while ensuring the continuity of the velocity and pressure in $x = y = 0$. It is noteworthy that the whole set of PDE is expressed as a function of the curvilinear abscissa $s$, instead of $x$. With a straightforward basis change, the $x$ and $y$ coordinates are related to $s$ with Eq. (19.a) and (19.b), which implies a fifth relation instead of four for Stephan's model. The following boundary conditions are associated to the system of equations,

The three first equations describe the contact with an origin in $(x, y) = (0, 0)$, and an initial contact angle $\theta_{\text{mic}}$. The fourth equation is for the initial integrated heat flux, which should be equal to zero, as it is defined as,

$$ Q(s) = \int_{s'=0}^{s} q \, ds'. \qquad (21) $$

Finally, an additional boundary condition has to be imposed, to ensure that at a sufficient distance of the contact line along the interface, the interface slope tends to a constant value,

$$ \left. \frac{\partial \theta}{\partial s} \right|_{s \to \infty} \to 0. \qquad (22) $$

Actually, this distance, which is not known a priori, corresponds to the spatial extension of the micro-region. From Eq. (19c), one can remark that the condition provided by Eq. (22) can be interpreted as a boundary condition for the pressure jump,

$$ \Delta P|_{s \to \infty} \to 0. \qquad (23) $$

which specifies that there is no curvature at the end of the micro-region. This assumption stems from the 2D solving of the model, which would consider a flat interface out of the micro-region. A more accurate point of view could be to impose an output boundary condition with a curvature corresponding to the bubble with a macroscopic size, and have $\Delta P|_{s \to \infty} = \frac{2\sigma}{r_{\text{eq}}}$. In practise, the influence of the output condition is low, and implementing the model for all the possible curvatures would increase the complexity significantly.

The total micro-region heat flux $Q_{\text{mic}}$, in Watt per meter of contact line, is switched into $q_{\text{mic}}$ by dividing it by the integration length of the micro-region $L_{\text{int}}$,

$$ q_{\text{mic}} = \frac{Q_{\text{mic}}}{L_{\text{int}}} \qquad (24) $$

In the rest of the paper, we will use the average value of the micro-region local heat flux $q_{\text{mic}}$ in W/m², since it is more convenient for future developments. Moreover, it can be emphasised that all the models developed in this section are purely 2D and neglect the true geometry involved in the triple line of a bubble. In all this study, the size of the micro-region is far smaller than the geometrical dimension of the bubble. The micro-region is defined on lengths
of the order of 100 nm whereas bubble are between four or five orders of magnitude larger, which justifies this local 2D approach, even for a coupling with axisymmetric or 3D computations.

### 3. Numerical Methods

We present in this section, the numerical methods developed for the coupling between a DNS Solver and the micro-region model presented previously. The DNS Solver is based on a sharp interface capturing method described by Gibou et al in [3], and extensively validated in [5, 6] against theoretical and experimental results within the framework of boiling flows and in [45, 46] for droplet vaporization and condensation. It will be briefly reminded in the next subsection. Thereafter, the solver for the micro-region model will be presented. Finally, a coupling strategy between both solvers will be detailed with specific emphasis on two important points. A subsection will focus on the coupling between the heat conduction in the solid wall and the micro-region heat flux in the vicinity of the contact line. Another subsection is devoted to the distribution of the boiling mass flow rate related to the micro-region on the computational grid. Indeed, the mass flux involved by vaporization in the micro-region takes much larger values than those computed in computational cells that do not contain the contact line. Consequently, it can appear as a singularity or a discontinuity on $\dot{m}$ in the tangential direction along the interface. This requires a specific treatment to avoid numerical instabilities due to this singular value of the micro-region vaporization mass flux.

#### 3.1. Solver for the Direct Numerical Simulation of Nucleate Boiling

The DNS solver is based on the Level Set method to capture the interface motion, the incompressible Navier-Stokes equations to describe the flow in both phases and the Ghost Fluid Method to impose jump conditions across the interface. The Level-Set-Method (LSM) proposed by Osher & Sethian in [47] consists in solving a convection equation for a $\phi$-function defined as the signed distance to the interface. The distance is set positive in the liquid phase and negative in the vapor phase. The interface motion can be computed by solving the following equation for the $\phi$-function,

$$ \frac{\partial \phi}{\partial t} + \mathbf{V}_{\Gamma} \cdot \nabla \phi = 0, \qquad (25) $$

where $\mathbf{V}_{\Gamma}$ is the interface velocity, defined as,

$$ \mathbf{V}_{\Gamma} = \mathbf{V}_{\text{liq}} + \frac{\dot{m}}{\rho_{\text{liq}}} \mathbf{n} = \mathbf{V}_{\text{vap}} + \frac{\dot{m}}{\rho_{\text{vap}}} \mathbf{n}, \qquad (26) $$

in configurations involving liquid-vapor phase change. This means that the interface velocity can be defined from two equivalent relations. In configurations involving droplet evaporation, the definition of the interface velocity based on the liquid will be preferred, whereas for configurations related to bubble growth and nucleate boiling, the definition based on the vapor velocity is more suited. A re-initialization algorithm, proposed by Sussman et al in [48], is solved after every time step to ensure that the $\phi$ function is maintained as a signed distance $d(\mathbf{r}, t)$ from the interface,

$$ \frac{\partial d}{\partial \tau} = \text{sign}(\phi)(1 - |\nabla d|). \qquad (27) $$

It is worth mentioning that this algorithm can move the zero level curve and thus induces some spurious interface displacements. However, this drawback has been corrected by applying a subcell resolution when solving Eq. (27), as proposed by Min & Gibou in [49]. The fundamental equations previously described lie on the concept of jump at the interface. This discontinuity is however discordant with classical numerical schemes which assume the continuity of the discretized variables. To preserve the continuity in the equations, the Ghost Fluid Method, proposed by Fedkiw & co-workers in [50, 51] suggests to populate ghost cells on each side of the interface. These ghost cells maintain the continuity across the interface of liquid and vapor phase variables which enables a numerical discretization on a smooth enough function. From the general guidelines of the Ghost Fluid Method, several numerical schemes have been proposed to account for pressure jump due to surface tension, density and viscosity jumps when solving Navier-Stokes equations for two-phase flows. The interested reader will find details on these points in the following references [52, 53, 54, 55, 56] for isothermal flows. Within the framework of boiling flows, a further jump condition on the normal velocity, given by Eq. (2), has to be imposed. This can be done by defining the following ghost velocities:
$$ \text{If } \phi > 0, \quad \mathbf{V}_{\text{liq}}^{\text{ghost}} = \mathbf{V}_{\text{vap}}^{n} - \dot{m} \left[ \frac{1}{\rho} \right]_{\Gamma} \mathbf{n}, \qquad (28a) $$

$$ \text{If } \phi < 0, \quad \mathbf{V}_{\text{vap}}^{\text{ghost}} = \mathbf{V}_{\text{liq}}^{n} + \dot{m} \left[ \frac{1}{\rho} \right]_{\Gamma} \mathbf{n}, \qquad (28b) $$

which have been coupled with an incompressible projection method in [51] to simulate the propagation of flame front, and used latter in [3, 5] to simulate the velocity jump condition for boiling flows. The phase change mass flow rate $\dot{m}$ is computed, from the jump condition on the energy conservation equation given by Eq. (11),

$$ \dot{m} = \frac{-\lambda_{\text{liq}} \nabla T_{\text{liq}} \cdot \mathbf{n} + \lambda_{\text{vap}} \nabla T_{\text{vap}} \cdot \mathbf{n}}{L_{\text{evap}}}, \qquad (29) $$

which requires to determine the thermal field. This last step is performed by solving the two following independent semi-implicit linear systems for the liquid and gas temperatures, respectively. Both systems are solved with the same prescribed Dirichlet boundary condition at the interface, $T_{\Gamma} = T_{\text{sat}}$, since the temperature is continuous across the interface,

$$ \text{If } \phi > 0, \quad \rho_{\text{liq}} C_{p,\text{liq}} T^{n+1} - \Delta t \nabla \cdot (\lambda_{\text{liq}} \nabla T^{n+1}) = \rho_{\text{liq}} C_{p,\text{liq}} (T^n - \Delta t \mathbf{V}_{\text{liq}} \cdot \nabla T_{\text{liq}}^n), \qquad (30a) $$

$$ \text{If } \phi < 0, \quad \rho_{\text{vap}} C_{p,\text{vap}} T^{n+1} - \Delta t \nabla \cdot (\lambda_{\text{vap}} \nabla T^{n+1}) = \rho_{\text{vap}} C_{p,\text{vap}} (T^n - \Delta t \mathbf{V}_{\text{vap}} \cdot \nabla T_{\text{vap}}^n). \qquad (30b) $$

The immersed Dirichlet boundary condition at the interface is imposed with the second order numerical scheme proposed by Gibou & Fedkiw in [57].

Similarly with the velocity ghost cells, ghost temperatures are required for gradient computations. In [58], Aslam developed a methodology to compute the field and preserve continuity of a given variable and its derivatives along the normal direction. He suggested the following set of equations which must be solved until a steady state with the fictitious time step $\tau$, to perform quadratic extrapolation,

$$ \frac{\partial T_{nn}}{\partial \tau} \pm H(\pm\phi)(\mathbf{n} \cdot \nabla T_{nn}) = 0 \qquad (31a) $$

$$ \frac{\partial T_n}{\partial \tau} \pm H(\pm\phi)(\mathbf{n} \cdot \nabla T_n - T_{nn}) = 0 \qquad (31b) $$

$$ \frac{\partial T}{\partial \tau} \pm H(\pm\phi)(\mathbf{n} \cdot \nabla T - T_n) = 0 \qquad (31c) $$

which maintains the continuity of the extrapolated field until the normal second order derivative, $H(\phi)$ being the Heaviside distribution. In Eqs. (31a), (31b), (31c), the + sign is kept for the liquid phase, and - sign for the vapor phase. Once the temperature field is properly defined, the mass flow rate is computed with Eq. (29). It is defined as the difference between the liquid heat flux crossing the interface and the vapor heat flux at the interface.

The mass flow rate is extended as well to ensure the continuity of the interface velocity in Eq. (26), before transporting the $\phi$ function with Eq. (25). The resolution of the linear systems is carried out with the Black-Box MultiGrid (BBMG) solver developed by Dendy in [59] to solve linear systems with strongly heterogeneous coefficients. It successfully has been benchmarked with several other linear system solvers in [60], and intensively used in [61] for primary atomization applications with an electro-hydrodynamic coupling. Convective derivatives in fundamental equations and in the redistance algorithm are computed with the WENO-Z scheme, proposed in [62], because of its attractive properties both for stability and accuracy. Other spatial derivatives are computed with standard second order finite volume centered schemes. The temporal integration is carried out with a second order Runge-Kutta scheme coupling the temporal evolution of all the physical variables. Time step for the overall solver accounts for the stability constraints for convective and capillary effects,

$$ \Delta t_{\text{conv}} = \frac{\Delta x}{\max \| \mathbf{V} \|} \qquad (32a) $$

$$ \Delta t_{\text{cap}} = \frac{1}{2} \sqrt{\frac{\rho_{\text{liq}} \Delta x^3}{\sigma}} \qquad (32b) $$
from which the global time step $\Delta t$ can be deduced with the following relation,

$$ \frac{1}{\Delta t} = \frac{1}{\Delta t_{\text{conv}}} + \frac{1}{\Delta t_{\text{cap}}} \qquad (33) $$

The temporal evolution of the viscosity term is carried out with an implicit temporal scheme, as proposed in [56], which removes the temporal constraint related to viscosity.

## 3.2. Numerical solver of the micro-region model

In the subsection 2.2, a micro-region model, known as the Mathieu's model, has been presented. The set of equations is composed of five ordinary differential equations. These equations are integrated with a fourth order Runge Kutta scheme with the boundary conditions given by Eq. (20). One of the difficulty to solve this system is that the boundary condition on the pressure jump $\Delta P|_{s=0}$ is not known at the contact line but at the outlet of the micro-region, instead. Indeed, we have to check that $\Delta P|_{s \to \infty} \to 0$, at a sufficient distance from the contact line, without knowing the value $\Delta P|_{s=0}$ to impose, since the latter will depend on the computed solution. As the integration process based on the Runge Kutta scheme starts at $s = 0$, a shooting method, as described in [63], is used to define an initial guess on $\Delta P|_{s=0}$. The algorithm starts with the initial guess on the pressure value, which is corrected after each iteration to reach the final condition. The computation requires multiple successive integrations on shorter ranges due to the high non-linearity of the equations. At each sub-step, the consistency of the result is verified, and the integration continued if so, until the boundary conditions $\Delta P|_{s \to \infty} \to 0$ is obtained in the far field. In particular as shown in Eq. (17), the following expression is obtained,

$$ \Delta P(s) = \frac{\rho_{\text{liq}} L_{\text{vap}} (T(s) - T_{\text{sat}})}{T_{\text{sat}}}, \qquad (34) $$

which, at the contact line, leads to,

$$ \Delta P|_{s=0} = \frac{\rho_{\text{liq}} L_{\text{vap}} (T_{\text{CL}} - T_{\text{sat}})}{T_{\text{sat}}}, \qquad (35) $$

where $T_{\text{CL}} = T_{\Gamma}(s = 0)$ is the interface temperature at the contact line, which is different from the wall temperature at the contact line, according to the following relation,

$$ q_{\text{CL}} = \frac{T_{\text{w}} - T_{\text{CL}}}{R_{\Gamma}} \qquad (36) $$

The temperature discontinuity at the contact line is consistent with the interface resistance defined in the model. We present here some preliminary results to illustrate the effect of the interface resistance $R_{\Gamma}$ on the contact line temperature discontinuity $T_{\text{w}} - T_{\text{CL}}$. For the sake of illustration, these results have been obtained from the integration of the micro-region model with the thermo-physical properties of water as this effect is more visible. The considered overall superheat $\Delta T = T_{\text{w}} - T_{\text{sat}}$ in this example is $\Delta T = 5$ K.

This temperature discontinuity between the contact line and the wall can be visualized in Fig. 3. As this local temperature discontinuity induces a singularity on the heat flux at the contact line, a thermal resistance is taken into account in the micro-region model which enables to obtain a finite heat flux despite the temperature jump at the contact line. The thermal resistance depends on an accommodation coefficient $f$ defined previously. The impact of the accommodation coefficient on the temperature difference $T_{\Gamma} - T_{\text{sat}}$ and on the contact line heat flux $q_{\text{CL}}$ is investigated. It is observed that the higher the accommodation coefficient $f$, the lower the thermal resistance and therefore the higher the contact line heat flux. The Kelvin effect related to the variations of the interface temperature with the local curvature is more important for high values of $f$, but is not sufficient to bring the interface temperature to the wall temperature. It is worth mentioning that other micro-region models as the one presented in [41] do not allow such a temperature discontinuity at the contact line. Indeed, the temperature continuity is imposed as a boundary condition at the contact line, despite their model also considers an interface resistance $R_{\Gamma}$ defined by Eq. (16). From this assumption results a zero heat flux at the contact line, unlike what is presented just prior with Mathieu's model. This issue remains an open question. Some new experimental results or molecular dynamics simulations would be required to assess definitely on the relevance of the different underlying assumptions between the existing micro-region models.

In the case of this study, the integration length of the micro-region is, for convenient reasons, chosen equal to the mesh size of the DNS solver, in order to be injected as a boundary condition in the DNS solver in the mesh cell
<figure><img src="image_5.png" /><figcaption>Fig. 3. Temperature and heat flux evolutions for three different values for the accommodation coefficient. Superheat $\Delta T = 5$ K, the liquid is water and the integration length $L_{int} = 1$ µm.</figcaption></figure>

containing the contact line. With sizes around micrometers, the computation for a given superheat $\Delta T$ can last about minutes. In practice, the longer the integration length, the higher the computational cost. The micro-region model is solved for a large range of the superheat $\Delta T$ on a specific length before performing the DNS, and the results are saved as a table integrated in the DNS solver to avoid prohibitive computations at each temporal iteration of the DNS solver.

Further investigations are performed with FC-72 at 537 mbar, the coolant used in the RUBI experiment that will be described in Section 4.1. For coolant, the accommodation coefficient is usual high. Therefore, we set $f$ to 1 for all further results involving FC-72. In Fig. 4 has been plotted the average heat flux and the macroscopic contact

<figure><img src="image_6.png" /><figcaption>Fig. 4. Comparison between Stephan's and Mathieu's model computed with FC-72. Integration length is 2 µm.</figcaption></figure>

angle given by the micro-region model computed with Mathieu's and Stephan's model. They provide almost identical results in terms of heat flux, but the macroscopic contact angle differs. However, in qualitative terms, similar trends are observed with strong contact angle variations with the wall superheat for low superheat values which tends to stagnate for higher superheat values. It is worth mentioning that Stephan's model assumes a zero microscopic contact angle, whereas Mathieu's model requires a non-zero microscopic contact angle value. However, for low microscopic contact angle ($\theta_{mic} = 8^\circ$ in the comparison presented), the solution obtained with Mathieu's model depends weakly
on the imposed microscopic contact angle.

### 3.3. Numerical resolution of the heat conduction in the superheated wall

Numerical simulations of nucleate boiling need to account for conjugate heat transfer between the fluid domain and the solid domain. In this subsection, we describe briefly the heat conduction coupling between the two domains. However, our analysis is restricted to the grid cells that do not contain the contact line. Indeed, the case of the grid cells containing the contact line will be thoroughly investigated in the next subsection. The unsteady heat equation is solved in the solid,

$$ \rho_w C_{p,w} \frac{\partial T}{\partial t} = \nabla \cdot (\lambda_w \nabla T) \quad (37) $$

where the subscript w refers to the solid wall physical properties. Spatial derivatives are solved with a standard second order numerical scheme, and the same temporal scheme as the fluid solver for the temporal derivatives, i.e. a second order Runge-Kutta scheme. It is coupled with the energy conservation equation in the fluid domain at the frontier between both domains through the continuity condition on the local heat flux which specifies that,

$$ -\lambda_w \nabla T|_w \cdot \mathbf{n} = -\lambda_f \nabla T|_f \cdot \mathbf{n} \quad (38) $$

with $\lambda_f$ the thermal conductivity associated to the fluid, which can be a liquid or a vapor within the framework of nucleate boiling simulations. For a plane solid frontier, this condition simply writes,

$$ -\lambda_w \frac{\partial T}{\partial y}\bigg|_w = -\lambda_f \frac{\partial T}{\partial y}\bigg|_f \quad (39) $$

From this condition on the local heat flux continuity and imposing the continuity of the temperature between the two domains, the following boundary conditions on the temperature $T_{f,i,0}$ at the bottom of the fluid domain, and $T_{w,i,N_y+1}$ at the top of the solid domain can be deduced,

$$ T_{f,i,0}^{n+1} = \frac{2\lambda_w}{\lambda_w + \lambda_f} (T_{w,i,N_y}^n - T_{f,i,1}^n) + T_{f,i,1}^n \quad (40a) $$

$$ T_{w,i,N_y}^{n+1} = \frac{2\lambda_f}{\lambda_w + \lambda_f} (T_{f,i,1}^n - T_{w,i,N_y}^n) + T_{w,i,N_y}^n \quad (40b) $$

where the index n indicates the temporal time step. As nucleate boiling simulations involve small time steps due to the capillary time step constraint, this explicit coupling between the boundary conditions of both domains at each temporal iteration gives a satisfactory accuracy in this framework. The linear system resulting of the implicit temporal discretization of Eq. (37) is solved with a classical Gauss-Seidel method which is efficient since the corresponding matrix is diagonally dominant.

### 3.4. Coupling the micro-region heat flux with the wall heat conduction

In the Section 3.3, the coupling between the heat conduction in the solid domain and the energy conservation in the fluid domain was presented for the grid cells that do not contain the contact line. Accounting for the micro-region heat flux in the cells containing the contact line is carefully addressed in the present subsection. In these cells, the following equality holds:

$$ q_{\text{mic}} = -\lambda_w \frac{\partial T}{\partial y}\bigg|_{\text{CL}} \quad (41) $$

where the micro-region heat flux $q_{\text{mic}}$ is computed from the microregion model, and thus depends on $T_{\text{sat}}$ and $T_w$, with $T_w$ the wall temperature in the mesh crossed by the contact line. However, as $T_w$ is not known a priori, since it is part of the solution, $q_{\text{mic}}$, the average micro-region local heat flux, is also unknown, and Eq. (41) can not be imposed as a simple Neumann boundary condition. In Fig. 5, the evolution of $q_{\text{mic}}$ with the local superheat $\Delta T = T_w - T_{\text{sat}}$ is plotted for FC-72, with an integration length of 2 µm with Mathieu's model and the fluid properties reported in Table 1. This figure highlights a quasi-linear relation between the local micro-region heat flux and the local superheat $\Delta T$ such as,

$$ q_{\text{mic}} = AT_w + B. \quad (42) $$
<table><caption>Table 1. Physical properties of N-perfluorohexane (FC-72) at 537 mbar, T<sub>sat</sub> = 39.5°C.</caption><thead><tr><th>ρ<sub>liq</sub> (kg.m<sup>-3</sup>)</th><th>ρ<sub>vap</sub> (kg.m<sup>-3</sup>)</th><th>C<sub>p,liq</sub> (J.kg<sup>-1</sup>.K<sup>-1</sup>)</th><th>σ (N.m<sup>-1</sup>)</th><th>L<sub>vap</sub> (J.kg<sup>-1</sup>)</th><th>λ<sub>liq</sub> (W.m<sup>-1</sup>.K<sup>-1</sup>)</th><th>μ<sub>liq</sub> (kg.m<sup>-1</sup>.s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1650.5</td><td>7.49</td><td>1072.0</td><td>9.75 × 10<sup>-3</sup></td><td>89.28 × 10<sup>3</sup></td><td>5.41 × 10<sup>-2</sup></td><td>5.38 × 10<sup>-4</sup></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Linear regression of the microregion heat flux according to Eq. (42). FC-72, 537 mbar. Integration length is 2 µm.</figcaption></figure>

Determining the coefficient A and B in Eq. (42) requires to perform pre-computations of the micro-region model to generate a data-base from which these coefficients can be extracted. For example, in Table 2, the coefficients obtained with Mathieu's micro-region model are given for computations carried out with water at ambient pressure, with a microscopic contact angle θ = 10°, over a superheat range of 0 to 3 K and. As q<sub>mic</sub> depends on the integration length L<sub>int</sub> (as it is averaged on this length), we report results in Table 2 for different values of L<sub>int</sub> which correspond to typical sizes of computational cells for Direct Numerical Simulations of nucleate boiling. This highlights that the coefficients A and B are mesh-dependent since the local heat flux is not homogeneous along the interface. Therefore, A and B must be adapted for each mesh size and fluid properties. Finally from Eqs. (42) and (41), the following Robin boundary condition is obtained on the temperature in cells containing the contact line:

$$AT_w + \lambda_w \left. \frac{\partial T}{\partial y} \right|_{\text{CL}} = -B. \qquad (43)$$

In Fig. 6, a schematic of the cell configuration at the frontier between the solid and the fluid domains is represented to show that this frontier is imposed at the point (i, N<sup>w</sup><sub>y</sub> + 1/2). If a wall cell is crossed by the interface that is, φ<sub>i,N<sup>w</sup><sub>y</sub>+1/2</sub> × φ<sub>i+1,N<sup>w</sup><sub>y</sub>+1/2</sub> < 0, then it contains the contact line and the wall temperature is refered as T<sub>w</sub>.

Then the following discrete approximation can be obtained for the conduction heat flux in the cell containing the contact line:

$$-\lambda_w \left. \frac{\partial T}{\partial y} \right|_w = -\lambda_w \frac{T_{w,i,N_y^w} - T_w}{\frac{\Delta y}{2}} \qquad (44)$$

<table><caption>Table 2. Values of coefficients A and B of Eq. (42) for different meshes. Calculated for FC-72, 537 mbar, with a microscopic contact angle of 10°, with Mathieu's model for superheat range of 0 to 3 K.</caption><thead><tr><th colspan="2">L<sub>int</sub> = 2 µm</th><th colspan="2">L<sub>int</sub> = 4 µm</th><th colspan="2">L<sub>int</sub> = 8 µm</th></tr><tr><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>B</th></tr></thead><tbody><tr><td>2.0001×10<sup>5</sup></td><td>-6.2535×10<sup>7</sup></td><td>1.1829×10<sup>5</sup></td><td>-3.6986×10<sup>7</sup></td><td>6.8893×10<sup>4</sup></td><td>-2.1540×10<sup>7</sup></td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>Fig. 6. Temperature boundary conditions for the solid domain. The interface is illustrated by the red dotted line.</figcaption></figure>

<table><caption>Table 3. Values of coefficients C and D of equation (50) for different configurations of mesh and solid thermal conductivity. Calculated with coefficient of table 2.</caption><thead><tr><th rowspan="2">λ<sub>w</sub></th><th colspan="2">L<sub>int</sub> = 2 µm</th><th colspan="2">L<sub>int</sub> = 4 µm</th><th colspan="2">L<sub>int</sub> = 8 µm</th></tr><tr><th>C</th><th>D</th><th>C</th><th>D</th><th>C</th><th>D</th></tr></thead><tbody><tr><td>1 W.m<sup>-1</sup>.K<sup>-1</sup></td><td>0.667</td><td>104.2</td><td>0.617</td><td>119.6</td><td>0.568</td><td>135.1</td></tr><tr><td>5 W.m<sup>-1</sup>.K<sup>-1</sup></td><td>0.923</td><td>24.05</td><td>0.910</td><td>28.25</td><td>0.896</td><td>32.66</td></tr><tr><td>50 W.m<sup>-1</sup>.K<sup>-1</sup></td><td>0.992</td><td>2.491</td><td>0.991</td><td>2.945</td><td>0.989</td><td>3.427</td></tr></tbody></table>

From this equality it follows that:

$$q_{\text{mic}} = AT_w + B = \lambda_w \frac{T_{w,i,N_y^w} - T_w}{\frac{\Delta y}{2}}. \quad (45)$$

Consequently, the following expression for the solid temperature at the contact line $T_w$ is obtained:

$$T_w = \frac{\frac{2\lambda_w}{\Delta y} T_{w,i,N_y^w} - B}{\frac{2\lambda_w}{\Delta y} + A} \quad (46)$$

The present solid thermal solver does not use $T_w$ for the boundary condition but a solid ghost cell value in the liquid domain instead, twice as far from the last solid mesh than the wall. By balancing the flux, the following equalities hold:

$$q_{\text{mic}} = \lambda_w \frac{T_{w,i,N_y^w} - T_w}{\frac{\Delta y}{2}} = \lambda_w \frac{T_w - T_{w,i,N_y^w+1}}{\frac{\Delta y}{2}} = \lambda_w \frac{T_{w,i,N_y^w} - T_{w,i,N_y^w+1}}{\Delta y}, \quad (47)$$

from which the solid ghost cell temperature $T_{w,i,N_y^w+1}$ can be deduced:

$$T_{w,i,N_y^w+1} = 2T_w - T_{w,i,N_y^w} = T_{w,i,N_y^w} \left( \frac{\frac{2\lambda_w}{\Delta y} - A}{\frac{2\lambda_w}{\Delta y} + A} \right) - \frac{2B}{\frac{2\lambda_w}{\Delta y} + A} \quad (48)$$

Eq. (48) can be written as:

$$T_{w,i,N_y^w+1} = CT_{w,i,N_y^w} + D \quad (49)$$

with

$$C = \left( \frac{\frac{2\lambda_w}{\Delta y} - A}{\frac{2\lambda_w}{\Delta y} + A} \right) \quad \text{and} \quad D = - \frac{2B}{\frac{2\lambda_w}{\Delta y} + A} \quad (50)$$

The coupling between the micro-region model and the solid heat conduction solver is taken into account through the coefficients C and D which both depend on the solid thermal conductivity $\lambda_w$ and on the coefficients A and B. In Table 3, some values of C and D are provided for different values of $L_{\text{int}}$ and $\lambda_w$. This demonstrates the strong influence of the solid thermal conductivity on the coefficients, and so on the ghost temperature $T_{w,i,N_y^w+1}$. In particular, it is noticeable that for $\lambda_w \to \infty$, C $\to 1$ and D $\to 0$ which means that $T_{w,i,N_y^w+1} \to T_{w,i,N_y^w}$. This corresponds to an
almost isothermal wall which is consistent with the asymptotic condition $\lambda_w \to \infty$. On the other hand, for a given local heat flux a lower value of the solid thermal conductivity will involve larger variations in the solid temperature field. From Eq. (50), the following expression of the temperature derivative in the solid domain at the contact line can be deduced:

$$ \left. \frac{\partial T_w}{\partial y} \right|_{i,N_y^w+1/2} = \frac{(C-1)T_{w,i,N_y^w} + D}{\Delta y}. \qquad (51) $$

As this expression of the temperature derivative in the cell containing the contact line is linear, it can be directly injected in the linear system resulting from the spatial and implicit temporal discretization of Eq. (37).

This modified boundary condition at the contact line enables an effective coupling between the micro-region model and the solid heat conduction solver. Afterward, to account for the model in the fluid domain, the heat flux computed previously with the temperature field in the wall below the contact line is converted into an evaporation mass flow rate injected in the liquid domain. To assure this contribution is not counted twice, the conduction heat flux from the wall to the fluid domain is set to zero in the cell containing the contact line, since the entirety of this heat flux related to the micro-region is absorbed by the liquid-vapor phase change. To do so, the ghost cell of the liquid mesh below the contact line is set to $T_{sat}$, the temperature of the interface. Investigations showed a very weak influence of this boundary condition on the vaporization flux at the interface.

### 3.5. Coupling between the micro-region model and the fluid solver

The previous section was devoted to the coupling between the micro-region model and the solid heat conduction solver. From this first part of the coupling, the local superheat at the contact line can be deduced, and so the micro-region heat flux and the apparent contact angle. These two outputs of the micro-region model must now be injected as inputs in the two-phase flows solver. Imposing the apparent contact angle is straightforward since the micro-region model provides a value which has just to be imposed with the following boundary condition on the $\phi$ function when solving the interface propagation equation, Eq. (25),

$$ \left. \frac{\partial \phi}{\partial y} \right|_{w} = -\cos \theta \qquad (52) $$

Compared to a nucleate boiling computation which does not account for a micro-region, this part of the coupling does not involve additional difficulties, except that the apparent contact angle can evolve with time. Temporal variations of the local superheat, and so of the apparent contact angle, can be computed due to the coupling with the solid heat conduction solver. A more complex issue is the coupling between the phase change mass flow rate of the micro-region with the fluid solver. Indeed, as the micro-region acts at scales much lower than the typical computational cell of a nucleate boiling simulation, a proper numerical strategy has to be defined to carry out this coupling.

#### 3.5.1. Micro-region mass flow rate in a single cell at the contact line

The most straightforward way to account for the contribution of the micro-region mass flux is to integrate the boiling mass flow rate only in cells containing the contact line, since the micro-region mass flux is integrated on a length equal to the mesh size. It simply consists in substituting the computed mass flow rate with the heat flux originating from the micro-region model converted in a vaporization mass flux. It is noticeable that the micro-region heat flux can be computed on a longer length than the micro-region length by using a simple heat conduction law outside of the micro-region. Therefore, an overall heat flux can be integrated on an arbitrary length which can be equal to the size of a computational mesh. The following expression of the boiling mass flow rate has to be imposed, in the cell containing the contact line,

$$ \dot{m}_{mic} = \frac{q_{mic}}{L_{vap}} \qquad (53) $$

Therefore, a mass flow rate spot is generated at the contact line. From this discontinuity, difficulties emerge to handle the contact line motion precisely. Especially, for high solid thermal conductivity with respect to the fluid conductivity, the magnitude of the mass flow rate related to the micro-region significantly exceeds the values outside of the micro-region, which leads to numerical instabilities at the contact line. This can be observed in Fig. 7 where the flow field in the vicinity of the contact line is displayed for two different values of the solid thermal conductivity. For low values,
<figure><img src="image_9.png" /><figcaption>Fig. 7. Mass flow rate field and its extension in the vicinity of the contact line for two different solid wall conductivity. Velocity scales are different in the two figures.</figcaption></figure>

the micro-region contribution is moderate compared to the overall vaporization rate, the spot at the contact line is of the same order of magnitude as the vaporization rate computed by the DNS. The velocity field at the contact line does not seem to be disturbed, as the flow resulting from the local vaporization at the contact line cell is consistent with a Stokes flow. At the opposite, for high values of the solid thermal conductivity, the method shows an erratic behavior with a contact line oscillating between adjacent cells and a contact angle not in accordance with the imposed one. The velocity field in this region is extremely unstable from an iteration to another. In Fig. 7(b), one can observe a thin liquid layer close to the micro-region spot, due to numerical flaws. This behavior is regularly observed with this method coupled with high solid/liquid thermal conductivity ratios. Thus, the method is not considered as valid as the singularity of the source term causes obvious inconsistencies and numerical instabilities. To alleviate these stability issues, other methods are investigated in what follows. This quite basic coupling method will be referred in the rest of the paper as DNS-LMR for DNS with Local Micro-Region model.

### 3.5.2. Micro-region mass flow rate integrated in the volumetric divergence term

Another possible method is to integrate and average the contribution of the micro-region mass flux term in the whole bubble, as proposed in [2]. To account for this, a source term is integrated in the divergence term of the mass balance equation. Indeed, the mass balance in the bubble allows us to write,

$$ \dot{M}_b = \frac{d}{dt} \iiint_{V_b} \rho \, dV = \iiint_{V_b} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{V}) \, d^3\mathbf{V} \quad \underbrace{=} \quad \rho_{\text{vap}} \iiint_{V_b} \nabla \cdot \mathbf{V} \, d^3\mathbf{V} \qquad (54) $$

incompressible

with $V_b$ the volume of the bubble. In order to account for the total mass flux related to the micro-region, $\dot{M}_{mic}$, defined as,

$$ \dot{M}_{mic} = \frac{2\pi r_{CL} \times Q_{mic}}{L_{vap}}, \qquad (55) $$

with $r_{CL}$ the contact line radius. We add an average source term in the velocity divergence term, which accounts for the mass gain due to the micro-region:

$$ \nabla \cdot \mathbf{V}_{|vap} = \frac{\dot{M}_{mic}}{\rho_{vap} V_b}. \qquad (56) $$

This contribution is added as a source term in the Poisson equation for the pressure field when solving the Navier-Stokes equations with a projection method. Despite its ease of implementation, this method suffers from an important shortage, since the local description of the vaporization process in the vicinity of the contact line is lost. Indeed, only
the average effect of this term is taken into account at the bubble scale. As a result, it could remove some physical flow structures in the vicinity of the contact line, and at the opposite, induce some artificial flow structures at the bubble scale. However, due to the integration of the micro-region mass flow rate at the bubble scale, some improvements are expected regarding the computation stability, with respect to the method described in the previous subsection.

This coupling method will be referred in the rest of the paper as DNS-IMR for DNS with Integrated Micro-Region model.

### 3.5.3. Spreading the micro-region mass flow rate along the interface

As both methods described in the previous subsections exhibit some deficiencies, i.e., stability issues for the first one, and the loss of local effect of the micro-region mass flow rate for the second one, we propose, in this subsection, an alternative method to overcome these shortcomings. As it has been identified that the strong spurious currents at the contact line originate from the discontinuity of the local mass flow rate along the interface, we propose to spread out the micro-region mass flow rate from the contact line to a given distance along the interface.

This can be achieved by carrying out the following steps:

* At first, a curvilinear abscissa $s$ is defined as the distance from the contact line along the liquid/vapor interface. This field is generated, starting at the contact line and his values is 0 in every cell not crossed by the interface. As, the microregion contribution is confined in a limited area, the curvilinear abscissa is not computed above a specific height. A typical field of $s$ can be visualized in Fig. 8(a) for a given interface shape.

* Afterward, an arbitrary function of the curvilinear abscissa $u(s)$ is defined in the cells crossed by the interface, as plotted in Fig. 8(b).

* The function is extended across the interface to ensure the continuity along the normal direction, as it can be observed in Fig. 8(c). This step is required for stability reasons, and is performed following the Aslam extensions [58].

<figure><img src="image_10.png" /><figcaption>Fig. 8. Example of curvilinear abscissa and spreading function, dimensionless here.</figcaption></figure>

We chose a Gaussian function to smooth $\dot{m}_{mic}(s)$ according to:

$$ \dot{m}_{mic}(s) = A_0 e^{-\frac{s^2}{a^2}}, \qquad (57) $$

since it provides a regular function which is maximum at the contact line, i.e. if $s = 0$, and tends quickly to zero when $s$ is increasing. $A_0$ is the maximum amplitude at the contact line and $a$ is a spreading coefficient which regulates the width of the distribution. In order to conserve the correct value of the total micro-region mass flux, the following condition on the integral of $\dot{m}_{mic}(s)$ must be fulfilled,

$$ \dot{M}_{mic} = \int_{0}^{\infty} \dot{m}_{mic}(s) ds = \frac{2\pi r_{CL} Q_{mic}}{L_{vap}}, \qquad (58) $$
which gives the following expression for the amplitude of the smoothing function,

$$A_0 = \frac{2\dot{M}_{\text{mic}}}{a\sqrt{\pi}}. \qquad (59)$$

Once the smoothing function has been defined the field $\dot{m}_{\text{mic}}(s)$ has to be defined in the computational domain. It consists in populating with appropriate values of $\dot{m}_{\text{mic}}(s)$ the computational nodes which are crossed by the interface. After this first step has been achieved, $\dot{m}_{\text{mic}}(s)$ is extended along the normal direction at the interface to build a continuous field in the whole computational domain. This extrapolation along the normal direction is carried out by solving the following PDE:

$$\frac{\partial \dot{m}_{\text{mic}}}{\partial \tau} + H(\phi)\mathbf{n} \cdot \nabla \dot{m}_{\text{mic}} = 0. \qquad (60)$$

<figure><img src="image_11.png" /><figcaption>Fig. 9. Implementation of the solid/fluid microregion coupling. The flux is here spread near the contact line. As the conduction heat flux in the cell below the contact line is entirely accounted for in the vaporization mass flow rate, no additional heat flux has to be added to the fluid in this cell.</figcaption></figure>

This extension is required on each side of the interface to build a smooth enough $\dot{m}$-field to define both a continuous interface velocity $V_{\Gamma}$ for the advection of the level set function in Eq. (25) and continuous extrapolation fields for the variables $V_{\text{vap}}^{\text{ghost}}$ and $V_{\text{liq}}^{\text{ghost}}$ in Eq. (28). Once the $\dot{m}$-field has been built in the whole domain, it is added to the classical boiling mass flow rate calculated by the DNS solver. Both fields can be observed in Fig.10.

<figure><img src="image_12.png" /><figcaption>(a) Smoothed micro-region mass flow rate field at the interface and its extension in the whole domain</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Total mass flow rate field at the interface and its extension in the whole domain</figcaption></figure>

<caption>Fig. 10. Boiling mass flow rate in the vicinity of the interface and its extension. In the present case, $a = 5r_{\text{CL}} \times dx$.</caption>

A five cells spreading length has been adopted as a good compromise between stability issues while keeping with the flow locality in the vicinity of the contact line. Indeed, stretching the micro-region on a too small number of cells
<figure><img src="image_14.png" /><figcaption>Fig. 11. RUBI numerical setup</figcaption></figure>

is similar to the first method described in part 3.5.1. On the other hand spreading lengths of 20 or 40 cell sizes seem unreasonable to correctly depict the local effects of the micro-region on the flow in the vicinity of the contact line. This coupling method will be referred in the rest of the paper as DNS-SMR for DNS with Smoothed Micro-Region model, and simulations which do not account for the micro-region will be referred as DNS-NMR for DNS with No Micro-Region model.

# 4. Numerical results

We present in this section the numerical simulations obtained from the coupling between the micro-region model at the contact line and the solver for the Direct Numerical Simulation of nucleate boiling. In the first test-case, we present comparisons between numerical simulations and experimental results in micro-gravity conditions. In the second benchmark, investigations on the impact of the solid wall conductivity on the micro-region heat flux are presented for bubbles in micro-gravity conditions. In the third example, the effect of the accommodation coefficient on the micro-region heat flux and so on the bubble growth is investigated both in micro-gravity and earth conditions. As all the previously mentioned benchmarks consider nucleate boiling in the contact line regime, we present in the last subsection some preliminary results of nucleate boiling in the micro-layer regime while accounting for the micro-region effect at the contact line.

## 4.1. Validation of DNS solver coupled with a micro-region model against RUBI experiment

The Multiscale Boiling experiment, called RUBI (Reference mUltiscale Boiling Investigation) described in [64], is a nucleate boiling experiment which has been performed in the international Space Station in order to benefit of long-time micro-gravity conditions, unlike drop tower or parabolic flight experiments. It provides an important number of configurations and includes large parameter scales. We present in this subsection comparisons between our numerical solver and an experiment of pool boiling with a low subcooling.

A schematic of the configuration is plotted in Fig. 11. The size of the computational domain is 2×4 mm², with a wall 5 mm thick. An adiabatic condition is imposed at the bottom of the wall. The fluid/wall interface is covered by a chromium thin coating (400 nm) where a volumetric heat flux $q'''$ is released due to the Joule effect. The initial conditions need to be carefully set up as well. Heating power is applied at the beginning of the experiment and after a specific time $t_{wait}$, a laser beam triggers the bubble nucleation on an artificial cavity. During this time, a thermal boundary layer is developed due to the liquid heating by the wall. This thermal initial condition is extracted
from a COMSOL simulation which computes the thermal field at the scale of the whole experiment. In the specific configuration considered here, the temperature has been increased so that the whole domain is superheated by 0.5 K. The volumetric heat flux $q'''$ is imposed inside the first computational cell of the solid domain so that it corresponds to a surface heat flux of 10000 W/m$^2$, and is applied 5 s before the bubble nucleation. The heat flux is maintained all along the experiment. The fluid is perfluorohexane at 500 mbar which properties are reported in Table 4.

<table><caption>Table 4. Physical properties used in the simulation. Fluid is FC-72, P=500 mbar, T<sub>sat</sub> = 310.85 K. Solid is sapphire.</caption><thead><tr><th></th><th>ρ (kg.m<sup>-3</sup>)</th><th>C<sub>p</sub> (J.kg<sup>-1</sup>.K<sup>-1</sup>)</th><th>σ (N.m<sup>-1</sup>)</th><th>L<sub>vap</sub> (J.kg<sup>-1</sup>)</th><th>λ (W.m<sup>-1</sup>.K<sup>-1</sup>)</th><th>μ (kg.m<sup>-1</sup>.s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Liquid</td><td>1653.9</td><td>1069.2</td><td>9.91×10<sup>-3</sup></td><td>89.77×10<sup>3</sup></td><td>5.44×10<sup>-2</sup></td><td>5.49×10<sup>-4</sup></td></tr><tr><td>Vapor</td><td>7.02</td><td>2034</td><td>/</td><td>/</td><td>0.024</td><td>1.228×10<sup>-5</sup></td></tr><tr><td>Solid</td><td>4890</td><td>410</td><td>/</td><td>/</td><td>10.9</td><td>/</td></tr></tbody></table>

The different methods introduced previously are now compared to the experimental results obtained with the RUBI experiment. For the sake of illustration, we show in Fig. 12 some snapshots of the bubble shape and thermal field at different times corresponding respectively to the beginning of the simulation, an intermediate time and the final time of the simulation. These snapshots have been obtained with the DNS-SMR solver, and they enable to visualize the bubble growth with time and the temperature field in the fluid and solid domains.

<figure><img src="image_15.png" /><figcaption>Fig. 12. Interface location and temperature field of the bubble at different times.</figcaption></figure>

In Fig. 13 has been plotted the temporal evolution of the equivalent bubble diameter $d_{eq}$ obtained with the different methods, i.e. DNS-NMR, DNS-SMR and DNS-IMR. Due to stability issues, the results with the DNS-LMR are not shown. One can observe that the results obtained with DNS-SMR matchs very well with the experimental results, even on long range of time. The DNS-NMR solver strongly underestimates the bubble growth, which demonstrates the strong impact of the micro-region heat flux on the overall bubble growth. On the opposite, the DNS-IMR overestimates the bubble growth rate which suggests that this coupling method suffers from a flaw to describe the flow around the bubble due to the loss of locality when accounting for the micro-region mass flux. To highlight this, the streamlines of the velocity field are plotted in the vicinity of the contact line at a given time in Fig. 14 for both DNS-SMR and DNS-IMR.

In the case of the DNS-SMR solver, as the mass flux is increased near the contact line, an aspiration of the liquid is observed in the corner formed between the liquid interface and the solid wall, where the phase change mass flow rate is stronger. On the other hand, no aspiration of the liquid is observed in the vicinity of the contact line with the DNS-IMR solver, since the streamlines show clearly that the liquid flows outward of the corner. This can be explained by the loss of locality of the DNS-IMR method since the micro-region mass flux is just integrated as an average volume source term spread out at any point of the bubble. For this reason, it cannot reproduce the effect of local liquid aspiration due to the strong vaporization in the vicinity of the contact line. Consequently, as it reproduces a better agreement with the experimental results and enables a better physical description of the flow in the vicinity
<figure><img src="image_16.png" /><figcaption>Fig. 13. Bubble equivalent diameter evolution. Comparison with the RUBI experiment. For both DNS-SMR end IMR, *f* is taken equal to 1.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Interface location, temperature field and streamlines near the contact line for the two methods at a given time.</figcaption></figure>
of the contact line, the DNS-SMR solver will be considered as more accurate for the coupling between a micro-region model and a nucleate boiling DNS solver, and it will be kept for the numerical investigations in the coming subsections.

<figure><img src="image_18.png" /><figcaption>Fig. 15. Convergence analysis for RUBI case.</figcaption></figure>

In Fig. 15, the grid convergence of simulations of the RUBI configuration, with DNS-SMR, is shown. Results are displayed for three consecutive meshes, containing 256 × 512, 512 × 1024 and 1024 × 2048 elements, which correspond respectively to $dx = 16\mu m$, $dx = 8\mu m$ and $dx = 4\mu m$. Due to significant computational cost, the 4 µm simulation is only shown on 0.4 s. This figure demonstrates a satisfactory space convergence, even with the coarsest grid, $dx = 16\mu m$.

## 4.2. Investigation on the influence of the wall conductivity for a water bubble in micro-gravity

In this subsection we take interest to the configuration of nucleate boiling in micro-gravity conditions presented by Urbano et al in [21]. It consists in a bubble in micro-gravity conditions which is immersed in a stratified thermal field so that the bottom of the bubble is vaporizing whereas the top of the bubble is condensing. As a consequence, an equilibrium radius is reached when the condensation at the top of the bubble, due to the subcooling, balances the evaporation at its bottom, due to the wall superheat. This is an interesting configuration, different from the Multiscale Boiling experiment where no steady state is reached due to the short-time of the experiment (10 s). In [21], this benchmark has been carried out with the DNS-NMR solver while accounting for the thermal conduction in the solid domain. It is reported that a steady state of the bubble volume is reached after a sufficient time, as it was expected. However, some convergence issues have also been reported for high conductivity ratio between the liquid and the solid, that is $r_λ > 10$, where $r_λ$ is defined as,

$$r_{\lambda} = \frac{\lambda_{w}}{\lambda_{liq}}. \qquad (61)$$

Indeed, the DNS-NMR solver imposes a cold spot at the contact line because the saturation temperature is imposed at the contact line by the liquid which tends to thermalize locally the solid wall. However, for high value of $r_λ > 10$, the cold spot region in the vicinity of the contact line is so tiny that space convergence can hardly be achieved, except
<figure><img src="image_19.png" /><figcaption>Fig. 16. Computations configuration used in [21]</figcaption></figure>

if much more refined grids are used in the solid domain than in the fluid domain. One of the objectives of the present benchmark is to investigate on the impact of the micro-region model on the overall bubble growth, and also to assess the effect of the micro-region model on the space convergence in a configuration with a high value of the thermal conductivity ratio, that is $r_λ = 74$. Fluid and solid properties are reported in Table 4.2, and the accommodation coefficient defined in the micro-region model is $f = 1$ which corresponds to a maximum effect of the micro-region model.

The configuration setup is described in Fig. 16. $T_{sub}$ and $T_{sup}$ are imposed all along the simulation, and initially the thermal field is adapted to ensure a constant flux in the wall and the liquid domain. In this study, we chose a Jakob ratio equal to 5, which lead to

$$ \frac{Ja_c}{Ja_e} = \frac{T_{sat} - T_{sub}}{T_w - T_{sat}} = 5 \qquad (62) $$

We chose to set $T_w = T_{sat} + 2$ K, which leads to $T_{sub} = T_{sat} - 10$ K. As we need to let $T_w$ free, we adapt the value of $T_{sup}$ by imposing conservation of heat flux in the wall, which lead to

$$ T_{sup} = T_w + \frac{L_w \Phi}{\lambda_w} \qquad (63) $$

With $L_w$ the thickness of the solid domain. The liquid domain is 3 × 3 mm², and the solid domain 1.5 mm thick. The bubble is initiated with a specific contact angle, and with a diameter of 500 µm.

<table><caption>Table 5. Physical properties used in the simulation. Fluid is water, $P_{atm}$, $T_{sat}$ = 373.12 K</caption><thead><tr><th></th><th>ρ (kg.m<sup>-3</sup>)</th><th>C<sub>p</sub> (J.kg<sup>-1</sup>.K<sup>-1</sup>)</th><th>σ (N.m<sup>-1</sup>)</th><th>L<sub>vap</sub> (J.kg<sup>-1</sup>)</th><th>λ (W.m<sup>-1</sup>.K<sup>-1</sup>)</th><th>μ (kg.m<sup>-1</sup>.s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Liquid</td><td>958</td><td>4216</td><td>5.8×10<sup>-2</sup></td><td>2.256×10<sup>6</sup></td><td>0.677</td><td>2.82×10<sup>-4</sup></td></tr><tr><td>Vapor</td><td>0.5974</td><td>2034</td><td>/</td><td>/</td><td>0.024</td><td>1.228×10<sup>-5</sup></td></tr><tr><td>Solid</td><td>1400</td><td>1500</td><td>/</td><td>/</td><td>50</td><td>/</td></tr></tbody></table>

In Fig. 17, snapshots of the thermal field and interface location can be visualized at the steady state for three simulations corresponding to three different microscopic angles. As the simulations do not consider gravity effect, the bubble shape is similar to a truncated sphere in all cases. In Fig. 18, the temporal evolution of the bubble equivalent diameter is plotted for both solvers with three consecutive meshes corresponding respectively to a cell size $dx$ = $dy$ = 12 µm, 6 µm and 3 µm and a mesh (width+(height of the fluid domain + height of the solid domain)) of
256 × (256 + 128), 512 × (512 + 256) and 1024 × (1024 + 512). With the DNS-SMR solver, a clear convergence is observed between the two finest grids, whereas with the DNS-NMR solver, the space convergence is not obvious since the gap between two consecutive meshes do not reduce when refining the grid. Such a convergence issue has already been reported in [21] for high values of $r_λ$. Thus, the coupling of the micro-region model with a DNS solver enables to solve a more realistic physical model which allows to compute a local superheat and an apparent angle that must be imposed as a boundary condition on the $\phi$ function, but it also strongly improves the space convergence of the overall solver. This can be explained by the fact that the coupling between the micro-region model and the DNS solver enables a better physical description of the heat flux in the vicinity of the contact line which is a critical issue for boiling simulations, especially in configurations involving a high value of $r_λ$.

<figure><img src="image_20.png" /><figcaption>Fig. 17. Snapshot of the bubble shapes and temperature fields with the DNS-SMR solver for three different microscopic contact angles.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 18. Mesh convergence analysis with and without microregion model for three different microscopic contact angles with $r_λ = 74$ and $f = 1$.</figcaption></figure>

To strengthen the present analysis on the effect of the thermal conductivity ratio on the DNS-SMR solver, the bubble evolution has been plotted in Fig. 19, for two different microscopic contact angles, that is $\theta_{mic} = 10^\circ$ and $\theta_{mic} = 50^\circ$ and two different solid thermal conductivities, $\lambda_w = 1 \text{ W.m}^{-1}.\text{K}^{-1}$ and $\lambda_w = 50 \text{ W.m}^{-1}.\text{K}^{-1}$ corresponding respectively to a low and a high $r_λ$ value. The local superheat at the contact line $\Delta T$ and the apparent contact angle $\theta_{app}$ are also reported in Fig. 19 for the four cases. For the sake of comparison, the results obtained with DNS-NMR solver are also reported in order to evaluate the impact of the micro-region on these configurations.

For the two cases involving a low value of $r_λ$, the influence of the micro-region on the overall bubble growth is weak, and may even be neglected for the low contact angle. The relative difference between the equivalent diameter with and without model for $10^\circ$ microscopic contact angle is 0.4%. For the higher angle, the relative difference between the solutions obtained with the two solvers is 4%. This weak effect of the micro-region in these two simulations is explained by the weak value of the local superheat at the contact line which has been computed by
<figure><img src="image_22.png" /><figcaption>Fig. 19. Bubble equivalent diameter evolution for different microscopic contact angle and wall conductivity. $e_r$ is the relative difference after 0.5 s. The accommodation coefficient is set to 1. Mesh size is 6 µm, 512 × (512 + 256) cells.</figcaption></figure>
the micro-region model. This is consistent with a low value of $r_\lambda$ which induces a low temperature discontinuity at the contact line. Low variations of the contact angle are also reported in both simulations, since it is correlated to the magnitude of the temperature discontinuity at the contact line. In the simulations involving a large $r_\lambda$ value, a more significant impact of the micro-region is observed with larger deviations between 17% and 20% on the steady equivalent diameter compared to simulations without model. It is noticeable that 20% on the diameter corresponds approximately to 70% on the bubble volume and so on the wall heat flux extracted. As expected, larger $\Delta T$ are computed in these configurations which is consistent with a larger micro-region heat flux. In these configurations, larger variations between the microscopic contact angle and the apparent contact angle are also observed for the low microscopic contact angle configuration with $\theta_{\text{app}} - \theta_{\text{mic}} = 3.75^\circ$, but remains moderate for the high contact angle configuration $\theta_{\text{app}} - \theta_{\text{mic}} = 0.6^\circ$. To summarize, the numerical results presented in this section highlight the strong coupling effect between the micro-region and the wall thermal conductivity. Indeed, it is reported that for low $r_\lambda$ values, we observe that the micro-region effect tends to be negligible, whereas it is a first order effect for large $r_\lambda$ values.

### 4.3. Sensitivity study of the numerical results to the accommodation coefficient

In this subsection, we take interest to the effect of the accommodation coefficient $f$ introduced in Eq. (16) to calculate the interface resistance in the micro-region. The correct value of the accommodation coefficient $f$ is not always perfectly known since its direct measurement is still challenging. However, it is defined between 0 and 1. A coefficient $f=0$ corresponds to an infinite thermal resistance, and thus a vanishing effect of the micro-region. In this case, $T_{\text{CL}}$ tends toward $T_{\text{sat}}$ and $\Delta P$ tends towards zero, which implies a constant interface slope and $\theta_{\text{mic}} = \theta_{\text{app}}$. On the other hand a coefficient $f=1$ corresponds to a maximum vaporization effect in the micro-region and so to a maximum heat flux at the contact line. In this case $T_{\text{CL}}$ increases and the $\Delta T$ is associated with a $\Delta P$ and consequent bending of the interface, leading to $\theta_{\text{app}} > \theta_{\text{mic}}$. Through two benchmarks involving high $r_\lambda$ values, two configurations have been simulated with different $f$ values ranging between [0.05 – 1]. The first configuration is identical to the one presented in the Sec. 4.2 in micro-gravity conditions with the microscopic contact angle $\theta_{\text{mic}} = 50^\circ$ and $r_\lambda = 74$. The second one is a classical nucleate boiling benchmark in terrestrial conditions which has already been presented in several anterior papers, as [2] or [6]. As the correct value of the accommodation coefficient is not known, the objective of this benchmark is to assess the sensitivity of the numerical simulations to this parameter. The numerical results obtained on the first configuration in micro-gravity conditions are presented in Fig. 20.

We observe that a higher value of the accommodation coefficient leads to an increase of the bubble growth rate and of the equilibrium diameter, as it was expected. The most important information is about the numerical simulations with low $f$ values which tend to the numerical results obtained with DNS-NMR solver. This demonstrates the continuity of the numerical solutions obtained with the DNS-SMR solver which tends to the results obtained with DNS-NMR solver if $f \to 0$ (or $r_\lambda \to 1$ as shown in Section 4.2). This behavior is consistent since a low value of $f$ entails a high interfacial resistance $R_f$, and therefore a low micro-region heat flux $Q_{\text{mic}}$.

The second configuration is a classical benchmark for nucleate pool boiling in terrestrial conditions. Nucleation occurs in water at saturation temperature, above a superheated wall. Solid and fluid properties are identical to [6]. A thermal boundary layer is considered in the liquid using usual Kays and Crawford correlation [65]. In this layer, we set a linear gradient, and the fluid is at saturation temperature above. In the solid, a linear thermal gradient is set as well, and the temperature boundary condition at the bottom imposes a 7 K superheat between the fluid and the solid. The Jakob number in this case is equal to 21, which is usual for nucleate pool boiling regime. The microscopic contact angle of the fluid is $50^\circ$. The fluid domain is a $4 \times 4 \text{ mm}^2$, the solid 1.33 mm thick. The simulations have been performed with a grid $512 \times (512+166)$ which corresponds to a 8 µm cell size, for which converged simulations have been already presented in [6] with the DNS-NMR solver. As the numerical results obtained with the DNS-NMR solver in [6] show a good agreement with the experimental results, this suggests that the micro-region may not have a significant impact in this configuration. The purpose of this benchmark is two-folds: (i) to check that the numerical results obtained with the DNS-SMR solver and $f \to 0$ tend to previous results obtained with the DNS-NMR solver, (ii) to investigate on the sensitivity of the numerical results to the $f$ coefficient in this configuration.

In Fig. 21, we observe that increasing $f$ leads to a decreasing departure frequency of almost 25-30% for $f=1$. This result was expected since a lower departure frequency is related to a higher bubble growth rate and so to a higher heat flux. On the other hand, the effect on the departure diameter is weak since the latter depends essentially on the buoyancy effect and on the contact angle which is weakly affected by the micro-region. We observe also in this configuration that the numerical results obtained with the DNS-SMR solver and a low accommodation coefficient tend
<figure><img src="image_23.png" /><figcaption>Fig. 20. Influence of the accommodation coefficient on the bubble diameter. The configuration is the same as presented in Fig. 16. Fluid is water, $P_{\text{atm}}$. Mesh size is 6 µm, 512 × (512 + 256) cells.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 21. Influence of the accommodation coefficient on the bubble diameter, normal gravity conditions. Mesh size is 8 µm.</figcaption></figure>
to the results obtained with the DNS-NMR solver. For water, experiments listed by Paul [66] report low *f* values, between 0.01 and 0.05. This could explain the accurate validation against experimental results reported in [6] with the DNS-NMR solver. On the opposite for FC-72 and other refrigerant liquids, computational researchers often agree for much greater values of *f* between 0.5 and 1, which strengthens the good agreement between our simulations with the DNS-SMR solver and *f* = 1 and the RUBI experiment.

These numerical results enable to demonstrate that numerical simulations with a low *f* value tend to the numerical results obtained with the DNS-NMR solver. From the numerical point of view, this important feature enables to show a smooth transition of the present DNS-SMR solver from configurations where the micro-region has a significant effect to configurations with a vanishing micro-region if *f* → 0 or *r*<sub>λ</sub> → 1. From the physical point of view, as *f* values are not always perfectly known, this provides a range of possible results for comparisons with experimental results which can be approximate in configurations where an accurate value of *f* is not known. Finally, another important information of this section is about the overall maximum impact of the micro-region which is about 20% on the bubble radius increase in the micro-gravity configuration, and 25-30% on the bubble departure frequency for the configuration in terrestrial conditions.

4.4. DNS-SMR solver in the micro-layer regime

As previously specified in the introduction, micro-region and micro-layer are two distinct phenomena which can both have a significant impact on the bubble growth. Although fundamentally different, these two phenomena are often confused which may cause some conceptual inconsistencies. The present paper is essentially focused on the heat flux at the contact line in the contact line regime. However, as a contact line also exists in the micro-layer regime (whereas it cannot be directly observed in the experiment since it is inside the micro-layer) and considering that both phenomena act at different length scales, one can assume they are spatially decoupled. This suggests that a micro-region could exist inside the micro-layer since the latter is much thicker than the micro-region. Therefore, our objective in this subsection is to show some preliminary simulations of nucleate boiling in the micro-layer regime accounting for the micro-region at the contact line inside the micro-layer.

The configuration of the simulations is described in Bureš [19], and correspond to an experiment carried out by Bucci [67]. This experiment shows some similarity with the Multiscale Boiling Experiment described in subsection 4.1, as the heat flux is imposed at the interface between the fluid and the wall. The solid is sapphire due to the IR-transparent properties, and the fluid is water, initially at saturation temperature. A Joule heat flux is set on in a titanium layer along the simulation, generating a thermal layer in both liquid and solid. After a short time (t = 88.1 ms), a bubble with a radius 10 µm is nucleated. This heat flux is variable in the r-coordinate, and is calculated so that it analytically matches the temperature values of the experiment. Our simulations are performed on a square fluid domain 1×1 mm², and a 0.5 mm thick sapphire wall of thermal conductivity λ<sub>w</sub> = 25.1 W.m<sup>-1</sup>.K<sup>-1</sup>. Due to the thin micro-layer thickness that has to be captured, simulations are performed with very refined uniform grid of 4096×(4096+2048) elements (2048 is for the solid domain), which corresponds to a cell size of 250 nm for which space convergence has been demonstrated for micro-layer simulations presented in [16]. Simulations have been performed with both DNS-SMR and DNS-NMR in order to evaluate the impact of the micro-region heat flux in such a configuration. The microscopic contact angle is taken at θ<sub>mic</sub> = 30°, and for DNS-SMR, the accommodation coefficient *f* is set to 0.005. In Fig. 22, are plotted snapshots of the overall bubble shape and of the temperature field on which a very thin liquid film between the bubble and the solid wall is computed. This liquid film is typical of the micro-layer regime. In Fig. 23, a zoom of the liquid film is plotted in order to visualize the interface shape, the temperature and velocity fields inside the micro-layer. This enables to highlight the intense vaporisation in the micro-layer, and especially in the vicinity of the contact line, along the liquid rim where the micro-region heat flux is spread.

In Fig. 24, the comparison between experimental data and two simulations is drawn in terms of bubble volume (left) and bubble lateral radius (right). One can observe that without micro-region model (DNS-NMR), the bubble volume is overestimated, whereas DNS-SMR shows a better agreement with the experimental results, and thus a lower bubble growth rate. It may seem contradictory that by adding a micro-region heat flux at the contact line, a decrease of the bubble growth rate is observed. However, the effect of the micro-region on the contact line tends to increase the contact line velocity and to shorten the micro-layer total length. Consequently, the evolution of the bubble radius is slowed down. Indeed, as in the micro-layer regime a large part of the vaporization mass flux is related to the micro-layer, if its length decreases the bubble growth rate also decreases. Although the micro-region heat flux has a significant impact on the contact line velocity, its contribution to the overall vaporization mass flux is less than 0.5%,
<figure><img src="image_25.png" /><figcaption>Fig. 22. Interface shape and temperature field at two given times with a grid 4096×(4096+2048), mesh size is 250 nm.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 23. Interface shape, velocity and temperature fields inside the micro-layer in the vicinity of the contact line, at t = 0.15 ms.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 24. Comparison between experimental results of Bucci [67], DNS-NMR and DNS-SMR, on the total bubble volume (left) and on the maximum bubble radius (right).</figcaption></figure>
whereas the micro-layer contribution is between 30 to 50%. Therefore, in such a regime the micro-region heat flux is negligible in comparison to the micro-layer heat flux, but it can have a significant impact on the micro-layer length, and so on the micro-layer heat flux. To strengthen this analysis, the temporal evolution of the micro-layer length and an instantaneous snapshot of the interface at a given time are plotted in Fig. 25 for both solvers. This highlights the strong effect of the micro-region on the micro-layer length.

<figure><img src="image_28.png" /><figcaption>Fig. 25. Temporal evolution of the micro-layer length (left), and a snapshot of the interface profile in the micro-layer with DNS-SMR (green dashed line) and DNS-NMR (blue line) at $t = 0.1$ ms (right).</figcaption></figure>

Such preliminary observations demonstrate the interest to account for the micro-region heat flux, even in configurations involving the formation of a micro-layer. Indeed, its formation, depletion and therefore heat flux are related to the contact line velocity. Development is still to be made to accurately describe the micro-region influence in such unsteady configurations. In particular, future works should focus on accurate comparisons between experiments and numerical simulations on the contact line velocity and on the micro-layer length.

## 5. Conclusion

We have presented in this paper a proper method to couple a micro-region model with a DNS solver for nucleate boiling. The main assets of our method are the following: (i)First, an accurate coupling between the micro-region heat flux and the wall thermal conduction solver is proposed. The important point is the coupling between the micro-region heat flux which is imposed as a boundary condition to the solid heat conduction solver. However, as the micro-region heat flux depends on the temperature at the contact line and thus on the temperature on the solid frontier, only an implicit coupling can be successfully achieved. This is carried out by proposing a suitable Robin boundary condition in the cells containing the contact line. Indeed, an explicit coupling based on a Neumann boundary condition imposing the heat flux computed from the previous iteration would not maintain a temperature discontinuity in the cell containing the contact line, which results in an artificial vanishing micro-region effect. (ii) The second important point is the coupling method between the boiling mass flow rate related to the micro-region to the overall numerical solver. As shown in the paper, a rough implementation of this mass flow rate in the cells containing the contact line leads systematically to numerical instabilities and local inconsistencies on the contact angle, especially for high micro-region heat flux. This method must be prohibited. Another coupling method based on a local smoothing of the boiling mass flow rate along the interface is proposed in the paper and shows satisfactory results both from the point of view of the numerical stability and of the accuracy. The proposed method, named DNS-SMR, is also compared to a coupling method integrating an average effect of the micro-region mass flow rate on the whole bubble. It is shown that this other method (DNS-IMR solver) does not describe correctly the local flow in the vicinity of the contact line and overpredicts the bubble growth, unlike DNS-SMR solver.

The interest of the proposed developments are highlighted on several benchmarks, including comparisons with experiments in micro-gravity conditions, detailed investigations on the solid thermal conductivity and accommodation
coefficient and first simulations in the micro-layer regime accounting for the micro-region heat flux at the contact line. The proposed benchmarks show a strong effect of the micro-region heat flux for high thermal conductivity and accommodation coefficient values around 1, whereas it tends to vanish for low values of the solid thermal conductivity and accommodation coefficient. Finally, in the last configuration involving the formation of a micro-layer, it is shown that, even in the micro-layer regime, the micro-region can have a significant impact on the overall bubble growth since it tends to shorten the micro-layer length.

## Acknowledgments

The authors gratefully acknowledge the French Research Agency for financial supporting in the framework of DROPSURF project ANR-20-CE05-0030. Moreover, this work has been realized thanks to the granted access to the HPC resources of CALMIP supercomputing centre under the allocation P23023.

## References

[1] Damir Juric and Grétar Tryggvason. Computations of Boiling Flows. *International Journal of Multiphase Flow*, 24(3):387-410, 1998.

[2] G. Son, V. K. Dhir, and N. Ramanujapu. Dynamics and Heat Transfer Associated With a Single Bubble During Nucleate Boiling on a Horizontal Surface. *Journal of Heat Transfer*, 121(3):623-631, 1999.

[3] Frédéric Gibou, Liguo Chen, Duc Nguyen, and Sanjoy Banerjee. A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change. *Journal of Computational Physics*, 222(2):536-555, 2007.

[4] Yohei Sato and Bojan Niceno. A sharp-interface phase change model for a mass-conservative interface tracking method. *Journal of Computational Physics*, 249:127-161, 2013.

[5] Sébastien Tanguy, Michaël Sagan, Benjamin Lalanne, Frédéric Couderc, and Catherine Colin. Benchmarks and numerical methods for the simulation of boiling flows. *Journal of Computational Physics*, 264:1-22, 2014.

[6] Grégory Huber, Sébastien Tanguy, Michael Sagan, and Catherine Colin. Direct numerical simulation of nucleate pool boiling at large microscopic contact angle and moderate Jakob number. *International Journal of Heat and Mass Transfer*, 113:662-682, 2017.

[7] B.M. Ningegowda, Z. Ge, G. Lupo, L. Brandt, and C. Duwig. A mass-preserving interface-correction level set/ghost fluid method for modeling of three-dimensional boiling flows. *International Journal of Heat and Mass Transfer*, 162:120382, 2020.

[8] Guillaume Sahut, Giovanni Ghigliotti, Guillaume Balarac, Manuel Bernard, Vincent Moureau, and Philippe Marty. Numerical simulation of boiling on unstructured grids. *Journal of Computational Physics*, 432:110161, 2021.

[9] L.C. Malan, A.G. Malan, S. Zaleski, and P.G. Rousseau. A geometric vof method for interface resolved phase change and conservative thermal energy advection. *Journal of Computational Physics*, 426:109920, 2021.

[10] Sebastian Fischer, Tatiana Gambaryan-Roisman, and Peter Stephan. On the development of a thin evaporating liquid film at a receding liquid/vapour-interface. *International Journal of Heat and Mass Transfer*, 88:346-356, 2015.

[11] M. G. Cooper and A. J. P. Lloyd. The microlayer in nucleate pool boiling. *International Journal of Heat and Mass Transfer*, 12(8):895-913, 1969.

[12] Yoshio Utaka, Yuki Kashiwabara, and Michio Ozaki. Microlayer structure in nucleate boiling of water and ethanol at atmospheric pressure. *International Journal of Heat and Mass Transfer*, 57(1):222-230, 2013.

[13] Satbioul Jung and Hyungdae Kim. Hydrodynamic formation of a microlayer underneath a boiling bubble. *International Journal of Heat and Mass Transfer*, 120:1229-1240, 2018.

[14] Yohei Sato and Bojan Niceno. A depletable micro-layer model for nucleate pool boiling. *Journal of Computational Physics*, 300:20-52, November 2015.

[15] Susann Hänsch and Simon Walker. The hydrodynamics of microlayer formation beneath vapour bubbles. *International Journal of Heat and Mass Transfer*, 2016:1282-1292, 1996.

[16] A. Urbano, S. Tanguy, G. Huber, and C. Colin. Direct numerical simulation of nucleate boiling in micro-layer regime. *International Journal of Heat and Mass Transfer*, 123:1128-1137, 2018.

[17] A Guion, S Afkhami, S Zaleski, and J Buongiorno. Simulations of microlayer formation in nucleate boiling. *International Journal of Heat and Mass Transfer*, 127:1271-1284, 2018.

[18] Lubomír Bureš and Yohei Sato. On the modelling of the transition between contact-line and microlayer evaporation regimes in nucleate boiling. *Journal of Fluid Mechanics*, 916:A53, 2021.

[19] Lubomír Bureš and Yohei Sato. Comprehensive simulations of boiling with a resolved microlayer: validation and sensitivity study. *Journal of Fluid Mechanics*, 933:A54, 2022.

[20] Sandipan Banerjee, Yang Liu, Mark Sussman, and Yongsheng Lian. Depletable micro-layer for nucleate boiling simulations in micro-gravity conditions: A new approach. *International Journal of Heat and Mass Transfer*, 190:122642, 2022.

[21] A. Urbano, S. Tanguy, and C. Colin. Direct numerical simulation of nucleate boiling in zero gravity conditions. *International Journal of Heat and Mass Transfer*, 143:118521, 2019.

[22] Robert W Schrage. *A theoretical study of interphase mass transfer*. Columbia University Press, New York, 1953.

[23] R. Raj, C. Kunkelmann, P. Stephan, J. Plawsky, and J. Kim. Contact line behavior for a highly wetting fluid under superheated conditions. *International Journal of Heat and Mass Transfer*, 55(9-10):2664-2675, 2012.

[24] P.C. Stephan and C.A. Busse. Analysis of the heat transfer coefficient of grooved heat pipe evaporator walls. *International Journal of Heat and Mass Transfer*, 35(2):383-391, 1992.
[25] Benoît Mathieu. *Études Physique, Expérimentale et Numérique des mécanismes de base intervenant dans les écoulement diphasiques en micro-fluidique*. PhD thesis, Université de Provence - CNRS - UMR 6595, 2003.

[26] Vadim S Nikolayev. Dynamics of the triple contact line on a nonisothermal heater at partial wetting. *Phys. Fluids*, page 13, 2010.

[27] Vladislav Janecek. *Evaporation à l'échelle microscopique et à haut Flux thermique*. PhD thesis, UPMC, 2012.

[28] A Ye Rednikov and P Colinet. Truncated versus Extended Microfilms at a Vapor-Liquid Contact Line on a Heated Substrate. *Langmuir*, page 12, 2011.

[29] Christian Kunkelmann and Peter Stephan. Numerical simulation of the transient heat transfer during nucleate boiling of refrigerant hfe-7100. *International Journal of Refrigeration*, 33(10):1221-1228, 2010.

[30] Vahab Mehdi, Mark Sussman, and Kourosh Shoele. Fluid-structure interaction of thin flexible bodies in multi-material multi-phase systems. *Journal of Computational Physics*, 429:110008, 2021.

[31] Grégory Huber, Sébastien Tanguy, Jean-Christophe Béra, and Bruno Gilles. A time splitting projection scheme for compressible two-phase flows. application to the interaction of bubbles with ultrasound waves. *Journal of Computational Physics*, 302:439-468, 2015.

[32] Daniel Fuster and Stéphane Popinet. An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension. *Journal Computational Physics*, 374:752-768, 2018.

[33] A Urbano, M Bibal, and S Tanguy. A semi implicit compressible solver for two-phase flows of real fluids. *J Comput Phys*, 456:111034, 2022.

[34] Youssef Saade, Detlef Lhose, and Daniel Fuster. A multigrid solver for the coupled pressure-temperature equations in an all-mach solver with vof. *Journal Computational Physics*, In press:In press, 2023.

[35] Alexis Dalmon, Kalyani Kentheswaran, Guillaume Mialhe, Benjamin Lalanne, and Sébastien Tanguy. Fluids-membrane interaction with a full Eulerian approach based on the level set method. *Journal of Computational Physics*, 406:109171, 2020.

[36] K Kentheswaran, N Dietrich, S Tanguy, and B Lalanne. Direct numerical simulation of gas-liquid mass transfer around a spherical contaminated bubble in the stagnant-cap regime. *International Journal of Heat and Mass Transfer*, 198:123325, 2022.

[37] Guillaume Mialhe, Sébastien Tanguy, Léo Tranier, Elena-Roxana Popescu, and Dominique Legendre. An extended model for the direct numerical simulation of droplet evaporation. influence of the marangoni convection on leidenfrost droplet. Submitted.

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

[39] Lucia Rueda Villegas, Romain Alis, Mathieu Lepilliez, and Sébastien Tanguy. A Ghost Fluid/Level Set Method for boiling flows and liquid evaporation: Application to the Leidenfrost effect. *Journal of Computational Physics*, 316:789-813, 2016.

[40] P. Wayner, Y.K. Kao, and L.V. LaCroix. The Interline Heat Transfer Coefficient of an Evaporating Wetting Film. *International Journal of Heat and Mass Transfer*, 19:487-492, 1976.

[41] V. Janeček and V. S. Nikolayev. Contact line singularity at partial wetting during evaporation driven by substrate heating. *EPL (Europhysics Letters)*, 100:14003, 2012.

[42] J. H. Lay and V. K. Dhir. Shape of a Vapor Stem During Nucleate Boiling of Saturated Liquids. *Journal of Heat Transfer*, 117(2):394-401, May 1995.

[43] H. K. Moffatt and B. R. Duffy. Local similarity solutions and their limitations. *Journal of Fluid Mechanics*, 96(2):299-313, January 1980.

[44] Stefan Batzdorf. *Heat transfer and evaporation during single drop impingement onto a superheated wall*. PhD thesis, Darmstadt, 2015.

[45] Romain Alis, Sébastien Tanguy, Kalyani Kentheswaran, Olivier Rouzaud, and Jean-Luc Estivalezes. Heat flux and forces acting on a vaporising droplet in a superheated vapor flow. *Int J Multiphase Flows*, In press, 2023.

[46] Annagrazia Orazzo and Sébastien Tanguy. Direct numerical simulations of droplet condensation. *International Journal of Heat and Mass Transfer*, 129:432-448, 2019.

[47] Stanley Osher and James A Sethian. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. *Journal of Computational Physics*, 79(1):12-49, 1988.

[48] Mark Sussman, Peter Smereka, and Stanley Osher. A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow. *Journal of Computational Physics*, 114:146-159, 1994.

[49] Chohong Min and Frederic Gibou. A second order accurate level set method on non-graded adaptive cartesian grids. *Journal Computational Physics*, 225:300-321, 2007.

[50] Ronald P Fedkiw, Tariq Aslam, Barry Merriman, and Stanley Osher. A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method). *Journal of Computational Physics*, 152(2):457-492, 1999.

[51] Duc Q. Nguyen, Ronald P. Fedkiw, and Myungjoo Kang. A Boundary Condition Capturing Method for Incompressible Flame Discontinuities. *Journal of Computational Physics*, 172(1):71-98, 2001.

[52] Xu-Dong Liu, Ronald P. Fedkiw, and Myungjoo Kang. A Boundary Condition Capturing Method for Poisson's Equation on Irregular Domains. *Journal of Computational Physics*, 160(1):151-178, 2000.

[53] Myungjoo Kang, Ronald P. Fedkiw, and Xu-Dong Liu. A Boundary Condition Capturing Method for Multiphase Incompressible Flow. *Journal of Scientific Computing*, 15(3):323-360, 2000.

[54] M. Sussman, K.M. Smith, M.Y. Hussaini, M. Ohta, and R. Zhi-Wei. A sharp interface method for incompressible two-phase flows. *Journal of Computational Physics*, 221(2):469-505, 2007.

[55] Benjamin Lalanne, Lucia Rueda Villegas, Sébastien Tanguy, and Frédéric Risso. On the computation of viscous terms for incompressible two-phase flows with Level Set/Ghost Fluid Method. *Journal of Computational Physics*, 301:289-307, 2015.

[56] Mathieu Lepilliez, Elena Roxana Popescu, Frederic Gibou, and Sébastien Tanguy. On two-phase flow solvers in irregular domains with contact line. *Journal of Computational Physics*, 321:1217-1251, 2016.

[57] Frederic Gibou, Ronald P. Fedkiw, Li-Tien Cheng, and Myungjoo Kang. A Second-Order-Accurate Symmetric Discretization of the Poisson Equation on Irregular Domains. *Journal of Computational Physics*, 176(1):205-227, 2002.

[58] Tariq D. Aslam. A partial differential equation approach to multidimensional extrapolation. *Journal of Computational Physics*, 193(1):349-355, 2004.

[59] J.E. Dendy. Black box multigrid. *Journal Computational Physics*, 48:366-386, 1982.

[60] S.P. MacLachlan, J.M. Tang, and C. Vuik. Fast and robust solvers for pressure-correction in bubbly flow problems. *Journal Computational Physics*, 227:9742-9761, 2008.
[61] B.P. Van Poppel, O. Desjardins, and J.W. Daily. A ghost fluid, level set methodology for simulating multiphase electrohydrodynamic flows with application to liquid fuel injection. *Journal Computational Physics*, 229:7977-7996, 2010.

[62] Rafael Borges, Monique Carmona, Bruno Costa, and Wai Sun Don. An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws. *Journal of Computational Physics*, 227(6):3191-3211, March 2008.

[63] Joab R Winkler. Numerical recipes in C: The art of scientific computing, second edition. *Endeavour*, 17(4):201, January 1993.

[64] A. Sielaff, D. Mangini, O. Kabov, M.Q. Raza, A.I. Garivalis, M. Zupančič, S. Dehaeck, S. Evgenidis, C. Jacobs, D. Van Hoof, O. Oikonomidou, X. Zabulis, P. Karamaounas, A. Bender, F. Ronshin, M. Schinnerl, J. Sebilleau, C. Colin, P. Di Marco, T. Karapantsios, I. Golobič, A. Rednikov, P. Colinet, P. Stephan, and L. Tadrist. The multiscale boiling investigation on-board the International Space Station: An overview. *Applied Thermal Engineering*, 205:117932, 2022.

[65] W.M Kays and M.E Crawford. *Convective Heat and Mass Transfer*. Mcgraw-hill edition, 1993.

[66] B. Paul. Compilation of Evaporation Coefficients. *ARS Journal*, 32(9):1321-1328, September 1962.

[67] Mattia Bucci. A theoritical and experimental study of vapor bubble dynamics in separate effect pool boiling conditions. September 2020.