# Effect of contact line evaporation during nucleate boiling studied by multiscale numerical simulation

Linkai Wei, Guillaume Bois, Vinod Pandey, Vadim Nikolayev

► To cite this version:

Linkai Wei, Guillaume Bois, Vinod Pandey, Vadim Nikolayev. Effect of contact line evaporation during nucleate boiling studied by multiscale numerical simulation. Applied Thermal Engineering, 2025, 275, pp.126690. 10.1016/j.applthermaleng.2025.126690. cea-05291094

HAL Id: cea-05291094

https://cea.hal.science/cea-05291094v1

Submitted on 30 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research Paper

# Effect of contact line evaporation during nucleate boiling studied by multiscale numerical simulation

Linkai Wei<sup>a ID</sup>, Guillaume Bois<sup>a ID,*</sup>, Vinod Pandey<sup>a</sup>, Vadim S. Nikolayev<sup>b ID</sup>

## ARTICLE INFO

Keywords:

Nucleate boiling

Microregion

Contact line

Phase change

TrioCFD

Multiscale simulation

## ABSTRACT

The present paper studies the bubble growth in saturated nucleate boiling considering its multiscale physical behaviour for the isothermal heater case. The effect of evaporation occurring in a close vicinity of the triple contact line has been included through a subgrid microregion model describing the partial wetting case. The microregion model is analysed and compared to other existing approaches. Numerical simulations have been performed using the open-source code TRUST/TrioCFD coupled with the microregion model. The coupled solver calculates the phase change taking into account the evaporation from both the microregion and the macroregion. The interface position is determined with a combined front tracking and volume-of-fluid algorithm. The microregion heat flux is applied locally, without smearing it over the bubble interface at the macroregion. First, the default phase change model in TrioCFD was validated using benchmark tests based on the Stefan and Scriven problems. Next, the multiscale coupling approach was validated against existing experimental data on nucleate boiling. Our simulations show that microregion evaporation can contribute between 15% and 26% of total bubble growth, with this contribution becoming increasingly significant as superheating intensifies. Neglecting microregion effect leads to an underestimation of the heat fluxes, dry spot sizes and bubble growth rates, thus causing longer departure times; ultimately it can lead to a failure in capturing the bubble dynamics observed in experimental studies.

## 1. Introduction

Vapour bubble nucleation and growth on the heating surface (that we call wall hereafter) favours high heat transfer to the liquid thus providing a high rate of overall heat exchange. The bubble growth modelling is a longstanding challenge [1] and is still debated in the literature. Because of the moving vapour–liquid interface, the bubble growth numerical modelling is not an easy task.

Moore and Mesler [2] are among the first researchers to advocate for the significance of evaporation at the bubble base which causes a localised drop in measured wall temperature. The bubble base is a special region because of the smallness of thickness of liquid layer that causes a small thermal resistance and thus high heat and mass fluxes. This theory was also supported by Zuber [3] while arguing in favour of the existence of two nucleate boiling regimes – one where the heat transfer is dominated by liquid agitation around isolated bubbles and another when it results from the transport of latent heat due to the evaporation at the bubble base.

The line of triple vapour–liquid–solid contact presents a separate issue because of its geometrical singularity that, if not properly treated,

causes a divergence of heat and mass fluxes in its vicinity [4,5]. Because of this, account of physical effects completely negligible at large scales but important microscopically is necessary in the CL vicinity (that we call microregion hereafter, to distinguish from macroregion, the remaining fluid part around a growing bubble, cf. Fig. 1). These effects relax the singularity thus leading to a coherent physical problem statement. In a numerical simulation, this coherence manifests itself as a grid independence, which is otherwise absent.

The bubble growth occurs thanks to evaporation from both the microregion and the macroregion [6,7], with an important contribution from microregion [8,9]. However, their relative contribution is still an open question. Nonetheless, it is extremely difficult to treat both macro and micro regions within the same numerical approach because of the large scale difference: the macroregion characteristic length scale is the bubble departure diameter that is usually millimetric. The microregion characteristic size is usually below 1 µm. For this reason, one needs a multiscale model, where the micro and macro regions are solved within different approaches and coupled. Their coupling however is a delicate issue.

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

<footnote><sup>b</sup> Service de Physique de l'Etat Condensé, CEA, CNRS, Université Paris-Saclay, 91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: guillaume.bois@cea.fr (G. Bois).</footnote>
## Nomenclature

### Acronyms

<table><tr><td>CL</td><td>contact line</td></tr><tr><td>DNS</td><td>direct numerical simulation</td></tr><tr><td>FT</td><td>front tracking</td></tr><tr><td>GFM</td><td>ghost fluid method</td></tr><tr><td>MR</td><td>microregion</td></tr><tr><td>NMR</td><td>no microregion</td></tr><tr><td>ODE</td><td>ordinary differential equation</td></tr><tr><td>VOF</td><td>volume of fluid</td></tr></table>

### Greek symbols

<table><tr><td>α</td><td>thermal diffusivity [m<sup>2</sup>/s]</td></tr><tr><td>β</td><td>volumetric coefficient of thermal expansion [K<sup>-1</sup>]</td></tr><tr><td>γ</td><td>Marangoni coefficient, [N/(m K)]</td></tr><tr><td>ΔT</td><td>wall superheating [K]</td></tr><tr><td>Δt</td><td>time step [s]</td></tr><tr><td>Δx, Δy, Δr, Δz</td><td>grid spacing [m]</td></tr><tr><td>δ</td><td>initial thickness of thermal boundary layer [m]</td></tr><tr><td>δ<sup>i</sup></td><td>Dirac delta function centred at the interface [m<sup>-1</sup>]</td></tr><tr><td>ε</td><td>relative density difference</td></tr><tr><td>ζ</td><td>length of circular arc [m]</td></tr><tr><td>η</td><td>relative microregion contribution to the bubble growth</td></tr><tr><td>θ</td><td>contact angle [rad]</td></tr><tr><td>κ</td><td>twice the mean interface curvature [m<sup>-1</sup>]</td></tr><tr><td>λ</td><td>thermal conductivity [W/(m K)]</td></tr><tr><td>μ</td><td>shear viscosity [Pa s]</td></tr><tr><td>ν</td><td>kinematic viscosity [m<sup>2</sup>/s]</td></tr><tr><td>ρ</td><td>density [kg/m<sup>3</sup>]</td></tr><tr><td>φ</td><td>radius in spherical coordinates [m]</td></tr><tr><td>σ</td><td>surface tension coefficient [N/m]</td></tr><tr><td>φ</td><td>constant of Scriven problem</td></tr><tr><td>φ</td><td>local interface slope [rad]</td></tr><tr><td>χ</td><td>phase indicator function</td></tr><tr><td>ω</td><td>constant of Stefan problem</td></tr></table>

### Latin

<table><tr><td>C<sub>p</sub></td><td>Specific heat at constant pressure [J/(kg K)]</td></tr><tr><td>C<sub>CFL</sub></td><td>Courant–Friedrichs–Lewy number</td></tr><tr><td>f</td><td>accommodation coefficient</td></tr><tr><td>ḡ</td><td>gravitational acceleration [m/s<sup>2</sup>]</td></tr><tr><td>h</td><td>local interface height (along z axis) [m]</td></tr><tr><td>K</td><td>local interface curvature in r-z plane [m<sup>-1</sup>]</td></tr><tr><td>L</td><td>latent heat of vaporisation [J/kg]</td></tr><tr><td>l</td><td>arc radius [m]</td></tr><tr><td>l<sub>c</sub></td><td>capillary length [m]</td></tr><tr><td>l<sub>s</sub></td><td>hydrodynamic slip length [m]</td></tr><tr><td>m̄</td><td>mass evaporation flux per interface area [kg/(m<sup>2</sup> s)]</td></tr></table>

<table><tr><td>n<sub>meso</sub></td><td>number of meso-region grid cell rows</td></tr><tr><td>n̄</td><td>unit normal vector to the interface (Fig. 1)</td></tr><tr><td>p</td><td>pressure [Pa]</td></tr><tr><td>Q</td><td>heat power per unit contact line length [W/m]</td></tr><tr><td>q</td><td>heat flux [W/m<sup>2</sup>]</td></tr><tr><td>R</td><td>bubble radius [m]</td></tr><tr><td>R<sup>i</sup></td><td>interfacial thermal resistance [K m<sup>2</sup>/W]</td></tr><tr><td>R<sub>v</sub></td><td>specific gas constant for vapour [J/(kg K)]</td></tr><tr><td>r</td><td>polar radius in x-y plane [m]</td></tr><tr><td>s</td><td>curvilinear coordinate along the interface contour [m]</td></tr><tr><td>T</td><td>temperature [K]</td></tr><tr><td>t</td><td>time [s]</td></tr><tr><td>ū</td><td>fluid velocity [m/s]</td></tr><tr><td>V<sub>v</sub></td><td>bubble volume [m<sup>3</sup>]</td></tr><tr><td>x, y, z</td><td>Cartesian coordinates [m]</td></tr><tr><td>x̄</td><td>local position vector [m]</td></tr></table>

### Subscript or superscript

<table><tr><td>0</td><td>initial value</td></tr><tr><td>app</td><td>apparent</td></tr><tr><td>cl</td><td>contact line</td></tr><tr><td>dep</td><td>at departure</td></tr><tr><td>i</td><td>interfacial value</td></tr><tr><td>l</td><td>liquid phase</td></tr><tr><td>meso</td><td>in the meso-region</td></tr><tr><td>micro</td><td>microscopic, microregion</td></tr><tr><td>p</td><td>at constant pressure</td></tr><tr><td>r</td><td>vapour recoil</td></tr><tr><td>ref</td><td>reference value</td></tr><tr><td>s</td><td>surface</td></tr><tr><td>sat</td><td>saturation property</td></tr><tr><td>ν</td><td>vapour phase or at constant volume</td></tr><tr><td>ω</td><td>wall</td></tr></table>

Probably the first work where such a multiscale modelling was applied is that by Son et al. [10] who modelled the bubble growth on an isothermal wall. They used the finite differences to solve the hydrodynamics problem and a level set approach for the description of the free vapour–liquid interface in the macroregion. The microregion approach of Wayner et al. [4] suitable for the complete wetting case was used. The input parameter for the model is the wall superheating. The micro and macro regions are coupled via two constant parameters, both coming from the microregion model. The first of them is the interface slope at the boundary separating the regions (the apparent contact angle). It is used as a boundary condition for the interface slope in the macroregion. The second coupled parameter is the evaporation flux integrated over the microregion; it was included into the macroregion fluid mass conservation equation, i.e smeared over the whole bubble interface.

A similar microregion model [5] suitable for the complete wetting case was used together with another coupling scheme by Kunkelmann and Stephan [11] who simulated the conjugate heat transfer problem with a non-isothermal wall with the hydrodynamic macroregion VOF approach implemented within OpenFOAM. Both the apparent angle and integrated flux that come from the microregion become variable and depend on the local superheating given by the macroregion model. The integrated heat flux does not contribute to the local fluid energy balance but rather creates a local heat sink at the CL position (for
<figure><img src="image_3.png" /><figcaption>Fig. 1. An illustration of a vapour bubble growing on a superheated wall. The inset shows an enlarged view of the microregion while the remaining part of the bubble comes under the macroregion.</figcaption></figure>

the heat conduction problem in the solid). The microregion was not considered as a special fluid region so the macroregion problem was solved in the whole fluid domain.

In practice however, the partial wetting case is much more important. Janeček and Nikolayev [12] were first to solve the truly multiscale bubble growth problem where the microregion and macroregion approaches were applied to the respective regions so the scales were entirely separated. Their microregion model suitable for the partial wetting case was studied in a series of preceding works [13–15] and validated later against experiment [16,17]. The conjugate problem of heat conduction in the liquid and in the solid was solved with the boundary element method similarly to Nikolayev et al. [18]. The curvilinear fluid wedge close to the CL was considered to be the microregion; the rest of the fluid was macroregion. The energy and mass conservation was enforced though the conditions at the microregion boundaries. However, the macroregion hydrodynamic flow was neglected so the bubble shape was quasi-static, which is not an adequate approximation in the most boiling problems.

Huber et al. [19] discussed a microregion model developed by Mathieu [20] for the partial wetting case and came to a conclusion that the contribution of microregion to the bubble growth was negligible for contact angles larger than 30°. For this reason, in this and some subsequent DNS of bubble growth, the microregion effects are neglected [21]. More recently, the same group considered smaller contact angles and coupled the above microregion model and their DNS VOF solver [22]. They highlighted a significant role of microregion for small contact angles. Their code was however unstable in the physically realistic situation of strong microregion evaporation concentrated at the computational cell adjacent to CL. They had to smear the microregion contribution over the bubble interface in the macroregion. However the smearing is likely to affect the flow in the CL vicinity and thus the overall bubble dynamics.

As it was underlined by several authors [19,23], it is essential to distinguish the microregion and a layer of liquid of a thickness ~5 µm called microlayer [24] that can be trapped between a nearly flat portion of bubble base and the wall, because the viscous effects in the layer oppose the bulk liquid inertia at high bubble growth rate.

In the recent microlayer regime simulations [25], instead of solving the microregion, the interfacial thermal resistance was used to limit the CL heat flux. An agreement with the experiment was obtained by adjusting the resistance value. Such an approach however neglects most microregion effects affecting the overall bubble growth.

In this work we present a novel physics-based multiscale coupling framework that explicitly considers a true scale separation between micro- and macro- regions. Implemented within the open-source TRUST/TrioCFD code, our approach incorporates a new microregion model that integrates recent advances in representing nanoscale effects and partial wetting. Focusing on the regime without microlayer, we simulate the growth of an isolated bubble on an isothermal wall as a first step towards the fully coupled problem. The numerical methodology adopted to resolve the macroregion comprising the governing equations and aspects of front tracking algorithm are presented in Section 2. Section 3 presents the microregion model [12] used here and compares it with that of Mathieu [20]. Section 4 details the

numerical multiscale coupling. The numerical setup for simulations is given in Section 5. Different verification tests of the model like grid independence and smoothness of multiscale coupling are presented in Section 6. These are followed by a detailed study of the bubble growth dynamics under different conditions of superheating and wettability.

## 2. Macroregion: numerical approach

A mixed numerical model of front tracking and volume of fluid algorithm with phase change implemented in TRUST/TrioCFD [26] is used for simulating the macroregion. In the FT framework, the interface is represented explicitly by a set of surface elements of a Lagrangian mesh, which moves through the fixed computational cells of an Eulerian mesh (where the volume variables are stored) [27]. The elements on the Lagrangian mesh consist of marker points which are moved based on the Eulerian velocities. This technique provides an accurate evaluation of volume fraction and results in a mass conserving methodology when combined with some elements of the VOF method.

### 2.1. Governing equations

In the present study, both liquid and vapour phases are considered as incompressible and the physical properties in the bulk phases are taken as constant. One-fluid formulation [28] considering the phase change at the interface has been used, resulting in the following form of the Navier–Stokes equations

$$ \nabla \cdot \vec{u} = -\dot{m} \llbracket 1/\rho \rrbracket \delta^l, \quad (1) $$

$$ \rho \left[ \frac{\partial \vec{u}}{\partial t} + \nabla \cdot (\vec{u}\vec{u}) \right] = -\nabla p + \nabla \cdot \left[ \mu \left( \nabla \vec{u} + \nabla^T \vec{u} - 2\dot{m} \llbracket 1/\rho \rrbracket \vec{n}\vec{n} \delta^l \right) \right] \\ + \rho \vec{g} + \sigma \kappa \vec{n} \delta^l, \quad (2) $$

$$ \rho C_p \left[ \frac{\partial T}{\partial t} + \nabla \cdot (T \vec{u}) \right] = \nabla \cdot (\lambda \nabla T) \quad (3) $$

where $\llbracket \rrbracket$ represents the interfacial jump in any variable $\phi$, such that $\llbracket \phi \rrbracket = \phi_l - \phi_v$. The temperature $T$ and pressure $p$ are both defined at the Eulerian cell centres while the normal components of the velocity vector $\vec{u}$ are defined at the faces of those cells. Interface is considered at saturation temperature, $T^l = T_{sat}$ and the temperature is continuous across the interface. The governing Eqs. (1)–(3) are the mass, momentum and energy conservation relations valid throughout the computational region. Vapour generation term $\dot{m}$ is of significance only at the interface. Due to phase change across the liquid–vapour interface, the following jump conditions [29,30] hold

$$ \llbracket \rho(\vec{u} - \vec{u}') \rrbracket = 0, \quad (4) $$

$$ \llbracket \lambda \nabla T \rrbracket \cdot \vec{n} = \mathcal{L} \dot{m}, \quad (5) $$

These jump conditions are satisfied through the introduction of source terms containing the phase change rate $\dot{m}$ in Eqs. (1)–(3). Inside the bulk phases, the continuity equation is always of the form $\nabla \cdot \vec{u} = 0$ since both liquid and vapour phases are assumed incompressible.

### 2.2. Front tracking and phase change description

The heat flux jump condition in Eq. (5) is used to determine the phase change rate $\dot{m}$ at the interface. In the present formulation, the heat flux from the vapour side is neglected with respect to that of the liquid side $q_l^l$ so Eq. (5) reduces to

$$ \dot{m} = \frac{q_l^l}{\mathcal{L}}, \quad (6) $$

with $q_l^l = \lambda_l (\vec{n} \cdot \nabla T_l)^l$. The phase change rate is used in the mass and momentum balance Eqs. (1) and (2) and is also used to construct continuous velocity fields for the interface transport and for the convection of the liquid temperature.
In the FT method, the sharp interface is explicitly defined using the moving marker mesh to ensure an accurate representation of its location, curvature, velocity, and the phase change rate. Some aspects of VOF and level set methods have also been incorporated in the present algorithm: the phase indicator function for accurate mass conservation and the signed distance to the interface, respectively. The latter distance is used to build a continuous extension of liquid temperature and velocity fields across the interface according to GFM. A continuous and continuously derivable liquid temperature extension to the vapour side is necessary to compute accurately the $T_l$ gradient for the assessment of the phase-change rate (6). The $T_l$ extension is also important to accurately enforce the saturation temperature boundary condition at the interface. Continuous velocity extension is required for interpolation of the fluid velocity (defined at the face nodes) to the positions of interface markers. This is needed for the evaluation of interfacial velocity and for the assessment of convective term in the energy Eq. (3). The phase change model of TRUST/TrioCFD with FT algorithm [27] has been successfully implemented and used in previous studies [20,31,32]. The key aspects of FT are briefly described below.

### 2.2.1. Interfacial mesh management and mass conservation

As the marker positions are updated with time, the phase indicator function also varies. Therefore, an accurate prediction of the marker positions with time is essential. A combination of FT and VOF approaches is adopted to transport the markers and to advect the phase indicator function (the vapour volume fraction) $\chi = (\rho_l - \rho)/[\rho]$.

The interfacial velocity $\vec{u}'$ is interpolated onto the marker positions $\vec{x}_j$ as $\vec{u}_j$ and the projected component of this velocity $P_n(\vec{u}_j)$ along the normal to the interface is then used to update the marker positions $\vec{x}_j^*$ after time step $\Delta t$ as

$$ \vec{x}_j^* = \vec{x}_j + \Delta t P_n(\vec{u}_j) \quad (7) $$

Additional operations on the markers are required to ensure their regular distribution [27]. They consist in the marker addition or deletion and tangential displacements. Limited normal displacements are used to increase curvature smoothness and interface stability. Subsequently, the positions are further improved through the volume corrections implemented and distributed throughout the interface. They are computed from the variation of $\chi$ with time calculated from the continuity equation as

$$ \frac{\partial \chi}{\partial t} - \frac{1}{[\rho]} \nabla \cdot (\rho \vec{u}) = 0. \quad (8) $$

The interface normal $\vec{n}$ is computed in the mixed (liquid and vapour) cells and is related to the phase indicator function as $\nabla \chi = -\vec{n}\delta'$. Here, $\vec{n}$ is physically defined at the interface and points towards the liquid phase. However, $\vec{n}$ is assessed at the centres of mixed elements before being numerically extended to the interface neighbourhood. This extension progresses iteratively from the layer of mixed cells to a determined thickness; a distance function and the liquid temperature extension are finally built from this evaluation. The total curvature $\kappa = \nabla_s \cdot \vec{n}$ is calculated as discussed in [20,27].

### 2.2.2. Interface velocity and its extension

The interface velocity $\vec{u}'$ can be determined from the mass conservation (4) by using the velocities of the pure liquid or pure vapour phases near the interface as

$$ \vec{u}' = \begin{cases} \vec{u}_l + \frac{\dot{m}}{\rho_l} \vec{n} & \text{liquid phase,} \\ \vec{u}_v + \frac{\dot{m}}{\rho_v} \vec{n} & \text{vapour phase,} \end{cases} \quad (9) $$

see Fig. 2. In this work, the Eulerian interfacial velocities are calculated based on the velocities of the pure liquid phase and the signed distance to the interface. These Eulerian interfacial velocities are interpolated to the marker positions of the Lagrangian mesh to transport the markers as described in Section 2.2.1. Note that the phase-change rate, the

<figure><img src="image_4.png" /><figcaption>Fig. 2. Velocities perpendicular to the interface in the 1D case for evaporation (i.e. positive vapour mass generation $\dot{m} > 0$). The ordinate axis is directed along $\vec{n}$.</figcaption></figure>

interface normal and the interface velocity are physically significant only over the surface of discontinuity. However, communication between the two meshes relies on interpolation. Consequently, variables in Eq. (9) must be defined and assessed in a small volume around the interface.

The volume source term in Eq. (1) is considered in the Poisson solver, scaled in proportion to the local interfacial area obtained from the intersection of the Lagrangian and Eulerian meshes. As the convection term in Eq. (3) is implemented in TrioCFD as $\nabla \cdot (T\vec{u})$ rather than $\vec{u} \cdot \nabla T$, this implementation requires a rigorously divergence-free velocity extension, even at the liquid-vapour interface, to ensure accurate advection of the temperature field. This velocity extension is achieved using a method similar to that proposed by Gennari et al. [33], where the expansion term in the projection step is shifted towards pure vapour cells near the interface. This approach not only guarantees a rigorously divergence-free velocity field at the liquid-vapour interface but also eliminates the need to solve an additional Poisson equation. It has demonstrated excellent performance in boiling simulations [33,34].

## 2.3. Numerical algorithm and time-stepping criteria

The solution algorithm follows the steps below:

1. Initialize the computational domain with $\vec{u} = 0$ and with a prescribed temperature field. Choose the vapour domain of a specified shape; obtain the indicator function $\chi$ over the domain and the initial front mesh.

2. Compute the mass transfer rate $\dot{m}$ in the mixed computational cells in the macroregion (with Eq. (6)), and in the micro- and meso- regions as discussed in Section 4.

3. Update the velocity and pressure fields by using the predictor-corrector approach (solution of Eqs. (1) and (2)). The surface tension term is calculated as a prerequisite using the local curvature obtained from the calculated normals at the interface elements.

4. Update the interface position by Eq. (7) using the extended interfacial velocity defined by Eq. (9). Correct the position to ensure mass preservation based on Eq. (8) and update the phase indicator function $\chi$.

5. Update the extended temperature field for each phase (details are described in [32]) and solve the conservative form of the energy Eq. (3) for each phase by using the continuous extension of the updated velocity field. The GFM is used to enforce the saturation temperature at the interface and to accurately assess the phase-change rate $\dot{m}$ in the step 2.

6. Repeat the steps 2–5 for the next time step.

The global time step restrictions during the simulations performed in the present study are imposed considering the constraint

$$ \Delta t \le \min \left[ \left( \Delta t_{\text{conv}}^{-1} + \Delta t_{\text{visc}}^{-1} \right)^{-1}, \Delta t_{\text{cap}} \right] \quad (10) $$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Schematics of the Stefan problem and initial $T_v$ distribution.</figcaption></figure>

<table><caption>**Table 1**
Properties of water at 1 bar.</caption><thead><tr><th></th><th>ρ,<br>kg/m<sup>3</sup></th><th>μ,<br>μPa s</th><th>λ,<br>W/(m K)</th><th>β,<br>1/K</th><th>C<sub>p</sub>,<br>kJ/(kg K)</th></tr></thead><tbody><tr><td>Liquid</td><td>958.37</td><td>280</td><td>0.679</td><td>0.0029</td><td>4.21</td></tr><tr><td>Vapour</td><td>0.597</td><td>12.27</td><td>0.025</td><td>0.00075</td><td>2.077</td></tr></tbody><tfoot><tr><td colspan="6">σ = 58.9 mN/m, L = 2.256 MJ/kg, T<sub>sat</sub> = 373 K</td></tr><tr><td colspan="6">γ = 0.192 mN/(m K), R<sub>v</sub> = 461.89 J/(kg K)</td></tr></tfoot></table>

where $\Delta t_{conv} = c_{CFL} \Delta x / u_{max}$ [35], $\Delta t_{visc} = \rho_v \Delta x^2 / \mu_v$ [36] and $\Delta t_{cap} = c_{cap} \sqrt{(\rho_v + \rho_l)/(8\sigma) \Delta x^{3/2}}$ [37,38], corresponding to the convective, viscous and capillary terms, respectively. It is observed that the capillary time step, $\Delta t_{cap}$, is the most restrictive, thereby limiting the simulation time step. In this study, the constants $c_{cap} \approx 0.1$ and $c_{CFL} = 1$ are employed.

## 2.4. Validation of the numerical approach

The verification of proper implementation of the numerical procedure in the phase change solver through various test cases has been performed and the numerical accuracy is checked by comparing the results with the corresponding analytical solutions.

### 2.4.1. Stefan test case

This is a classical 1D test case considered by various researchers [39–41] for verification of the phase change numerical algorithm.

The schematic of the problem is shown in Fig. 3. The interface position is denoted by $x^l(t)$ at time $t$, with the initial interface position $x_0^l = 0.05$ mm. The vapour phase (depicted in blue on the left side of the tube of length $L = 1$ mm) is heated by an isothermal wall at $x = 0$, where $T_v = T_{sat} + \Delta T$ while the liquid phase remains at the saturation temperature $T_{sat}$. Heat flux from the vapour to the interface induces liquid evaporation, causing the interface to move towards the liquid. As the interface displaces, a thermal distribution develops in the vapour phase, while the liquid remains isothermal at $T_{sat}$.

For the 1D numerical test, the outflow boundary condition (constant pressure and zero normal derivatives at the boundary for all other quantities) was imposed at $x = L$. The symmetry conditions (zero normal derivatives and zero normal velocity) was set to the y-boundaries (top and bottom). To verify our phase change numerical algorithm, we compared our computational results with the analytical solutions (Appendix A) using the properties of water at 1 bar (Table 1). The wall superheating is characterised by the dimensionless Stefan number, $Ste = C_{p\nu} \Delta T / L$. For the reference test case, the wall superheating is set to $\Delta T = 10$ K, which corresponds to $Ste_{ref} = 9.2065 \times 10^{-3}$.

Once $Ste$ is chosen, one can determine the time $t_0$ by solving Eq. (A.1) as the time $t = t_0$ corresponding to the interface position $x^l = x_0^l$. For instance, for $Ste = Ste_{ref}$, $t_0 = 6.77$ ms. This instant is taken as initial for the numerical simulation. A linear temperature profile is

<figure><img src="image_6.png" /><figcaption>Fig. 4. Interface position evolution for Stefan problem calculated for $Ste = Ste_{ref}$ using three different grid spacings of $\Delta x = 50$ µm, 25 µm, and 12.5 µm.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Spatial distribution of vapour temperature at $t \approx 74t_0$ for different $Ste$ values. The temperature $T_v$ and the coordinate $x$ are scaled with $\Delta T$ and $x^l(t = 74t_0)$, respectively. The grid spacing $\Delta x = 12.5$ µm is used in all the cases.</figcaption></figure>

imposed in the vapour phase and both phases are initially in a quiescent state, which is an extremely accurate approximation of the theory.

The numerically calculated interface position and its theoretical value are both presented in Fig. 4, showing excellent agreement across all the considered grid spacings.

To further verify the accuracy of calculations, temperature distributions are compared for three cases: $Ste = Ste_{ref}$, $5000Ste_{ref}$ (an extremely high value), and an intermediate value of $5000Ste_{ref}$. To compare the temperature evolution, $T_v$ is scaled by $\Delta T$. Due to the different Stefan numbers, each case has a different time $t_0$, determined as discussed above. Simulations are conducted for each case over the time range $t \in (t_0, 74t_0)$.

The resulting temperature distributions for the different Stefan numbers are shown in Fig. 5. The distributions closely match the analytical results at realistic $Ste$ values, with very slight deviations observed for the extreme $Ste$ value. These results demonstrate the overall precision of TrioCFD in capturing the temperature profiles.

### 2.4.2. Scriven test case

We consider the growth of a spherically symmetric bubble in the initially quiescent and superheated infinite liquid. This is a well known problem [42] which was used as a benchmark test case for validating several phase change solvers [41,43]. As compared to the 1D Stefan test case, this problem is more challenging since it accounts for the interface curvature and is sensitive to the thermal boundary layer development around the bubble in the liquid phase (as shown through the schematics in Fig. 6).
<figure><img src="image_8.png" /><figcaption>Fig. 6. Schematics for the Scriven test case.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Final ($t = 4t_0$) liquid temperature distributions in the radial direction for the Scriven test case computed with a grid spacing of 25 µm at a superheating of $\Delta T = 1$ K.</figcaption></figure>

Since the liquid is incompressible, the mass balance results in an explicit expression for the liquid velocity $u_l$. The liquid velocity and temperature are governed by the equations

$$u_l = \epsilon \frac{R^2 \dot{R}}{\rho^2}, \\ \frac{\partial T_l}{\partial t} + u_l \frac{\partial T_l}{\partial \rho} = \alpha_l \left( \frac{\partial^2 T_l}{\partial \rho^2} + \frac{2}{\rho} \frac{\partial T_l}{\partial \rho} \right), \quad (11)$$

where $\epsilon = 1 - \rho_v/\rho_l$. The solution should satisfy the initial and boundary conditions

$$u_l(t=0) = 0, \quad (12)$$

$$T_l(\rho = R) = T_{sat}, \quad (13)$$

$$T_l(\rho \to \infty) = T_{sat} + \Delta T_{\infty}. \quad (14)$$

It turns out that such a problem has a self-similar solution $T_l = T_l[\rho/R(t)]$ presented in Appendix B.

To validate our DNS approach, the simulation has been performed within a 2D axisymmetric numerical approach for water as the working fluid (cf. Table 1 for its parameters) and $\Delta T_\infty = 1$ K. It corresponds to $J_a \approx 3$, where the Jacob number is defined as

$$J_a = \frac{\rho_l C_{pl} \Delta T_\infty}{\rho_v \mathcal{L}}. \quad (15)$$

The outflow boundary conditions, along with $\Delta T_\infty$, are imposed at the domain boundaries located at $r = 4$ mm and $z = 4$ mm. A numerical

<figure><img src="image_10.png" /><figcaption>Fig. 8. Bubble radius evolution for the Scriven test case at a superheating of $\Delta T = 1$ K using three different grid spacings of $\Delta r = \Delta z = 100$ µm, 50 µm, and 25 µm.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Relative error in the final radius $R$ at $t = 4t_0$. The slanted line shows the slope -1.</figcaption></figure>

value of the Scriven constant $\varphi \approx 3.32$ is obtained by solving Eq. (B.8). Since one cannot start the simulation from $R=0$ like in the theory, the initial bubble radius $R_0 = 1$ mm is imposed at $t = t_0 = 0.1346$ s (obtained from Eq. (B.1)). The initial temperature distribution in the liquid phase is determined using the analytical solution (B.7) for $t = t_0$. The simulation was performed until $t = 4t_0$, at which $R = 2R_0$. The same configuration can be found in previous studies for phase change model validation [34,43].

The mesh independence analysis is presented in Fig. 8. The relative error $|R - R_{theo}|/R_{theo}$ in the simulated radius $R$ at $t = 4t_0$ with respect to its theoretical value $R_{theo}$ is shown in 9. The first-order mesh convergence, i.e. the error inversely proportional to the total grid cell number, is observed. The temperature profile at the end of simulation agrees with the analytical solution as shown in Fig. 7. The initial and final radius can be seen as points of deviation of temperature from $T_{sat}$. As expected, $T_l(\rho)$ distribution shows a propagating boundary layer near the bubble interface, the size of which is proportional to that of the bubble.

### 3. Microregion model

Some macroregion effects can be neglected in the microregion. They are linked to the small thickness of the liquid. As a consequence, the inertial effects can be ignored so the Stokes equations hold. One can then use the generalised lubrication approach [12,17,20] valid for the interface slope $\phi < 90^\circ$. The normal viscous stresses at the vapour-liquid interface can be neglected with respect to the surface tension
too so the conventional normal stress balance reduces to the Laplace equation $\llbracket p \rrbracket = -\kappa\sigma$ that defines the local interfacial pressure jump. The interfacial curvature can be written as $\kappa = K + \sin\phi/r$, where K is the bubble contour line curvature in the $r-z$ plane. As we will see in the following, the curvature itself does not matter but its spatial derivative does. Since the CL radius $r_{cl}$ is typically much larger than the microregion size, the term $\sin\phi/r$ varies much weaker than K and can be approximated by a constant contribution $\sin\theta_{app}/r_{cl}$ absorbed by $p_l$ so it becomes $p_l^*$. The microregion model can then be described in the 2D approximation in the $r-z$ plane and the pressure balance becomes

$$ \Delta p = K\sigma, \quad (16) $$

where $\Delta p = p_v - p_l^*$. However several more effects are important in the microregion and need to be accounted for. Here we give only a short description. More detailed review is given by Nikolayev [44].

### 3.1. Relevant microscale phenomena

First of such effects is the vapour recoil. Every fluid molecule evaporated from the liquid interface causes a recoil force analogous to that created by the gas emitted by a rocket engine. It pushes the interface towards the liquid side in the normal direction. This force appears because the fluid necessarily expands while transforming from liquid to gas phase. The vapour recoil is important wherever the evaporation rate is high, which is the case of microregion. From the interfacial momentum balance and the mass balance (4), one can obtain [45] the vapour recoil force $\vec{p}_r$ per unit interface area

$$ \vec{p}_r = -\vec{n}\dot{m}^2 \llbracket 1/\rho \rrbracket. \quad (17) $$

Note that $\vec{p}_r$ is directed towards the liquid both for evaporation ($\dot{m} > 0$) and condensation ($\dot{m} < 0$) cases. The vapour recoil pressure enters the normal pressure balance equation as [45]

$$ K\sigma = p_r + \Delta p, \quad (18) $$

where $p_r = |\vec{p}_r|$. One notes that the vapour recoil term is present in Eq. (2) but is not necessarily important in the macroregion.

As for the thermal part of the problem, it can be simplified by neglecting the thermal inertia of the thin liquid layer formed near the CL, which means that one can consider the stationary problem. For a straight wedge of the opening angle $\phi$, $T_l$ varies linearly [46] along the wedge arc of the length $\zeta$, which causes

$$ q_l^i = \lambda_l \frac{T_w - T^i}{\zeta}. \quad (19) $$

For a curved wedge, one can assume that at each point A of the interface the temperature distribution is defined by the same Eq. (19); the equivalent “local” straight wedge is defined by the tangent to the interface at the point A (Fig. 1, dotted line) with the tangent slope $\phi$. One can introduce the length $l$ from the point A to the apex of the local wedge, so

$$ \zeta = l\phi. \quad (20) $$

The local saturation temperature (i.e. the temperature corresponding to the local equilibrium occurring at a given point s at the interface) is different from the system saturation temperature $T_{sat}$ (defined by the system pressure $p_v$) because of the Kelvin effect. In addition to the Kelvin term that involves the interfacial pressure jump [47], a vapour recoil term appears in the interface temperature expression [48]

$$ T_{sat,local} = T_{sat} \left( 1 + \frac{\Delta p}{\mathcal{L}\rho_l} - \frac{\dot{m}^2}{2\mathcal{L}} \llbracket \rho^{-2} \rrbracket \right). \quad (21) $$

Another effect related to the strong evaporation is the kinetic interfacial resistance [5]

$$ R^i = \frac{2-f}{2f} \frac{T_{sat} \sqrt{2\pi R_v T_{sat}} \llbracket \rho \rrbracket}{2\mathcal{L}^2 \rho_l \rho_v}. \quad (22) $$

that limits the evaporation rate when it becomes to be comparable to the thermal velocity of vapour molecules. As a consequence, the interfacial temperature $T^i$ begins to deviate from the local saturation temperature $T_{sat,local}$ at high heat flux:

$$ T^i = T_{sat,local} + R^i q_l^i. \quad (23) $$

At relatively low heat fluxes observed in the macroregion the second term is negligible because $R^i$ value is usually quite small. In Eq. (22), $f$ represents the accommodation coefficient. A unit accommodation coefficient ($f = 1$) is used in the present work, which results in $R^i \approx 6.38 \times 10^{-8} \text{ K m}^2/\text{W}$.

Because of the interfacial temperature variation in the CL vicinity, the Marangoni effect needs to be accounted for. The surface tension variation along the surface reads

$$ \frac{\partial \sigma}{\partial s} = -\gamma \frac{\partial T^i}{\partial s}, \quad (24) $$

where $\gamma = -d\sigma/dT$ is assumed to be constant. Finally, one can account for the hydrodynamic slip that can occur in the CL vicinity because of the high hydrodynamic stress. It is characterised by the slip length $l_s$ so that the no-slip condition applied in the macroregion at the wall $z = 0$ is replaced by the Navier slip condition

$$ u_x = l_s \frac{\partial u_x}{\partial z}. \quad (25) $$

Throughout this work we use the value $l_s = 10 \text{ nm}$ which is a reasonable value for the partial wetting case [17]. As for the surface forces that can be accounted through the disjoining pressure [4], their effect is negligible for the partial wetting case [14].

Apart from large local evaporation rates in the microregion mentioned above, the most notable effect of the microregion on the macroscopic phenomena is the emergence of the apparent $\theta_{app}$ contact angle that can be much larger than the true (microscopic) value $\theta_{micro}$. This effect has been discovered by Wayner et al. [4] and studied by many other researchers, e.g. [5]. It was however erroneously attributed to the impact of disjoining pressure instead of the vapour recoil and the viscous pressure drop in the wedge [44]. In other words, a large $\theta_{app}$ appears even if the disjoining pressure is neglected.

### 3.2. Microregion theory

The conventional lubrication theory [47,49] has a severe limitation imposed on the interface slope that in practice should be smaller than 30° to obtain a reasonable accuracy. It happens that there is a mean to extend the range of the admissible slopes by describing the interface in a parametric form $h = h(s)$ where $s$ is the coordinate varying along the interface; $s = 0$ at the CL (Fig. 1). One can write

$$ \frac{\partial h}{\partial s} = \sin \phi. \quad (26) $$

Such a generalisation of the lubrication approach [12,17,20] for high interface slopes results in the following equation:

$$ \frac{\partial}{\partial s} \left\{ \frac{1}{G(\phi)} \left[ \frac{\partial \sigma}{\partial s} \left( \frac{\zeta^2}{2} + \zeta l_s \right) + \left( \frac{\zeta^3}{3} + \zeta^2 l_s \right) \frac{\partial \Delta p}{\partial s} \right] \right\} = -\mu \frac{\dot{m}}{\rho_l}, \quad (27) $$

where

$$ G(\phi) = \frac{4\phi^3}{3[\cos\phi \sin\phi - \phi \cos(2\phi)]} $$

is the correction factor to the conventional lubrication theory that corresponds to the limit $\phi \to 0$ where $G = 1$. Eq. (18) becomes

$$ \sigma \frac{\partial \phi}{\partial s} = \Delta p - \dot{m}^2 \llbracket 1/\rho \rrbracket \quad (28) $$

since within such a parametrisation, the bubble contour line curvature writes

$$ K = \frac{\partial \phi}{\partial s}. \quad (29) $$
The height h is related to the arc length ζ introduced above as

$$h = \zeta \frac{\sin \phi}{\phi}. \qquad (30)$$

Evidently, the formulation becomes complete if one accounts for the Eqs. (6), (19)–(23).

The set of Eqs. (26)–(28) is of the fourth order. To solve it, one thus needs four boundary conditions. Two of them are geometrical constraints defined at the CL,

$$\begin{aligned} \zeta(s \to 0) &= 0, \\ \phi(s \to 0) &= \theta_{\text{micro}}. \end{aligned} \qquad (31)$$

The third boundary condition is given by the matching of the micro- and macro- regions. The interfacial pressure jump in microregion is expected to be much smaller than in microregion because of much smaller K and m. For this reason, the condition

$$\Delta p = 0 \qquad (32)$$

valid at $s \to \infty$ is applied. In a numerical treatment, the condition (32) is imposed at a finite microregion boundary point M (cf. Fig. 1) with the coordinate $s = s_M$.

The fourth boundary condition comes from the regularity of all the variables at the CL, in particular from the temperature continuity

$$T^i(s \to 0) = T_{\text{iv}}(r \to r_c). \qquad (33)$$

The authors of the works [19,22] admitted a temperature discontinuity at the CL by arguing that it could be related to the interfacial thermal resistance, the effect of which is introduced in Eq. (23). One mentions however that its combination with Eq. (19) yields the expression

$$T^i = T_{\text{sat,local}} + R^i \frac{\zeta}{\lambda_l} + R^i. \qquad (34)$$

In the limit $ζ \to 0$ valid at the CL, the above expression yields Eq. (33). The introduction of interfacial resistance cannot thus lead to a temperature discontinuity at the CL. Moreover, the above reasoning shows that Eqs. (19) and (23) imply the temperature continuity at the CL. Once it is established, the asymptotic argument [13,17] shows that

$$q_l^i(s \to 0) = 0, \qquad (35)$$

so the interfacial mass flux $m(s \to 0)$ is zero too. By using Eqs. (21), (23), one comes to the fourth boundary condition

$$\Delta p(s \to 0) = \frac{\mathcal{L} \rho_l \Delta T}{T_{\text{sat}}}. \qquad (36a)$$

Alternatively, one can use its equivalent form [47]

$$\left. \frac{\partial \Delta p}{\partial s} \right|_{s \to 0} = 0. \qquad (36b)$$

Note the role of the Kelvin effect that is responsible for the explicit form of the fourth boundary condition. It can be proven mathematically [15] that a unique regular solution of a similar (slightly reduced) problem exists.

The above model is close to that of Mathieu [20] (used by Huber et al. [19],Torres et al. [22]) except for the fourth boundary condition that they write as

$$\lim_{s \to 0} \int_0^s q_l^i(s) ds = 0. \qquad (37)$$

Such a condition admits a nonzero $q_l^i$ value at the CL as obtained numerically by Huber et al. [19]. We argue hereafter that this boundary condition is insufficient to close their system of equations and render their formulation well posed. Since the less constraining condition (37) is compatible with (35), its adoption leads to two alternatives. Either it results in the same solution as (35) or the solution based on (37) is not unique. The latter alternative would invalidate such an approach so one must recognise the temperature continuity at the CL as the only

viable choice. From the numerical perspective, the conditions (35) and (37) are not equivalent. Since the latter condition admits a nonzero $q_l^i$ value at the CL, it can result in important numerical errors as shown below by comparing the results of both approaches.

To compute the interface point abscissas (r coordinates), one more equation, decoupled from the others, should be solved

$$\frac{\partial r}{\partial s} = \cos \phi \qquad (38)$$

with the boundary condition $r(s \to 0) = r_{cl}$.

### 3.3. Microregion results

We present here a comparison with the results of Huber et al. [19] as it is the only published microregion results (other than ours) on large angles. They used the generalised lubrication formulation of Mathieu [20]. The same parameters (Table 1) are used here. The spatial variation of several important quantities along the vapour-liquid interface is presented in Figs. 10. to make a viable comparison of the microregion behaviour with Huber et al., we use the same $\theta_{micro} = 49.5^\circ$ although the resulting $\theta_{app}$ are different in both models.

The pressure jump variation (Fig. 10(a)) is the most important because $\Delta p$ is the main contribution to the slope variation defining the apparent contact angle $\theta_{app}$, cf. Eq. (28). One can see that $\Delta p$ coming from our model smoothly approaches the asymptotic value (36a) at $s \to 0$.

Fig. 10(b) demonstrates that the interface slope $\phi$ saturates with increasing $s$. Therefore, one can introduce the microregion size $s_M$ as an $s$ value beyond which the $\phi$ variation (i.e. the curvature, cf. Eq. (29)) is much weaker than near the CL. In this section we use $s_M \approx 1 \mu m$ to comply with Huber et al. [19]. However the value $\approx 5 \mu m$ will be used in the rest of the paper.

The pressure jump calculated within the present model is larger than that of Huber et al. The difference probably comes from the imprecise boundary condition (37) used by them as discussed above. This causes a quite strong difference in the interface slope. While Huber et al. [19] obtained $\theta_{app} - \theta_{micro} \approx 0.5^\circ$, our model results in an order of value larger difference $\theta_{app} - \theta_{micro} \approx 5^\circ$ (Fig. 10(b)). This result thus invalidates a conclusion of Huber et al. [19] about a negligible effect of microregion for moderate contact angles.

One can compare the $q_l^i$ microregion behaviour with the heat flux that would occur if the macroregion expression (6) was applied:

$$q_l^i = \frac{\Delta T \lambda_l}{s \theta_{app}}. \qquad (39)$$

In Fig. 10(c) the corresponding curve is labelled “NMR” (no microregion). It describes well the flux for $s \ge 0.4 \mu m$ but is divergent as $s \to 0$. It means that if no microregion model is applied, the flux obtained with the macroregion approach is limited by the size of the grid cell to which the CL belongs. This provides an effective cutoff, which leads to the grid dependence and, ultimately, in an erroneous bubble growth rate.

Huber et al. [19] found that their microregion heat flux agrees with the model that they called “without microregion” where the straight wedge with the opening angle $\theta_{app}$ is assumed and only the interfacial thermal resistance $R^i$ is accounted for. Note that such an name is misleading as it suggests that the macroregion expression (39) is applied to the microregion. For this reason, we call such an approach “R^i only”, as illustrated in Fig. 10(c). By posing $T_{sat,local} = T_{sat}$ in Eq. (23) and coupling it with Eq. (19), one gets

$$q_l^i = \frac{\Delta T \lambda_l}{s \theta_{app} + R^i \lambda_l}. \qquad (40)$$

To obtain an apparent angle $\theta_{app} = 50^\circ$ within our model, we had to choose $\theta_{micro} = 45.2^\circ$. The respective flux curve (Fig. 10(c)) coincides with both “R^i only” and NMR curves for large $r - r_{cl}$ while the flux for $\theta_{micro} = 49.5^\circ$ is lower. Such a value of $\theta_{micro}$ results in a larger $\theta_{app} = 54.2^\circ$ (Fig. 10(b)) so the liquid layer is thicker.
<figure><img src="image_12.png" /><figcaption>Fig. 10. Variation of several variables along the wedge interface in the microregion calculated for $\Delta T = 7$ K.</figcaption></figure>

The “R' only” curve coincides with that of Huber et al. However this curve differs from the behaviour (35) of our approach. Note that because of the imprecise boundary condition used by Huber et al., their flux does not go to zero at the CL contrary to the mathematical analysis of the microregion model [13]. Also, note that for the large $\theta_{micro}$ case, as this difference concerns only a very narrow vicinity of CL, it should

<figure><img src="image_13.png" /><figcaption>(b) $Q_{micro}$ calculated for $s_M = 1$ µm (the value used in Huber et al. [19]).</figcaption></figure>

<caption>Fig. 11. $\theta_{app}$ and $Q_{micro}$ computed for $\Delta T = 7$ K as functions of $\theta_{micro}$.</caption>

not contribute to the integrated heat flux (i.e., the heat power per unit perimeter of the CL)

$$Q_{micro} = \int_{0}^{s_M} q_{i}^{l}(s) ds, \qquad (41)$$

which is one of the main parameters characterising the microregion. One can compare the integrated result with the case that neglects fluid flow and nanoscale effects (“R' only” case), which is obtained by integrating the flux (40) from zero to $s_M$:

$$Q_{micro} = \frac{\Delta T \lambda_{j}}{\theta_{app}} \ln \left( \frac{s_M \theta_{app}}{R' \lambda_{j}} + 1 \right). \qquad (42)$$

In investigations performed for a given solid-fluid couple, $\theta_{micro}$ is imposed by their properties while $\theta_{app}$ remains a priori unknown, limiting the predictive capacity of Eq. (42). One can however determine $\theta_{app}$ for fixed $\theta_{micro}$ and $\Delta T$ from the microregion simulation (Fig. 11(a)). When the deduced $\theta_{app}$ is substituted into Eq. (42), the resulting $Q_{micro}$ can be compared with Eq. (41). Fig. 11(b) shows a good agreement for $\theta_{micro} > 20^{\circ}$, which means that Eq. (42) can be used for the estimation of $Q_{micro}$ at comparison with the experiment where $\theta_{app}$ is fixed. It should be however applied with caution for small angles as it underestimates $Q_{micro}$. This conclusion on the flux agrees with Huber et al. One can see however the difference between both approaches on the contact angle. For $\theta_{app} = 50^{\circ}$, Huber et al. [19] compute $\theta_{micro} = 49.5^{\circ}$ (black dot in Fig. 11(a)) while the present model results in $\theta_{micro} = 45.2^{\circ}$ (red dot).

The variation of $\theta_{app}$ and $Q_{micro}$ with $\Delta T$ are displayed in Fig. 12 for $\theta_{micro} = 15^{\circ}$ used below in bubble growth simulations.

The final remark on the microregion model concerns the nonphysically small length scales that can be found in Figs. 10. They are shown to illustrate the asymptotic behaviour of the variables.

## 4. Multiscale coupling

Fig. 13 illustrates the coupling strategy applied in this study. The point M marks the end of microregion. Depending on the CL position, the microregion may be located either inside the macroregion grid cell
<figure><img src="image_14.png" /><figcaption>Fig. 12. Variation of $\theta_{app}$ (black solid lines) and $Q_{micro}$ (computed for $s_M = 5$ µm, dashed blue lines) with $\Delta T$ for $\theta_{micro} = 15^{\circ}$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. Schematic of the multi-scale coupling. Three regions are illustrated: the microregion (not to scale, in red), the meso-region (purple), and the macro-region (white). The point M denotes the boundary between micro- and meso- regions.</figcaption></figure>

containing the CL or span both the CL cell and an adjacent cell. In the example of Fig. 13, the microregion spans two cells necessitating a redistribution of the heat flux $Q_{micro}$ between them. A straightforward approach is a flux distribution proportional to the curvilinear length $s$ of the vapour-liquid interface in each of the cells. However, such an approach is physically incorrect because the flux distribution is in reality far from linear (cf. Fig. 10(c)). By following the long-range flux behaviour (39) that results in a logarithmic tendency of its integral, we distribute $Q_{micro}$ between the cell CL (that contains the CL) and the cell M (that contains the point M) proportionally to the logarithm of the interface height:

$$Q_{micro,CL} = Q_{micro} \frac{\ln s(L)}{\ln s_M} \approx Q_{micro} \frac{\ln z(L)}{\ln z(M)} \quad (43)$$

where $Q_{micro,CL}$ is a part of $Q_{micro}$ applied to the cell CL and the point L is the rightmost interfacial point belonging to it (Fig. 13). In the classical GFM, at least one neighbouring cell (better, several cells) containing pure liquid (i.e. not containing a portion of interface) in the normal to the interface direction is needed to determine the temperature gradient. Therefore, uncertainties arise in calculating the temperature gradient (and, consequently, the heat flux) within macroscopic cells near CL where the liquid layer is thinner (cf. Fig. 13). To address this issue, we introduce an intermediate “meso-region” to increase the precision of computation of heat flux. To achieve the continuity between the micro- and meso- regions, the heat flux $q'_l$ should be treated there similarly to the microregion, i.e. to comply with Eq. (19). In contrast to the microregion, the microscale effects (Section 3.1) can be neglected in the meso-region so $T^l = T_{sat}$; Eq. (19) reduces to

$$q'_l = \frac{\lambda_l \Delta T}{\phi_l} \quad (44)$$

To identify all the cells defining the meso-region, we search for all liquid and vapour interface cells whose height from the wall is less than a specified number of cells, $n_{meso}$. The ending interfacial point of the meso-region (Q in Fig. 13) is always located at the right face of a grid cell of the row $n_{meso}$. In this example ($n_{meso} = 2$), in the meso-region, there are three grid cells containing the interface. For the

cell containing the interface part (N,P), the total heat power per unit perimeter length (situating in the x-y plane) can be calculated as

$$Q_{meso}(N, P) = \int_{N}^{P} q'_{l}(s) ds = \frac{\lambda_{l} \Delta T}{\phi(N, P)} \ln \left[ \frac{s(P)}{s(N)} \right] \\ = \frac{\lambda_{l} \Delta T}{\phi(N, P)} \ln \left[ \frac{z(P)}{z(N)} \right], \quad (45)$$

where the interfacial slope $\phi(N, P)$ is constant within the cell. This analytical integration provides a more accurate assessment of heat flux compared to the default finite difference algorithm of TrioCFD.

The integration of the microregion model into the macro-scale CFD framework involves three steps:

1. Enforcement of the apparent contact angle: The apparent contact angle $\theta_{app}$ predicted by the microregion model is enforced as a boundary condition for the interface at the cell containing CL.

2. Adjustment of wall heat flux in micro- and meso- regions: An appropriate heat flux is imposed at the bottom face of the cells containing these regions. For the cell ‘M’ containing both the micro- and meso-regions (i.e. containing the point M), the heat flux contributions from both regions are added. The heat flux $q_M$ at the bottom face of the cell M is

$$q_M = \frac{2\pi r_{cl} Q_{micro,M} + \pi [r(M) + r(N)] Q_{meso}(M, N)}{A_M}, \quad (46)$$

where $Q_{micro,M} = Q_{micro} - Q_{micro,CL}$ (the last term is defined by Eq. (43)) is a part of $Q_{micro}$ inside the cell M, N is the ending interfacial point in the cell M, and $A_M$ is the area of its bottom face.

3. Adjustment of evaporation rate in micro- and meso-regions: The heat fluxes in the micro- and meso-regions are also used to evaluate the evaporation rate in the respective cells. For instance, the mass transfer rate $\dot{m}_M$ for the cell M is determined based on the combined heat flux contributions from the micro- and meso-regions:

$$\dot{m}_M = \frac{q_M A_M}{\Gamma_M L}, \quad (47)$$

where $\Gamma_M$ represents the area of the liquid-vapour interface in the cell M.

It is essential to note that the microregion model is driven by the wall superheating imposed in the present problem setup. Therefore, for each simulation, $Q_{micro}$ and $\theta_{app}$ are fixed values shown in Table 2 below. They are derived from the microregion model and vary with the wall superheating $\Delta T$ and microscopic angle $\theta_{micro}$ (Section 3). These values serve as inputs for macroregion calculations.

Although Fig. 13 shows the meso-region spanning only the first two rows of grid cells, it can extend over multiple rows to provide a smooth transition between the micro- and macro- regions. In this study $n_{meso} = 4$ is found to be optimal, cf. Appendix C.

## 5. Simulated system and boundary conditions

A single bubble is studied during its growth over a heated wall using a 2D axisymmetric simulation. Free outflow boundary conditions are applied at $r = r_{max} = 2$ mm and $z = z_{max} = 4$ mm (Fig. 14). The simulation initiates with a spherical cap-shaped bubble seed, with radius $R_0$ and contact angle $\theta_{app}$, positioned over a wall with the imposed temperature $T_{iv} = T_{sat} + \Delta T$.

To incorporate natural convection effects, a linear thermal boundary layer along $z$ (cf. Fig. 14) is initialised above the wall. The initial thickness $\delta$ of this layer, representing the stored thermal energy in the superheated liquid, is estimated via the Kays and Crawford correlation [50]:

$$\delta = 7.14 \frac{v_l \alpha_l}{g \beta_l \Delta T}.$$
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Schematics of the computational domain (not to scale) and of the initial and boundary conditions for 2D axisymmetric simulation.</figcaption></figure>

<table><caption>**Table 2**
Various simulated cases and the corresponding microregion parameters.</caption><thead><tr><th>Case</th><th>δ<br>[mm]</th><th>Δr<br>[μm]</th><th>ΔT<br>[K]</th><th>θ<sub>app</sub><br>[°]</th><th>θ<sub>micro</sub><br>[°]</th><th>s<sub>M</sub><br>[μm]</th><th>h<sub>M</sub><br>[μm]</th><th>Q<sub>micro</sub><br>[W/m]</th><th>n<sub>meso</sub><br>[-]</th></tr></thead><tbody><tr><td>S0</td><td>0.733</td><td>10</td><td>6.2</td><td>38</td><td>38</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>S1a</td><td></td><td>5</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>S1b</td><td></td><td>10</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>S1c</td><td>0.733</td><td>10</td><td>6.2</td><td>38</td><td>32.8</td><td>5</td><td>3.8</td><td>27.7</td><td>4</td></tr><tr><td>S1d</td><td></td><td>20</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>40</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>S1e</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>2</td></tr><tr><td>S1f</td><td>0.733</td><td>10</td><td>6.2</td><td>38</td><td>32.8</td><td>5</td><td>3.8</td><td>27.7</td><td>6</td></tr><tr><td>S1g</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>8</td></tr><tr><td>S2a</td><td>0.6588</td><td>10</td><td>8.5</td><td>50</td><td>44.0</td><td>5</td><td>3.8</td><td>30.5</td><td>4</td></tr><tr><td>S2b</td><td></td><td></td><td></td><td></td><td>50</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>S3</td><td>0.9323</td><td>10</td><td>3</td><td>19.63</td><td>15</td><td>5</td><td>1.6</td><td>22.5</td><td>4</td></tr><tr><td>S4</td><td>0.7399</td><td>10</td><td>6</td><td>23.16</td><td>15</td><td>5</td><td>1.9</td><td>40.3</td><td>4</td></tr><tr><td>S5</td><td>0.6464</td><td>10</td><td>9</td><td>26.21</td><td>15</td><td>5</td><td>2.1</td><td>55.3</td><td>4</td></tr><tr><td>S6a</td><td>0.7029</td><td>10</td><td>7</td><td>50</td><td>45.2</td><td>5</td><td>3.785</td><td>25.032</td><td>4</td></tr><tr><td>S6b</td><td></td><td></td><td></td><td></td><td>50</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

The rest of the domain is set to the saturation temperature, and the vapour temperature remains constant at $T_{sat}$ during the simulation, i.e. the energy Eq. (3) is not solved in the vapour (anyway, the heat fluxes in the vapour are negligibly small). Such boundary and initial conditions are commonly applied in nucleate boiling simulations [10, 19].

The vapour bubble grows over time due to heat transfer from the superheated liquid in the thermal layer and the evaporative heat flux from the microregion $Q_{micro}$. Both $Q_{micro}$ and the apparent contact angle $\theta_{app}$ are derived from the microregion model and vary with the wall superheating $\Delta T$ and microscopic angle $\theta_{micro}$ (Section 3). These values serve as inputs for macroregion calculations.

Simulations are conducted with pure water at atmospheric pressure, using properties listed in Table 1. A regular square grid is used, i.e., $\Delta r = \Delta z$. The combinations of input parameters used for the simulations are summarised in Table 2.

## 6. Results and discussion

The present study is focused on the simulation of the nucleate boiling with a single bubble, on the influence of microregion on bubble dynamics and on the effect of superheating. However, the validation checks should be performed first. In the comparisons with experiment, the values of $\theta_{app}$ and $\Delta T$ are fixed since these parameters are measurable. The value of $\theta_{micro}$ is determined from these parameters with the microregion model. When the superheating effect is studied, $\theta_{micro}$ is

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Variation of $R_{eq}$ with time using various grid spacings computed for $\Delta T = 6.2$ K and $\theta_{app} = 38^{\circ}$ (cases S1a–S1d in Table 2).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Bubble shapes just before detachment using various grid spacings (for the cases S1a–S1d of Table 2).</figcaption></figure>

fixed as it is imposed by the wettability of the wall; $\theta_{app}$ is determined for each $\Delta T$.

### 6.1. Grid independence and validation

The grid independence study was conducted with grids of different spacings (cases S1a–S1d in Table 2). The conditions of Son et al. [10] were adopted, namely the apparent contact angle $\theta_{app} = 38^{\circ}$ and the superheating $\Delta T = 6.2$ K. The equivalent radius of bubble, $R_{eq}$, defined as a radius of a sphere with the same volume as the bubble, was calculated during the bubble growth phase. The results obtained for different spacings are presented in Fig. 15, and the final bubble shapes prior to the detachment are shown in Fig. 16.

All simulations predict the same departure size and shape regardless of the grid spacing, confirming that the contact angle is accurately modelled and demonstrating that the code correctly captures the buoyancy effect. However, the detachment time varies: bubbles simulated with denser meshes detach earlier. This occurs because of an initial temperature discontinuity across the bubble interface, which requires denser meshes (and smaller time steps) to resolve the transient behaviour. As grid spacing decreases, the difference in departure time between different grids becomes smaller. Between the grid spacings 5 µm and 10 µm the difference in bubble growth time is less than 6%. One can conclude that the choice 10 µm provides a trade-off between a sufficient accuracy and extensive computational resources. The simulations were conducted on the Rome partition of the TGCC supercomputer, which is equipped with the dual 64-core AMD Rome CPUs operating at 2.6 GHz and featuring AVX2 instruction set. The computational time and resources required for each case are summarised in Table 3. A dramatic CPU time increase is observed with mesh refinement. This
<table><caption>**Table 3**
Summary of physical bubble growth (departure) time $t_{\text{dep}}$, time step $\Delta t$, CPU number $n_{\text{CPU}}$, computational time $t_{\text{exe}}$, and CPU time $t_{\text{CPU}}$ used for each case.</caption><thead><tr><th>Case</th><th>t<sub>dep</sub>, ms</th><th>Δt, μs</th><th>n<sub>CPU</sub></th><th>t<sub>exe</sub>, hours</th><th>t<sub>CPU</sub>, hours</th></tr></thead><tbody><tr><td>S1a</td><td>31</td><td>0.05</td><td>78</td><td>17.54</td><td>1368</td></tr><tr><td>S1b</td><td>33</td><td>0.075</td><td>18</td><td>9.42</td><td>170</td></tr><tr><td>S1c</td><td>37</td><td>0.3</td><td>6</td><td>1.40</td><td>8.41</td></tr><tr><td>S1d</td><td>42</td><td>0.5</td><td>2</td><td>0.66</td><td>1.32</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Variation of $R_{\text{eq}}$ with time for 10 µm grid spacing and comparison with the experimental [51] and numerical [10] results corresponding to the case S1b.</figcaption></figure>

increase in computational demand is attributed not only to a higher number of grid points but, more significantly, to a substantially smaller time step $\Delta t$ required for simulations with finer meshes.

A comparison with the numerical results [10] and the experimental data [51] is presented in Fig. 17. Two sets of black markers represent the observed variation range in the experiment for different bubbles showing the upper and lower bounds. The growth rate obtained from the present simulation agrees with both the experimental and numerical results.

Moreover, the current study successfully reproduces the stages of the CL receding and advancing without numerical artifacts, such as thermal layer distortion near CL caused by parasite flows typically present in the bubble growth simulations of other groups. The fluid superheating and interface evolution at different stages of the boiling cycle are shown in Fig. 18. Initially, rapid bubble growth occurs due to intense evaporation at both the microregion and the liquid-vapour interface. While the bubble remains small, evaporation dominates buoyancy, causing the CL receding (to the liquid side). As the bubble grows, buoyancy becomes significant and eventually prevails, creating a vortex that reverses the flow (Fig. 18) and causing the CL advancing that occurs from $t \approx 17$ ms. CL keeps advancing until it reaches the symmetry axis, at which point the bubble detaches.

A simulation is also performed with $\theta_{\text{app}} = 50^\circ$ and superheating $\Delta T = 8.5$ K (case S2a of Table 2). Our simulations are compared with the experimental data of Siegel and Keshock [51] and numerical results of Son et al. [10] obtained for the same values of $\theta_{\text{app}}$ and $\Delta T$. The bubble profile prior to departure is compared in Fig. 19(a). A comparison of bubble growth rate, based on equivalent radius, is provided in Fig. 19(b). Our simulation results show excellent agreement with those from [10]. While a slight difference in bubble cap shape is observed, the neck region aligns closely with the experimental profile. The equivalent radius evolution agrees too. The numerical approach can thus be considered as validated.

## 6.2. Verification of multiscale coupling

To determine the influence of microregion on the bubble dynamics, the results with the microregion model labelled MR (case S1b of Table 2) are compared with those without the microregion model (case S0 of Table 2) labelled NMR. The latter is obtained with $Q_{\text{micro}} = 0$ and the heat flux in the meso-region calculated by TrioCFD. The same apparent contact angle is used for both cases as the macroregion slope at the CL. It is important to note that, in reality, the microregion model describes the slope variation from $\theta_{\text{micro}}$ to $\theta_{\text{app}}$. Consequently, to achieve the same apparent contact angle in both cases, the microscopic angles differ.

Fig. 20 compares the bubble growth and CL dynamics under different conditions. The simulation without microregion modelling slightly underestimates the bubble departure size, as shown in the zoomed view in Fig. 21(a). Given that the contact angle and wall superheating (and thus buoyancy force) remain constant, the bubble departure diameter is very similar between cases. However, a significant difference appears in the bubble growth rates: the simulation without microregion modelling predicts a slower growth rate, demonstrated in Fig. 21(b) for bubble shapes at the same instant. This reduced growth rate results in the case without microregion modelling failing to capture experimental bubble dynamics accurately, underscoring the essential role of microregion modelling in simulating nucleate boiling phenomena. Such a behaviour is consistently observed for different superheating values (see Appendix D for additional results).

At $t = 24$ ms, the CL and interface profiles near the wall are aligned for both cases (see Fig. 21(b)). As the bubble interfaces coincide, one can compare the $\dot{m}$ distributions along the interface contour. The phase change rate along the interface at $t = 24$ ms for both cases is presented in Fig. 22. In the micro- and meso- regions, the phase change rate in the MR case matches the analytical solution (40). In contrast, the NMR solution deviates from the analytical. This means that the phase change rate is significantly underestimated by the NMR solution leading to an underestimation of the bubble growth rate. In the macroregion, both MR and NMR solutions coincide. Such a graph clearly shows the cross-over from the microregion to the macroregion behaviour. Note that a smooth transition of mass change rate is obtained thanks to the meso-region introduction.

## 6.3. Multi-cycle behaviour

By artificially introducing a newly nucleated seed after the bubble detachment from its nucleation site, multi-cycle simulations become feasible. Using such an approach, Huber et al. [19] observed that the first two bubble departures are transient, with both the bubble departure diameter and frequency stabilising by the fourth cycle. We have obtained similar quasi-repetitive cycles from the second bubble departure using the same configuration. Besides, we find that the bubble departure diameter remains consistent across cycles, and $t_{\text{dep}}$ becomes smaller in repetitive cycles than for the first bubble.

In simulations without microregion modelling (Fig. 23(a)), the simulated bubble departure time is always longer than the experimental one, indicating the importance of microregion. Conversely, when the microregion model is included (Fig. 23(b)), the bubble departure time for repetitive cycles becomes shorter than the experimental value while the bubble volume is higher. These discrepancies are expected, as the current study does not account for thermal coupling with the solid. The latter should lead to a wall temperature reduction near the CL [52–54]. In addition, the current multi-cycle configuration (imposed by the constant wall overheating assumption) where the bubble starts nucleating immediately after the bubble departure does not model the waiting period that in reality occurs. During this period, the liquid above the nucleation site becomes cooler (due to the heat exchange with the surrounding) than in our simulation. These two factors lead to an overestimation of the bubble growth rate and thus to an underestimation of $t_{\text{dep}}$ that are seen in Fig. 23. Given these deficiencies of the simulation without the wall thermal coupling and accounting for a large simulation time, our subsequent analysis focuses on the dynamics of the first bubble to emphasise the effect of microregion.
<figure><img src="image_20.png" /><figcaption>Fig. 18. Dynamics of bubble interface, fluid superheating (colour) and streamlines computed for the case S1b of Table 2.</figcaption></figure>

## 6.4. Contribution of microregion to bubble growth

To estimate the microregion contribution, the instantaneous microregion evaporation rate $\dot{m}_{micro}$ and the instantaneous total evaporation rate $\dot{m}_{tot}$ inside the bubble can be defined as

$$ \dot{m}_{micro} = \frac{2\pi r_{cl} Q_{micro}}{L}, \qquad (48) $$

$$ \dot{m}_{tot} = \rho_{v} \frac{dV_{v}}{dt}. \qquad (49) $$

The time variations of both evaporation rates during the bubble growth and of their ratio are shown in Fig. 24. These two phase change rates are of the same order of magnitude, underscoring the critical role of microregion modelling in accurately capturing bubble growth dynamics.

The overall relative contribution of microregion to the bubble growth can be calculated as

$$ \eta = \frac{\int_{0}^{t_{dep}} \dot{m}_{micro} dt}{\int_{0}^{t_{dep}} \dot{m}_{tot} dt} \qquad (50) $$

where $t_{dep}$ is the bubble departure time. For the case S1b of Table 2, $\eta \approx 26\%$, which highlights the importance of microregion modelling in nucleate boiling.

## 6.5. Effect of wall superheating

The investigation continues to analyse the influence of wall superheating on the bubble growth characteristics. It is to be noted that the comparison is not performed at a constant $\theta_{app}$. The apparent contact angle increases with the superheating $\Delta T$ while the microscopic contact angle $\theta_{micro}$ remains constant. The study is performed for $\theta_{micro} = 15^{\circ}$. Three superheating values are considered: 3 K ($Ja \sim 9$), case S3 of Table 2), 6 K ($Ja \sim 18$, case S4) and 9 K ($Ja \sim 27$, case S5). The microregion model results in increasingly large values of $\theta_{app}$ and $Q_{micro}$ (cf. Fig. 12) that are displayed in Table 2.

The bubble shapes just before the departure, the zoomed in views of the near-wall region at the maximum dry spot size, and the variations of $R_{eq}$ and $r_{cl}$ with time are shown in Figs. 25 to illustrate the effect of superheating. From the results it can be observed that the bubble departure volume increases with the superheating. This can be explained by a higher apparent contact angle at higher superheating. The bubble departure occurs at a radius close to the Fritz radius, $R_F$, which represents the maximum size an attached bubble can achieve before detaching, as determined by the bubble static force balance. It can be shown that this balance involves the apparent, and not the microscopic contact angle [55], The Fritz radius is thus defined as [56]:

$$ R_F = 0.0104\theta_{app} \sqrt{\frac{\sigma}{g(\rho_l - \rho_v)}} \qquad (51) $$
<figure><img src="image_21.png" /></figure>

(a) Bubble profile just before departure superposed on the experimental bubble image. Blue line: our simulation; black line: simulation of Son et al. [10]. The axes are scaled with $l_c = \sqrt{\sigma/g(\rho_l - \rho_v)} = 2.5$ mm.

<figure><img src="image_22.png" /><figcaption>**Fig. 19.** Comparison with the numerical results of Son et al. [10] and the experimental data from Siegel and Keshock [51] for the case S2a of Table 2.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 20.** Effect of microregion: variation of equivalent radius $R_{eq}$ and dry spot extension $r_{cl}$ with time. The microregion model is activated for curves labelled with MR (case S1b of Table 2), but not for NMR (case S0).</figcaption></figure>

**Table 4**

Summary of departure time $t_{dep}$, departure radius $R_{eq,max}$, Fritz radius $R_F$ and maximum dry spot size $r_{cl,max}$ under different superheating conditions.

<table><thead><tr><th>Case</th><th>ΔT, K</th><th>t<sub>dep</sub>, ms</th><th>R<sub>eq,max</sub>, mm</th><th>R<sub>F</sub>, mm</th><th>r<sub>cl,max</sub>, mm</th></tr></thead><tbody><tr><td>S3</td><td>3</td><td>15.4</td><td>0.635</td><td>0.688</td><td>0.13</td></tr><tr><td>S4</td><td>6</td><td>12.8</td><td>0.802</td><td>0.812</td><td>0.21</td></tr><tr><td>S5</td><td>9</td><td>12.6</td><td>0.932</td><td>0.919</td><td>0.29</td></tr></tbody></table>

It is shown in Table 4 together with some other characteristic parameters relevant to the bubble departure.

The departure time decreases, while the maximum dry spot size increases with superheating. This can be explained by the increase of $\theta_{app}$ too.

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

<figure><img src="image_25.png" /><figcaption>**Fig. 21.** Effect of microregion: bubble profiles in MR and NMR cases (S1b and S0 of Table 2, respectively). (a) Bubble profiles just before departure ($t_{dep}^{MR}$ = 33.2 ms and $t_{dep}^{NMR}$ = 47.4 ms). (b) Bubble profiles at the same instant $t = 24$ ms.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 22.** $\dot{m}$ variation along the interface close to the CL for the enabled (MR) and disabled (NMR) microregion modelling (corresponding to the cases S1b and S0 of Table 2, respectively). The plots correspond to $t = 24$ ms.</figcaption></figure>

The effect of superheating on the bubble shape evolution is shown in Figs. 26. The bubble is more spherical at departure in the case of the low superheating (e.g. at 3 K). As the superheating increases, the bubble stretches vertically and becomes narrower close to the neck.

Fig. 27 illustrates the contribution (50) of microregion to the bubble growth for different $\Delta T$. The microregion contribution is significant in these simulations, accounting for over 15% of total bubble growth and almost linearly increases with superheating.
<figure><img src="image_27.png" /><figcaption>(a) $\Delta T = 7$ K: NMR case (S6b of Table 2)</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(b) $\Delta T = 7$ K: MR case (S6a of Table 2)</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Fig. 24.** Instantaneous phase change rate due to the microregion, $\dot{m}_{micro}$, and the instantaneous phase change rate of the overall bubble, $\dot{m}_{tot}$, for the case S1b of Table 2.</figcaption></figure>

## 7. Conclusions

In the present work, a multiscale simulation approach where a microregion model is coupled with the existing phase change solver of TRUST/TrioCFD is developed and validated to simulate the bubble growth during saturated nucleate boiling of water at atmospheric pressure for the wall superheating values up to 10 K and high contact angles. To our best knowledge, it is the first paper where the scales are rigorously separated, i.e. the microregion heat flux is applied in the microregion and not smeared over the macroregion. Both the flow and temperature distributions do not show parasite rolls. An intermediate meso-region is introduced to provide a smooth interfacial heat flux variation between the micro- and macro- regions. The microregion model of Janeček and Nikolayev was compared to that of Mathieu used by several authors for boiling simulations. The microregion models are carefully analysed. While being quite similar, the Mathieu model

<figure><img src="image_30.png" /><figcaption>(a) Bubble shape just before departure.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(b) Zoom in of the near-wall region corresponding to the maximum dry spot size (for 8.6 ms at $\Delta T = 3$ K, for 5.6 ms at $\Delta T = 6$ K, for 5.2 ms at $\Delta T = 9$ K).</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(c) Variation of equivalent radius $R_{eq}$ and dry spot radius $r_{el}$ over time.</figcaption></figure>

<caption>**Fig. 25.** Bubble behaviour at different superheating conditions (cases S3, S4 and S5 of Table 2).</caption>

uses an imprecise boundary condition that can cause an important numerical error in the evaluation of the apparent contact angle.

By comparing simulations with and without the microregion model, it is found that the microregion evaporation can contribute up to 25% of total bubble growth, which shows the microregion importance even for moderately high contact angles. For a given apparent contact angle, neglecting microregion causes an underestimation of both the heat flux and the evaporation rate at the contact line, leading to a smaller dry spot size, a slower bubble growth, a longer departure time and, ultimately, a failure to capture bubble dynamics observed in experiments. Therefore, the microregion account is essential for accurate simulations of the bubble growth in boiling and for the correct heat transfer evaluation.

A parametric investigation on the wall superheating is performed. It shows that the contribution of microregion evaporation becomes increasingly significant as the superheating rises. The dry spot size and the bubble departure volume both increase while the departure time decreases with superheating, which is explained mostly by the increase of the apparent contact angle.

The microregion effect is evaluated here for the case of isothermal heater. It is well known however that the wall temperature distribution reflecting the thermal response of the heating substrate can have a strong effect on the bubble growth. Although this influence is limited for highly conductive substrates, accounting for the thermal coupling is
<figure><img src="image_33.png" /><figcaption>Fig. 26. Bubble shape evolution for different superheating values corresponding to the cases S3–S5 of Table 2. In each case, the largest time is close to departure.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Fig. 27. Relative contribution of the microregion to the total vapour volume until bubble detachment as a function of the wall superheating. Line is a linear fit.</figcaption></figure>

essential to explore more general heater cases and, more importantly, to accurately simulate multi-cycle bubble growth. The thermal response of the heater will be considered in a forthcoming paper.

## Declaration of competing interest

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

## Acknowledgements

This work was provided with computer and storage resources by GENCI at TGCC thanks to the grant A0162b7712 on the supercomputer Joliot Curie's ROME partition.

# Appendix A. Analytical solution to Stefan problem

The following expressions [57] define the interface position $x^l(t)$, velocity $u^l(t)$

$$x^l(t) = 2\omega\sqrt{\alpha_v t}, \qquad (A.1)$$

$$u^l(t) = \omega \sqrt{\frac{\alpha_v}{t}}, \qquad (A.2)$$

and vapour temperature profile $T_v(x, t)$ for $x \in [0, x^l(t)]$,

$$T_v(x,t) = T_w - \frac{\Delta T}{\operatorname{erf}(\omega)} \operatorname{erf}\left(\frac{x}{2\sqrt{\alpha_v t}}\right), \qquad (A.3)$$

where $\operatorname{erf}()$ is the error function and $\omega$ is a constant given by a solution of the following transcendental equation,

$$\omega \exp(\omega^2) \operatorname{erf}(\omega) = \frac{\operatorname{Ste}}{\sqrt{\pi}}, \qquad (A.4)$$

where the Stefan number is defined as $\operatorname{Ste} = C_{pv} \Delta T / \mathcal{L}$.

# Appendix B. Analytical solutions to Scriven problem

Scriven [42] expressed the bubble radius as

$$R = 2\varphi\sqrt{\alpha_v t}, \qquad (B.1)$$

where $\varphi$ is a constant that should be determined as a part of the problem. Evidently, such a hypothesis is reasonable and is compliant with the boundary condition (12). Then Eq. (11) can be reduced to an ODE

$$T_l'' = 2T_l'(\epsilon\varphi^3\psi^{-2} - \psi - \psi^{-1}), \qquad (B.2)$$

where $T_l$ depends on

$$\psi = \frac{\varphi}{2\sqrt{\alpha_v t}}. \qquad (B.3)$$

Note that such a choice is quite reasonable as it describes a boundary layer solution. By expressing Eq. (B.2) as a first order ODE with respect to the function $T_l'$ and separating the variables, it can be easily integrated, which results in

$$T_l' = A\psi^{-2} \exp(-\psi^2 - 2\epsilon\varphi^3\psi^{-1}). \qquad (B.4)$$

To determine the integration constant $A$, one needs an extra boundary condition

$$\lambda_l \left. \frac{\partial T_l}{\partial \varphi} \right|_{\varphi=R} = \rho_v \dot{R} \mathcal{L}, \qquad (B.5)$$

similar to Eq. (6). Scriven used a slightly different equation

$$\lambda_l \left. \frac{\partial T_l}{\partial \varphi} \right|_{\varphi=R} = \rho_v \dot{R} [\mathcal{L} + (C_{pl} - C_{pv})\Delta T_{\infty}].$$

While the second term in the square brackets is not justified, it usually remains to be much smaller than $\mathcal{L}$ and thus hardly affects the solution. The resulting value is

$$A = \frac{2\varphi^3 \Delta T_{\infty}}{Ja} \exp[\varphi^2(1 + 2\epsilon)],$$

where the modified Jacob number is given by Eq. (15). $T_l$ can be expressed in an integral form by using the boundary condition (14):

$$T_l = T_{\infty} - A \int_{\psi}^{\infty} \psi^{-2} \exp(-\psi^2 - 2\epsilon\varphi^3\psi^{-1}) d\psi. \qquad (B.6)$$

The substitution $\tau = 1 - \varphi/\psi$ leads to a form more convenient for $T_l$ numerical computation:

$$T_l(\varphi/R) = T_{\infty} - \frac{2\varphi^2}{Ja} \Delta T_{\infty} \int_{R/\varphi}^{1} \exp\{-\varphi^2 [(1-\tau)^{-2} - 2\epsilon\tau - 1]\} d\tau, \qquad (B.7)$$
<figure><img src="image_35.png" /><figcaption>Fig. C.28. Comparison of bubble shapes prior to departure for different $n_{meso}$ values. (Cases S1b, S1e-S1g of Table 2).</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Fig. C.29. Variation of $R_{eq}$ and $r_{cl}$ with time for different $n_{meso}$ values. (S1b, S1e-S1g of Table 2).</figcaption></figure>

where the equality $\varphi/\psi = R/\rho$ coming from Eqs. (B.1), (B.3)) was used. The quantity $\varphi$ can be obtained by solving numerically the equation

$$Ja = 2\varphi^2 \int_0^1 \exp \{-\varphi^2 [(1-\tau)^{-2} - 2\epsilon\tau - 1]\} d\tau. \quad (B.8)$$

obtained from Eq. (B.7) by applying the boundary condition (13).

## Appendix C. Validation of meso-region model

As discussed in Section 6.4, a significant heat flux is applied in the microregion, typically covering one or two cells near the CL. To ensure a smooth transition to the macroregion, a meso-region spanning several ($n_{meso}$) cells following microregion is introduced, where the interfacial heat flux is modelled with Eq. (40). To make the $n_{meso}$ choice, we include an analysis varying $n_{meso}$ from 2 to 8 (Cases S1b, S1e-S1g of Table 2).

The number of meso-cells has a minimal impact on bubble dynamics, as indicated by same bubble shapes prior to departure (cf. Fig. C.28), nearly identical bubble growth rates, and CL dynamics (cf. Fig. C.29). The only notable difference is a slightly longer departure time for large $n_{meso} = 6, 8$ (cases S1f, S1g). To explain this, we consider below the heat flux variation along the interface.

Figs. C.30 illustrate the heat flux spatial variation along the interface near the CL for different $n_{meso}$ values, comparing with both the analytical solution (40) and the NMR case (S0 of Table 2) as it was done in Section 6.2. Specific time instants are chosen to ensure the superposition of the two interface contours near the CL (similarly to Fig. 21) so the fluxes become comparable (they become more and more different for larger $s$ values). In all the MR cases, regardless of the $n_{meso}$ value, the predicted heat flux closely matches the analytical

<figure><img src="image_37.png" /><figcaption>(a) Cases S1e, $n_{meso} = 2$, and S0 (NMR)</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(b) Case S1b, $n_{meso} = 4$, and S0 (NMR)</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>(c) Case S1f, $n_{meso} = 6$, and S0 (NMR)</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>(d) Case S1g, $n_{meso} = 8$, and S0 (NMR)</figcaption></figure>

<caption>Fig. C.30. Variation of the interfacial heat flux $q'_i$ along the interface for different $n_{meso}$ for the cases S1b, S1e-S1g compared to the case S0 of Table 2.</caption>

solution within the meso-region, demonstrating a successful model implementation.

However, $n_{meso}$ should be chosen carefully. On one hand, it should not be too large, as the meso-region model assumes negligible thermal advection relative to diffusion in the reference frame of CL, which is valid only in a narrow zone near the CL. Beyond this zone, the
<figure><img src="image_41.png" /><figcaption>**Fig. D.31.** Variation of the equivalent radius, $R_{eq}$, and dry spot extension, $r_{cl}$, over time for different cases considered in the literature, with and without the account of microregion.</figcaption></figure>

analytical solution (44) underestimates the heat flux as seen from the comparison of the analytical and numerical solutions (in particular, NMR solution) in Figs. C.30. This means that if $n_{meso}$ is too large, there is a sharp increase at the boundary between the meso- and macro-regions seen for $n_{meso} = 8$ in Fig. C.30(d). Conversely, for small $s$, TrioCFD algorithms underestimate the heat flux close to the wall as seen from the behaviour of the NMR solution. Too small $n_{meso}$ will thus cause a sharp flux decrease at the boundary between the meso- and macro-regions seen for $n_{meso} = 2$ in Fig. C.30(a). To have a smooth flux variation that produces a good-quality simulation, one needs to choose an intermediate value. A smooth transition between the meso- and macro-regions is achieved for $n_{meso} = 4$. (cf. Fig. C.30(b), which is the choice applied to all other cases considered in this paper (cf. Table 2).

## Appendix D. Additional simulation results showing the microregion contribution

In Fig. D.31 we present comparisons with the experimental [51] and numerical [19] data obtained for two different superheating values but the same apparent contact angle shown in Table 2. We also show the effect of the microregion model by comparing the simulation with (MR) and without its account (NMR). These graphs additionally validate our numerical approach by showing an agreement with the data of other authors. On the other hand, they show the importance of the microregion account. Indeed, without it, the simulation data fall far away from the realistic curves.

## Data availability

Data will be made available on request.

## References

[1] V.K. Dhir, Boiling heat transfer, Ann. Rev. Fluid Mech. 30 (1998) 365-401.

[2] F.D. Moore, R.B. Mesler, The measurement of rapid surface temperature fluctuations during nucleate boiling of water, AIChE J. 7 (1961) 620-624.

[3] N. Zuber, Nucleate boiling. The region of isolated bubbles and the similarity with natural convection, Int. J. Heat Mass Transfer 6 (1963) 53-78.

[4] P.C. Wayner, Y.K. Kao, L.V. LaCroix, The interline heat-transfer coefficient of an evaporating wetting film, Int. J. Heat Mass Transfer 19 (1976) 487-492.

[5] P. Stephan, J. Hammer, A new model for nucleate boiling heat transfer, Heat Mass Transf. 30 (1994) 119-125.

[6] J. Kim, Review of nucleate pool boiling bubble heat transfer mechanisms, Int. J. Multiph. Flow. 35 (2009) 1067-1076.

[7] V.V. Yagov, Generic features and puzzles of nucleate boiling, Int. J. Heat Mass Transfer 52 (2009) 5241-5249.

[8] C. Sodtke, J. Kern, N. Schweizer, P. Stephan, High resolution measurements of wall temperature distribution underneath a single vapour bubble under low gravity conditions, Int. J. Heat Mass Transfer 49 (2006) 1100-1106.

[9] Z. Chen, Y. Utaka, On heat transfer and evaporation characteristics in the growth process of a bubble with microlayer structure during nucleate boiling, Int. J. Heat Mass Transfer 81 (2015) 750-759.

[10] G. Son, V.K. Dhir, N. Ramanujapu, Dynamics and heat transfer associated with a single bubble during nucleate boiling on a horizontal surface, J. Heat Transf. 121 (1999) 623-631.

[11] C. Kunkelmann, P. Stephan, CFD simulation of boiling flows using the volume-of-fluid method within OpenFOAM, Numer. Heat Transf. A 56 (2009) 631-646.

[12] V. Janeček, V.S. Nikolayev, Triggering the boiling crisis: a study of the dry spot spreading mechanism, Interfacial Phenom. Heat Transf. 2 (2014) 363-383.

[13] V. Janeček, V.S. Nikolayev, Contact line singularity at partial wetting during evaporation driven by substrate heating, Europhys. Lett. 100 (2012) 14003.

[14] V. Janeček, V.S. Nikolayev, Apparent-contact-angle model at partial wetting and evaporation: impact of surface forces, Phys. Rev. E 87 (2013) 012404.

[15] V. Janeček, B. Andreotti, D. Pražák, T. Bárta, V.S. Nikolayev, Moving contact line of a volatile fluid, Phys. Rev. E 88 (2013) 060404(R).

[16] L. Fourgeaud, E. Ercolani, J. Duplat, P. Gully, V.S. Nikolayev, Evaporation-driven dewetting of a liquid film, Phys. Rev. Fluids 1 (2016) 041901(R).

[17] X. Zhang, V.S. Nikolayev, Dewetting acceleration by evaporation, J. Fluid Mech. 948 (2022) A49.

[18] V.S. Nikolayev, D.A. Beysens, G.-L. Lagier, J. Hegseth, Growth of a dry spot under a vapor bubble at high heat flux and high pressure, Int. J. Heat Mass Transfer 44 (2001) 3499-3511.

[19] G. Huber, S. Tanguy, M. Sagan, C. Colin, Direct numerical simulation of nucleate pool boiling at large microscopic contact angle and moderate Jakob number, Int. J. Heat Mass Transfer 113 (2017) 662-682.

[20] B. Mathieu, Etudes Physiques, Expérimentale Et Numérique Des Mécanismes De Base Intervenant Dans Les Écoulements Diphasiques (Phd. Thesis), Polytech Marseille, Université de Provence, 2003.

[21] A. Urbano, S. Tanguy, G. Huber, C. Colin, Direct numerical simulation of nucleate boiling in micro-layer regime, Int. J. Heat Mass Transfer 123 (2018) 1128 - 1137.

[22] L. Torres, A. Urbano, C. Colin, S. Tanguy, On the coupling between direct numerical simulation of nucleate boiling and a micro-region model at the contact line, J. Comput. Phys. 497 (2024) 112602.

[23] L. Bureš, Y. Sato, Direct numerical simulation of evaporation and condensation with the geometric VOF method and a sharp-interface phase-change model, Int. J. Heat Mass Transfer 173 (2021) 121233.

[24] Y. Sato, B. Ničeno, A depletable micro-layer model for nucleate pool boiling, J. Comput. Phys. 300 (2015) 20-52.

[25] L. Bureš, Y. Sato, Comprehensive simulations of boiling with a resolved microlayer: validation and sensitivity study, J. Fluid Mech. 933 (2022) A54.

[26] CEA/DES TrioCFD, 2023, URL: http://triocfd.cea.fr/.

[27] A. Toutant, B. Mathieu, O. Lebaigne, Volume-conserving mesh smoothing for front-tracking methods, Comput. Fluids 67 (2012) 16 - 25.

[28] G. Tryggvason, B. Bunner, A. Esmaeeli, N. Al-Rawahi, Computations of multiphase flows, Adv. Appl. Mech. 39 (2003) 81-120.

[29] J.-M. Delhaye, Jump conditions and entropy sources in two-phase systems. Local instant formulation, Int. J. Multiph. Flow 1 (1974).

[30] I. Kataoka, Local instant formulation of two-phase flow, Int. J. Multiph. Flow. 12 (1986) 745-758.

[31] G. Bois, Transferts De Masse Et D'énergie Aux Interfaces Liquide/vapeur Avec Changement De Phase: Proposition De Modélisation Aux Grandes Échelles Des Interfaces (Ph.D. thesis), Grenoble, 2011.

[32] M. Grosso, G. Bois, A. Toutant, Thermal boundary layer modelling for heat flux prediction of bubbles at saturation: a priori analysis based on fully-resolved simulations, Int. J. Heat Mass Transf vol. 222 (2024) 124980.

[33] G. Gennari, R. Jefferson-Loveday, S.J. Pickering, A phase-change model for diffusion-driven mass transfer problems in incompressible two-phase flows, Chem. Eng. Sci. 259 (2022) 117791.
[34] E. Cipriano, A.E. Saufi, A. Frassoldati, T. Faravelli, S. Popinet, A. Cuoci, Multicomponent droplet evaporation in a geometric volume-of-fluid framework, J. Comput. Phys. 507 (2024) 112955.

[35] R. Courant, K. Friedrichs, H. Lewy, On the partial difference equations of mathematical physics, IBM J. Res. Dev. 11 (1967) 215-234.

[36] M. Kang, R.P. Fedkiw, X.-D. Liu, A boundary condition capturing method for multiphase incompressible flow, J. Sci. Comput. 15 (2000) 323-360.

[37] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335-354.

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

[39] S.W.J. Welch, J. Wilson, A volume of fluid based method for fluid flows with phase change, J. Comput. Phys. 160 (2000) 662-682.

[40] S. Hardt, F. Wondra, Evaporation model for interfacial flows based on a continuum-field representation of the source terms, J. Comput. Phys. 227 (2008) 5871-5895.

[41] Y. Sato, B. Ničeno, A sharp-interface phase change model for a mass-conservative interface tracking method, J. Comput. Phys. 249 (2013) 127-161.

[42] L.E. Scriven, On the dynamics of phase growth, Chem. Eng. Sci. 10 (1959) 1-13.

[43] S. Tanguy, M. Sagan, B. Lalanne, F. Couderc, C. Colin, Benchmarks and numerical methods for the simulation of boiling flows, J. Comput. Phys. 264 (2014) 1-22.

[44] V.S. Nikolayev, Evaporation effect on the contact angle and contact line dynamics, in: M. Marengo, J. De Coninck (Eds.), The Surface Wettability Effect on Phase Change, Springer, 2022, pp. 133-187.

[45] V.S. Nikolayev, D.A. Beysens, Boiling crisis and non-equilibrium drying transition, Europhys. Lett. 47 (1999) 345-351.

[46] D.M. Anderson, S.H. Davis, Local fluid and heat flow near contact lines, J. Fluid Mech. 268 (1994) 231-265.

[47] V.S. Nikolayev, Dynamics of the triple contact line on a nonisothermal heater at partial wetting, Phys. Fluids 22 (2010) 082105.

[48] D.M. Anderson, P. Cermelli, E. Fried, M.E. Gurtin, G.B. McFadden, General dynamical sharp-interface conditions for phase transformations in viscous heat-conducting fluids, J. Fluid Mech. 581 (2007) 323-370.

[49] X. Zhang, V.S. Nikolayev, Liquid film dynamics with immobile contact line during meniscus oscillation, J. Fluid Mech. 923 (2021) A4.

[50] W.M. Kays, M.E. Crawford, B. Weigand, Convective Heat and Mass Transfer, McGraw-Hill, 2005.

[54] L. Zhang, Z.-D. Li, K. Li, H.-X. Li, J.-F. Zhao, Influence of heater thermal capacity on bubble dynamics and heat transfer in nucleate pool boiling, Appl. Therm. Eng. 88 (2015) 118-126.

[55] V.S. Nikolayev, V. Janeček, Impact of the apparent contact angle on the bubble departure at boiling, Int. J. Heat Mass Transfer 55 (2012) 7352-7354.

[56] W. Fritz, Maximum volume of vapor bubbles, Phys. Zeitschr. 36 (1935) 379-384.

[57] S.C. Gupta, The Classical Stefan Problem: Basic Concepts, Modelling and Analysis with Quasi-Analytical Solutions and Methods, Elsevier, 2017.