# Direct numerical simulation of bubble cavitation at a wall in micro-gravity

Maxence Deferrez, Sébastien Tanguy, Catherine Colin, Annafederica Urbano

► To cite this version:

Maxence Deferrez, Sébastien Tanguy, Catherine Colin, Annafederica Urbano. Direct numerical simulation of bubble cavitation at a wall in micro-gravity. International Journal of Heat and Mass Transfer, 2026, 254, pp.127612. 10.1016/j.ijheatmasstransfer.2025.127612. hal-05322300

HAL Id: hal-05322300

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

Submitted on 20 Oct 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>

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

# Direct numerical simulation of bubble cavitation at a wall in micro-gravity

Maxence Deferrez<sup>a</sup>, Sébastien Tanguy<sup>b</sup>, Catherine Colin<sup>b ID</sup>, Annafederica Urbano<sup>a ID</sup>,*

## ARTICLE INFO

### Keywords:

- Compressible two-phase flows

- Pool cavitation

- Flash boiling

- Cryogens

## ABSTRACT

Cryogenic tanks used in space applications face challenges due to propellant management, among which is phase change induced by a depressurization. Pool cavitation at the wall in pure methane under micro-gravity conditions is investigated, analyzing its impact on vapor creation and wall heat transfer. A compressible two-phase flow solver, including phase change and conjugate heat transfer, is used. After a validation with available experimental data for single bubble cavitation at the wall in micro-gravity, numerical simulations varying the depressurization amplitude and rate are carried out. Three distinct growth phases are put into evidence: an *expansion-controlled* phase, occurring within the depressurization time frame, and a *transition* phase leading to a *diffusion-controlled* growth where the bubble behaves as in a superheated liquid at constant pressure. Results allow to show the impact that compressible effects have on the phenomenon, and to demonstrate the importance of accounting for a nano-region at the contact line, strongly affecting the wall heat transfer.

## 1. Introduction

Cryogenic fluids are substances that are in a vapor state for standard conditions but are stored in a liquid state at extremely low temperatures. They are used in various technologies such as scientific equipment, energy storage systems and in space applications. Specifically for space propulsion, liquid hydrogen, methane and oxygen are used because of the high specific impulse they can provide [1]. Moreover, the concept of a space depot, a storage facility in space that can refuel and resupply spacecrafts, has gained traction in recent years as a potential solution to the challenges of deep space exploration [2]. Both launcher upper stages or in space depots must be able to store and manage cryogenic fluids for long periods. However, the transport and handling of cryogenic fluids pose challenges due to their tendency to undergo phase change [3]. This can have significant consequences, including the loss of fuel (that need to be vented when the pressure increases in the tank), reduced mission duration, and increased risk of accidents during propellant transfer [4]. Phase change can be induced either by a pressure or temperature variation [5]. For instance, boiling occurs if a liquid is brought to a superheated temperature, that is a temperature higher than the saturation temperature at a constant pressure. Nucleate boiling refers to the process in which small bubbles nucleated in specific sites over a superheated wall grow under isobaric conditions. On the other hand, cavitation refers to phase change induced by pressure variation: it can occur when the pressure drops

in a liquid and becomes lower than the saturation pressure [6] (see Fig. 1). If the pressure drop is induced by a flow (like for instance in a convergent nozzle) the phenomena is referred as *hydrodynamic cavitation* [7,8]. On the other hand, if the depressurization is induced in a quiescent environment, like for instance in a tank when decreasing the pressure, the phenomena is referred as *pool cavitation*. Nucleation of the bubbles can occur either in the liquid core, for instance at the location of impurities, or at the wall. A drop of pressure causing the pool cavitation at the wall in a tank can happen due to a variety of factors, including leakage, or other system failures but particularly during the depressurization process required to cool the propellant down before use [9–11]. Despite its importance, the phenomenon of pool cavitation at the wall in microgravity under controlled depressurization has been investigated only in a few experimental works [4]. We can mention in particular SOURCE I and II sounding rocket experiments [5,12] using HFE7000 with and without non-condensable gas and more recently a pure methane experiment carried out in the ZARM drop tower [13].

To develop suitable pressure and control systems for the tank it is necessary to master the phenomena and this starts with the understanding of fundamental mechanisms at the bubble scale. The problem is complicated by the nature of cryogenic fluids, which are thermosensitive [14]. Cavitation of cryogenic fluids cannot be assimilated with an isothermal process and thermal effects have to be taken into account. Existing models in the literature for cavitation typically focus either on

<footnote><sup>a</sup> Fédération ENAC ISAE-SUPAERO ONERA, Toulouse, France</footnote>

<footnote><sup>b</sup> IMFT, Université de Toulouse, CNRS, INP, Toulouse, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: maxence.deferrez@isae-supaero.fr (M. Deferrez), tanguy@imft.fr (S. Tanguy), colin@imft.fr (C. Colin), annafederica.urbano@isae-supaero.fr (A. Urbano).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Pressure-Temperature diagram for CH₄, liquid vapor equilibrium curve, cavitation and boiling.</figcaption></figure>

an isothermal assumption, like the Rayleigh Plesset equation [6], or
on instantaneous depressurization, neglecting the effects of a gradual
decrease in pressure [15–17] or neglect the very early stages of the bub-
ble growth [18–20]. These last models for the bubble growth typically
tend towards the behavior of a bubble in a superheated liquid, that is
with a diameter proportional to the square root of time in agreement
with Scriven or Plesset and Zwick theories [21–25]. The number of
experimental works dealing with phase change in cryogenic fluids in
micro-gravity conditions is quite limited in the literature especially
if focusing on variable pressure environment. Recently, Weber and
Dreyer [13] proposed an extensive experimental work investigating
wall pool cavitation in pure methane in micro-gravity conditions (drop
tower). These studies have provided valuable insights into the dynamics
of bubble growth induced by a depressurization in micro-gravity. They
could provide images showing the evolution of the bubble in time.
Their results show that after the depressurization phase is over, during
the constant pressure phase, the bubble behavior tends towards the one
of a bubble growing in a superheated liquid. However, their results
do not allow to conclude about the behavior in the depressurization
phase. Moreover, no data are available in terms of wall heat transfer
nearby the bubble, whereas such information would be required for a
better understanding of the overall phenomenon. For instance, it is no
possible to differentiate between contributions arising from compress-
ibility effects (expansion, variable density) and phase change effects
over the bubble growth. Moreover, available experimental diagnostics
do not allow to capture some small scale phenomena, in particular for
cryogenic fluids, like the possible existence of a nano-region at the
contact line strongly affecting the bubble behavior. The existence of
a nano-metric region at the contact line where non-equilibrium ther-
modynamic holds has been theoretically proposed by several authors
since the 90s [26,27]. The main idea is that, at the liquid/vapor/solid
interface, a small region exists where, in non-isothermal conditions,
the interface slope varies from a microscopic contact angle (the one
provided by Young theory, which is almost zero for cryogens), further
associated to the nano-region, to an apparent contact angle which is
measurable. In the small region the heat flux from the wall is large and
induces strong phase change mass flow rate that can strongly contribute
to the bubble growth. Simulating the nano-region is not realistically
feasible with a direct numerical simulation, and recently approaches
to include the effects of the nano-region as a subgrid model in a DNS
simulation have been proposed [28]. It has been demonstrated that the

effect of the nano-region is particularly important if the conductivity of
the wall material underneath the bubble is high. Moreover, important
differences between microscopic and apparent contact angles are to be
expected for wettable fluids, like cryogens. The present work focuses
on cavitation in pure methane, which is a cryogenic fluid. As such,
the existence of a nano-region is expected to play a major role on
the outline of the bubble and on the associated wall heat transfer.
Therefore, the present work will enable to observe if the nano-region
plays indeed an important role for wall cavitation.

In summary, to the authors knowledge, there are no available stud-
ies in the literature analyzing in detail the interaction between phase
change, variable density and expansion over the growth of bubble at
the wall in a variable pressure environment, and quantifying the impact
that these bubbles have on wall heat transfer. Moreover, as already
stated, it is not clear if a nano-region at the wall exists and plays an
important role for pool cavitation, and wall heat transfer. The objective
of the present work is to answer the above research questions carrying
out direct numerical simulations of single bubble cavitation at the wall
in micro-gravity conditions.

To address this problem with numerical simulation, it is mandatory
to have a solver able to accurately simulate two-phase flows with phase
change, including compressible effects and an accurate description of
the thermodynamics of both phases and of the saturation conditions
at the interface. Consequently, usual approaches for the simulation of
boiling flows based on incompressible solvers and constant densities, as
the following [29–33], are not suitable to describe the present problem.

Different methodologies for the simulation of compressible two-
phase flows have been proposed during the last twenty years (see
for instance [34] for a review). In particular, semi-implicit projection
approaches have been developed [35,36]. These algorithms are asymp-
totically preserving (that is they are accurate at low Mach number)
and they allow to avoid the acoustic time step constraint (implicit
algorithm) [34]. In particular, these last points are very beneficial to
carry out numerical simulations in configurations where compressible
effects are important, while not related to high speed flows, such as
depressurization cases. These methods have been extended to two-
phase isentropic flows [37,38] and to two-phase flows including viscous
terms, heat conduction [34,39–41], and phase change [42–45]. The
solver used in the present paper is based on a semi-implicit projection
scheme coupled with a level set/ghost fluid method to handle the inter-
face, as detailed in [41], and includes phase change between a liquid
and its vapor following the framework proposed in [45]. The solver has
been integrated in the DIVA code [30,46–50] yet developed for in-
compressible and variable density flows. It has been recently extended,
in [41,45], to compressible two-phase flows driven by heat conduction,
acoustic interface interaction and phase change phenomena induced by
pressure and temperature variations. The compressible solver for two-
phase flows is based on a generalized cubic Equation of State (EoS) to
describe the vapor phase, the liquid phase and saturation conditions
at the interface in a consistent way. For the purpose of this paper,
the solver has been further developed to include contact lines (liquid,
vapor, solid interface) and conjugate heat transfer. These additional
developments will enable to validate the solver against experimental
data to demonstrate its ability in describing bubble growing at the wall
induced by phase change caused either by temperature (boiling) or
pressure (cavitation) variations.

In the following, the governing equations and numerical methods
used to investigate the problem of a bubble cavitation at the wall under
micro-gravity conditions are first introduced. Attention is given to the
description of the nano-region model used at the contact line. Then
the solver is validated against experimental data with the simulation
of an experiment of cavitation in pure methane under micro-gravity
conditions. The importance of the initial thermal field in the liquid
and of the nano-region at the contact line is discussed in terms of
impact on the bubble growth rate and on the apparent contact angle.
The validated solver is then used to study in detail the cavitation
phenomena and to understand how it is affected by the amplitude of the depressurization and its rate in time. The importance of compressibility effects in this configuration is described. Finally, the impact of the bubble on the wall heat transfer is studied.

## 2. Numerical methods

Describing cavitation phenomena requires a numerical solver to simulate compressible two-phase flows with phase change, including contact angles and conjugate heat transfer with a wall. In this section the governing equations describing the flow are introduced followed by the description of the employed numerical method.

### 2.1. Governing equations

The following system of the governing equations for mass, momentum and energy in each phase, written in terms of primitive variables, i.e. the density $ρ$, the velocity vector $\vec{u}$ and the pressure $p$, is considered:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{u}) = 0, \quad (1) $$

$$ \frac{\partial \vec{u}}{\partial t} + \vec{u} \cdot \nabla \vec{u} + \frac{\nabla p}{\rho} = \frac{1}{\rho} \nabla \cdot \tau + \rho \vec{g}, \quad (2) $$

$$ \frac{\partial p}{\partial t} + \vec{u} \cdot \nabla p + \rho c^2 \nabla \cdot \vec{u} = \frac{c^2 \alpha}{c_p} \left[ \tau \otimes \nabla \vec{u} - \nabla \cdot \vec{q} \right], \quad (3) $$

where $c^2$ is the speed of sound, $\alpha$ the local isobaric expansion coefficient, $c_p$ is the specific heat at constant pressure, $\vec{q} = -k\nabla T$ is the local heat flux (with $T$ temperature) and $\tau$ is the shear stress tensor defined as:

$$ \tau = \mu (\nabla \vec{u} + \nabla \vec{u}^T) - \frac{2}{3} \mu \nabla \cdot \vec{u} \mathbf{I}, \quad (4) $$

with $\mu$ dynamic viscosity and $\mathbf{I}$ is the identity matrix. The following set of jump conditions at the interface $\Gamma$, that account for phase change, must be verified in order to maintain the conservation of mass, momentum and energy across the interface [51]:

$$ [p]_\Gamma = \sigma\kappa + 2 \left[ \mu \frac{\partial u_n}{\partial n} \right]_\Gamma - \frac{2}{3} [\mu \nabla \cdot \vec{u}]_\Gamma - \hat{m}^2 \left[ \frac{1}{\rho} \right]_\Gamma, \quad (5) $$

$$ [\vec{u}]_\Gamma = \hat{m} \left[ \frac{1}{\rho} \right]_\Gamma \vec{n}, \quad (6) $$

$$ [-k \nabla T \cdot \vec{n}]_\Gamma = \hat{m} h_{lv}, \quad (7) $$

where $\sigma$ is the surface tension, $h_{lv}$ is the latent heat of vaporization, $\kappa$ denotes the local interface curvature and $u_n$ is the velocity component in the $\vec{n}$ direction, with $\vec{n}$ normal vector pointing towards the liquid phase [52]. The operator $[\cdot]_\Gamma$ indicates the jump across the interface $\Gamma$. It is defined as $[f]_\Gamma = f_v - f_l$ with the subscripts $l$ and $v$ referring respectively to the liquid and vapor phases. The energy jump condition, Eq. (7), assumes that the interface is at the local saturation temperature ($T_\Gamma = T_{sat}$). Furthermore, thermodynamic equilibrium and zero entropy production at the interface are assumed, implying that the temperature is continuous across the interface ($[T]_\Gamma = 0$) [53]. Small temperature variations are expected in this study and therefore constant values for $\mu$, $k$ and $\sigma$ will be considered. The governing equations in the fluid are coupled to the resolution of the energy equation in the solid wall:

$$ \rho_w c_{p,w} \frac{\partial T}{\partial t} = \nabla \cdot (k_w \nabla T), \quad (8) $$

assuming an equality of the fluxes and a continuity of the temperature between the solid and fluid domain at the solid/fluid interface. The solid material properties, density $\rho_w$, conductivity $k_w$ and specific heat $c_{p,w}$ are considered constant.

<table><caption>**Table 1**
Properties for methane CH$_4$ needed for the evaluation of the SRK EoS coefficients.</caption><tr><td>T<sub>c</sub></td><td>190.564 K</td></tr><tr><td>ρ<sub>c</sub></td><td>45.992 bar</td></tr><tr><td>v<sub>c</sub></td><td>9.8629 10<sup>-5</sup> m<sup>3</sup> mol<sup>-1</sup></td></tr><tr><td>ω</td><td>1.142 10<sup>-2</sup></td></tr><tr><td>M</td><td>1.60425 10<sup>-2</sup> kg mol<sup>-1</sup></td></tr></table>

### 2.2. Equation of state

An unified description of both the liquid and the vapor phases, together with saturation conditions at the interface, is achieved using a cubic Equation of State (EoS) that is required to close the system of Eqs. (1)–(6) as done in [41]. Among the different cubic EoS, the Soave Redlich–Kwong (SRK) is considered for the present study because of its accuracy for the description of methane in the thermodynamic conditions of interest for the present temperature range (see Appendix A). The EoS is expressed as:

$$ p = \frac{\rho RT}{1 - \rho b} - \frac{a \rho^2}{(1 + b_1 \rho)(1 + b_2 \rho)}, \quad (9) $$

where $R = R/M$ is the gas constant, $R = 8.31446 \text{ J K}^{-1} \text{ mol}^{-1}$ the universal gas constant and $M$ the molecular weight.

The coefficients $a$, $b$, $b_1$ and $b_2$, that can be found in [45,54], are functions of the critical parameters, which for methane are reported in Table 1. These coefficients need the definition of the acentric factor $\omega$ which is also reported in this table. Latent heat of vaporization $h_{lv}$, specific heat at constant pressure $c_p$, sound speed $c$ and expansion coefficients are locally computed using the EoS as described in [45].

### 2.3. Contact line and nano-region

The nano-region at the contact line can have a significant impact on nucleate boiling on the wall heat transfer, the contact line dynamics and the bubble’s growth. Note that the nano-region is most of the time referred as micro-region in the literature. However, since it is related with phenomena occurring at the nanometer scale and to avoid confusion with the microlayer [31] we will refer to it with the term nano-region. The description of the nano-region can be done with following system of ordinary equations Eq. (10), derived from the work of Mathieu [55], which is itself a modification of the original model proposed by Stephan and Busse [26] modified for partially wetting fluids:

$$ \left\{ \begin{aligned} \frac{\partial x}{\partial s} &= \cos(\theta), \\ \frac{\partial y}{\partial s} &= \sin(\theta), \\ \frac{\partial \theta}{\partial s} &= \frac{\Delta P}{\sigma}, \\ \frac{\partial \Delta P}{\partial s} &= \frac{8Q\mu_l}{\rho_l h_{lv} r^3 (2\theta \cos(2\theta) - \sin(2\theta))}, \\ \frac{\partial Q}{\partial s} &= \frac{T_w - T_{sat} \left(1 + \frac{\Delta P}{\rho_l h_{lv}}\right)}{R_\Gamma + r\theta/k_l} \end{aligned} \right. \quad (10) $$

where $s$ is a curvilinear abscissa along the bubble interface, $s = 0$ being the contact line, whereas $x$ and $y$ are the cartesian coordinates (Fig. 2). The local interface slope is $\theta$ and $r$ is a modified abscissa, including a slip length $l_s = 0.1$ nm, defined by $r = (y + l_s)/\sin(\theta)$. These equations come from the resolution of Stokes equations in the liquid wedge. They allow relating the pressure difference between the liquid and vapor phases $\Delta p$ with the total heat flux $Q$ (in W m$^{-1}$) in the wedge between the contact line and $s$. The local wall heat flux is $q = \frac{\partial Q}{\partial s}$ and depends on the thermal resistance of the interface $R_\Gamma$ defined as:

$$ R_\Gamma = \frac{2 - f}{2f} \frac{T_{sat} \sqrt{2\pi R_\Gamma T_{sat}}}{\rho_v h_{lv}^2}, \quad (11) $$
<figure><img src="image_5.png" /><figcaption>Fig. 2. Schematic of the nano-region for a partially-wetting fluid.</figcaption></figure>

where f is an accommodation coefficient that can vary between 0 and 1.

The following boundary conditions are considered:

$$ \left\{ \begin{array}{l} x|_{s=0} = 0, \\ y|_{s=0} = 0, \\ \theta|_{s=0} = \theta_{\text{mic}}, \\ \Delta p|_{s=0} = \Delta p_0, \\ Q|_{s=0} = 0, \end{array} \right. \qquad (12) $$

where $\theta_{\text{mic}}$ is the microscopic contact angle associated to the nano-region. The $\Delta p$ at the contact line is not known a priori and the system is solved with a shooting technique requiring that far from the contact line the slope $\theta$ does not vary anymore which is translated in the following condition on the pressure: $\Delta p|_{s \to \infty} \to 0$.

Therefore to summarize the system of Eqs. (10) with the initial conditions Eqs. (12) can be solved for a specific fluid and specific thermodynamic conditions once the following parameters are provided: $\theta_{\text{mic}}$, the wall temperature $T_w$ and the accommodation coefficient f. The outputs of the model are the apparent contact angle at the outlet of the nano-region, $\theta_{\text{app}} = \theta_{s=L_{\text{int}}}$ and the total nano-region heat flux $Q_{\text{nr}} = Q_{s=L_{\text{int}}}$ where $L_{\text{int}}$ is the integration length for the nano-region. The corresponding nano-region heat flux (in W.m$^{-2}$) is $q_{\text{nr}} = Q_{\text{nr}}/L_{\text{int}}$.

## 2.4. Numerical schemes

The system of Eqs. (1)–(3) with the jump conditions Eqs. (5)–(7) are solved with an implicit projection method included in the DIVA code, presented in [41,45]. The solver is based on a Level Set method [56] to capture interface motion. The Level Set function $\phi$ is defined as the signed distance from the interface, $\phi = 0$ representing the interface, and its evolution follows the transport equation:

$$ \frac{\partial \phi}{\partial t} + \vec{u}_T \cdot \nabla \phi = 0, \qquad (13) $$

where $\vec{u}_T$ is the interface local velocity deduced from the jump condition Eq. (6). A re-initialization algorithm ensures that the $\phi$ function keeps its distance property.

The Ghost-Fluid method [57–59] allows to handle the jumps across the sharp interface. In particular, second order extrapolations, proposed in [60], enable to build a ghost field for the temperature field T on each side of the interface while maintaining the continuity of its normal derivative and second order normal derivative, as detailed in [30]. Such extrapolations are critical for an accurate computation of the phase change mass flow rate $\dot{m}$. Regarding the pressure and density fields, $\rho$ and $\rho$, constant extrapolations have shown to be sufficient. The phase change mass flow rate $\dot{m}$ is computed according to the energy jump equation Eq. (7) and then used to extrapolate the velocity ghost field in order to impose the jump condition on the velocity Eq. (6) [61]. The pressure velocity correction algorithm presented in [41,45] brings to the resolution of two linear systems for the pressure field that are solved with a Black-Box Multi-grid solver [62]. Convective terms are solved using 5th order Weno Z schemes [63], diffusive terms using second-order finite differences and a second-order Runge–Kutta scheme is used

for the temporal integration. The solver has been proved to be able to describe a large variety of compressible two-phase flows including natural convection, acoustic interface interaction and cavitation in the fluid bulk, induced either by a depressurization in a tank or by a flow in a convergent duct. Its capability to account for phase change phenomena induced by temperature and pressure variations has been demonstrated. Numerical results have been compared with analytical results of typical test cases like the Scriven case [21] describing the growth of a bubble in a superheated liquid.

The present work proposes novel developments of this numerical solver in order to simulate contact line problems. This consists in three major points. The first point is to impose a Neumann boundary condition for the Level-Set function to impose the contact angle $\theta_{\text{app}}$ at the liquid/vapor/solid interface, such as,

$$ \vec{n} \cdot \nabla(\phi)|_{w} = -\cos(\theta_{\text{app}}). \qquad (14) $$

The second point is the coupling of the fluid resolution with the heat transfer in the solid domain, thus solving Eq. (8) with appropriate interface conditions at the solid/fluid interface. A linear system results from the implicit temporal discretization of Eq. (8). This linear system is solved with a Gauss–Seidel method. The time step is chosen to be the same as the fluid solver time step. A heat flux continuity condition is imposed at the boundary between the solid and fluid domains to ensure energy conservation. Similar implementations have already been performed in previous works with the incompressible and variable density solvers in DIVA, as presented in [32,64] for example, and therefore are not reported in more details here. The third point is related to the coupling between a subgrid model for the nano-region and the compressible two-phase flow solver. A recent methodology to include the nano-region as a subgrid model in an incompressible solver has been presented in [28] and implemented in DIVA solver. This coupling has demonstrated its interest to accurately capture the bubble growth and heat transfer for sub-cooled nuclease boiling in micro-gravity conditions [28]. In particular, accurate comparisons between simulations and experiments have highlighted the significant impact of the nano-region on the overall bubble growth rate. In the present work, the nano-region model described previously is implemented in the compressible two-phase flow solver. As it is shown in [28], that heat conduction in the solid substrate has a strong effect on the nano-region heat flux, a coupled resolution of the energy conservation equations in the solid and fluid domains is mandatory to ensure a correct prediction of the nano-region heat flux. The underlying assumption is that the nano-region heat flux corresponds to the heat flux coming from the wall and is computed at the cell crossed by the contact line. The integration length $L_{\text{int}}$ for the nano-region model is taken equal to the cell size. The heat flux $q_{\text{nr}}$ and wall temperature $T_w$ obtained are imposed via a Robin boundary condition for the temperature linear system resolution in the solid which ensures continuity of the heat flux in the cell crossed by the contact line, such as,

$$ \lambda \frac{\partial T}{\partial y}|_{w} = q_{\text{nr}}(T_w) = AT_w + B, \qquad (15) $$

where the nano-region heat flux $q_{\text{nr}}$ has been linearized with respect to $T_w$, A and B being the coefficients. Once it has been computed when solving the thermal field, the nano-region mass flow rate can then be injected into the simulation and spread out on 5 cells along the interface of the bubble for stability reasons. All details about the coupling between the nano-region model and the solid solver are provided in [28].

The correct implementation of these new developments will be validated with a comparison against experimental data presented in this paper.

## 3. Validation against experimental data

In this section, the experiment of the cavitation of a single bubble at the wall in micro-gravity conditions carried out by Weber and Dreyer [13] is introduced in a first time. Next, comparisons between
<figure><img src="image_6.png" /><figcaption>Fig. 3. Experimental pressure profile imposed at the top of the numerical domain.</figcaption></figure>

experimental data and our numerical simulations are presented. Thor-
ough interpretations of the numerical results will enable to estimate the
impact of the nano-region.

3.1. Experimental setup and data

The goal of the experiment led by Weber and Dreyer [13] is to
study the behavior of a single bubble cavitation of pure methane
during a depressurization event under micro-gravity conditions. The
experimental setup involves a cryogenic test facility designed for use
in the compensated gravity conditions of the ZARM drop tower at the
University of Bremen. A simplified view of the test cell is shown in Fig.
4.

During the experiment in micro-gravity conditions, a depressuriza-
tion through a valve is performed in a vessel containing pure liquid
methane in equilibrium with its vapor at saturation conditions: p₀ =
1.404 bar and Tₐt,₀ = Tₐt(p₀) = 116.06 K. The pressure profile in
the experiment is reported in Fig. 3. Before the depressurization, a
bubble is nucleated with a heater in an artificial cavity having a
radius of Rₑ = 266.5 µm. During the depressurization, the saturation
temperature decreases and the fluid's temperature becomes higher than
its saturation temperature: Tᵢ = Tᵢ > Tₐt(p(t)), with p₀ > p(t) > p_f and
p_f = 1.1 bar. As a consequence the liquid becomes superheated and
the bubble grows.

Five experiments are performed. The available experimental data
in [13] are the evolution of the bubble radius *R* in time and the images
of the bubbles. The available images are used in order to extract the
apparent contact angle θ<sub>app</sub> versus time. Details on the algorithm used
to extract the angles are provided in the Appendix B.

Other available experimental data are temperature measurements
in the liquid at the glass cylinder, (i.e. about 40 mm away from the
nucleation site), that remain equal to the initial saturation temperature
T<sub>sat</sub>(p<sub>0</sub>) during the test. For what concerns the region around the
nucleation site, a platinum Resistance Temperature Detector (RTD)
PT1000 is used to generate heat for the nucleation. The RTD surface
is 7 mm<sup>2</sup> and the power varies in the range 5–150 mW (corresponding
to a heat flux in the range 0.7–21 kW m<sup>-2</sup>), during the nucleation phase
before the micro-gravity phase, and is kept constant at 5 mW during the
micro-gravity phase (corresponding to an heat flux of 0.7 kW m<sup>-2</sup>).

3.2. Computational domain and boundary conditions

The computational domain, schematically reported in Fig. 5, is a
cylinder with dimensions of l<sub>x</sub> = l<sub>y</sub> = 8mm filled with liquid methane.
Axisymmetric simulations are carried out, x = 0 being the symmetry
axis. The cavity support, which is a tube in the experiment, is modeled
as a flat wall, with a height of l<sub>w</sub> = 500 µm, and the solid material is

<figure><img src="image_7.png" /><figcaption>Fig. 4. Simplified scheme of the experimental set-up used in [13], as acquired through private correspondence with N. Weber. The computational domain used in this study is reported (around the location of the nucleation site).</figcaption></figure>

stainless steel 304. The diameter of the stainless steel tube is D<sub>tube</sub> = 6 mm. The pressure variation is applied at the top of the domain, following the profile illustrated in Fig. 3, which is directly derived from experimental data. A no-slip and adiabatic condition is applied to the external wall of the cylinder. The initial conditions are set at saturation both in the liquid and in the vapor phase, with p<sub>0</sub> = 1.404 bar and T<sub>i</sub> = T<sub>sat</sub>(p<sub>0</sub>) = T<sub>sat,0</sub> = 116.06 K, as determined by the SRK EoS (see Section 2.2). The initial saturation conditions and the corresponding thermophysical properties (kept constant) are reported in Tables 1–3. Moreover, in order to investigate the possible impact of the heat input from the PT1000 prior to the zero gravity phase (i.e in normal gravity, see Section 3.1) a superheat can be considered in the wall initialized at T<sub>w,0</sub> with ΔT = T<sub>w,0</sub> – T<sub>sat,0</sub>. A thermal boundary layer for the temperature is considered nearby the wall. The thickness δ<sub>tbl</sub> of the thermal boundary layer is computed from Eq. (34):

$$
\delta_{tbl} = \frac{D_{tube}}{Nu_{tbl}}, \tag{16}
$$

where the Nusselt number Nu<sub>tbl</sub> is computed from a natural convection
correlation provided by the work of Kuehn and Goldstein [65] for
a horizontal infinitely long cylinder. Natural convection correlation
under normal gravity conditions has been chosen because most the
injected heat inducing this thermal boundary layer happens before the
drop, so in normal gravity conditions. See Appendix C for details about
the correlation and the values used.

A linear profile is imposed in the boundary layer between T<sub>sat</sub> + ΔT at the wall and T<sub>sat</sub> at y = δ<sub>tbl</sub>. Therefore, for axisymmetrical simulations, the corresponding heat flux from the wall to the liquid can be computed from :

$$
q_{tbl} = \frac{k_f \Delta T}{\delta_{tbl}}. \tag{17}
$$
<figure><img src="image_8.png" /><figcaption>Fig. 5. Computational domain and boundary conditions used in the numerical simulations.</figcaption></figure>

<table><caption>**Table 2**

Thermodynamical and thermophysical properties.</caption><thead><tr><th>Property</th><th>Liquid (l)</th><th>Vapor (v)</th></tr></thead><tbody><tr><td>T<sub>sat,0</sub> [K]</td><td>116.06</td><td>116.06</td></tr><tr><td>p<sub>sat,0</sub> [bar]</td><td>1.404</td><td>1.404</td></tr><tr><td>k [mW m<sup>-1</sup> K<sup>-1</sup>]</td><td>177.57</td><td>12.03</td></tr><tr><td>c<sub>0</sub> [kJ g<sup>-1</sup> K<sup>-1</sup>]</td><td>2.047</td><td>2.283</td></tr><tr><td>μ [μPa s]</td><td>106.93</td><td>4.4820</td></tr><tr><td>σ [mN m<sup>-1</sup>]</td><td>12.0656</td><td></td></tr></tbody></table>

<table><caption>**Table 3**

Stainless Steel 304 tube's properties at T = 116 K.</caption><thead><tr><th>Property</th><th>Wall (w)</th></tr></thead><tbody><tr><td>k [W m<sup>-1</sup> K<sup>-1</sup>]</td><td>10.2</td></tr><tr><td>ρ [kg m<sup>-3</sup>]</td><td>7900</td></tr><tr><td>c<sub>p</sub> [J kg<sup>-1</sup> K<sup>-1</sup>]</td><td>282</td></tr></tbody></table>

The superheat, $\Delta T$, will be varied between 0 (no thermal boundary layer, and wall initialized at liquid saturation conditions) to 1.5 K. This corresponds to a $q_{tbl}$ between 0 and 50 W m$^{-2}$, and a power input $q_{tbl} * S$, with $S = \pi l_r^2$, between 0 and 100 mW. These fluxes are in the range of the experimental values since, as described in Section 3.1, the power input of the PT1000 varies between 5 and 3000 mW, and the heat flux between 0.7–21 kW m$^{-2}$.

The level set function is initialized to obtain an initial vapor bubble of radius $R_0$ at the wall (see Fig. 5). The initial radius of the bubble is defined as the radius of the cavity $R_0 = R_c = 266.5$ µm. The initial apparent contact angle is set to a minimal value of $\theta_{app,0} = 1^\circ$ equal to the micro-scopic contact angle considered $\theta_{mic}$. The apparent contact angle $\theta_{app}$ that will be imposed along the simulation after initialization will be predicted by the nano-region model.

Cartesian and uniform grids are used, with cell sizes ranging from 4 µm to 16 µm depending on the case, corresponding to 2048 × 2048 to 256 × 256 cells in the fluid domain and 2048 × 128 to 256 × 16 in the solid domain. The impact of the grid resolution over the results will be discussed in Section 3.3.

### 3.3. Simulation of the experiment and analysis

Numerical simulations, accounting for conjugate heat transfer and nano-region, have been performed with the solver introduced in Sec-

tion 2, and are presented in this section. In a first time, a base simulation with an accommodation factor set to $f = 1$ (see Eq. (11)) and a $\Delta T = 0$ has been carried out. A grid sensitivity is presented and analyzed in this preliminary configuration. Then, the influence of $f$ and $\Delta T$ will be investigated.

#### 3.3.1. Base simulation and convergence analysis

As explained above, the base simulation is obtained with the accommodation factor $f = 1$ and the superheat in the thermal boundary layer $\Delta T = 0$. This base case is performed for different grid resolutions and the resulting equivalent radius of the bubble versus time is shown in Fig. 7. The results are compared with the experimental data on the same figure. A numerical convergence of the results is clearly demonstrated since the results obtained with the different grids are almost superimposed. For $t < \Delta t$, during the depressurization time frame, there is an expansion phase, where bubble growth is driven both by the phase change and by the expansion. For $t > \Delta t$ there is a transition towards a diffusive regime, where the bubble growth is solely due to phase change. Hence, we define these two different phases as expansion-controlled growth and diffusion-controlled growth which will be further discussed in Section 4.1. A change in the slope of the bubble's growth is clearly observable in the numerical results as well as in the experimental data at $t = \Delta t$. While the simulation accurately captures the transition to the diffusion-controlled growth's slope, it underestimates the expansion-controlled growth's slope, resulting in a 20% shift in the final radius.

#### 3.3.2. Impact of the nano-region

To evaluate the impact of the nano-region on the results, the base simulation ($\Delta T = 0$) was repeated while varying the accommodation coefficient $f$ from 0 to 1 (representing no effect from the nano-region to maximum effect). The impact on the equivalent radius of the bubble versus time can be observed in Fig. 7. The nano-region has a very low impact at the beginning of the expansion no matter the value of the accommodation coefficient. Indeed, in the initial condition there is no temperature difference since $T_{w,0} = T_{sat,0}$, and therefore the nano-region does not exist. During the expansion-controlled growth the pressure decreases (from $p_0$ to $p_f$, see Fig. 3) reducing the saturation temperature while, overall, the solid temperature tends to stick to its initial temperature. Consequently, the temperature difference increases reaching its maximum at the end of the expansion-controlled growth when the wall temperature $T_{sat}(p_0) > T_w > T_{sat}(p_f) = T_f$, and this induces an increasing effect of the nano-region. The phase change mass flow rate associated to the nano-region becomes more important as this temperature difference increases, increasing the total phase change mass flow rate, and consequently the bubble radius increases as well compared to a simulation without nano-region. A difference up to 9% is therefore observed on the final radius ($f = 0$ against $f = 1$ cases). Overall, the impact of $f$ over the bubble growth is very small for $f > 0.25$ as shown by the negligible differences between $f = 0.25$ and $f = 1$.

On the other hand, $f$ has a non negligible impact on the outline of the bubble, especially at the contact line as shown in Fig. 6. This figure shows the temperature fields and the interface location at the end of the depressurization phase ($t = \Delta t$) obtained with the different $f$ reported. For $f = 0$ the apparent contact angle is very close to the microscopic one: $\theta_{app} = 4^\circ$, resulting in an almost spherical bubble. The wall appears to be slightly cooled down only in the tiny region around the contact line. It is worth mentioning that $\theta_{app}$ increases when increasing the value of $f$, with typical values between 18° and 25°. The vapor at the bubble foot appears to be warmer and the cooled region in the wall increases in size as the solid-vapor contact surface increases. The impact of $f$ on the apparent contact angle is represented in Fig. 8 where $\theta_{app}$ is reported versus time, and compared with the experimental angle extracted from the images (see Section 3.1 and Appendix B). At the beginning of the simulation, $\theta_{app} = \theta_{mic}$. As discussed above, the
<figure><img src="image_9.png" /><figcaption>Fig. 6. Comparison of the temperature field for different nano-region accommodation factor with a focus on the apparent contact angle.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Comparison of the growth of the bubble between the experimental data and the numerical simulation with the accommodation factor set to 1 for different mesh sizes along equivalent radius of the bubbles versus time obtained with different simulations varying the accommodation coefficient *f* in the nano-region for a constant mesh size *dx* = 8 µm.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Evolution of the apparent contact angle for different nano-region accommodation factor and comparison to the experimental data.</figcaption></figure>

apparent angle increases over time as a consequence of the increase on the superheat at the contact line. The numerical results exhibit an increase in the contact angle during the *expansion-controlled* growth. After the *expansion-controlled* growth, the pressure in the vessel rises, and so does the saturation temperature *T*<sub>sat</sub>. At the same time, the wall cools down, inducing a local decrease in the solid wall temperature. Overall, this induces a decrease of the wall superheat and consequently a decrease on the nano-region heat flux resulting in a slight decrease of θ<sub>app</sub>.

Based on the comparison between experimental and numerical results, it appears that the measured θ<sub>app</sub> is compatible with the existence of a nano-region. Indeed, without a nano-region, in micro-gravity conditions the apparent contact angle should be equal to θ<sub>mic</sub> and thus very far from the experimental observation (f = 0 in Fig. 8). On the other hand, it is worth noting that the uncertainty in the measurement of θ<sub>app</sub> is very high, and this does not allow to select a specific f based on this comparison, and any value between 0.5 and 1 would be suitable (also considering the negligible impact it has over the diameter, Fig. 7). For the remainder of this paper it has been decided to consider a value of f = 0.75.

### 3.3.3. Impact of the thermal boundary layer

As previously observed, numerical simulations starting in an isothermal liquid underestimate the temporal evolution of the bubble radius. A possible explanation of this discrepancy could be the initial temperature field in the liquid which cannot be perfectly known and controlled in the experiments. In particular, in Fig. 7, two sequences of experimental results are plotted which correspond to two different shots. Deviations are observed between both sequences, which suggest that the experimental results are sensitive to variations on the initial conditions. In particular, such deviations could be attributed to slightly different temperature initial conditions in the liquid. Therefore, we present in this section a sensitivity study on the initial condition of the temperature field, to highlight the relevance of this parameter. For this purpose, simulations have been performed with ΔT = [1.0; 1.5] K respectively giving δ<sub>tbl</sub> = [4.30; 4.09] mm. These simulations are compared with the simulation without thermal boundary layer, i.e. ΔT = 0 K as well as with the experimental data. For all the simulations, the nano-region model is activated with f = 0.75 as discussed in Section 3.3.2. The results are shown in Fig. 9 in terms of equivalent radius versus time. For every simulation, results obtained with two different grid resolutions are shown (with a cell size 8 and 16 µm) demonstrating the numerical independence of the numerical results with respect to the grid. The superheat has an important impact over the results and the experimental data are retrieved with a ΔT between 1 and 1.5 K. A very good agreement is found in the *diffusion-controlled* phase while the slope of the curve continue to be very slightly underestimated in the *expansion-controlled* growth. This underestimation might be due to two different hypothesis that have been made: the thickness of the thermal boundary layer might not exactly follow a natural convection correlation, and, in reality, the distribution might not be linear. A higher gradient at the wall would help for a faster growth when the bubble is small in the *expansion-controlled* growth while still preserving the good agreement in the *expansion-controlled* growth slope.

The evolution in time of θ<sub>app</sub> is shown in Fig. 10. As expected, for ΔT > 0, the nano-region has an impact since the beginning of the simulation, providing a θ<sub>app</sub> > θ<sub>mic</sub>. In the *diffusion-controlled* growth a good agreement between experimental and numerical data is obtained and the impact of ΔT appears to be limited.

Fig. 11 illustrates the temperature field and bubble outline for the simulation with ΔT = 1.5 K at three given times *t* (0 s, 0.14 s, and 0.5 s), corresponding to the initialization, end of depressurization *t* = Δ*t*, and end of simulation, respectively. During the expansion the temperature in the bubble decreases following the evolution of the saturation temperature. After the depressurization, the vapor temperature increases as a consequence of *T*<sub>sat</sub> which slightly increases and of the heat coming from the solid wall through the dry surface. The impact of the bubble on the wall temperature which decreases in time under the bubble is visible.

In conclusion, the comparison obtained between the numerical results reproduce conveniently the experimental data both in terms of bubble growth rate in time: in *expansion-controlled* and *diffusion-controlled* regions, as well as in terms of θ<sub>app</sub>. In the rest of the paper, the solver will be used to investigate in detail the pool cavitation phenomena.
<figure><img src="image_12.png" /><figcaption>Fig. 9. Equivalent radius of the bubble versus time, from simulations with different ΔT and different meshes.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Apparent contact angle evolution versus time. Comparison between simulations obtained with different ΔT and the experimental data with a coarse and a refined meshes.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Temperature field of the simulation with nano-region accommodation factor set to 0.75 and ΔT = 1.5 K.</figcaption></figure>

## 4. Analysis of the growth of the bubble

In Section 3, it was observed that the growth of the bubble is divided in an expansion-controlled and diffusion-controlled growth.

In the diffusion-controlled growth, the bubble growth is induced by vaporization in a superheated liquid. Consequently, it is expected that in this case the growth of the bubble is driven by thermal diffusion only and should follow a slope proportional to $\sqrt{t}$, as demonstrated by analytical theories and experimental observations [15,25]. In Fig. 12, the temporal derivative of the ratio $R/\sqrt{t}$ is reported versus time. It appears indeed that after the end of the expansion, the growth rate tends towards a diffusion-controlled behavior. It appears also that

<figure><img src="image_15.png" /><figcaption>Fig. 12. Example of the evolution of the derivative of the growth divided by the square root of time for $\Delta t = 65$ ms and $\epsilon_p = 20\%$ along its corresponding equivalent radius.</figcaption></figure>

the growth rate is very far from a diffusion-controlled growth in the expansion-controlled phase. After the expansion some time is required before reaching a growth rate driven by diffusion. An intermediate zone is observed, we will call it the transition phase, for $t < \Delta t < t_{diff}$ where bubble growth rate changes to eventually reach the $\sqrt{t}$ behavior when $t = t_{diff}$ the time where the diffusion-controlled growth is reached. These observations motivate the rest of the paper that aims at characterizing the expansion-controlled growth and the transition towards a diffusion-controlled growth. To address this, investigations on the two main parameters of the problem, that is the depressurization amplitude $\Delta p$, and the depressurization rate $\Delta t$, have to be carried out. An ideal case is considered, simplifying the experimental one. First, liquid and solid wall temperatures are initially assumed uniform (no thermal boundary layer). Second, a linear profile of the temporal pressure variations during the expansion-controlled growth is considered and the pressure is kept constant after the end of the depressurization, so as to keep the $T_{sat}$ constant for $t > \Delta t$:

$$p(t) = \begin{cases} p_0 - \frac{\Delta p}{\Delta t} t & \text{for } 0 < t < \Delta t, \\ p_0 - \Delta p & \text{for } \Delta t < t. \end{cases} \qquad (18)$$

Four simulations are carried out with $\epsilon_p = \frac{\Delta p}{p_0} = [5\%; 25\%]$ and $\Delta t = [50; 200]$ ms. The domain dimensions, initial conditions and boundary conditions are the same as the one described in Section 3.2 except that no thermal boundary layer is considered. The cell size is $dx = 8$ µm and the initial radius of the bubble is $R_0 = 80$ µm. This initial radius as been chosen from preliminary simulations in order to ensure that there are no effects of the initial radius on the final diameter.

The results are shown in Figs. 13–19 and Table 4.

### 4.1. Depressurization rate and amplitude

The temperature field and the bubble shape are plotted in Fig. 13 for the four configurations and at two given times, $t = \Delta t$ and $t = 0.5$ s which correspond respectively to the end of the expansion-controlled growth and a longer time after the depressurization. The evolution of the equivalent radius versus time is shown in Fig. 14. Furthermore, eight numbered tags are displayed across the figures in this section to help the comparison between the figures for the different simulations and displayed time.

The first observation is that increasing the $\Delta p$ the bubble grows more. This is expected, since increasing the $\Delta p$, the liquid superheat increases. Indeed, Fig. 13 shows that the temperature in the bubble (which is driven by the $T_{sat}$ at the interface) is lower for higher $\Delta p$. The second observation is that at the end of the expansion-controlled growth, the bubble radius is greater for higher $\Delta t$ cases (see Fig.
<figure><img src="image_16.png" /><figcaption>Fig. 13. Temperature fields and bubble shape from four simulations with different Δp and Δt at two different time instants t: (a) for t = Δt; (b) t = 0.5 s.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Equivalent radius from four simulations with different Δp and Δt against time.</figcaption></figure>

14). The phase change process is driven by the diffusion time scales, which are slower in comparison to the considered depressurization rate. Meaning, a higher depressurization rate Δt allows the bubble more time for phase change to happen. A complementary information can be found evaluating the thermal boundary layer around the bubble δb,tbl at different times (see Appendix D). The δb,tbl is a measure of the intensity of the phase change mass flow rate (see Eq. (7)), which is shown in Fig. 15. At t = Δt, δb,tbl is thicker for slower depressurization rates: the amplitude of the phase change mass flow rate decreases for higher Δt but has more time to act and therefore results in a higher bubble's diameter D. On the other hand, if we observe during the diffusion-controlled phase and compare the bubbles for a same time t = 0.5 s (i.e enough time after the depressurization ends) we observe that the bubble size tends to be slightly bigger for a faster depressurization. Indeed, since the Δt is smaller, for a same t overall the bubble spends more time in the maximum superheated liquid conditions, inducing the maximum of mass flow rate sooner. This explains why the radius is greater. However, as observed in Fig. 14, for t → ∞, the impact of the Δt on the radius become negligible as the relative difference tends to zero. It is understandable since the same quantity of energy is involved, simply not on the same time-scale. Thus, for a time big enough, the final size of the bubble is almost only driven by the parameter Δp.

Note that it is possible to evaluate the diffusion time $t_{diff}$: it is defined as the end of the transition region when the growth of the

<figure><img src="image_18.png" /><figcaption>Fig. 15. Thickness δb,tbl of the bubble's thermal boundary layer for different Δp and Δt against time.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. Reach of the diffusion-controlled growth.</figcaption></figure>

bubble becomes proportional to $\sqrt{t}$ meaning that the diffusion-controlled growth has been reached and that the influence of the initial expansion-controlled growth starts to be negligible. The time $t_{diff}$, defined in this work as the time t when $\frac{\partial R/\sqrt{t}}{\partial t} = 0$ for t > Δt, is shown in Fig. 16 for the four simulations (see Appendix E for the exact evaluation method used). For all the cases the diffusion-controlled regime is reached before the end of the simulation and appears to be only slightly dependent on Δp. On the other hand $t_{diff}$ increases with Δt but the ratio $t_{diff}/\Delta t$ is reduced for higher Δt. Therefore, both Δp and Δt are important parameters in the following section that aims at extracting models for the diameter evolution in time.

## 4.2. Impact of compressible effects on the expansion-controlled growth

In the expansion-controlled growth, the bubble grows for two reasons: the density decrease because of the depressurization and because of phase change is occurring. It is of interest to quantify the relative importance of the density decrease and phase change over the bubble growth in this phase: this will provide an indication on the importance of accounting for compressible effects in the simulations of bubble cavitation at the wall. If the flow is incompressible, the bubble volume variation $\Delta v_{I,PC}$ induced by the depressurization is solely induced by phase change:

$$ \Delta v_{I,PC} = \frac{\Delta m}{\rho_0} = \frac{m_f - m_0}{\rho_0}, \qquad (19) $$

with

$$ \Delta m = m_f - m_0 = \int_{0}^{\Delta t} \dot{m} dt, \qquad (20) $$
where vapor density is constant and equal to the initial saturation conditions $ρ_0$, and $m_0$ and $m_f$ are the vapor masses at the beginning and at the end of the depressurization.

On the other hand, for a compressible flow, the volume variation $Δv_C$ will be induced both by the density variation and by the vapor mass variation:

$$Δv_C = \frac{m_f}{ρ_f} - \frac{m_0}{ρ_0}, \qquad (21)$$

$Δv_C$ can be further decomposed in three distinct parts as proposed in Eq. (22)–(23):

$$\begin{aligned} Δv_C = & \frac{Δm}{ρ_0} \\ & + Δm \left( \frac{1}{ρ_f} - \frac{1}{ρ_0} \right) \\ & + m_0 \left( \frac{1}{ρ_f} - \frac{1}{ρ_0} \right) \end{aligned} \qquad (22)$$

$$Δv_C = Δv_{I,PC} + Δv_{C,PC} + Δv_{C,E} \qquad (23)$$

Defining $Δv_{C,E}$ and $Δv_{C,PC}$ as:

$$Δv_{C,PC} = Δm \left( \frac{1}{ρ_f} - \frac{1}{ρ_0} \right), \qquad (24)$$

$$Δv_{C,E} = m_0 \left( \frac{1}{ρ_f} - \frac{1}{ρ_0} \right). \qquad (25)$$

This shows that the $Δv_C$ can be decomposed in several terms. A variation of volume induced by incompressible phase change $Δv_{I,PC}$ (that means for a constant density of the vapor) plus two contributions associated to compressible effects. The first one is $Δv_{C,E}$, the volume variation induced by the expansion i.e by the density variation induced by the variation of pressure. The second one is $Δv_{C,PC}$, the gain of volume induced by compressible phase change, i.e phase change while $ρ$ is changing. Note that $Δv_{C,PC}$ and $Δv_{C,E}$ are both strictly positive values for an expansion. This implies that the volume variation accounting for compressible effects is larger than the one for an incompressible flow $Δv_C > Δv_{I,PC}$.

The bar chart Fig. 17 illustrates the three contributions $Δv_{C,E}$, $Δv_{I,PC}$ and $Δv_{C,PC}$ for the four simulations under analysis, for the expansion-controlled growth (corresponding to Fig. (13.a)). It appears that $Δv_C$ increases with $Δρ$ and $Δt$ as previously discussed. Looking into the different contributions, it appears that $Δv_{C,E}$ only depends on $Δρ$, which is consistent with the fact that this contribution only depends on the variation of the density which follows instantaneously the variation of pressure in the vessel. On the other hand, the volume variation induced by phase change increases with $Δt$. In particular it is interesting to note that $Δv_{C,PC}$, the contribution associated with the compressible phase change, becomes more important for higher $Δρ$. Quantitatively, the gain of volume resulting from the compressible effects, $Δv_{C,E} + Δv_{C,PC}$ at the end of depressurization accounts for approximately 10 to 20% of the total volume variation $Δv_C$, indicating that this contribution is indeed significant and has to be taken into account. An important result of this analysis is that 99% of total volume variation $Δv_C$ is attributed to phase change $Δv_{I,PC} + Δv_{C,PC}$, and the expansion only has a minor contribution in the bubble growth $Δv_{C,E}$ (around 1%). This highlights the importance of accurately describing phase change when modeling cavitation in cryogenic fluids.

### 4.3. Wall heat transfer

Pool cavitation at the wall is similar to some extent to nucleate boiling in a superheated liquid, at least in the diffusion-controlled growth. The simulations of the experimental test case as well as the four simulations under investigation in this section show that, even for a wetting fluid, an apparent contact angle around 25° is established because of the effect of the nano-region. It was qualitatively observed that this has an impact on the solid wall temperature. The aim of the

<figure><img src="image_20.png" /><figcaption>Fig. 17. Gain of volume at the end of the expansion split into the three contributions.</figcaption></figure>

present section is to quantify the impact of the bubble over the wall heat transfer.

The integral heat flux $Q$ is defined by the following equation:

$$Q = 2\pi \int_{0}^{l_{x}} q(r) r dr, \qquad (26)$$

where $q(r)$ is the local radial specific wall heat flux (in W/m²). The absolute value of $Q$ is reported in Fig. 18 for the four simulations as a function of $t$. It starts from a zero value for $t = 0$ s as expected since the fluid and the wall are at the same temperature in the initial condition. During the expansion, $T_{sat}$ decreases and therefore the wall is cooled down by the bubble. It appears that the higher the $Δρ$, the higher the $Q$, in agreement with the observed behavior of the bubble growth which is greater for higher $Δρ$. The influence of the $Δt$ is also visible in the expansion-controlled growth. For longer expansions, the bubble grows more and so does the $Q$ extracted from the wall.

In order to separate the effects of the bubble growth and the effect of the heat transfer efficiency, we can define a Nusselt number with the following equation:

$$Nu = \frac{q_w}{q_{w,0}}. \qquad (27)$$

The numerator $q_w$ is the equivalent specific heat flux over a reference circular surface of radius $R$ considered as the influence surface of the bubble:

$$q_w = \frac{Q}{\pi R^2}. \qquad (28)$$

The reference heat flux $q_{w,0}$, i.e. the denominator in Eq. (27), can be interpreted as the specific heat flux we would have over the same surface $S$ by pure conduction because of a temperature varying between $T_{sat,0}$ at the wall and $T_{sat,f}$ in the liquid out of a boundary layer of thickness $R$:

$$q_{w,0} = \frac{k_l(T_{sat,0} - T_{sat,f})}{R}. \qquad (29)$$

Combining Eqs. (28) and (29) the following Nu number definition is obtained:

$$Nu = \frac{Q}{\pi R k_l (T_{sat,0} - T_{sat,f})}. \qquad (30)$$

Note that if there is no bubble, $Q = 0$ and therefore the definition of the denominator of Eq. (27) cannot be the heat flux we would have without a bubble (as usually done for the definition of the Nusselt number for nucleate boiling, see for instance [32]). The evolution of Nu versus time is shown in Fig. 19 for the four simulations.

Nu = 0 for $t = 0$ s because $Q = 0$. Then, as the expansion starts, Nu increases because there is a decrease of $T_{sat}$ and the wall is cooled down. Also the bubble radius increases, and eventually Nu reaches a maximum value before the end of the expansion. Then, Nu tends
<figure><img src="image_21.png" /><figcaption>Fig. 18. Evolution of the integrated heat flux over time.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 19. Evolution of the Nusselt number over time.</figcaption></figure>

to a constant value in the diffusion-controlled region where both Q (numerator) and R (denominator) are proportional to $\sqrt{t}$.

It is of interest to compare the energy extracted from the wall $E_Q$ (the time integrated Q), with respect to the energy used for liquid vapor phase change $E_{PC}$:

$$E_{PC} = \int_{t=0.0s}^{t=0.6s} \frac{\dot{m}(t)}{h_{lv}} dt, \qquad (31)$$

$$E_Q = \int_{t=0.0s}^{t=0.6s} Q(t) dt. \qquad (32)$$

The ratios $E_r = E_Q/E_{PC}$ obtained for the four simulations are reported in Table 4. In agreement with the observations made in Section 4.1, $E_{PC}$ is higher for a higher depressurization amplitude, and is only slightly affected by the depressurization rate. The energy ratio $E_r$ varies from 18% to 27%. These ratios are quite important and show that the interaction between the wall heat transfer and bubble cavitation is not negligible. It also appears that $E_r$ decreases for higher depressurization amplitudes. This can be explained by the fact that a higher depressurization implies a larger bubble and, therefore, more surface exchange between the vapor and the liquid for a similar dry surface.

## 5. Conclusion

A direct numerical solver for compressible flows, including phase change, has been employed to simulate a single bubble at the wall during a cavitation process in microgravity, i.e pool cavitation. This solver

<table><caption>Table 4
Comparison between the energy used for phase change $E_{PC}$ and energy associated with the heat transferred from the solid wall to the fluid $E_Q$.</caption><thead><tr><th>ε<sub>p</sub>[%]</th><th>Δt [ms]</th><th>E<sub>Q</sub> [mJ]</th><th>E<sub>PC</sub> [mJ]</th><th>E<sub>r</sub> [%]</th></tr></thead><tbody><tr><td rowspan="2">5</td><td>50</td><td>3.31</td><td>12.5</td><td>26</td></tr><tr><td>200</td><td>2.88</td><td>10.3</td><td>27</td></tr><tr><td rowspan="2">25</td><td>50</td><td>22.1</td><td>116</td><td>18</td></tr><tr><td>200</td><td>18.2</td><td>93.4</td><td>19</td></tr></tbody></table>

incorporates a nano-region model at the contact line and conjugate heat transfer at the wall. It has been validated against experimental data, demonstrating its reliability. It was known that the nano-region plays a critical role in the nucleate boiling process, and the present work shows that it is also the case in the cavitation process for wettable fluids. Indeed the nano-region significantly modifies the apparent contact angle and therefore the dry surface underneath the bubble. Consequently, this has a substantial impact on wall heat transfer. The present results emphasize that this effect cannot be neglected in accurate modeling.

The study of bubble's growth during the process has identified three distinct growth phases: an expansion-controlled phase, a diffusion-controlled growth phase and a transition phase between the two. The bubble growths as the square root of time, typical of boiling in a superheated liquid, during the diffusion controlled phase. On the other hand, the bubble growth and associated wall heat transfer are very different in the first two phases. In particular, it has been shown that during the expansion phase, compressible and variable density effects cannot be neglected, contributing up 20% on the total volume change. The results also demonstrated that wall heat transfer plays an important role in bubble cavitation at the wall, in particular for low depressurization amplitudes.

Future works will aim at extending the present work to different cryogenic fluids and different operating conditions in order to eventually deduce generalized models for the bubble growth and associated wall heat transfer.

## CRediT authorship contribution statement

**Maxence Deferrez:** Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Sébastien Tanguy:** Writing – review & editing, Supervision, Software, Methodology, Formal analysis. **Catherine Colin:** Writing – review & editing, Investigation, Conceptualization. **Annafederica Urbano:** Writing – review & editing, Writing – original draft, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

## Declaration of competing interest

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

## Acknowledgments

The authors acknowledge N. Weber and M. Dreyer for giving access to the experimental data, ESA for co-funding the PhD thesis of M. Deferrez (M2CRYO project), and S. Vincent-Bonnieu for his advices. This work was supported by the Chair SaCLaB2, resulting from the partnership between Airbus Defence and Space, Ariane Group and ISAE-SUPAERO, that co-fund the PhD thesis of M. Deferrez. This work was granted access to the CALMIP HPC resources (P23052). The authors acknowledge CNES for supporting the study in the framework of the GDR "Microgravité fondamentale et appliquée" (W-EEXP/10-02-01-02). All authors approved the version of the manuscript to be published.
<table><caption>**Table 5**
Deviation between thermodynamics properties at saturation conditions for T between 100 and 120 K computed with different cubic EoS and the NIST.</caption><thead><tr><th>[%]</th><th>PTV</th><th>VDW</th><th>PR</th><th>SRK</th></tr></thead><tbody><tr><td>r ∝ ρ<sub>v</sub><sup>-1/3</sup></td><td>3.5</td><td>-33.5</td><td>-0.1</td><td>1.0</td></tr><tr><td>ρ<sub>t</sub></td><td>6.8</td><td>-31.4</td><td>12.4</td><td>-0.4</td></tr><tr><td>ρ<sub>sat</sub></td><td>-9.1</td><td>220</td><td>0.7</td><td>-2.4</td></tr><tr><td>h<sub>L,V</sub></td><td>2.9</td><td>-41.1</td><td>-1.4</td><td>1.1</td></tr></tbody></table>

<figure><img src="image_23.png" /><figcaption>**Fig. 20.** Example of bubble detection with the presented algorithm.</figcaption></figure>

## Appendix A. Cubic EoS benchmark for pure methane

Among the different cubic EoS, the Soave Ridge Kwong (SRK) has been selected to conduct the present study because of its accuracy in reproducing the thermodynamics properties of methane in the conditions of the experiments conducted in [13]. In Table 5 a comparison between SRK and the NIST database [66] is shown for four properties at saturation conditions on a temperature range [100:120] K: $ρ_v$, $ρ_t$, $ρ_{sat}$ and $h_{L,V}$. The differences for all the properties are smaller than 2.5%. Other cubic EoS are also shown for comparison: Patel-Teja-Valderrama (PTV), Van der Waals (VDW) and Peng and Robinson (PR) [45].

## Appendix B. Apparent contact angle evaluation from experimental data

In the present work, experimental data available in [13] have been used to extract the apparent contact angle $θ_{app}$ versus time. This has been done using OpenCV image recognition library using the images available in [13]. An example of the bubble detection is given in Fig. 20 and an example of the $θ_{app}$ is given in Fig. 21. The bubble is considered a perfect circle of radius R and therefore detected using HoughCircles. The tube is detected with a combination between Canny, to detect edges, and HoughLinesP, to detect the lines within the edges. The contact angle is determined by computing the distance $l_c$ between the two points of the circle that cross the tube line:

$$θ_{app} = \sin^{-1} \left( \frac{l_c}{2R} \right). \qquad (33)$$

Note that only about 40% of the frames were usable by the algorithm over the different experiments, notably all frames before 0.16s are considered unusable since the bubble is too small.

## Appendix C. Nusselt number correlation for natural convection over a horizontal cylinder.

The following Nusselt number correlation from Kuehn and Goldstein [65] is considered:

$$Nu_{tbl} = \frac{2}{\ln \left[ 1 + \frac{2}{Nu_{tbl,0}} \right]} \qquad (34)$$

<figure><img src="image_24.png" /><figcaption>**Fig. 21.** Measured contact angle by the algorithm during growth of the bubble from the available experiments.</figcaption></figure>

$$Nu_{tbl,0} = \left[ Nu_{tbl,a}^{15} + Nu_{tbl,b}^{15} \right]^{1/15} \qquad (35)$$

$$Nu_{tbl,a} = 0.518Ra^{1/4} \left[ 1 + \left( \frac{0.559}{Pr} \right)^{3/5} \right]^{-5/12} \qquad (36)$$

$$Nu_{tbl,0b} = 0.1Ra^{1/3} \qquad (37)$$

with Rayleigh, Ra, and Prandtl, Pr, numbers defined as:

$$Ra = \frac{g\beta_l \Delta T D_{\text{tube}}}{v_l \alpha_l}, \qquad (38)$$

$$Pr = \frac{v_l}{\alpha_l}. \qquad (39)$$

with $g = 9.81 \text{ m/s}^2$ the gravity acceleration at sea level, $v_l$ the kinematic viscosity, $\alpha_l$ the thermal diffusivity, and $\beta_l$ thermal expansion coefficient. When using this correlation in Section 3 the thermodynamic and thermophysical properties are evaluated in saturation conditions at the initial pressure $p_{sat,0} = 1.404 \text{ bar}$. Therefore, the following values for the properties are obtained using the SRK EoS and the thermophysical properties reported as follow: $c_{p,l}(p_{sat,0}) = 3522.8 \text{ J g}^{-1} \text{ K}^{-1}$, $ρ_l = 413.7 \text{ kg m}^{-3}$ and $β_l = 3.33 \times 10^{-9} \text{ K}^{-1}$.

## Appendix D. Bubble's thermal boundary layer thickness

In order to compute the thermal boundary layer thickness around the bubble $δ_b$, the radial evolution of the temperature $T(x)$ across the center of the bubble is extracted. The dimensionless temperature profile $\bar{T}(r)$ is computed as:

$$\bar{T}(x) = \frac{T(x) - \min T}{\max T - \min T} \qquad (40)$$

with $\max T = T_{sat,0}$. $δ_b$ is defined as the distance between the interface and the radial position where the temperature reaches the initial $T_{sat,0}$ using the following criteria: $\bar{T} > (1 - \epsilon)$ with $\epsilon = 10^{-3}$ (see an example in Fig. 22 for two different depressurization rate and amplitude).

## Appendix E. Evaluation of the starting time for the diffusion controlled phase $t_{diff}$

In order to compute the diffusion time $t_{diff}$, the following dimensionless parameter is defined:

$$\bar{R}_{\text{diff}} = \frac{R_{\text{diff}} - \min R_{\text{diff}}}{\max R_{\text{diff}} - \min R_{\text{diff}}} \qquad (41)$$

with $\bar{R}_{\text{diff}}$ of $R_{\text{diff}} = \frac{\partial R/\sqrt{t}}{\partial t}$. For $t > \Delta t$, the starting time $t_{diff}$ for the diffusion-controlled regime is defined as the time $t$ for which $\bar{R}_{\text{diff}} < \epsilon$ with $\epsilon = 10^{-3}$.
<figure><img src="image_25.png" /><figcaption>Fig. 22. Dimensionless temperature profile for four different simulations for a same time $t = 0.5$ s.</figcaption></figure>

## Data availability

Data will be made available on request.

## References

[1] NASA, NASA technology roadmaps TA 14: Thermal management systems, 2015.

[2] NASA, Cryogenic fluid management technologies enabling for the artemis program and beyond, 2020.

[3] K. Petersen, S. Rahbarimanesh, J. Brinkerhoff, Progress in physical modelling and numerical simulation of phase transitions in cryogenic pool boiling and cavitation, Applied Math Modeling 116 (2023) 327-349, http://dx.doi.org/10.1016/j.apm.2022.11.028.

[4] A. Simonini, M. Dreyer, A. Urbano, F. Sanfedino, T. Himeno, P. Behruzi, M. Avila, J. Pinho, L. Peveroni, J.-B. Gouriet, Cryogenic propellant management in space: open challenges and perspectives, Npj Microgravity 10 (34) (2024) 1-13, http://dx.doi.org/10.1038/s41526-024-00377-5.

[5] C. Colin, O. Kannengieser, W. Bergez, M. Lebon, J. Sebilleau, M. Sagan, S. Tanguy, Nucleate pool boiling in microgravity: Recent progress and future prospects, Comptes Rendus Mecanique (2017).

[6] C.E. Brennen, Cavitation and Bubble Dynamics, vol. 44, Oxford University Press, 1995, http://dx.doi.org/10.1017/CBO9781107338760.

[7] J.-P. Franc, J.-P. Michel, Fundamentals of Cavitation, Kluwer Academic Publishers, 2004.

[8] A. Ferrari, Fluid dynamics of acoustic and hydrodynamic cavitation in hydraulic power systems, Proc. R. Soc. Lond. A 473 (2017) 20160345.

[9] J. Bae, J. Yoo, L. Jin, S. Jeong, Experimental investigation of passive thermodynamic vent system (TVS) with liquid nitrogen, Cryogenic 89 (2018) 147-156, http://dx.doi.org/10.1016/j.cryogenics.2017.11.001.

[10] Z. Liu, Y. Lia, G. Lei, S. Xia, C. Li, Experimental study on refrigeration performance and fluid thermal stratification of thermodynamic vent experimental study on refrigeration performance and fluid thermal stratification of thermodynamic vent, Int. J. Refrig. 88 (2018) 496-505, http://dx.doi.org/10.1016/j.ijrefrig.2018.02.022.

[11] J. Hartwig, Propellant management devices for low-gravity fluid management: Past, present, and future applications, J. Spacecrafts and Rockets 54 (4) (2017) 808-824, http://dx.doi.org/10.2514/1.A33750.

[12] O. Kannengieser, C. Colin, W. Bergez, Pool boiling with non-condensable gas in microgravity: Results of a sounding rocket experiment, Microgravity Sci. Technol. 22 (3) (2010) 447-454, http://dx.doi.org/10.1007/s12217-010-9211-9.

[13] N.H. Weber, M.E. Dreyer, Depressurization induced vapor bubble growth in liquid methane during microgravity, Cryogenics 134 (2023) 103716, http://dx.doi.org/10.1016/j.cryogenics.2023.103716.

[14] M. Rodio, M.D. Giorgi, A. Ficarella, Influence of convective heat transfer modeling on the estimation of thermal effects in cryogenic cavitating flows, Int. J. Heat and Mass Transfer 55 (2012) 6538-6554, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2012.06.060.

[15] B.B. Mikic, W.M. Rohsenow, P. Griffith, On bubble growth rates, Int. J. Heat Mass Transfer 13 (4) (1970) 657-666, http://dx.doi.org/10.1016/0017-9310(70)90040-2.

[16] Z. Wang, S. Bankoff, Bubble growth on a solid wall in a rapidly-depressurizing liquid pool, Int. J. Multiph. Flow 17 (4) (1991) 425-437, http://dx.doi.org/10.1016/0301-9322(91)90040-A.

[17] O.C. Jones, N. Zuber, Bubble growth in variable pressure fields, J. Heat Transf. 100 (3) (1978) 453-459, http://dx.doi.org/10.1115/1.3450830.

[18] A. Prosperetti, M.S. Plesset, Vapour-bubble growth in a superheated liquid, J. Fluid Mech. 85 (2) (1978) 349-368, http://dx.doi.org/10.1017/S0022112078000671.

[19] Y.S. Cha, R.E. Henry, Bubble growth during decompression of a liquid, J. Heat Transf. 103 (1) (1981) 56-60, http://dx.doi.org/10.1115/1.3244430.

[20] M.-H. Chun, A model for vapor bubble growth under variable liquid pressure fields, Int. Commun. Heat Mass Transfer 15 (6) (1988) 731-739, http://dx.doi.org/10.1016/0735-1933(88)90016-4.

[21] L.E. Scriven, On the dynamics of phase growth, Chem. Eng. Sci. 10 (1) (1959) 1-13, http://dx.doi.org/10.1016/0009-2509(59)80019-1.

[22] M.S. Plesset, S.A. Zwick, The growth of vapor bubbles in superheated liquids, J. Appl. Phys. 25 (4) (1954) 493-500, http://dx.doi.org/10.1063/1.1721668.

[23] Y. Huang, X. Qin, Growth behavior of bubbles containing non-condensable gas in superheated cryogenic liquids, Cryogenics 146 (2025) 104016, http://dx.doi.org/10.1016/j.cryogenics.2024.104016.

[24] E. Valero, I.E. Parra, The role of thermal disequilibrium in critical two-phase flow, Int. J. Multiph. Flow 28 (1) (2002) 21-50, http://dx.doi.org/10.1016/S0301-9322(01)00064-7.

[25] A. Prosperetti, Vapour bubbles, Ann. Rev. Fluid Mech. 49 (2017) 221-248, http://dx.doi.org/10.1146/annurev-fluid-010816-060221.

[26] P.C. Stephan, C.A. Busse, Analysis of the heat transfer coefficient of grooved heat pipe evaporator walls, Int. J. Heat Mass Transfer 35 (2) (1992) 383-391, http://dx.doi.org/10.1016/0017-9310(92)90276-X.

[27] V. Janeček, V.S. Nikolayev, Apparent-contact-angle model at partial wetting and evaporation: Impact of surface forces, Phys. Rev. E 87 (1) (2013) 012404, http://dx.doi.org/10.1103/PhysRevE.87.012404.

[28] 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) http://dx.doi.org/10.1016/j.jcp.2023.112602.

[29] Y. Sato, B. Ničeno, A sharp-interface phase change model for a mass-concervative interface tracking method, J. Comput. Phys. 249 (2013) 127-161, http://dx.doi.org/10.1016/j.jcp.2013.04.035.

[30] 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.

[31] A. Urbano, S. Tanguy, G. Huber, C. Colin, Direct numerical simulation of nucleate boiling in micro-layer regime, Int. J. Heat and Mass Transfer 123 (2018) 1128-1137, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.02.104.

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

[33] T. Long, J. Pan, S. Zaleski, An edge-based interface tracking (EBIT) method for multiphase flows with phase change, J. Comput. Phys. 513 (2024) 113159, http://dx.doi.org/10.1016/j.jcp.2024.113159.

[34] M. Jemison, M. Sussman, M. Arienti, Compressible, multiphase semi-implicit method with moment of fluid interface representation, J. Comput. Phys. 279 (2014) 182-217.

[35] F. Xiao, R. Akoh, S. Ii, Unified formulation for compressible and incompressible flows by using multi-integrated moments II: Multi-dimensional version for compressible and incompressible flows, J. Comput. Phys. 213 (2006) 31-56, http://dx.doi.org/10.1016/j.jcp.2003.10.014.

[36] N. Kwater, J. Su, J.T. Grétarsson, R. Fedkiw, A method for avoiding the acoustic time step restriction in compressible flow, J. Comput. Phys. 228 (2009) 4146-4161.

[37] G. Huber, S. Tanguy, J.-C. Béra, B. Gilles, A time splitting projection scheme for compressible two-phase flows, application to the interaction of bubbles with ultrasound waves, J. Comput. Phys. 302 (2015) 439-468.

[38] D. Fuster, S. Popinet, An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension, J. Comput. Phys. 374 (2018) 752-768, http://dx.doi.org/10.1016/j.ijmultiphaseflow.2016.09.017.

[39] M. Arienti, M. Sussman, A numerical study of the thermal transient in high-pressure diesel injection, Int. J. Multiphase Flow 88 (2017) 205-221, http://dx.doi.org/10.1016/j.ijmultiphaseflow.2016.09.017.

[40] M.B. Kuhn, O. Desjardins, An all-mach, low-dissipation strategy for simulating multiphase flows, J. Comput. Phys. 445 (2021) 110602.

[41] A. Urbano, M. Bibal, S. Tanguy, A semi implicit compressible solver for two-phase flows of real fluids., J. Comput. Phys. 456 (2022) 111034, http://dx.doi.org/10.1016/j.jcp.2022.111034.

[42] B. Duret, R. Canu, J. Reveillon, F. Demoulin, A pressure based method for vaporizing compressible two-phase flows with interface capturing approach, Int. J. Multiphase Flow 108 (2019) 42-50, http://dx.doi.org/10.1016/j.ijmultiphaseflow.2018.06.022.

[43] S. Cho, G. Son, A level-set method for bubble growth in acoustic droplet vaporization, Int. Commun. Heat Mass Transfer 93 (2018) 83-92, http://dx.doi.org/10.1016/j.icheatmasstransfer.2018.02.016.

[44] S. Cho, G. Son, Numerical simulation of acoustic droplet vaporization near a wall, Int. Commun. Heat Mass Transfer 99 (2018) 7-17, http://dx.doi.org/10.1016/j.icheatmasstransfer.2018.10.007.

[45] M. Bibal, M. Deferrez, S. Tanguy, A. Urbano, A compressible solver for two-phase flows with phase change for bubble cavitation, J. Comput. Phys. 500 (2024) 112750, http://dx.doi.org/10.1016/j.jcp.2023.112750.
[46] S. Tanguy, T. Ménard, A. Berlemont, J.-L. Estivalezes, F. Couderc, Level set method for interface tracking: development and applications, Houille Blanche 92 (2006) 23-31.

[47] S. Tanguy, T. Ménard, A. Berlemont, A level set method for vaporizing two-phase flows, J. Comput. Phys. 221 (2007) 837-853, http://dx.doi.org/10.1016/j.jcp.2006.07.003.

[48] B. Lalanne, L.R. Villegas, S. Tanguy, F. Risso, On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method, J. Comput. Phys. 301 (2015) 289-307, http://dx.doi.org/10.1016/j.jcp.2015.08.036.

[49] A. Dalmon, M. Lepilliez, S. Tanguy, A. Pedrono, B. Busset, H. Bavestrello, J. Mignot, Direct numerical simulation of a bubble motion in a spherical tank under external forces and microgravity conditions, J. Fluid Mech. 849 (2018) 467-497.

[50] G. Mialhe, S. Tanguy, L. Tranier, E.-R. Popescu, D. Legendre, An extended model for the direct numerical simulation of droplet evaporation. Influence of the marangoni convection on leidenfrost droplet, J. Comput. Phys. (2023) 112366.

[51] R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Ann. Rev. Fluid Mech. 31 (1) (1999) 567-603, http://dx.doi.org/10.1146/annurev.fluid.31.1.567.

[52] F. Gibou, L. Chen, D.Q. Nguyen, S. Banerjee, A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change, J. Comput. Phys. 222 (2007) 536-555, http://dx.doi.org/10.1016/j.jcp.2006.07.035.

[53] M. Ishii, T. Hibiki, Thermo-Fluid Dynamics of Two-Phase Flow, second ed., Springer-Verlag New York, 2011.

[54] J.O. Valderrama, A generalized patel-teja equation of state for polar and nonpolar fluids and their mixtures, J. Chem. Eng. Japan 23 (1) (1990) 87-91, http://dx.doi.org/10.1252/jcej.23.87.

[55] B. Mathieu, Etudes physique, expérimentale et numérique des mécanismes de base intervenant dans les écoulements diphasiques (Ph.D. thesis), Aix-Marseille 1, 2003, URL https://theses.fr/2003AIX11039.

[56] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1) (1994-09-01) 146-159, http://dx.doi.org/10.1006/jcph.1994.1155.

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

[58] X.-D. Liu, R.P. Fedkiw, M. Kang, A boundary condition capturing method for Poisson's equation on irregular domains, J. Comput. Phys. 160 (1) (2000) 151-178.

[59] M. Sussman, K.M. Smith, M.Y. Hussaini, M. Ohta, R. Zhi-Wei, A sharp interface method for incompressible two-phase flows, J. Comput. Phys. 221 (2) (2007-02-10) 469-505, http://dx.doi.org/10.1016/j.jcp.2006.06.020.

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

[61] L. Rueda Villegas, R. Alis, M. Lepilliez, S. Tanguy, A ghost fluid/level set method for boiling flows and liquid evaporation: Application to the leidenfrost effect, J. Comput. Phys. 316 (2016-07-01) 789-813, http://dx.doi.org/10.1016/j.jcp.2016.04.031.

[62] J.E. Dendy, Black box multigrid, J. Comput. Phys. 48 (3) (1982) 366-386, http://dx.doi.org/10.1016/0021-9991(82)90057-2.

[63] R. Borges, M. Carmona, B. Costa, W.S. Don, An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws, J. Comput. Phys. 227 (6) (2008) 3191-3211, http://dx.doi.org/10.1016/j.jcp.2007.11.038.

[64] 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, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.05.083.

[65] T.H. Kuehn, R.J. Goldstein, Correlating equations for natural convection heat transfer between horizontal circular cylinders, Int. J. Heat Mass Transfer 19 (10) (1976) 1127-1134, http://dx.doi.org/10.1016/0017-9310(76)90145-9.

[66] N.S.R. Data, Thermophysical properties of fluid systems, 2025, URL https://webbook.nist.gov/chemistry/fluid. [Online] (Accessed 9 April 2025).