1. Introduction
Nucleate boiling is recognised as one of the most efficient heat-transfer processes due to the significant latent heat of evaporation (Chen et al. Reference Chen, Yu, Lu, Wang, Sun, Jiao, Zhang and Tao2024). Consequently, it plays a critical role in various industrial applications, such as nuclear reactors (Manglik Reference Manglik2006), electronics cooling (Cheng & Xia Reference Cheng and Xia2017) and thermal management subsystems in aerospace engineering (Dhruv et al. Reference Dhruv, Balaras, Riaz and Kim2019). In nucleate boiling, a very thin liquid layer, known as the microlayer, may form between the heated wall and the liquid–vapour interface during bubble growth (Hänsch & Walker Reference Hänsch and Walker2019). This bubble growth regime is referred to as the microlayer regime (Urbano et al. Reference Urbano, Tanguy, Huber and Colin2018). As depicted in figure 1, the microlayer, which is merely a few microns thick, extends radially up to a few millimetres. The small thickness of the microlayer leads to very high heat flux (Zupančič et al. Reference Zupančič, Gregorčič, Bucci, Wang, Aguiar and Bucci2022) and its considerable radial extent over the heated surface results in a significant contribution to the overall heat transfer (Yabuki & Nakabeppu Reference Yabuki and Nakabeppu2014). In addition, the drying of the microlayer leads to the spreading of dry spots, making its dynamics crucial for understanding critical heat flux (CHF) (Zhao, Masuoka & Tsuruta Reference Zhao, Masuoka and Tsuruta2002), which is the maximum heat flux that can be transferred through nucleate boiling. These features of the microlayer have motivated intense scientific investigation into it.
In recent decades, various modern high-resolution experimental techniques have been developed, significantly advancing the understanding of the dynamics and heat-transfer characteristics of the microlayer (Utaka, Kashiwabara & Ozaki Reference Utaka, Kashiwabara and Ozaki2013; Jung & Kim Reference Jung and Kim2014; Bucci et al. Reference Bucci, Richenderfer, Su, McKrell and Buongiorno2016). Despite these advances, experimental approaches still face several limitations and shortcomings. For example, direct measurement of microlayer thickness using laser interferometry (Jung & Kim Reference Jung and Kim2018; Narayan & Srivastava Reference Narayan and Srivastava2021) inevitably introduces systematic errors due to the loss of the first fringes near the contact line, where the interface slope exceeds the observable limit (Tecchio et al. Reference Tecchio, Zhang, Cariteau, Zalczer, Bulkin, Charliac, Vassant and Nikolayev2024). In addition, the inherent uncertainties in experiments require careful evaluation (Kim et al. Reference Kim, Sergis, Kim and Hardalupas2020). The range of controllable parameters in experiments is also limited by achievable laboratory conditions (Hänsch & Walker Reference Hänsch and Walker2019), which may affect the generality of experimental findings.
Compared with experimental approaches, numerical simulations are more flexible in changing operating conditions and can provide comprehensive physical information (Bureš & Sato Reference Bureš and Sato2022). Nowadays, with rapid advances in high-performance computing, direct numerical simulation (DNS) has become a powerful tool for studying nucleate boiling and has gained increasing interest in academia (Chen et al. Reference Chen, Yu, Lu, Wang, Sun, Jiao, Zhang and Tao2024). Despite remarkable progress in recent years, DNS of nucleate boiling in the microlayer regime remains a significant challenge due to the inherent multiscale nature of the problem. To capture the entire bubble, the domain must extend to several millimetres, while grid sizes below one micron are required to explicitly resolve the microlayer (Urbano et al. Reference Urbano, Tanguy, Huber and Colin2018). Although a number of microlayer models have been developed to conduct multiscale simulations of nucleate boiling (Sato & Niceno Reference Sato and Niceno2015; Chen et al. Reference Chen, Jin, Yu, Ling, Sun, Zhang, Jiao and Tao2023), their application to flow boiling, which involves bubble slipping and coalescence, is not straightforward (Chen et al. Reference Chen, Yu, Lu, Wang, Sun, Jiao, Zhang and Tao2024). In contrast, once a DNS code is validated, it should be applicable across different scenarios. Without accounting for heat transfer, Hänsch & Walker (2016) first performed a pioneering DNS study on the hydrodynamics of microlayer formation in nucleate boiling. In their simulations, the bubble-growth rate was calculated from the analytical solution of Scriven (Reference Scriven1959). Subsequently, several researchers have studied the hydrodynamics of the microlayer using the same modelling strategy, except that the bubble-growth rate may be deduced from other models or experimental results. Guion et al. (Reference Guion, Afkhami, Zaleski and Buongiorno2018) performed simulations over a broad range of conditions and identified the minimum set of dimensionless parameters governing the hydrodynamics of microlayer formation. Giustini et al. (Reference Giustini, Kim, Issa and Bluck2020) simulated microlayer formation in the boiling of industrially relevant fluids, whose properties differ significantly from typically studied fluids such as water and ethanol. Following this, Giustini (Reference Giustini2024) studied the formation of a dewetting ridge near the contact line during bubble growth. Most recently, Saini et al. (Reference Saini, Chen, Zaleski and Fuster2024) investigated microlayer formation during heterogeneous bubble nucleation triggered by a sudden drop in ambient pressure.
Schematic of nucleate boiling with a microlayer (not to scale).
In contrast to previous studies, Urbano et al. (Reference Urbano, Tanguy, Huber and Colin2018) were the first to successfully simulate nucleate boiling in the microlayer regime with resolved thermal effects. A rigorous parametric study was carried out to determine the criterion for microlayer formation. Zhang et al. (Reference Zhang, Rafique, Ding, Bolotnov and Hampel2023a ) developed an unstructured-mesh-based solver to simulate microlayer formation and evaporation driven by the local temperature gradient. This solver was later extended to study nucleate boiling on surfaces with micro-pillars (Zhang et al. Reference Zhang, Rafique, Ding, Bolotnov and Hampel2023b ), where the microlayer may be disturbed or disrupted by the pillars. It should be noted that in these works (Urbano et al. Reference Urbano, Tanguy, Huber and Colin2018; Zhang et al. Reference Zhang, Rafique, Ding, Bolotnov and Hampel2023a , Reference Zhang, Rafique, Ding, Bolotnov and Hampelb ), only heat transfer in the fluid is considered, while conjugate heat transfer between the fluid and solid, which has been shown to be important in microlayer evaporation (Sato & Niceno Reference Sato and Niceno2015; Ding, Krepper & Hampel Reference Ding, Krepper and Hampel2018), is neglected. Hänsch & Walker (Reference Hänsch and Walker2019) extended their previous work (Hänsch & Walker 2016) to simulate the depletion of the microlayer, where the evaporation of the microlayer is determined based on the solid temperature obtained by solving the conjugate heat transfer. However, in this work, the bubble growth rate is still approximated using the Scriven solution (Scriven Reference Scriven1959) instead of solving the heat transfer in the liquid domain. Subsequently, significant progress was made by Bureš & Sato (Reference Bureš and Sato2022), who were the first to perform detailed DNS studies of nucleate boiling in the microlayer regime by explicitly resolving heat and mass transfer in the fluid, along with conjugate heat transfer. They simulated an experiment conducted at Massachusetts Institute of Technology (MIT) (Bucci Reference Bucci2020), and their numerical results were in good agreement with experimental measurements. The same experiment was later adopted by Torres et al. (Reference Torres, Urbano, Colin and Tanguy2024) and El Mellas et al. (Reference El Mellas, Samkhaniani, Falsetti, Stroh, Icardi and Magnini2024) to validate their computational models with resolved conjugate heat transfer. However, the microlayer was not the focus of Torres et al. (Reference Torres, Urbano, Colin and Tanguy2024) and El Mellas et al. (Reference El Mellas, Samkhaniani, Falsetti, Stroh, Icardi and Magnini2024), and it was not studied in detail.
Although the above-mentioned works demonstrated the capabilities of DNS solvers in modelling the microlayer in nucleate boiling, the computational cost of DNS still limits its application. As discussed earlier, the multiscale nature of nucleate boiling requires a large computational domain with a small grid size, leading to significant computational costs. Despite the efforts by Bureš & Sato (Reference Bureš and Sato2022) to stretch the grid and minimise the computational burden, the simulations remained highly time-consuming. For a simulation at the finest resolution, 336 cores were used, requiring approximately 400 000 CPU hours to advance the physical time to 0.5 ms (Bureš & Sato Reference Bureš and Sato2022). By contrast, a complete bubble cycle from nucleation to detachment is typically of the order of 10 ms. In fact, to the best of our knowledge, DNS of a complete bubble cycle with a resolved microlayer and conjugate heat transfer has not been achieved before due to the high computational cost. Employing a quad/octree-based AMR method, the free software Basilisk (Popinet Reference Popinet2009, Reference Popinet2015) provides a highly efficient framework for the DNS of multiphase flows (Wang et al. Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023; Pan et al. Reference Pan, Long, Chirco, Scardovelli, Popinet and Zaleski2024). Several phase-change models (Zhao, Zhang & Ni Reference Zhao, Zhang and Ni2022; Long, Pan & Zaleski Reference Long, Pan and Zaleski2024; Cipriano et al. Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024) have been developed using Basilisk, demonstrating its potential for studying complex boiling problems. In this work, we aim to extend the previous open-source phase-change model developed by Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024) to include conjugate heat transfer and to simulate nucleate boiling with a microlayer. In addition to incorporating conjugate heat transfer, the original phase-change model requires careful modifications to correctly capture the physics within the microlayer. This is due to the approximate projection method required by AMR, which leads to unphysical oscillations in phase-change problems as the velocity is discontinuous across the interface (Zhao et al. Reference Zhao, Zhang and Ni2022; Long et al. Reference Long, Pan and Zaleski2024). Although the previous model of Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024) appears to perform well in benchmark tests, the problem of interest in this paper is much more challenging due to the intense heat and mass transfer within the microlayer. We will demonstrate that the unphysical oscillations induced by the approximate projection method lead to incorrect microlayer dynamics, resulting in inaccurate predictions of heat and mass transfer. We adopt the ghost fluid method (Nguyen, Fedkiw & Kang Reference Nguyen, Fedkiw and Kang2001; Tanguy et al. Reference Tanguy, Sagan, Lalanne, Couderc and Colin2014) to address this issue, in which the singularity is removed by setting the ghost velocity according to the jump condition. By doing so, we achieve highly efficient and accurate DNS of nucleate boiling in the microlayer regime. With the present solver, we have successfully simulated the experiment conducted at MIT (Bucci Reference Bucci2020). The results of the present solver not only agree well with the numerical results of Bureš & Sato (Reference Bureš and Sato2022), but also achieve a remarkable reduction in computational cost by three orders of magnitude, while significantly enhancing stability. The influence of the contact angle is further investigated by analysing the contributions of different regions of the bubble to evaporation for various contact angles. Moreover, a complete bubble cycle for nucleate boiling in the microlayer regime is directly simulated with all effects explicitly resolved, which, to our knowledge, is the first such study reported in the open literature.
The remainder of this paper is structured as follows. In § 2, we briefly review the original phase-change model of Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024). Then, we introduce the implementation of the ghost fluid method and the treatment of the solid. The extended model is verified in § 3, followed by the numerical results and related discussions presented in § 4. Finally, concluding remarks are provided in § 5.
2. Numerical method
2.1. Governing equations
In nucleate boiling simulations, the fluid domain is occupied by the liquid and vapour phases, which are separated by a zero-thickness interface,
$\varGamma$
. Both phases are assumed to be incompressible and monocomponent. With phase change considered, the governing equations for fluids are as follows:
\begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{u} &= S_{\textit{pc}}, \end{align}
\begin{align} \rho \left ( \frac {\partial \boldsymbol{u}}{\partial t} + \left (\boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{\nabla }\right )\boldsymbol{u} \right ) &= -\boldsymbol{\nabla} \! p + \boldsymbol{\nabla } \boldsymbol{\cdot } \big ( \mu \big( \boldsymbol{\nabla} \kern-1pt \boldsymbol{u} + \boldsymbol{\nabla} \kern-1pt \boldsymbol{u}^T\big) \big ) + \rho \boldsymbol{g} + \boldsymbol{f}_{\kern-2pt \sigma} , \end{align}
\begin{align} \rho C_{\kern-1.5pt p} \left ( \frac {\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{\nabla} \kern-1pt T\right ) & = \boldsymbol{\nabla } \boldsymbol{\cdot } (\lambda \boldsymbol{\nabla} \kern-1pt T). \end{align}
Here,
$\boldsymbol{u}$
,
$\rho$
,
$\mu$
,
$\boldsymbol{g}$
,
$T$
,
$C_{\kern-1.5pt p}$
and
$\lambda$
represent the fluid velocity, density, dynamic viscosity, gravitational acceleration, temperature, specific heat and thermal conductivity, respectively. The source term
$S_{\textit{pc}}$
, introduced due to phase change, will be elaborated on later. The surface tension force is
$\boldsymbol{f}_{\kern-2pt \sigma} = \sigma \kappa \delta _s \boldsymbol{n}$
, where
$ \sigma$
is the surface tension,
$\kappa$
the interface curvature,
$\boldsymbol{n}$
the interface normal vector pointing from the liquid phase to the vapour phase and
$\delta _s$
the Dirac delta function concentrated at the interface. Following previous studies (Bureš & Sato Reference Bureš and Sato2022; Torres et al. Reference Torres, Urbano, Colin and Tanguy2024), the surface tension
$\sigma$
is assumed to be constant, implying that Marangoni convection effects are neglected. To resolve conjugate heat transfer, the following energy equation in the solid domain is considered:
\begin{equation} \rho _s C_{\kern-1.5pt p,s} \frac {\partial T_s}{\partial t} = \boldsymbol{\nabla } \boldsymbol{\cdot } (\lambda _s \boldsymbol{\nabla} \kern-1pt T_s) + Q_h, \end{equation}
where the subscript
$s$
denotes the prosperity of solid and
$Q_h$
is a volumetric power term due to electrical resistance heating.
2.2. Jump conditions
The main challenge in the simulations of boiling flows is the accurate implementation of jump conditions at the interface. Let
$H$
represent the Heaviside function (Tanguy et al. Reference Tanguy, Sagan, Lalanne, Couderc and Colin2014), defined as
$1$
in the liquid phase and
$0$
in the vapour phase. Accordingly, the jump of a given fluid property
$\phi$
, such as density and viscosity, across the interface can be expressed as
\begin{equation} \phi = \phi _l H + \phi _v (1 - H), \end{equation}
where the subscripts
$l$
and
$v$
indicate the physical properties of the liquid and vapour, respectively. When phase change occurs, the mass balance at the interface gives
\begin{equation} \dot {m} = \rho _l (\boldsymbol{u}_l - \boldsymbol{u}_{\kern-1.5pt \varGamma} )\boldsymbol{\cdot } \boldsymbol{n} = \rho _v (\boldsymbol{u}_v - \boldsymbol{u}_{\kern-1.5pt \varGamma} )\boldsymbol{\cdot } \boldsymbol{n}, \end{equation}
where
$\dot {m}$
represents the mass flux. Introducing the jump operator
$[\phi ]_\varGamma = \phi _l - \phi _v$
in (2.6), the velocity jump can be formulated as
\begin{equation} [\boldsymbol{u}]_\varGamma = \dot {m}\left [ \frac {1}{\rho } \right ]_\varGamma \boldsymbol{n}. \end{equation}
Moreover, the stress jump condition across the interface can be obtained by integrating the momentum (2.2) and is usually given in the form of the pressure jump as
\begin{equation} [p]_\varGamma = \sigma \kern-1pt \kappa + 2\left [\mu \frac {\partial u_n}{\partial n} \right ]_\varGamma - \dot {m}^2 \left [ \frac {1}{\rho } \right ]_\varGamma , \end{equation}
where
${\partial u_n}/{\partial n}$
is the normal derivative of the normal velocity component
$u_n = \boldsymbol{u}\boldsymbol{\cdot } \boldsymbol{n}$
. Note that the above-mentioned jump conditions are also validate in the absence of phase change, i.e.
$\dot {m} = 0$
.
Schematic for the concept of equivalent conductive resistance: (a) IHTR located at the liquid–vapour interface; (b) IHTR located at the fluid–solid boundary.
For the temperature field, energy conservation across the liquid–vapour interface leads to the jump condition
\begin{equation} [q]_\varGamma = \dot {m}h_{lv}, \end{equation}
where
$q = \lambda \boldsymbol{\nabla} \kern-1pt T \boldsymbol{\cdot } \boldsymbol{n}$
is the heat flux and
$h_{lv}$
is the latent heat. It should be noted that thermodynamic equilibrium,
\begin{equation} T_{\varGamma ,l} = T_{\varGamma ,v} = T_{\textit{sat}}, \end{equation}
is assumed in the derivation of (2.9) (Tanguy et al. Reference Tanguy, Sagan, Lalanne, Couderc and Colin2014), where
$T_{\textit{sat}}$
denotes the saturation temperature at the ambient system pressure. However, this assumption leads to large deviations in thermodynamic characteristics from experimental observations in the presence of a microlayer (Giustini et al. Reference Giustini, Jung, Kim and Walker2016). In fact, the degree of thermodynamic non-equilibrium within the microlayer remains an open question. Currently, a widely used modelling strategy to account for non-equilibrium effects is to introduce an additional interfacial heat-transfer resistance (IHTR) (Hänsch & Walker Reference Hänsch and Walker2019; Bureš & Sato Reference Bureš and Sato2022; Giustini et al. Reference Giustini, Jung, Kim and Walker2016). As shown in figure 2(a), IHTR models assume that the vapour temperature at the liquid–vapour interface remains the saturation temperature, while the liquid temperature is higher. In the present study, a simplified IHTR model proposed by Bureš & Sato (Reference Bureš and Sato2022) is employed. Since IHTR is usually assumed to be localised in the microlayer (Giustini et al. Reference Giustini, Jung, Kim and Walker2016), Bureš & Sato (Reference Bureš and Sato2022) demonstrated that its implementation could be simplified by introducing a numerically equivalent contact heat-transfer resistance at the fluid–solid boundary. In the presence of conjugate heat transfer, the heat flux is balanced between the solid and fluid domains (Torres et al. Reference Torres, Urbano, Colin and Tanguy2024), as expressed in the following equation:
\begin{equation} \lambda _s \boldsymbol{\nabla} \kern-1pt T_s \boldsymbol{\cdot } \boldsymbol{n}_s + \delta _q = \lambda _{\!f} \boldsymbol{\nabla} \kern-1pt T_{\!f} \boldsymbol{\cdot } \boldsymbol{n}_s\quad (f=l\ \text{or}\ v), \end{equation}
where
$\delta _q\ (\mathrm{W\,m^{-2}})$
denotes a Dirac source term representing a heat source localised at the solid–fluid boundary and
$\boldsymbol{n}_s$
is the normal vector at the fluid–solid boundary, pointing from the solid domain to the fluid domain. As shown in figure 2(a), with the IHTR, the heat flux from solid to vapour is
\begin{equation} j_q = \frac {T_{b,s} - T_{\textit{sat}}}{d/\lambda _l + R_\varGamma }, \end{equation}
where
$T_{b,s}$
is the solid temperature at the fluid–solid boundary,
$d$
is the microlayer thickness and
$R_\varGamma$
is the IHTR factor. Assuming the existence of a contact heat-transfer resistance
$R_c$
at the fluid–solid boundary, as shown in figure 2(b), the following equation is obtained:
\begin{equation} j_q = \frac {T_{b,s} - T_{\textit{sat}}}{d/\lambda _l + R_c} = \frac {T_{b,s} - T_{\textit{sat}}}{d/\lambda _l + R_\varGamma }, \end{equation}
provided that
$R_c = R_\varGamma$
. In this way, the temperature discontinuity due to heat-transfer resistance is shifted from the liquid–vapour interface to the liquid–solid interface, while the overall heat flux remains unchanged. With this simplified model, special care is only needed for the implementation of conjugate heat transfer, while the original phase-change model in the fluid domain (Cipriano et al. Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024), which assumes that temperature is continuous across the interface, can be applied directly, as (2.10) still holds.
2.3. Numerical scheme
2.3.1. One-fluid method
In the phase-change model of Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024), a geometric VOF method with piecewise linear interface construction (PLIC) was adopted to capture the liquid–vapour interface. The VOF function
$f$
, defined as the volume fraction of the reference phase (which is the liquid phase in the present study) in the control cell, is updated by solving
\begin{equation} \frac {\partial f}{\partial t}+\boldsymbol{u}_{\varGamma } \boldsymbol{\cdot } \boldsymbol{\nabla} \! f=0 \end{equation}
with a directional split advection method (Weymouth & Yue Reference Weymouth and Yue2010; Zhao et al. Reference Zhao, Zhang and Ni2022). By considering the mass balance (2.6), the interfacial velocity is calculated as
\begin{equation} \boldsymbol{u}_{\kern-1.5pt \varGamma} = \boldsymbol{u}_l - \frac {\dot {m}}{\rho _l} \boldsymbol{n} = \boldsymbol{u}_v - \frac {\dot {m}}{\rho _v} \boldsymbol{n}. \end{equation}
With the VOF function, the jumps of physical properties (2.5) can be approximated by
\begin{equation} \phi = \phi _l f + \phi _v (1 - f), \end{equation}
which gives the so-called one-fluid method (Tanguy et al. Reference Tanguy, Sagan, Lalanne, Couderc and Colin2014). Accordingly, the mass equation for the one-fluid velocity
$\boldsymbol{u}$
can be obtained by considering the divergence-free condition in the bulk region of each phase and the velocity jump at the interface (2.7), and is given by
\begin{equation} \boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{u} = \dot {m} \delta _s\left [ \frac {1}{\rho } \right ]_\varGamma . \end{equation}
Here,
$\delta _s$
is the Dirac delta function concentrated at the interface, which is approximated as
$\delta _s = S_\varGamma /V_c$
in the discretised formulation, where
$V_c$
is the volume of the computational cell and
$S_{\varGamma }$
denotes the area of the interface within it.
2.3.2. Solving the mass and momentum equations
The phase-change model of Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024) is developed in the free software Basilisk (Popinet Reference Popinet2009, Reference Popinet2015), which employs a quad/octree grid with AMR, where the velocity and pressure are collocated at cell centres. The incompressible Navier–Stokes equations are solved using a time-staggered approximate projection method, leading to the following discretisation:
\begin{align} &\rho _c^{n+\frac {1}{2}}\left (\frac {\boldsymbol{u}^*-\boldsymbol{u}^n}{\Delta t}+\left (\boldsymbol{u}^{n+\frac {1}{2}} \boldsymbol{\cdot } \boldsymbol{\nabla }\right ) \boldsymbol{u}^{n+\frac {1}{2}}\right )_c=\boldsymbol{\nabla }_{\! c} \boldsymbol{\cdot }\left [\mu _{\!f}^{n+\frac {1}{2}}\big (\boldsymbol{\nabla} \kern-1pt \boldsymbol{u}+\boldsymbol{\nabla} \kern-1pt \boldsymbol{u}^T\big )^*\right ]\nonumber\\&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad +\left [\left (\sigma \kappa \delta _s \boldsymbol{n}\right )^{n-\frac {1}{2}}- \boldsymbol{\nabla} \! p^{n-\frac {1}{2}}\right ]_{f \rightarrow c}, \end{align}
\begin{align} &\qquad\qquad\quad\boldsymbol{u}_c^{* *}=\boldsymbol{u}_c^*-\frac {\Delta t}{\rho _c^{n+\frac {1}{2}}}\left [\left (\sigma \kappa \delta _s \boldsymbol{n}\right )^{n-\frac {1}{2}}-\boldsymbol{\nabla} \! p^{n-\frac {1}{2}}\right ]_{f \rightarrow c}, \end{align}
\begin{align} &\qquad\qquad\qquad\quad \boldsymbol{u}_{\!f}^{*}=\boldsymbol{u}_{c \rightarrow f}^{* *} + \frac {\Delta t}{\rho _{\!f}^{n+\frac {1}{2}}}\left (\sigma \kappa \delta _s \boldsymbol{n}\right )^{n+\frac {1}{2}}, \end{align}
\begin{align} &\qquad\qquad\qquad\qquad\boldsymbol{u}_{\!f}^{n+1}=\boldsymbol{u}_{\!f}^{*} - \frac {\Delta t}{\rho _{\!f}^{n+\frac {1}{2}}}\boldsymbol{\nabla} \! p^{n+\frac {1}{2}}, \end{align}
\begin{align} &\qquad\qquad\,\,\, \boldsymbol{u}_c^{n+1}=\boldsymbol{u}_c^{* *}+\frac {\Delta t}{\rho _c^{n+\frac {1}{2}}}\left [\left (\sigma \kappa \delta _s \boldsymbol{n}\right )^{n+\frac {1}{2}}-\boldsymbol{\nabla} \! p^{n+\frac {1}{2}}\right ]_{f \rightarrow c}, \end{align}
where the superscripts
$n- ({1}/{2})$
,
$n+ ({1}/{2})$
,
$n$
and
$n+1$
represent the states at different time steps. The subscripts
$c$
and
$f$
denote the cell-centred and face-centred variables, respectively, and the conversion between them is achieved by a second-order accurate interpolation scheme (Zhao et al. Reference Zhao, Zhang and Ni2022), denoted by the symbol
$c \rightarrow f$
and
$f \rightarrow c$
. At each time step, the physical variables at
$n- ({1}/{2})$
and
$n$
are known, and their values at
$n + ({1}/{2})$
and
$n+1$
are obtained by solving the Navier–Stokes equations.
To solve the above-mentioned equations, the advection term is discretised using the Bell–Colella–Glaz (BCG) scheme (Bell, Colella & Glaz Reference Bell, Colella and Glaz1989), and the viscous term is discretised using the implicit Crank–Nicolson scheme. Moreover, the pressure at time
$n + ({1}/{2})$
in (2.21) is obtained by solving the Poisson equation,
\begin{equation} \boldsymbol{\nabla }_{\! c} \boldsymbol{\cdot }\left (\frac {\Delta t}{\rho ^{n+\frac {1}{2}}_{\!f}} \boldsymbol{\nabla} \! p^{n+\frac {1}{2}}\right )=\boldsymbol{\nabla }_{\! c} \boldsymbol{\cdot } \boldsymbol{u}_{\!f}^*-\boldsymbol{\nabla }_{\! c} \boldsymbol{\cdot } \boldsymbol{u}^{n+1}_{\!f}, \end{equation}
where the second term on the right-hand side is determined according to the mass (2.17). This is the so-called projection step, in which the intermediate velocity
$\boldsymbol{u}_{\!f}^*$
is projected onto a velocity field fulfilling the mass equation,
$\boldsymbol{u}_{\!f}^{n+1}$
. However, as
$\boldsymbol{u}_c^{n+1}$
is interpolated using (2.22), it is only approximately projected (Zhao et al. Reference Zhao, Zhang and Ni2022). The approximate projection method is adopted to facilitate the implementation of AMR, as the cell-centred velocity is difficult to project exactly due to the spatial decoupling of the stencils used for the relaxation operator (Popinet Reference Popinet2003). It has been shown that the approximate projection method introduces unphysical oscillations in the presence of phase change (Zhao et al. Reference Zhao, Zhang and Ni2022; Long et al. Reference Long, Pan and Zaleski2024). As shown in Appendix A, although the original model of Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024) performs well in many benchmark tests, it fails to correctly capture the physics of the microlayer during bubble growth. In the present work, a ghost fluid method is adopted and will be introduced later, enabling highly efficient and accurate simulations of boiling flows.
2.3.3. Solving the energy equation
To solve the mass and momentum equations in the presence of phase change, the mass flux
$\dot {m}$
is needed to determine the source term in (2.17). As the latent heat is usually given a priori, the mass flux can be calculated using (2.9) based on the difference between the heat fluxes across the interface. To compute the heat flux on both sides of the interface with a Dirichlet boundary condition (2.10) imposed, the normal temperature gradient
$\boldsymbol{\nabla} \kern-1pt T \boldsymbol{\cdot } \boldsymbol{n}$
is calculated using an embedded boundary method (Zhao et al. Reference Zhao, Zhang and Ni2022; Cipriano et al. Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024).
Schematic for the discretisation of the diffusion terms in the energy equation.
To resolve the temperature evolution, the energy equation in the fluid domain is solved using the finite volume method. In the phase-change model of Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024), two temperature fields, one for the liquid phase and one for the vapour phase, are solved separately using an operator splitting approach. This method handles the convective and diffusive parts of the energy equation independently (Cipriano et al. Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024). For brevity, the volume fractions of the liquid and vapour phases are defined as
$\theta _l = f$
and
$\theta _v = 1 - f$
. First, the advection term is solved as
\begin{equation} \left (\frac {(\theta _k T_k)^* - (\theta _k T_k)^{n-\frac {1}{2}}}{\Delta t}\right )_c = -\frac {1}{V_c} \sum _{\textit{cell}\ \textit{faces}}T_{k,f}^nF_{\theta _k,f}, \end{equation}
where
$k = l\ \text{or}\ v$
indicates the phase of interest,
$V_c$
denotes the cell volume,
$T^n_{k,f}$
represents the face-centred temperature obtained using the BCG scheme (Bell et al. Reference Bell, Colella and Glaz1989) and
$F_{\theta _k,f}$
is the volume flux computed during the geometrical advection of the VOF function (Cipriano et al. Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024). In this way, the energy is advected as a tracer associated with the VOF advection, avoiding numerical diffusion across the interface.
Once we obtain the intermediate temperature field, which takes into account the advection effect, we can solve the diffusion term using finite volume discretisation:
\begin{align} \left (\rho _k C_{p,k}\theta _k^{n+\frac {1}{2}}\frac {T_k^{n+\frac {1}{2}} - T_k^*}{\Delta t}\right )_c &= -\frac {1}{V_c} \left [\sum _{\textit{cell}\ \textit{faces}} \lambda _k \boldsymbol{\nabla} \kern-1pt T_{k,f}^{n+\frac {1}{2}}\boldsymbol{\cdot } \boldsymbol{S}_{k,f} \right . \nonumber\\&\quad \left .+ \left (\lambda _k \frac {\partial T}{\partial n }\bigg |_{\varGamma ,k} \boldsymbol{n} \boldsymbol{\cdot } \boldsymbol{S}_{\varGamma ,k} \right )^{n - \frac {1}{2}}\right ]\!, \end{align}
where
$\boldsymbol{S}_{k,f}$
and
$\boldsymbol{S}_{\varGamma ,k}$
denote the area vectors of the cell face and the interface segment, respectively, as illustrated in figure 3. The face-centred temperature gradient
$\boldsymbol{\nabla} \kern-1pt T_{k,f}$
is computed using a second-order accurate central difference method, and the normal temperature gradient at the interface
$({\partial T}/{\partial n })|_{\varGamma ,k}$
is calculated using the embedded boundary method (Zhao et al. Reference Zhao, Zhang and Ni2022; Cipriano et al. Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024). Note that the second term on the right-hand side of (2.25) is introduced only in the cells cut by the interface.
2.3.4. Ghost fluid method
In previous sections, the original phase-change model of Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024) is briefly reviewed. In this section, the modifications and extensions for simulating nucleate boiling with a resolved microlayer and conjugate heat transfer are introduced. As mentioned earlier, the source term in the mass equation induced by phase change (2.17) is singular and non-zero only in the interfacial cells. As a result, numerical oscillations may arise near the interface, especially for the cell-centred velocity, which is only approximately projected (Zhao et al. Reference Zhao, Zhang and Ni2022; Long et al. Reference Long, Pan and Zaleski2024). However, the approximate projection method is required for the implementation of AMR, which can significantly improve computational efficiency. To achieve a highly efficient and accurate phase-change model, the original model of Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024) is improved by adopting the ghost fluid method (Nguyen et al. Reference Nguyen, Fedkiw and Kang2001; Tanguy et al. Reference Tanguy, Sagan, Lalanne, Couderc and Colin2014). As shown in Appendix A, this modification is crucial for obtaining an accurate prediction of heat transfer within the microlayer. The principle of the ghost fluid method is to solve two separate velocity fields for each phase. Let
$C$
denote the colour function, which is computed by
\begin{equation} C = \left \{ \begin{aligned} &1 \quad \text{if $f \geqslant 0.5$}, \\ &0 \quad \text{if $f \lt 0.5$}. \\ \end{aligned} \right . \end{equation}
Accordingly, the velocity fields can be expressed as
\begin{equation} \begin{aligned} \boldsymbol{u}_l & = C\boldsymbol{u}_l + (1 - C) \boldsymbol{u}_l^g, \\ \boldsymbol{u}_v & = (1 - C)\boldsymbol{u}_v + C \boldsymbol{u}_v^g, \\ \end{aligned} \end{equation}
where the ghost velocities
$\boldsymbol{u}_l^g$
and
$\boldsymbol{u}_v^g$
are calculated as
\begin{equation} \begin{aligned} \boldsymbol{u}_{v}^{g} &= \boldsymbol{u}_{l} - \dot {m}\left [ \frac {1}{\rho } \right ]_\varGamma \boldsymbol{n}, \\ \boldsymbol{u}_{l}^{g} &= \boldsymbol{u}_{v} + \dot {m}\left [ \frac {1}{\rho } \right ]_\varGamma \boldsymbol{n}. \\ \end{aligned} \end{equation}
Since the above-mentioned equation is obtained according to the velocity jump condition (2.7), the singularity at the interface is removed (Nguyen et al. Reference Nguyen, Fedkiw and Kang2001; Tanguy et al. Reference Tanguy, Sagan, Lalanne, Couderc and Colin2014), leading to the new mass equations:
\begin{equation} \boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{u}_l = \boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{u}_v = 0. \end{equation}
Consequently, during the projection step, (2.23) becomes
\begin{equation} \boldsymbol{\nabla }_{\! c} \boldsymbol{\cdot } \left (\frac {\Delta t}{\rho ^{n+\frac {1}{2}}_{\!f}}\boldsymbol{\nabla} \! p^{n+1} \right ) = \left \{ \begin{aligned} &\boldsymbol{\nabla }_{\! c} \boldsymbol{\cdot } \boldsymbol{u}_{f,l}^* \quad \text{if $C = 1$}, \\ &\boldsymbol{\nabla }_{\! c} \boldsymbol{\cdot } \boldsymbol{u}_{f,v}^* \quad \text{if $C = 0$}, \\ \end{aligned} \right . \end{equation}
where the singular source term due to phase change on the right-hand side is eliminated. In practice, there is no need to populate ghost velocities throughout the domain. Equation (2.28) is applied only within a narrow band near the interface, defined by whether a cell or any of its neighbours is cut by the interface.
2.3.5. Treatment of solid
Finally, the solution of heat conduction in the solid domain and the implementation of conjugate heat transfer between the fluid and solid domains are elaborated. Applying the finite volume discretisation to (2.4) leads to
\begin{equation} \begin{aligned} \left (\rho _s C_{\kern-1.5pt p,s}\frac {T_s^{n+\frac {1}{2}} - T_s^{n-\frac {1}{2}}}{\Delta t}\right )_c = -\frac {1}{V_c} \left (\sum _{\textit{cell}\ \textit{faces}} \lambda _s \boldsymbol{\nabla} \kern-1pt T_{s,f}^{n+\frac {1}{2}}\boldsymbol{\cdot } \boldsymbol{S}_{s,f} \right ) + Q_h^{n-\frac {1}{2}}. \\ \end{aligned} \end{equation}
It can be seen from (2.25) and (2.31) that the heat conduction terms are discretised using an implicit scheme in both the fluid and solid domains, eliminating the strict time step constraint imposed by these terms. During simulations, these two equations are solved simultaneously with the associated boundary conditions. For fluid and solid cells near the boundary, the discretisation stencil is incomplete when the face-centred temperature gradient is computed with a central scheme. As shown in figure 4, this issue is addressed by employing a one-sided difference scheme, relying on two associated boundary temperatures,
$T_{b,s}$
and
$T_{b,f}$
. With the contact heat-transfer resistance and the continuity of heat flux, the following equation is obtained:
\begin{equation} 2\lambda _s \frac {T_s - T_{b,s}}{\varDelta } = \frac {T_{b,s} - T_{b,f}}{R_c}= 2\lambda _{\!f} \frac {T_{b,f} - T_{\!f}}{\varDelta }, \end{equation}
where
$\varDelta$
is the grid size. This equation can be simplified by introducing
$R_s = \varDelta / (2 \lambda _s)$
and
$R_{\!f} = \varDelta / (2 \lambda _{\!f})$
. The solution of (2.32) yields the boundary values as follows:
\begin{equation} \begin{aligned} T_{b,s} = & \frac {(R_c + R_{\!f})T_s + R_s T_{\!f}}{R_c + R_{\!f} + R_s},\\ T_{b,f} = &\frac {(R_c + R_s)T_{\!f} + R_{\!f} T_s }{R_c + R_{\!f} + R_s}. \end{aligned} \end{equation}
The temperature gradient at the fluid–solid boundary can then be calculated using (2.32). Note that the continuity of temperature (Huber et al. Reference Huber, Tanguy, Sagan and Colin2017) across the fluid–solid boundary is automatically recovered when
$R_c = 0$
. In addition, all the above-mentioned procedures can be directly applied when the solid consists of two different materials. Following Bureš & Sato (Reference Bureš and Sato2022), the properties of the solid are mixed based on the solid volume fractions. In this study, binary mixed solids are considered, and the cell-specific values of
$\lambda _s$
and
$C_{\kern-1.5pt p,s}$
are computed by
\begin{equation} \begin{aligned} \lambda _s &= \lambda _{1}f_{s,1} + \lambda _{2}(1 - f_{s,1}),\\ C_{\kern-1.5pt p,s} &= C_{\! p,1}f_{s,1} + C_{\! p,2}(1 - f_{s,1}),\\ \end{aligned} \end{equation}
where
$f_{s,1}$
represent the volume fraction of material 1.
Schematic for the implicit discretisation of the heat diffusion terms at the fluid–solid boundary.
2.3.6. Summary
For clarity, the numerical procedures for each time step are summarised as follows.
-
(i) Solve the advection of the VOF function using the split scheme (2.14).
-
(ii) Solve the advection term of the energy equation in the fluid domain (2.24).
-
(iii) Solve heat conduction in the fluid and solid domains simultaneously using an implicitly coupled approach ((2.25) and (2.31)).
-
(iv) Compute the mass flux (2.9) and set the ghost velocities (2.28).
3. Verification
The phase-change model of Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024) has been verified for the liquid–vapour system using a wide range of benchmark tests, with all codes available on the Basilisk website (Cipriano Reference Cipriano2023). Here, we focus primarily on the validation of the conjugate heat transfer extension. In this section, we consider a film evaporating in the presence of conjugate heat transfer, as illustrated in figure 5(a). This benchmark case, proposed by Burěs (Reference Burěs2021), is specifically designed for code verification. The theoretical solution can be derived using the method of manufactured solutions (MMS) (Roache Reference Roache1998). Before detailing the case set-up, some notation is introduced for clarity. In the following derivations,
$\phi '(x,t)$
and
$\dot {\phi }(x,t)$
denote the derivatives of the variable
$\phi$
with respect to
$x$
and
$t$
, respectively. Additional primes and dots indicate corresponding higher-order derivatives.
Film evaporation with conjugate heat transfer. (a) Schematic of the one-dimensional film evaporation with conjugate heat transfer. (b) Temperature distribution at
$ t = 0.5$
s. (c) Time history of the interface position. (d) Relative error of the interface position on different grid resolutions. Grid levels 6–8 correspond to effective grid resolutions ranging from
$ 64 \times 1$
cells to
$ 256 \times 1$
cells, resulting in minimum grid sizes from
$0.05$
to
$0.0125$
.
In this case, the interface position is specified as
\begin{equation} x_{\varGamma }(t) = x_0 - Z t^2, \end{equation}
where
$x_0$
is the initial position and
$Z$
is a user-defined parameter. The vapour temperature is maintained at the saturation temperature
$T_{\textit{sat}}$
throughout the process, while the liquid temperature is modelled with an exponential profile:
\begin{equation} T_l(x, t)=T_{\textit{sat}} E(x, t)=T_{\textit{sat}} \exp \left [M \dot {x}_\varGamma (t)\left (x-x_\varGamma (t)\right )\right]\!. \end{equation}
Here, the factor
$M$
is defined as
\begin{equation} M = \frac {h_{lv} \rho _l}{T_{\textit{sat}}\lambda _l}, \end{equation}
ensuring the Stefan condition
\begin{equation} -\frac {\lambda _l }{h_{lv} \rho _l}T_l(x_\varGamma , t) = -\dot {x}_\varGamma (t) \end{equation}
is satisfied. To enforce the flux and temperature jump conditions at the fluid–solid boundary, the solid temperature distribution is given by
\begin{equation} T_{s}(x, t)=T_s\left [\left (\frac {\lambda _l}{\lambda _{s}}+D(t)\right ) E(x, t)+C(t) E(0, t)\right ]\!, \end{equation}
where
$D(t)$
and
$C(t)$
are defined as
\begin{equation} D(t)=\frac {\delta _q(t)}{\lambda _{s} T_s} \frac {1}{E^{\prime }(0, t)} \end{equation}
and
\begin{equation} C(t)=-\left [\lambda _l R_{c} M \dot {x}_\varGamma (t)+D(t)+\left (\frac {\lambda _l}{\lambda _{s}}-1\right )\right ]\!, \end{equation}
respectively. Here, the wall source term
$\delta _q(t)$
is chosen as
$\delta _q(t)=-\psi \dot {x}_\varGamma (t)$
, with
$\psi$
being a given control parameter. To achieve the temperature evolution described here, additional source terms must be imposed in both the liquid and solid domains during the simulation, which are given by
\begin{align} S_l(x, t) &= C_{\kern-1.5pt p, l} \dot {T}_l(x, t)-\lambda _l T_l^{\prime \prime }(x, t)=T_s\left [C_{\kern-1.5pt p, l} \dot {E}(x, t)-\lambda _l E^{\prime \prime }(x, t)\right]\!, \end{align}
\begin{align} S_{s}(x, t) &= C_{\kern-1.5pt p,s} \dot {T}_s(x, t)-\lambda _{s} T_s^{\prime \prime }(x, t) \nonumber \\ &= T_s\left \{C_{\kern-1.5pt p,s}\left [\left (\frac {\lambda _l}{\lambda _{s}}+D(t)\right ) \dot {E}(x, t)+\dot {D}(t) E(x, t)+\dot {C}(t) E(0, t)+C(t) \dot {E}(0, t)\right ] \right . \nonumber \\ &\quad\left .-\lambda _{s}\left (\frac {\lambda _l}{\lambda _{s}}+D(t)\right ) E^{\prime \prime }(x, t)\right \} \! . \end{align}
In the simulation, the fluid–solid boundary is placed at the origin, with the lengths of the fluid and solid domains chosen as
$2.2$
and
$1$
, respectively. Following Burěs (Reference Burěs2021), with viscosity neglected, the physical properties are set as
\begin{equation} \left \{\begin{aligned} \rho _{l}&=1, \lambda _{l}=1,C_{\kern-1pt p,l} =2, \rho _{v}=1, \lambda _{v}=1,C_{\kern-1.5pt p,v} =2,\\ \rho _{s}&=4, \lambda _{s}=7,C_{\kern-1.5pt p,s} =5, T_{s}=1, h_{lv}=1, R_c = 2.3,\\ x_0&=1.03, Z=1, \psi = 1.65.\\ \end{aligned}\right . \end{equation}
Note that all the above-mentioned parameters are dimensionless. The simulation is performed up to
$t = 0.5$
with three grid levels ranging from
$6$
to
$8$
, resulting in grid sizes from
$0.05$
to
$0.0125$
. Note that for a given grid level
$L$
, the corresponding number of cells is
$N = 2^{L}$
. In figure 5(b), the temperature distributions at the final time for different grid refinements are compared with the theoretical solution. It can be observed that the temperature discontinuity at the fluid–solid boundary is accurately and sharply captured by the present method. Moreover, as shown in figure 5(c), the time evolutions of the interface positions for all grid levels are in good agreement with the theoretical solution. For a quantitative comparison, the relative error of the final interface position is computed by
\begin{equation} E(x_\varGamma (0.5)) = \frac {|x_\varGamma (0.5) - x_\varGamma (0.5)^{\textit{num}}|}{x_\varGamma (0.5)}, \end{equation}
where the superscript
$\textit{num}$
represents the numerical solution. As can be observed from figure 5(d), the present method exhibits a second-order convergence rate, demonstrating its efficacy for phase-change problems involving conjugate heat transfer.
4. Results and discussion
After validating the phase-change model with conjugate heat transfer, it is applied to simulate the experiments of Bucci (Reference Bucci2020). In these experiments, the bubble growth in pool boiling of water under atmospheric conditions was investigated. Note that all codes used in this study are available from the Basilisk website (Long Reference Long2024).
4.1. Simulation set-up
As mentioned, the IHTR is considered in the simulations by imposing a temperature contact discontinuity at the fluid–solid interface. The contact heat-transfer resistance factor
$R_c$
is set equal to the IHTR factor
$R_\varGamma$
to preserve the overall heat flux. In the work of Bureš & Sato (Reference Bureš and Sato2022), the Hertz–Knudsen relation is combined with the Clausius–Clapeyron relation to model
$R_\varGamma$
, leading to the following expression:
\begin{equation} R_\varGamma = \frac {1}{\omega }\frac {1}{\rho _v h_{lv}^2}\sqrt {\frac {2\pi R_g T_{\textit{sat}}^3}{M_v}}, \end{equation}
where
$M_v$
is the molar mass,
$R_g$
is the universal gas constant and
$\omega$
is the so-called accommodation coefficient, which is a priori unknown. By fitting experimental data, Bureš & Sato (Reference Bureš and Sato2022) evaluated the bounds of
$\omega$
and selected two values for their simulations:
$\omega = 0.0345$
and
$\omega = 0.0460$
. In the present work, the same two accommodation coefficients are adopted, referring to the cases with
$\omega = 0.0345$
and
$\omega = 0.0460$
as case A and case B, respectively.
Physical properties used in the simulations of the experiment of Bucci (Reference Bucci2020).
Schematic of the computational domain used for the simulations of the experiment of Bucci (Reference Bucci2020) (not to scale).
In addition to accommodation coefficients, the other physical properties considered in the simulations are listed in table 1. The simulations are performed with a two-dimensional (2-D) axisymmetric configuration, which accurately represents the real experimental conditions, where perfect axial symmetry of the growing bubble has been observed (Bucci Reference Bucci2020). As shown in figure 6, a rectangular domain of
$ [ 0\ \mathrm{mm}, 1.3 \ \mathrm{mm} ] \times [ -1.0\ \mathrm{mm}, 1.3 \ \mathrm{mm} ]$
is employed, with the fluid–solid boundary placed at
$z = 0\ \mathrm{mm}$
. The solid phase consists of a
$1\ \mathrm{mm}$
-thick sapphire substrate with a
$500\ \mathrm{nm}$
-thick titanium heater. The heater is modelled using a volumetric source term (
$Q_h$
in (2.4)), and, following Bureš & Sato (Reference Bureš and Sato2022), we consider the case with an applied heat flux of
$425\ \mathrm{kW}\,\mathrm{m}^{-2}$
. With the left boundary as the axis of symmetry, outflow boundary conditions are applied to the right and top boundaries of the fluid domain, while the fluid–solid boundary is treated as a no-slip wall. Note that since the VOF function is advected using the discretised velocity located half a grid spacing above the wall, an implicit slip condition is introduced for the interface motion (Afkhami & Bussmann Reference Afkhami and Bussmann2008). Generally, this implicit slip length of half a grid spacing can lead to mesh-dependent contact line behaviour (Afkhami, Zaleski & Bussmann Reference Afkhami, Zaleski and Bussmann2009) for which more advanced contact line models are typically required. However, several studies (Han et al. Reference Han, Zhang, Tan and Ni2021; Huang et al. Reference Huang, Han, Zhang and Ni2025) have shown that, in certain contexts, accurate contact line dynamics can still be captured using this approach, provided that the grid resolution is sufficiently high. In the present study, as shown in figure 13, we observe the convergence of both contact line mobility and the microlayer profile with increasing mesh resolution, indicating that this simplified treatment is adequate for the nucleate boiling problems considered. In the solid domain, a Neumann boundary condition (zero heat flux) is imposed on the right and bottom boundaries. At the fluid–solid boundary, a three-phase contact line forms, requiring the specification of a contact angle
$\theta _C$
. In the work of Bureš & Sato (Reference Bureš and Sato2022), a dynamic contact angle model is adopted, yielding an angle of less than
$1^\circ$
throughout the simulations. However, as will be shown in this paper, the influence of the contact angle is minimal when it is less than
$10^\circ$
. In this section, following El Mellas et al. (Reference El Mellas, Samkhaniani, Falsetti, Stroh, Icardi and Magnini2024), a static contact angle model with
$\theta _C = 5^\circ$
is used. The static contact angle is imposed using the height function approach of Afkhami & Bussmann (Reference Afkhami and Bussmann2008), which is already implemented in Basilisk.
With the above-mentioned set-up, simulations have been conducted at increasing grid levels from
$10$
to
$12$
, with corresponding minimum grid sizes ranging from
$2.246$
to
$0.562\ \mathrm{\unicode{x03BC} m}$
. The number of CPU cores and the total CPU time are compared with those reported by Bureš & Sato (Reference Bureš and Sato2022) in table 2. It is observed that the number of CPU cores and CPU hours required for the simulations are significantly reduced using the present model. For the finest resolution, only approximately 680 core-hours are needed with 48 cores, compared with approximately
$400\ 000$
core-hours with 336 cores reported by Bureš & Sato (Reference Bureš and Sato2022). The significant improvement in computational efficiency is attributed to two factors. First, for the same effective resolution, the use of AMR greatly reduces the total number of grid cells in the simulation. Second, the present model is more robust and stable, allowing for a larger time step. During the simulation, as required by numerical stability (Bureš & Sato Reference Bureš and Sato2022), the time step
$\Delta t$
is determined by
\begin{equation} \Delta t = \mathrm{min}\left (C_{\textit{ad}v}\frac {\varDelta }{u_{\textit{max}}}, C_\sigma \sqrt {\frac {(\rho _v + \rho _l)\varDelta ^ 3}{\sigma }}\right)\!, \end{equation}
where
$\varDelta$
is the minimum grid size used and
$u_{\textit{max}}$
is the maximum velocity component in all directions. Here,
$C_{\textit{ad}v}$
and
$C_\sigma$
are the limiting coefficients for the advection and capillary terms, which are set to
$C_{\textit{ad}v} = 0.5$
and
$C_\sigma = 0.282$
in the present study. In the work of Bureš & Sato (Reference Bureš and Sato2022), the time step limits are much more stringent (
$C_{\textit{ad}v} = 0.02$
and
$C_{\sigma } = 0.063$
) due to the oscillations induced by the strong heat transfer involved in this problem.
Comparison of computational efficiency between the present work and the study of Bureš & Sato (Reference Bureš and Sato2022).
4.2. Initial condition
In the experiment, the system starts from saturated and stagnant conditions. After a period of heating, the first bubble nucleates in a cavity on the heated surface and then grows rapidly. Within the sharp-interface framework, an initial bubble seed is required to trigger the phase-change process. To model the nucleation process, following Bureš & Sato (Reference Bureš and Sato2022), a transient heating problem without the vapour phase is solved until the temperature at the nucleation site reaches the nucleation temperature of
$112.55\,^\circ$
C. Given the nucleation superheat of
$\Delta T = 12.55\ \mathrm{K}$
, this problem can be characterised by the Jakob number:
\begin{equation} \textit{Ja} = \frac {\rho _{l}C_{\kern-1pt p,l} \Delta T}{\rho _{v} h_{lg}}, \end{equation}
yielding a value of 37.5. Note that although the bubble growth results reported by Bucci (Reference Bucci2020) are axisymmetric, the heater geometry is rectangular. As a result, a uniform heat flux distribution under the axisymmetric numerical configuration leads to deviations in the calculated surface temperature distribution compared with experimental measurements. To address this issue, Bureš & Sato (Reference Bureš and Sato2022) adopted a modified heat flux distribution for the initial temperature field calculation, given as
\begin{equation} j_{\textit{ini}} = \mathrm{max}\left (\frac {4 - r}{4}\times 481, 0\right ) (\mathrm{kW}\,\mathrm{m}^{-2}). \end{equation}
This linear power input is used to solve the liquid–solid heat transfer problem. To compute the initial temperature field, the domain is extended to
$ [ 0\ \mathrm{mm}, 7.0 \ \mathrm{mm} ] \times [ -1.0\ \mathrm{mm}, 6.0 \ \mathrm{mm} ]$
. The simulations are performed with three grid levels ranging from 10 to 12, ensuring grid convergence (not shown here). In figure 7, the superheat distribution on the solid surface at different time instants, obtained at grid level 12, is compared with experimental data (Bucci Reference Bucci2020) and the numerical results of Bureš & Sato (Reference Bureš and Sato2022), showing excellent agreement. The waiting time (the time from the beginning of heating to nucleation) in the experiment (
$88.1\ \mathrm{ms}$
) is also well reproduced numerically. Since the macroscopic sharp-interface method is employed in the present study, the actual nucleation process can only be taken into account indirectly, i.e. by artificially introducing a bubble embryo at the nucleation site (Bureš & Sato Reference Bureš and Sato2022) when its temperature reaches a prescribed nucleation temperature. The same strategy has been adopted by Cai et al. (Reference Cai, Wang, Jiang and Tan2024) to simulate the experiment of Jung & Kim (Reference Jung and Kim2014), where multiple bubble cycles were successfully replicated, demonstrating the effectiveness of this simple approach. More advanced models could be introduced to compute the nucleation time by considering liquid metastability and the heat release rate. However, since our objective is to replicate the bubble growth stage in the experiment of Bucci (Reference Bucci2020) and the nucleation process is not the main focus of this study, the incorporation of these models is left for future work. It should also be emphasised that accurately simulating the nucleation process is challenging even with more detailed microscopic methods, such as molecular dynamics (MD) (Zhang, Xu & Lei Reference Zhang, Xu and Lei2018), phase field (PF) Zhang et al. (Reference Zhang, Li, Liang and Shu2025) or the lattice Boltzmann method (LBM) (Mu et al. Reference Mu, Chen, He, Kang and Tao2017), which can simulate nucleation without artificially introducing a bubble embryo. To the authors’ knowledge, the nucleation temperatures obtained using these methods still depend on the potential models employed and therefore require careful validation against experimental results (Chen et al. Reference Chen, Yu, Lu, Wang, Sun, Jiao, Zhang and Tao2024).
Evolution of the heater surface superheat before the onset of nucleation. The results of the present study are compared with the experimental data of Bucci (Reference Bucci2020) and the numerical results of Bureš & Sato (Reference Bureš and Sato2022).
4.3. Numerical results
Using the temperature distribution at
$t = 88.1\ \mathrm{ms}$
as the initial condition, we position a bubble seed with a radius of
$R_0 = 20\ \mathrm{\unicode{x03BC} m}$
at the origin to initiate bubble growth. The simulations are performed up to
$t = 0.5\ \mathrm{ms}$
, during which a uniform heat flux of
$425\ \mathrm{kW}\,\mathrm{m}^{-2}$
is applied to the heater (Bureš & Sato Reference Bureš and Sato2022). In this section, the numerical results obtained by the present model are validated against experimental data and previous numerical studies. Before the quantitative comparison, the superheat distributions within the computational domain at two time instants obtained at grid level 12 are shown in figure 8 for case A (
$\omega = 0.0345$
). As in previous studies (Bureš & Sato Reference Bureš and Sato2022; Torres et al. Reference Torres, Urbano, Colin and Tanguy2024), the thermal boundary layer near the bubble surface becomes progressively thinner as the bubble grows. Additionally, the solid surface beneath the bubble is cooled due to vigorous evaporation within the microlayer. This phenomenon is more clearly observed in the 2-D distribution of the superheat on the solid surface in figure 9, which is generated by rotating the axisymmetric data. Notably, the highest superheat on the solid surface is located at the outer edge of the microlayer rather than at the lateral edge of the bubble.
Superheat distributions at (a)
$t = 0.2\ \mathrm{ms}$
and (b)
$t = 0.4\ \mathrm{ms}$
for case A (
$\omega = 0.0345$
), obtained at grid level 12. The white dashed line and the black solid line represent the liquid–vapour interface and the fluid–solid boundary, respectively. The results are mirrored about
$r = 0$
for better visualisation.
Three-dimensional (3-D) bubble interfaces and the superheat distributions on the solid surface at different time instants for case A (
$\omega = 0.0345$
), obtained at grid level 12. The plots are generated by rotating the axisymmetric results.
We then turn to a quantitative comparison of growth characteristics. In figures 10 and 11, the volume and lateral radius of the bubble are plotted against time for all grid resolutions, demonstrating grid convergence. It should be noted that grid convergence was not achieved in the results reported by Bureš & Sato (Reference Bureš and Sato2022). This peculiar convergence behaviour is probably caused by the one-fluid method used in their simulations, which, as shown in the Appendix, can lead to numerical oscillations within the microlayer due to strong heat transfer. In general, the convergent results obtained with the present method align well with the experimental data of Bucci (Reference Bucci2020) and the numerical results of Bureš & Sato (Reference Bureš and Sato2022) at their finest grid resolution with the smallest cell size of
$0.475\ \mathrm{\unicode{x03BC} m}$
. For case A, the bubble volume and lateral radius are slightly overestimated before 0.3 ms and underestimated thereafter, compared with the experimental results. In contrast, for case B, the results show improved agreement in later stages. This is expected due to the higher accommodation coefficient in case B, which leads to lower IHTR and thus larger heat flux within the microlayer for the same superheat on the solid surface, according to (4.1). The same trend is also observed in the results of Bureš & Sato (Reference Bureš and Sato2022). In their simulations, for case A, the bubble volume and lateral radius are underestimated after 0.2 ms, whereas for case B, better agreement is observed in later stages. It is evident from figures 10 and 11 that the accommodation coefficient significantly influences the bubble growth. The discrepancies between the convergent numerical results and the experimental data may be attributed to the constant accommodation coefficient used throughout the simulations. Constant accommodation coefficients are commonly used in numerical simulations of microlayers (Bureš & Sato Reference Bureš and Sato2022; Hänsch & Walker Reference Hänsch and Walker2019; El Mellas et al. Reference El Mellas, Samkhaniani, Falsetti, Stroh, Icardi and Magnini2024). Regarding experimental studies, Giustini et al. (Reference Giustini, Jung, Kim and Walker2016) also reported a constant accommodation coefficient of 0.02 based on their analysis of evaporation kinetics and data. However, it is known that the accommodation coefficient can be influenced by various factors, such as vapour–liquid relative velocity (Nie et al. Reference Nie, Chandra, Liang and Keblinski2019), local pressure (Marek & Straub Reference Marek and Straub2001) and temperature (Vaartstra et al. Reference Vaartstra, Lu, Lienhard and Wang2022), and can vary significantly from 0.01 to 1. In the work of Cai et al. (Reference Cai, Wang, Jiang and Tan2024), the experiment of Jung & Kim (Reference Jung and Kim2014) is successfully replicated using a more advanced IHTR model, where the accommodation coefficient is determined by the free energy and the local interfacial temperature. Depending on the local conditions, the accommodation varies in both space and time, with the time-varying feature being the dominant effect observed. It is reasonable to deduce that a more advanced model for predicting the accommodation coefficient instantaneously, based on local conditions, can help accurately simulate real situations. However, to our knowledge, while several theories exist (Nathanson et al. Reference Nathanson, Davidovits, Worsnop and Kolb1996; Persad & Ward Reference Persad and Ward2016), accurately modelling the accommodation coefficient remains an open question in the literature. Moreover, it should also be emphasised that although the Schrage model, based on the accommodation coefficient (Schrage Reference Schrage1953), is chosen by many researchers (Giustini et al. Reference Giustini, Jung, Kim and Walker2016; Hänsch & Walker Reference Hänsch and Walker2019; Bureš & Sato Reference Bureš and Sato2022) to model the non-equilibrium effects at the fluid interface, it might not be the most appropriate model for computing the phase change mass flow rate in the microlayer. The assessment of the performance of other models (Liu et al. Reference Liu, Tang, Sun, Mo and Xie2020), such as the Lee model (Lee Reference Lee1980) and the Tanasawa model (Tanasawa Reference Tanasawa1991), within the same numerical framework will be left for future work.
Time histories of the bubble volume for (a) case A (
$\omega = 0.0345$
) and (b) case B (
$\omega = 0.0460$
) at different grid levels, compared with the experimental data of Bucci (Reference Bucci2020) and the numerical results of Bureš & Sato (Reference Bureš and Sato2022). Note that the results for case B are not reported in the work of Bureš & Sato (Reference Bureš and Sato2022).
Time histories of the bubble lateral radius for (a) case A (
$\omega = 0.0345$
) and (b) case B (
$\omega = 0.0460$
) at different grid levels, compared with the experimental data of Bucci (Reference Bucci2020) and the numerical results of Bureš & Sato (Reference Bureš and Sato2022).
Bubble interfaces at
$t = 0.31\ \mathrm{ms}$
for case A (
$\omega = 0.0345$
) at different grid levels, compared with the experimental data of Bucci (Reference Bucci2020) and the numerical results of Bureš & Sato (Reference Bureš and Sato2022). The results are mirrored about
$r = 0$
for better visualisation.
A more detailed comparison of the bubble shape at
$t = 0.31\ \mathrm{ms}$
is presented in figure 12. For case A, the bubble interfaces at different grid levels are compared with the experimental data of Bucci (Reference Bucci2020) and the numerical result of Bureš & Sato (Reference Bureš and Sato2022). In our results, the interfaces at grid levels 11 and 12 nearly overlap and agree well with the experimental data, demonstrating good grid convergence. Moreover, our results show a smoother bubble shape compared with the result of Bureš & Sato (Reference Bureš and Sato2022), which is attributed to the reduced numerical oscillations achieved by the ghost fluid method. The bubble shape obtained by Bureš & Sato (Reference Bureš and Sato2022) is more flattened compared with both our results and the experimental data. This explains why their results are in good agreement with the experiment regarding the lateral radius, while the bubble volume is significantly underestimated. In addition to macroscopic bubble shapes, the microlayer profiles for case A at different time instants and grid levels are illustrated in figure 13. The microlayer profile can be measured using laser interferometry (Jung & Kim Reference Jung and Kim2018; Narayan & Srivastava Reference Narayan and Srivastava2021), though this technique was not used in Bucci’s experiments (Bucci Reference Bucci2020). Thus, we only compare our results with the numerical results of Bureš & Sato (Reference Bureš and Sato2022), showing good agreement. The step-like profiles observed are numerical artefacts associated with the finite grid resolution, as also noted in previous studies (Bureš & Sato Reference Bureš and Sato2021, Reference Bureš and Sato2022). The result at grid level 10 for
$t = 0.21\ \mathrm{ms}$
exhibits oscillations due to the coarse grid size (
$2.246\ \mathrm{\unicode{x03BC} m}$
), which is insufficient to resolve the heat transfer within the microlayer. It is emphasised that the current set-up is numerically challenging due to the involvement of multiple spatial scales and intense mass and heat transfer. In the work of Bureš & Sato (Reference Bureš and Sato2022), strong spurious waves along the interface were reported. To prevent simulation crashes, they imposed very strict time step limits (
$C_{\textit{ad}v} = 0.02$
and
$C_{\sigma } = 0.063$
) and applied an artificial averaging procedure for the phase-change rate during the initial stage of bubble growth. In the present method, numerical stability is improved by the ghost fluid method, allowing the use of a larger time step and eliminating the artificial averaging procedure.
Microlayer profiles at different time instants and grid levels for case A (
$\omega = 0.0345$
), compared with the numerical results of Bureš & Sato (Reference Bureš and Sato2022).
Surface superheat distribution at different time instants and grid levels, compared with the experimental data of Bucci (Reference Bucci2020) and the numerical results of Bureš & Sato (Reference Bureš and Sato2022). (a,b) Results of case A (
$\omega = 0.0345$
); (c,d) results of case B (
$\omega = 0.0460$
).
Surface heat flux distribution at different time instants and grid levels, compared with the experimental data of Bucci (Reference Bucci2020) and the numerical results of Bureš & Sato (Reference Bureš and Sato2022). (a,b) Results of case A (
$\omega = 0.0345$
); (c,d) results of case B (
$\omega = 0.0460$
).
Finally, the thermodynamic characteristics for cases A and B are quantitatively validated. The surface superheat distributions at various time instants and grid levels are presented in figure 14. It is observed that as the grid resolution increases, our results progressively converge to the experimental data of Bucci (Reference Bucci2020), showing excellent agreement with the numerical results of Bureš & Sato (Reference Bureš and Sato2022). The extent of the microlayer can be identified from surface temperature variations: the temperature decreases from the origin to the contact line, increases along the microlayer and then declines beyond the microlayer front. These results demonstrate that our method effectively captures the evolution of the microlayer. As in the work of Bureš & Sato (Reference Bureš and Sato2022), larger deviations between numerical and experimental results are observed within the dry patch region (from the origin to the contact line). Note that detailed measurement uncertainties are not reported by Bucci (Reference Bucci2020), and this discrepancy may be attributed to less accurate measurements on surfaces not covered by liquid (Bureš & Sato Reference Bureš and Sato2022). Furthermore, the heat flux distributions at different times and resolutions are also measured and shown in figure 15. It can be seen that the heat flux peaks at the contact line and then decreases to zero along the radial extent of the microlayer. Notably, the heat flux within the microlayer region significantly exceeds the electrical power input. Before the bubble nucleation, the energy released from electrical resistance is distributed between both the liquid and solid. Given that the thermal diffusivity of the solid is considerably higher than that of the liquid, a greater proportion of energy is transferred to the solid. The subsequent cooling of the liquid due to boiling induces energy release from the solid, resulting in a substantially higher heat flux on the solid surface. This highlights the importance of including conjugate heat transfer between the fluid and solid in pool boiling simulations. It can be concluded from figure 15 that the experimental measurements are well bounded by the simulation results of cases A and B, as observed by Bureš & Sato (Reference Bureš and Sato2022). With a smaller accommodation coefficient (case A), better agreement is achieved in the early stages (
$t = 0.21\ \mathrm{ms}$
), while a larger accommodation coefficient (case B) yields improved results in the later stages (
$t = 0.42\ \mathrm{ms}$
). In future work, we plan to improve the current model by incorporating more advanced models to predict the accommodation coefficient based on local conditions, as demonstrated by the work of Cai et al. (Reference Cai, Wang, Jiang and Tan2024).
4.4. Effect of contact angle
In the work of Bureš & Sato (Reference Bureš and Sato2022), a dynamic contact angle model was proposed, where the predicted values of the contact angle throughout the simulation did not exceed
$1^\circ$
. As shown in § 4.3, our numerical results obtained with a static contact angle of
$5^\circ$
are in very good agreement with those obtained using the dynamic contact angle model (Bureš & Sato Reference Bureš and Sato2022). To further explore the influence of the contact angle, the same problem is simulated with varying contact angles
$\theta _{C}$
. Since it has been validated that grid level 11 is sufficient to capture the main physics of this problem, all simulations are performed at this grid level, with an accommodation coefficient of
$\omega = 0.0460$
. Seven contact angles, equally distributed within
$[5^\circ , 35^\circ ]$
, are considered.
Time histories of (a) the bubble volume, (b) the bubble lateral radius, (c) the position of the contact line and (d) the microlayer length for various contact angles. The results are obtained with
$\omega = 0.0460$
at grid level 11.
Microlayer profiles at different time instants for various contact angles. The results are obtained with
$\omega = 0.0460$
at grid level 11.
Schematic of region partitioning. Region 1 extends from the contact line to the microlayer front. Region 2 spans from the microlayer front to the bubble nose, where the maximum width occurs. Region 3 extends from the bubble nose to the bubble apex.
We begin with a quantitative analysis of growth characteristics, focusing on the effect of the contact angle on microlayer development. In addition to the volume and lateral radius of the bubble, the contact line position
$x_{\textit{CL}}$
and the microlayer length
$l_{\textit{ML}}$
are also measured. The microlayer length is determined as the distance from the contact line to the microlayer front
$x_{\textit{MF}}$
(Urbano et al. Reference Urbano, Tanguy, Huber and Colin2018), which is defined as the position where the interface slope exceeds
$5^\circ$
(see figure 18). The results for various contact angles are presented in figure 16. It is demonstrated that the influence of the contact angle on bubble growth diminishes as the angle decreases, particularly when it is less than
$15^\circ$
. As illustrated in figure 16(c), after a short transition stage, the contact line moves at a constant velocity that is positively correlated with the contact angle. In particular, the results for
$\theta _{C} = 35^\circ$
differ from others in the later stages due to the complete depletion of the microlayer, as indicated in figure 16(d). Furthermore, a decrease in the growth velocity of the bubble lateral radius is observed over time, which corresponds to the bubble growth within the diffusion-controlled regime (Guion et al. Reference Guion, Afkhami, Zaleski and Buongiorno2018). As the bubble growth velocity decreases, the contact line velocity remains constant and the extent of the microlayer is influenced by the competition between these two velocities. Simultaneously, evaporation contributes to microlayer depletion. At smaller contact angles, evaporation effects become more pronounced due to reduced contact-line mobility (Bureš & Sato Reference Bureš and Sato2022), whereas dryout driven by hydrodynamic effects is more significant at larger contact angles (Urbano et al. Reference Urbano, Tanguy, Huber and Colin2018). The microlayer profiles at various time instants for different contact angles are presented in figure 17. The transition from the microlayer regime to the contact line regime is observed for
$\theta _{C} = 35^\circ$
. Notably, at large contact angles, a dewetting ridge forms near the contact line, with its height increasing as the contact angle increases. This phenomenon results from mass conservation in the presence of enhanced contact-line mobility associated with larger contact angles (Giustini Reference Giustini2024). Additionally, the central portions of the microlayer overlap for different contact angles, consistent with findings from previous studies (Guion et al. Reference Guion, Afkhami, Zaleski and Buongiorno2018; Giustini Reference Giustini2024). In conclusion, different contact angles primarily influence interface motion near the contact line.
Time histories of the evaporation rate over (a) region 1, (b) region 2, (c) region 3 and (d) the entire bubble for various contact angles. The results are obtained with
$\omega = 0.0460$
at level 11.
Subsequently, to thoroughly investigate the influence of the contact angle on thermodynamic characteristics, the bubble interface is divided into three regions, as illustrated in figure 18. Region 1, the microlayer region, extends from the contact line to the microlayer front. Region 2 spans from the microlayer front to the bubble nose, which is the farthest point along the radial direction. The remaining area, from the bubble nose to the bubble apex, is labelled as Region 3. It is important to note that the following analyses primarily focus on bubble growth in the microlayer regime. Consequently, the case of
$\theta _C = 35^\circ$
is excluded, as it involves the transition from the microlayer regime to the contact line regime.
We begin the quantitative comparison by evaluating the evaporation rate, which is calculated as follows:
\begin{equation} r_e = \iint \dot {m} \,\text{d}A. \end{equation}
The time histories of the evaporation rate within different regions are presented in figure 19. It is observed that, generally, larger contact angles result in smaller evaporation rates, a trend validated across all three regions. In region 1, the results vary more significantly with different contact angles, while in regions 2 and 3, the results are more consistent across various contact angles. This phenomenon is expected since the evaporation rate is related to the integral area, and the microlayer extent is more sensitive to variations in the contact angle compared with the bubble volume and lateral radius. Following the trend of the microlayer length shown in figure 16(d), the evaporation rate in region 1 increases at a decreasing rate over time. Notably, with a contact angle of
$30^\circ$
, the evaporation rate in region 1 begins to decrease after
$t = 0.3$
ms, as the velocity of the microlayer front falls below that of the contact line. In region 2, the evaporation rate first increases and then stabilises at
$t = 0.2$
ms. The situation in region 3 is noticeably different: the evaporation rate initially increases, peaks at
$t = 0.13$
ms and then decreases rapidly. The different trends in the three regions arise because regions 1 and 2 are close to the solid surface and remain within the thermal boundary layer. In contrast, region 3, the upper surface region of the bubble, moves above the thermal boundary layer due to its axial growth, as observed in figure 8.
Time histories of the average mass flux over (a) region 1, (b) region 2, (c) region 3 and (d) the entire bubble for various contact angles. The results are obtained with
$\omega = 0.0460$
at grid level 11.
(a) Total evaporated mass in different regions for different contact angles. (b) Percentage of the evaporated mass in each region relative to the total evaporated mass for different various angles. (c) Time-averaged surface area for different regions and contact angles. (d) Total average mass flux in different regions for different contact angles. The results are obtained with
$\omega = 0.0460$
at grid level 11.
As the evaporation rate is strongly related to the integral area, the average mass flux, defined as
\begin{equation} \widetilde {\dot {m}} = \frac {\iint \dot {m} \,\text{d}A}{\iint \,\text{d}A}, \end{equation}
is considered to further investigate the influence of the contact angle. The time histories of the average mass flux within different regions are presented in figure 20. It is observed that, for all regions, the average mass flux decreases at a decelerating rate, consistent with the theoretical analysis (Scriven Reference Scriven1959) of bubble growth driven by heat diffusion. Notably, as shown in figure 20, the average mass flux for different contact angles matches closely, especially in regions 2 and 3. This suggests that the observed differences in evaporation rates are primarily due to variations in integral area.
Additionally, beyond the temporal variations of thermal characteristics, the cumulative behaviour is investigated. The total evaporated mass
$M_e$
is calculated by integrating the evaporation rate over time:
\begin{equation} M_e = \int _0^{t_{\!f}}\left (\iint \dot {m} \,\text{d}A\right ) \,\text{d}t, \end{equation}
where
$t_{\!f} = 0.5\ \mathrm{ms}$
is the final time of the simulations. For different contact angles, the total evaporated mass in different regions is computed and shown in figure 21(a). It can be seen that as the contact angle increases, the evaporated mass in the microlayer region is most affected. The proportions of the evaporated mass in each region relative to the total evaporated mass for different contact angles are detailed in figure 21(b). As the contact angle increases from
$5^\circ$
to
$30^\circ$
, the contribution from the microlayer region decreases from 31 % to 25 %. In figure 21(c), the time-averaged surface area in different regions,
\begin{equation} \overline {S} = \frac {1}{t_{\!f}}\int _0^{t_{\!f}}\left (\iint \,\text{d}A\right ) \,\text{d}t, \end{equation}
is evaluated for different contact angles. The correlation between the trends of the average surface area and the evaporated mass relative to the contact angle is evident. To specifically investigate the influence of the contact angle on global thermal effects, we define the total average mass flux as
\begin{equation} \overline {\dot {m}} = \frac {M_e}{\overline {S} t_{\!f}} = \frac {\int _0^{t_{\!f}}\left (\iint \dot {m} \,\text{d}A\right ) \,\text{d}t}{\int _0^{t_{\!f}}\left (\iint \,\text{d}A\right ) \,\text{d}t}. \end{equation}
The results for different regions are plotted against the contact angle in figure 21(d). It is observed that, although the evaporated mass and surface area vary with different contact angles, the total average mass flux remains consistent. In particular, in regions 2 and 3, the total average mass fluxes are almost constant, regardless of the contact angle. As the contact angle increases from
$5^\circ$
to
$30^\circ$
, the total average mass flux over the entire bubble decreases slightly (from
$0.47$
to
$0.44\ \mathrm{kg}\,\mathrm{(m^2\,s)}^{-1}$
), mainly due to minor changes in the microlayer thickness. Overall, it can be concluded that the hydrodynamic effect is the dominant factor influencing bubble growth over different contact angles, while the thermal effect remains consistent regardless of the contact angle. A larger contact angle negatively affects the evaporation process because the increased mobility of the contact line results in a smaller surface area within the microlayer region. Over time, the reduction of area in the microlayer region, region 1, slows down the bubble expansion due to the reduced amount of vapour evaporated from this region, yielding a smaller area for regions 2 and 3. The reduction of the evaporation rate in regions 2 and 3, in turn, further decelerates bubble growth. This explains why the differences in results for different contact angles increase with time, as shown in figure 16.
4.5. Complete bubble cycle
In previous DNS studies of nucleate boiling in the microlayer regime (Bureš & Sato Reference Bureš and Sato2022; Urbano et al. Reference Urbano, Tanguy, Huber and Colin2018), the total physical time simulated was generally quite short (less than 2 ms) due to high computational costs. This duration is notably brief compared with the entire bubble cycle from nucleation to detachment. In this section, using the present, more efficient and stable solver, we perform DNS for the entire cycle of bubble growth, which, to our knowledge, is the first reported in the open literature. It should be noted that, in the present study, no additional physical model is included for the bubble detachment process. The topological change associated with detachment occurs numerically when the neck connecting the bubble to the heated wall becomes thinner than the cell size. However, as shown by Eggers (Reference Eggers1997), the cylinder pinching is a finite time singularity which will happen in very similar ways and similar time at sufficiently high grid resolutions. Moreover, as reported in previous studies (Rana, Das & Das Reference Rana, Das and Das2017; Allred, Weibel & Garimella Reference Allred, Weibel and Garimella2021), this simplified treatment allows the interface evolution up to the point of detachment to be correctly captured, provided that the mesh resolution is sufficiently high. Since the detailed dynamics of detachment is not the primary focus of this work, a more advanced detachment model is left for future development.
The set-up remains consistent with the previous configuration, as shown in figure 6, except that the computational domain is extended to
$ [ 0\ \mathrm{mm}, 5.5 \ \mathrm{mm} ] \times [ -1.0\ \mathrm{mm}, 4.5 \ \mathrm{mm} ]$
to capture the larger bubble encountered during the simulation. Three grid levels, from 11 to 13, are adopted, resulting in minimum grid sizes ranging from
$2.69$
to
$0.67\ \mathrm{\unicode{x03BC} m}$
. The computational costs for simulations up to the physical time of 18 ms are detailed in table 3 for different grid resolutions. The promising computational efficiency demonstrates the potential applications of the present solver in studies of more challenging nucleate boiling problems, such as flow boiling.
Computational cost for the DNS of a complete bubble cycle.
In the experiment of Bucci (Reference Bucci2020), after the depletion of the microlayer, contact angle hysteresis is observed. The contact angle hysteresis effect is that a contact line can recede (i.e. expansion of the dry region) only for contact angles lower than a critical value
$\theta _{\textit{rec}}$
and advance (i.e. shrinkage of the dry region) for contact angles higher than another critical value
$\theta _{\textit{ad}v}$
. These two critical angles are referred to as the receding contact angle and the advancing contact angle. Due to the presence of the microlayer, the receding contact angle cannot be measured, while the advancing contact angle is measured at approximately
$50^\circ{-}55^\circ$
. In the present work, a contact-angle hysteresis model (Bureš et al. Reference Bureš, Bucci, Sato and Bucci2024; Fang et al. Reference Fang, Hidrovo, Wang, Eaton and Goodson2008) is adopted, and the contact angle is updated from time step
$n$
to
$n+1$
by
\begin{equation} \theta ^{n+1}=\max \left [\theta _{\textit{rec}}, \min \left (\theta _{\textit{ad}v}, \theta ^n-2 \sin ^2\left (\theta ^n\right ) \frac {H_l^{n+1}-H_l^n}{\varDelta }\right )\right ]\!, \end{equation}
where
$H_l$
represents the total wall-parallel liquid height in the contact-line cell row. This formula is derived based on the assumption that, for a contact angle between the receding angle
$\theta _{\textit{rec}}$
and the advancing angle
$\theta _{\textit{ad}v}$
, any change in the liquid content within the contact-line cell row corresponds to a rotation of the interface rather than a shift in the contact-line position. Following Bureš et al. (Reference Bureš, Bucci, Sato and Bucci2024), we use
$\theta _{\textit{rec}} = 5^\circ$
and
$\theta _{\textit{ad}v} = 55^\circ$
in our simulations.
Subsequently, the numerical results are compared with the experimental data. With
$\omega = 0.0345$
, the time histories of the bubble volume and the bubble lateral radius at different grid levels are presented in figures 22(a) and 22(b). It can be seen that the bubble growth rate is underestimated after
$t = 2.0$
ms, although good agreement is observed in the early stages. The oscillations of the lateral radius in the later stages are caused by deformations of the bubble after its detachment. In addition to
$\omega = 0.0345$
, the other accommodation coefficient used by Bureš & Sato (Reference Bureš and Sato2022),
$\omega = 0.0460$
, still leads to a significant underestimation of the bubble growth rate, and the corresponding results are not shown here for brevity. It is important to note that the modelling of the accommodation coefficient remains an open question in the literature (Cai et al. Reference Cai, Wang, Jiang and Tan2024). Our goal here is not to solve this problem, but to demonstrate the effect of the accommodation coefficient on the DNS of complete bubble cycles in nucleate boiling. In previous studies (Bureš & Sato Reference Bureš and Sato2022; El Mellas et al. Reference El Mellas, Samkhaniani, Falsetti, Stroh, Icardi and Magnini2024; Torres et al. Reference Torres, Urbano, Colin and Tanguy2024), the physical time of the simulations was not large enough to observe these differences. In the present work, the perfect evaporation case (
$\omega = 1.0$
) is also simulated, and the results are shown in figures 22(c) and 22(d). It can be observed that, with perfect evaporation within the microlayer, the growth of the bubble volume and the bubble lateral radius is in better agreement with the experimental measurements. This is expected because diminishing interfacial heat-transfer resistance, corresponding to an increasing accommodation coefficient, leads to a higher evaporation rate.
Time histories of the bubble volume and the bubble lateral radius at different grid levels. The simulation results obtained with (a,b)
$\omega = 0.0345$
and (c,d)
$\omega = 1.0$
are compared with the experimental data of Bucci (Reference Bucci2020).
Superheat distributions at grid level 13 for different time instants. The simulation results obtained with (a)
$\omega = 0.0345$
and (b)
$\omega = 1.0$
are compared with the experimental data of Bucci (Reference Bucci2020).
Bubble interfaces at grid level 13 for different time instants. The simulation results obtained with (a)
$\omega = 0.0345$
and (b)
$\omega = 1.0$
are compared with the experimental data of Bucci (Reference Bucci2020). The results are mirrored about
$r = 0$
for better visualisation.
In figure 23, the superheat distributions along the wall at different time instants are compared with the experimental data. With
$\omega = 0.0345$
, at
$t = 3$
ms and
$t = 9$
ms, despite the deviations in the positions of the contact line and the microlayer front (caused by the underestimated bubble growth rate), the predicted lowest superheat on the wall is in better agreement with the experiment compared with that in the perfect evaporation case. A possible explanation for these deviations is that the rectangular heater is modelled as an axisymmetric configuration in numerical simulations. Consequently, the energy input is lower than in actual experiments due to the missing heater area. This accumulated effect might not have been identified in previous studies due to the limited simulation time. In the perfect evaporation case, this energy loss is compensated by a lower IHTR, allowing more energy transfer from the solid to the fluid. However, in the current case, the selected accommodation coefficient overcompensates for this loss, leading to an underestimation of the superheat on the solid wall. These findings highlight the necessity of 3-D simulations to better assess the cumulative effect of heater geometry. The highly efficient solver presented in this paper provides a promising platform for tackling such challenging simulations in future work. At
$t = 15$
ms, the superheat distribution with
$\omega = 0.0345$
deviates more from the experimental data than that with
$\omega = 1.0$
. This deviation is due to the advanced detachment of the bubble, as shown in figure 24, where the bubble interfaces at different time instants are presented. It can be seen that at
$t = 3$
ms and
$t = 9$
ms, the interfaces obtained with
$\omega = 1.0$
are in excellent agreement with the experimental measurements. A larger deviation is observed at
$t = 15.6$
ms, which could be attributed to the contact angle model. During the detachment stage, according to (4.10), the contact angle in the simulations is fixed at
$\theta _C = \theta _{\textit{ad}v} = 55^\circ$
. The time-dependent contact angle observed during the detachment stage in the experiment (Bucci Reference Bucci2020) is not modelled in the present study. Hence, it can be concluded that a better dynamic contact angle model is needed for accurately predicting the detachment at the end of the bubble cycle. This has not been discussed in previous DNS studies (Bureš & Sato Reference Bureš and Sato2022; Urbano et al. Reference Urbano, Tanguy, Huber and Colin2018; Torres et al. Reference Torres, Urbano, Colin and Tanguy2024; El Mellas et al. Reference El Mellas, Samkhaniani, Falsetti, Stroh, Icardi and Magnini2024), as the simulation time was not long enough to reach the detachment stage.
Microlayer profiles and velocity vectors in the liquid phase at different time instants, obtained at grid level 13.
Finally, to better visualise the depletion of the microlayer during the long-time simulation, the microlayer profiles along with the velocity vectors in the liquid phase at different times are shown in figure 25. It is evident that both the movement of the contact line and the evaporation of the microlayer contribute to the dryout process, consistent with previous studies (Urbano et al. Reference Urbano, Tanguy, Huber and Colin2018; Bureš & Sato Reference Bureš and Sato2022). Throughout the simulation, the liquid velocity remains small, with the flow mainly driven by the contact line movement directed radially outward. Parasitic currents are observed and become pronounced as the microlayer depletes, particularly when capillary waves are transported along the interface and the microlayer becomes increasingly thin (Bureš & Sato Reference Bureš and Sato2022), as shown in figures 25(b) and 25(c). As the microlayer thickness eventually approaches the size of a computational cell, these parasitic currents cannot be fully resolved by simply increasing mesh resolution. Instead, they may be addressed through the incorporation of depletion models (Sato & Niceno Reference Sato and Niceno2015). However, the implementation of such models is beyond the scope of the present study.
5. Conclusion
We have extended the open-source phase-change model developed by Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024) to simulate nucleate boiling in the microlayer regime with resolved conjugate heat transfer. Following Bureš & Sato (Reference Bureš and Sato2022), heat transfer in the fluid and solid is coupled in a fully implicit manner, with a temperature jump condition accounting for the IHTR. The current work is based on the free software Basilisk (Popinet Reference Popinet2009, Reference Popinet2015), in which the quad/octree-based AMR technique is employed to improve computational efficiency. To facilitate the implementation of AMR, a cell-centred velocity is adopted, though it can only be approximately projected (2.22). Based on the approximate projection method, the original model of Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024) works well for numerous benchmark tests (Cipriano et al. Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024). However, we have shown in the Appendix that the original model of Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024) introduces significant numerical oscillations within the microlayer (Zhao et al. Reference Zhao, Zhang and Ni2022; Long et al. Reference Long, Pan and Zaleski2024), thereby failing to accurately predict its development. The ghost fluid method, which removes the singularity at the interface by setting the ghost velocity, is employed and effectively suppresses these oscillations. With the ghost fluid method, we have successfully replicated the pool boiling experiments conducted at MIT (Bucci Reference Bucci2020). The numerical results are in good agreement with the experimental data of Bucci (Reference Bucci2020) and the previous numerical results of Bureš & Sato (Reference Bureš and Sato2022). With the help of AMR, computational efficiency is significantly improved and the required CPU hours are reduced by three orders of magnitude. We thus believe the current work makes the present model more applicable for complex phase-change problems with high fidelity. The codes and configurations for all simulations are released in the Basilisk sandbox (Long Reference Long2024).
Subsequently, with the present model, the influence of the contact angle on nucleate boiling in the microlayer regime is investigated. We have shown that the value of the contact angle influences the results in a decelerating manner, and very small contact angles, such as those predicted by the dynamic contact angle model proposed by Bureš & Sato (Reference Bureš and Sato2022), are not necessary for the current scenario. Moreover, by dividing the bubble surface into three regions, we have shown that the influence of the contact angle is mainly confined to the microlayer region. It is demonstrated that thermal effects exhibit similarity across different contact angles, while hydrodynamic effects predominantly influence bubble growth. As the contact angle increases, the growing contact line mobility leads to a smaller surface area, while the total average mass flux remains approximately constant. Moreover, a complete bubble cycle from nucleation to detachment has been directly simulated with a resolved microlayer and conjugate heat transfer. The predicted bubble shapes show a good agreement with the experimental data, and the influence of dynamic contact angle models and accommodation coefficient on the long-term behaviour of the bubble are discussed. These aspects were previously obscured by the challenges posed by the high computational burden. To the best of our knowledge, the present work represents the first such effort reported in the open literature. We believe the present study has effectively demonstrated the capability of the DNS solver for nucleate boiling problems. Several improvements are considered for our future work, including more advanced IHTR models and dynamic contact angle models in the context of nucleate boiling. Furthermore, the bubble growth rate is governed by the interplay between inertial and heat diffusion effects, with the maximum expansion rate occurring during the initial inertial-controlled stage (
${\sim}10$
–100
$\unicode{x03BC}$
s), where compressible effects may become important. In the present set-up, characterised by atmospheric pressure and low superheat, the observed bubble expansion speed is relatively low (
${\sim}1$
m s−1), which justifies the use of the incompressible assumption. This assumption is also supported by previous studies on microlayer formation in nucleate boiling using an incompressible model (Guion et al. Reference Guion, Afkhami, Zaleski and Buongiorno2018). However, under more extreme conditions, such as high pressure and superheat, rapid bubble expansion in the early inertial-controlled stage can generate acoustic waves in the liquid, potentially leading to coupling effects between different nucleation sites when multiple nucleating bubbles are present (Mallozzi, Judd & Balakrishnan Reference Mallozzi, Judd and Balakrishnan2000). In such cases, compressible phase change solvers are desirable to accurately capture these dynamics. Recent developments in this area (Bibal et al. Reference Bibal, Deferrez, Tanguy and Urbano2024; Coseru et al. Reference Coseru, Tanguy, Freton, Gonzalez, Urbano, Bibal and Bourdon2025) may be considered for future extensions of the present model.
Funding
This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement number 883849). We are grateful to GENCI for generous allocation on Adastra supercomputers (grant agreement number A0152B14629) and their teams for assistance and the use of Irene-Rome at TGCC.
Declaration of interests
The authors report no conflict of interest.
Appendix. Ghost fluid method versus one-fluid method
In the present study, it is evident that the computational cost can be significantly reduced using AMR. However, it is emphasised that the use of AMR necessitates a compromise: the cell-centred velocity field can only be approximately projected, leading to unphysical oscillations in the presence of phase change. Based on the well-balanced VOF framework (Popinet Reference Popinet2009), the original model of Cipriano et al. (Reference Cipriano, Frassoldati, Faravelli, Popinet and Cuoci2024) performs well in a number of benchmark tests, even with the approximate projection method. However, as shown here, the original model fails when applied to nucleate boiling with a microlayer, mainly due to the intense heat and mass transfer within this thin layer. Here, the one-fluid method is compared with the ghost fluid method. Using the same set-up as given in § 4.1, case A (
$\omega = 0.0345$
) is simulated with the one-fluid method across grid levels 10–12. The time histories of the volume and lateral radius of the bubble at these grid levels are presented in figure 26. It is shown that the results at grid level 10 obtained with the one-fluid method agree better with the experimental data of Bucci (Reference Bucci2020). However, this better agreement is merely a numerical artefact, as increasing grid resolution leads to significant deviations rather than grid convergence. The microlayer profiles at various time instants and grid levels are depicted in figure 27. It can be seen that the microlayer obtained with the one-fluid method exhibits more oscillations and is thinner than that obtained with the ghost fluid method (figure 13). As indicated in figures 28 and 29, compared with the experimental data, larger heat fluxes and smaller surface temperatures are obtained using the one-fluid method, attributed to the reduced microlayer thickness. Moreover, these oscillations in the microlayer also lead to erroneous surface temperature and heat flux distributions.
Time histories of (a) the bubble volume and (b) the bubble lateral radius for case A (
$\omega$
= 0.0345) at different grid levels. The results are obtained with the one-fluid method and are compared with the experimental data of Bucci (Reference Bucci2020) and the numerical results of Bureš & Sato (Reference Bureš and Sato2022).
Microlayer profiles at different time instants and grid levels for case A (
$\omega = 0.0345$
). The results are obtained with the one-fluid method and are compared with the numerical results of Bureš & Sato (Reference Bureš and Sato2022).
Surface superheat distribution at different time instants and grid levels for case A (
$\omega = 0.0345$
). The results are obtained with the one-fluid method and are compared with the experimental data of Bucci (Reference Bucci2020) and the numerical results of Bureš & Sato (Reference Bureš and Sato2022).
Surface heat flux distribution at different time instants and grid levels for case A (
$\omega = 0.0345$
). The results are obtained with the one-fluid method and are compared with the experimental data of Bucci (Reference Bucci2020) and the numerical results of Bureš & Sato (Reference Bureš and Sato2022).
Magnitudes of the cell-centered velocity
$\boldsymbol{u}_c$
obtained with (a) the one-fluid method and (b) the ghost fluid method at
$t = 0.21\ \mathrm{ms}$
and level 12. The blue line represents the liquid–vapour interface, while the red line represents the fluid–solid boundary. (c,d) Close-ups of the selected region, with the velocity vectors plotted.
For a better comparison between the one-fluid method and the ghost fluid method, the magnitudes of the cell-centred velocity
$\boldsymbol{u}_c$
obtained with the two methods at
$t = 0.21\ \mathrm{ms}$
and grid level 12 are given in figure 30. It can be observed from figures 30(a) and 30(b) that more pronounced numerical oscillations are introduced by the one-fluid method. When examining a selected region and comparing the velocity fields in figures 30(c) and 30(d), we observe that the fluid within the microlayer is almost stagnant when using the ghost fluid method, consistent with findings in previous studies (Urbano et al. Reference Urbano, Tanguy, Huber and Colin2018; Bureš & Sato Reference Bureš and Sato2022). In contrast, significant oscillations are observed with the one-fluid method, as indicated by the erroneous vertical component of the velocity. This vertical component consequently results in an incorrect microlayer thickness.
Open access
