2.1. Phase conservation and interface tracking

Two phases, liquid and vapour, are included in PSI-BOIL. Their mass densities are considered to be constant, $\rho _l$ and $\rho _v$ (${\rm kg}~{\rm m}^{-3}$); thus, to track the phasic distribution, it is sufficient to solve the transport equation for the liquid volume fraction, $\phi$, defined for a computational cell with volume $\Delta V$ (m$^3$) as

(2.1)\begin{equation} \phi = \frac{V_l}{\Delta V}, \end{equation}

where $V_l$ is the volume of the cell occupied by liquid. The vapour volume fraction is then $(1-\phi )$.

The transport equation for $\phi$ can be derived from general two-phase integral balance considerations and reads (Bureš & Sato Reference Bureš and Sato2021b)

(2.2)\begin{equation} \frac{\partial\phi}{\partial t} ={-}\boldsymbol{\nabla} \boldsymbol{\cdot} (\phi \boldsymbol{u}_l) - \frac{{\dot{m}}'''}{\rho_l} ={-} \frac{1}{\Delta V} \sum_i F_{l,i} - \frac{{\dot{m}}'''}{\rho_l}. \end{equation}

Here, $\boldsymbol {u}_l$ (${\rm m}~{\rm s}^{-1}$) is the (divergence-free) liquid velocity field, satisfying the incompressibility condition, and computed from a modified version of the extrapolation algorithm of Malan et al. (Reference Malan, Malan, Zaleski and Rousseau2021), without which consistent advection of $\phi$ in the presence of phase change would not be possible. Furthermore, $\dot {m}'''$ is the evaporative mass source density $({\rm kg}~({\rm m}^{3}~{\rm s})^{-1})$. The advection term is treated using the GVOF method of Weymouth & Yue (Reference Weymouth and Yue2010), expressed in terms of the liquid surface flows $F_{l,i}$ $({\rm m}^3~{\rm s}^{-1})$ defined for each bounding surface $i$ of the computational cell. Equation (2.2) is discretised in explicit form as follows:

(2.3)\begin{gather} \frac{\phi^\star{-}\phi^{n}}{\Delta t} ={-} \frac{\dot{m}'''}{\rho_l}, \end{gather}

(2.4)\begin{gather}\frac{\phi^{n+1}-\phi^\star}{\Delta t} ={-} \frac{1}{\Delta V}\sum_i F_{l,i}, \end{gather}

where $\Delta t$ (s) is the discrete time step, $\phi ^\star$ is an intermediate value of the volume fraction, while $n$ and $n+1$ indicate the previous and new time steps, respectively. The evaporative mass source density needed for the phase-change step (2.3) is calculated according to

(2.5)\begin{equation} \dot{m}''' = a_\gamma \frac{1}{L} (\lambda_l\boldsymbol{\nabla} T|_{\gamma,l} \boldsymbol{\cdot}\boldsymbol{n}-\lambda_v\boldsymbol{\nabla} T|_{\gamma,v} \boldsymbol{\cdot}\boldsymbol{n}). \end{equation}

Here, $a_\gamma$ (m$^2$ m$^{-3}$) is the interfacial area density, computed by means of the marching cubes algorithm (Lorensen & Cline Reference Lorensen and Cline1987), and based on the local $\phi$ stencil, adapted for use in axisymmetric geometry; $L$ (J kg$^{-1}$) is the latent heat of vaporisation; and $\lambda _l$ and $\lambda _v$ are the liquid and vapour thermal conductivities $({\rm W}~({\rm m}~{\rm K})^{-1})$, respectively. Furthermore, $\boldsymbol {\nabla } T|_{\gamma,l}$ and $\boldsymbol {\nabla } T|_{\gamma,v}$ are the temperature gradients on the liquid and vapour sides of the interface $({\rm K}~{\rm m}^{-1})$. They are calculated using fourth-order-accurate upwind differences for non-uniform grids, evaluated at the interface; details appear in Bureš & Sato (Reference Bureš and Sato2021b). The unit normal vector to the interface pointing from the vapour phase to the liquid phase, $\boldsymbol {n}$, is evaluated using the ELVIRA algorithm of Pilliod & Puckett (Reference Pilliod and Puckett2004). Note that the interfacial jump in kinetic energy is neglected, being in general small in comparison to the latent heat of phase transition, as can be deduced from a scaling analysis.

The second step in the solution of the transport equation (2.4) is performed using the directional-split, bounded-conservative GVOF algorithm of Weymouth & Yue (Reference Weymouth and Yue2010), also adapted for use in axisymmetric geometry.

2.2. Momentum conservation

The two-phase momentum conservation equations, written in the single-fluid approximation, read as follows:

(2.6)\begin{equation} \rho \frac{\text{D}\boldsymbol{u}}{\text{D}t} ={-}\boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot} [\mu(\boldsymbol{\nabla} \boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^{\rm T})] + \rho \boldsymbol{g} + \,\boldsymbol{f}_\sigma. \end{equation}

Note that the recoil pressure due to phase change has not been considered, since it is generally negligible within the present application area, as we concluded using a scaling analysis. In this equation, $\boldsymbol {u}$ is the single-fluid mechanical velocity, $p$ is the pressure (Pa), $\boldsymbol {g}$ is the gravitational acceleration (${\rm m}~{\rm s}^{-2}$), $\,\boldsymbol {f}_\sigma$ is the surface-tension force density $({\rm N}~{\rm m}^{-3})$ and $\rho$ is the mixture density, given by

(2.7)\begin{equation} \rho = \phi\rho_l + (1-\phi)\rho_v. \end{equation}

In the single-fluid approximation, $\boldsymbol {u}$ is taken to be continuous across the interface and the stress balance at the interface is not represented exactly. Instead, a heuristic harmonic mixing law is applied for the computation of mixture dynamic viscosity $\mu$ $({\rm Pa}~{\rm s})$ in the vicinity of the interface. In the simulations shown below, we use the harmonic mixing law (Prosperetti Reference Prosperetti2002):

(2.8)\begin{equation} \frac{1}{\mu} = \frac{\phi}{\mu_l} + \frac{1-\phi}{\mu_v}. \end{equation}

We have also conducted simulations of the validation problem presented in § 4 using other heuristic viscosity averaging schemes; no discernible differences in the overall dynamics have been observed. Note that, since the near-interface stress balance is anyway treated only approximately, we neglect the $(\boldsymbol {\nabla }\boldsymbol {u})^{\rm T}$ term in (2.6) for reasons of simplicity and robustness of the numerical implementation.

In PSI-BOIL, (2.6) is solved in the finite-volume formulation:

(2.9)\begin{align} &\rho\left(\frac{\partial\boldsymbol{u}}{\partial t}\,\Delta V + \sum_i (\boldsymbol{u}\otimes\boldsymbol{u}) \boldsymbol{\cdot}\Delta\boldsymbol{S}_i - (\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u})\boldsymbol{u}\,\Delta V \right) \nonumber\\ & \quad ={-} \boldsymbol{\nabla} p\, \Delta V + \sum_i \mu\boldsymbol{\nabla}\boldsymbol{u}\boldsymbol{\cdot} \Delta\boldsymbol{S}_i + \rho\boldsymbol{g}\,\Delta V + \,\boldsymbol{f}_\sigma \,\Delta V. \end{align}

The discretised variables are arranged using the staggered-grid approach of Harlow & Welch (Reference Harlow and Welch1965). For the spatial discretisation, a second-order-accurate, central-difference scheme and a second-order flux-limiting, total-variation-diminishing (TVD) scheme (Roe Reference Roe1986) are used for the diffusion and advection terms, respectively. For the time discretisation, backward and forward Euler methods are used to treat the diffusion and advection terms, respectively. The Boussinesq approximation (Gray & Giorgini Reference Gray and Giorgini1976) for density is adopted: i.e. the variability of densities with temperature is only taken into account in the calculation of the buoyancy force. The surface-tension force density is given by

(2.10)\begin{equation} \boldsymbol{f}_\sigma = \sigma \kappa a_\gamma \boldsymbol{n}, \end{equation}

with $\sigma$ $({\rm N}~{\rm m}^{-1})$ being the surface tension (assumed constant), and $\kappa$ $({\rm m}^{-1})$ the interfacial curvature, which is taken to be positive if the centre of curvature is in the liquid phase and negative otherwise. Note that Marangoni stresses (Faghri & Zhang Reference Faghri and Zhang2006) have been neglected based on our analysis of the heat-flux conditions within the microlayer for the experiment discussed in § 4 (see § 3 for the description of the relation between heat flux and interfacial temperature). Additionally, disjoining pressure (Faghri & Zhang Reference Faghri and Zhang2006; Janeček & Nikolayev Reference Janeček and Nikolayev2013) has not been considered, since this phenomenon becomes relevant at length scales well below the grid resolution employed in the simulations described in this paper.

As is common practice, (2.10) is discretised using the continuum surface force approach of Brackbill, Kothe & Zemach (Reference Brackbill, Kothe and Zemach1992):

(2.11)\begin{equation} \boldsymbol{f}_\sigma = \sigma\kappa \frac{2\rho}{\rho_l+\rho_v}\, \frac{\boldsymbol{\nabla}\rho}{\rho_l-\rho_v}. \end{equation}

In our approach, the curvature $\kappa$ is calculated by means of the height function method (Popinet Reference Popinet2018), with a $3 \times 7$ stencil utilising the local-topology adaptation approach of López et al. (Reference López, Zanzi, Gómez, Zamora, Faura and Hernández2009), again modified for use in axisymmetric geometry.

2.3. Mass conservation

As the densities of the phases are here considered constant, the mass continuity equation may be conveniently replaced by the equivalent volume-conservation condition, resulting from the application of the two-phase divergence theorem, as

(2.12)\begin{equation} \sum_i \boldsymbol{u}\boldsymbol{\cdot}\Delta \boldsymbol{S}_i = \dot{m}'''\left(\frac{1}{\rho_v}-\frac{1}{\rho_l}\right)\Delta V. \end{equation}

Here, $\Delta \boldsymbol {S}_i$ is the (outward-oriented) area of the bounding surface $i$ within the given computational cell. This condition is enforced at each time step by the projection algorithm of Chorin (Reference Chorin1968), which consists of solving the discrete Poisson equation for the pressure $p$ at the advanced time step $n+1$,

(2.13)\begin{equation} \sum_i \frac{\Delta t}{\rho}\,\boldsymbol{\nabla} p^{n+1} \boldsymbol{\cdot} \Delta\boldsymbol{S}_i = \sum_i \boldsymbol{u}^\star \boldsymbol{\cdot} \Delta\boldsymbol{S}_i - \dot{m}^{\prime\prime\prime,n}\left(\frac{1}{\rho_v}-\frac{1}{\rho_l}\right) \Delta V, \end{equation}

followed by the projection step,

(2.14)\begin{equation} \boldsymbol{u}^{n+1} = \boldsymbol{u}^\star{-} \frac{\Delta t}{\rho} \, \boldsymbol{\nabla} p^{n+1}. \end{equation}

In these equations, $\boldsymbol {u}^\star$ is a tentative estimate of the velocity field, obtained from the numerical solution of the following equation:

(2.15)\begin{align} &\frac{\rho^n\boldsymbol{u}^\star}{\Delta t}\,\Delta V - \sum_i \mu\boldsymbol{\nabla}\boldsymbol{u}^{{\star}}\boldsymbol{\cdot} \Delta\boldsymbol{S}_i \nonumber\\ & \quad = \rho^n\left(\frac{\boldsymbol{u}^n}{\Delta t}\,\Delta V -\sum_i (\boldsymbol{u}^n\otimes\boldsymbol{u}^n) \boldsymbol{\cdot} \Delta\boldsymbol{S}_i + (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^n) \boldsymbol{u}^n\, \Delta V \right) + \rho^n\boldsymbol{g}\,\Delta V + \,\boldsymbol{f}^n_\sigma\,\Delta V. \end{align}

Note that, by assuming the validity of (2.12), we are effectively identifying the single-fluid velocity $\boldsymbol {u}$ with the volume-averaged velocity.

2.4. Energy conservation

The energy conservation equation is the only balance equation in PSI-BOIL not recast in integral form, but rather solved in a finite-difference approximation to its differential form:

(2.16)\begin{gather} C_{p,l} \left(\frac{\partial T}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(T\boldsymbol{u}_l) -T(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_l)\right)= \lambda_l\nabla^2 T \quad \text{in the liquid}, \end{gather}

(2.17)\begin{gather}C_{p,v} \left(\frac{\partial T}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(T\boldsymbol{u}_v) -T(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_v)\right)= \lambda_v\nabla^2 T\quad \text{in the vapour}, \end{gather}

(2.18)\begin{gather}T = T_\gamma\quad \text{at the interface}. \end{gather}

The equation is discretised separately in the liquid and vapour regions of the domain, with the phasic interface being treated as an internal boundary, unlike for the momentum equations. Since the finite-volume scheme cannot be easily employed in the vicinity of the phasic interface, where the computational cells are cut, the finite-difference approach is employed here. In the above equations, $C_{p,l}=\rho _l c_{p,l}$ and $C_{p,v}=\rho _v c_{p,v}$ are the volumetric heat capacities $({\rm J}~({\rm m}^3~{\rm K})^{-1})$ for the liquid and vapour, respectively; and $c_{p,l}$ and $c_{p,v}$ are the respective specific heat capacities $({\rm J}~({\rm kg}~{\rm K})^{-1})$.

Note that the form of the energy conservation equation presented here could be derived from the internal energy balance equation, assuming zero internal heating and neglecting heat generation due to viscous stress work and compressibility effects, all justified within the context of the present application. A sharp-interface, finite-difference algorithm is employed for its spatial discretisation: the diffusion term is treated implicitly, with the interfacial position resolved with subgrid accuracy, utilising the information provided by the GVOF interface-tracking method referred to earlier. A three-point, central-difference scheme is employed, of second-order accuracy for steady-state conditions. For problems depending on the gradients of the solution, such as the bubble-growth problem of Scriven (Reference Scriven1959), first-order accuracy is retrieved (Gibou et al. Reference Gibou, Fedkiw, Cheng and Kang2002). (We thank one of the referees for clarifying this point during the review process.) The advection terms are treated explicitly, employing a second-order TVD scheme. The phasic velocities, $\boldsymbol {u}_l$ and $\boldsymbol {u}_v$, are obtained from a variation of the extrapolation algorithm of Malan et al. (Reference Malan, Malan, Zaleski and Rousseau2021) and ‘ghost’ values of temperature are calculated by means of linear extrapolation across the interface, assuming continuity of temperature at the interface itself, calculated at the appropriate interface position. Note that the temperature gradient is not continuous at the interface, due to the step change in the thermal conductivities occurring there, and the presence of a non-zero phase-change rate (2.5).

The interfacial temperature $T_\gamma$ is assumed to be constant and equal to the saturation temperature at the prevailing system pressure: see § 3 for a discussion of the role and implementation of the IHTR.

2.5. Treatment of the solid

In the solid part of the domain, only heat conduction needs to be considered:

(2.19)\begin{equation} C_{p,s}\frac{\partial T}{\partial t} = \sum_i\lambda_s\boldsymbol{\nabla} T|_i \boldsymbol{\cdot}\frac{\Delta\boldsymbol{S}_i}{\Delta V} + q, \end{equation}

where $q$ $({\rm W}~{\rm m}^{-3}$) is the volumetric heat source, e.g. due to electrical resistance heating. Note that a finite-volume formulation, rather than a finite-difference one, is employed, since it was found to have better convergence characteristics in the presence of a heat source in preliminary 1-D test studies. With the solid assumed to be a single phase, the finite-volume formulation can be employed without issues. If it is composed of several layered materials, simple averaging based on volume fractions is used to obtain the cell-specific values of $C_{p,s}$ and $\lambda _s$. This treatment is exact for $C_{p,s}$ but only approximate for $\lambda _s$, the degree of approximation depending on the grid resolution.

2.5.1. Wall boundary conditions for heat transfer

At the solid–fluid boundary, continuity of heat flux is assumed. The existence of a temperature contact discontinuity is permitted, and modelled in terms of a contact heat-transfer resistance $\mathcal {R}_c$ $({\rm K~m}^2~{\rm W}^{-1})$. As a result, the solid-side and fluid-side boundary temperatures can be derived from

(2.20)\begin{gather} T_{{b},s} = \frac{(\mathcal{R}_c+\mathcal{R}_f) T_s+ \mathcal{R}_s T_f}{\mathcal{R}_c+\mathcal{R}_f+\mathcal{R}_s } = \frac{(\mathcal{R}_c+\Delta z_f/\lambda_f) T_s+ (\Delta z_s/\lambda_s) T_f}{\mathcal{R}_c+\Delta z_f/\lambda_f +\Delta z_s/\lambda_s }, \end{gather}

(2.21)\begin{gather}T_{{b},f} = \frac{\mathcal{R}_f T_s+ (\mathcal{R}_c+\mathcal{R}_s) T_f}{\mathcal{R}_c+\mathcal{R}_f+\mathcal{R}_s } = \frac{(\Delta z_f/\lambda_f) T_s+ (\mathcal{R}_c+\Delta z_s/\lambda_s) T_f}{\mathcal{R}_c+\Delta z_f/\lambda_f +\Delta z_s/\lambda_s}, \end{gather}

where $T_s$ is the temperature in the solid at distance $\Delta z_s$ from the boundary, and $T_f$ is the temperature in the fluid at distance $\Delta z_f$ from the boundary; this could be either the distance to the nearest mesh point or that to the phasic interface.

For the purpose of implicit discretisation of the Laplacian operator in (2.16) and (2.17), we start with a generic, three-point central-difference formula:

(2.22)\begin{equation} \frac{\partial^2 T}{\partial z^2}\approx K_P (T_P-T_C) + K_B (T_b-T_C). \end{equation}

Here, $T_C$ is the temperature at the point where the difference is calculated; for illustration, we choose the solid-adjacent fluid cell, shown as the central cell in figure 2. The coefficients $K_P$ and $K_B$ are functions of the distances to the solid–fluid boundary $\Delta z_p$ and $\Delta z_b$, also illustrated in figure 2. Assuming that $T_b$ has been computed from $T_M$ and $T_C$ using (2.21), we can then develop the second term on the right-hand side of (2.22) as follows:

(2.23)\begin{align} K_B (T_b - T_C) &= K_B\left(\frac{\mathcal{R}_M T_C + \mathcal{R}_C T_M}{\mathcal{R}_M+\mathcal{R}_C}-T_C\right) = K_B\left(\frac{\mathcal{R}_C T_M-\mathcal{R}_C T_C}{\mathcal{R}_B+\mathcal{R}_C}\right)\nonumber\\ &= K_B\frac{\mathcal{R}_C}{\mathcal{R}_M+\mathcal{R}_C} (T_M - T_C)\equiv K_M (T_M - T_C). \end{align}

Figure 2.

Schematic representation of the stencil used for the implicit discretisation of the heat diffusion term of the energy transport equation in the vicinity of a solid–fluid boundary, indicated here by the thick vertical line: $T_{b,s}$ and $T_{b,f}$ are the solid-side and liquid-side boundary temperatures, respectively.

Similar estimations are used within the solid finite-volume formulation on the solid side of the boundary. Thus, the conjugate heat-transfer problem between the solid and the fluid can be solved by means of a fully coupled algorithm, with heat conduction treated implicitly on both the fluid and solid sides of the boundary. Note that, in our application domain, the physical solid–liquid contact resistance is assumed to be zero; we use $\mathcal {R}_c$, however, for reduced-order modelling of IHTR. This is discussed further in § 3.5.

2.5.2. Wall boundary conditions for interface reconstruction

Our subgrid-accurate implementation of phase change using (2.5) necessitates the calculation of the interfacial area density $a_\gamma$ in all wall-adjacent cells. The nodal values of $\phi$ required by the marching cubes algorithm (Lorensen & Cline Reference Lorensen and Cline1987) necessitate the use of ghost values in the wall cells. These are obtained using the following iterative extrapolation procedure:

1. (i) Perform the GVOF reconstruction of the interfacial geometry in wall-adjacent cells, utilising only information in the fluid cells themselves. Here, GVOF reconstruction means computing normal vectors and line constants defining the GVOF interfacial orientation based on the $\phi$ stencil described in Bureš et al. (Reference Bureš, Sato and Pautz2021).

2. (ii) Extend the linear interface from wall-adjacent cells to the ghost cells to obtain $\phi ^{{ghost}}$.

3. (iii) Perform the GVOF reconstruction of the interfacial geometry in wall-adjacent cells again, utilising the full stencil this time, and including the ghost-cell estimates. Note that the normal vector calculated at this step could be different from the one obtained in step  (i).

4. (iv) Compare the new normal vector field with the previous estimate. If the change is outside the specified tolerance, go back to step  (ii) and repeat.

In practical terms, the method converges very rapidly, requiring only one or two iterations to reach negligible variation of the normal vector field, as stipulated in the predefined tolerances. A schematic representation of the result of the wall iterative extrapolation procedure is presented in figure 3.

Figure 3.

Schematic representation of the wall iterative extrapolation procedure. Note that the shape of the interface in the fluid domain is slightly perturbed in the process due to realignment of the normal vector. Note also that the linear extension of the interface is limited only to the adjacent ‘ghost’ cell.

2.6. Overall solution algorithm

In summary, the following algorithm is employed to solve the governing equations in PSI-BOIL:

1. (i) Calculate the interfacial area density $a_\gamma$ using the marching cubes algorithm needed for (2.5).

2. (ii) Calculate the mass transfer $\dot {m}'''$ between the phases directly using (2.5).

3. (iii) Calculate the curvature $\kappa$ using the height-function method, and use it to estimate the surface tension force $\,\boldsymbol {f}_\sigma$ via (2.11).

4. (iv) Solve the momentum conservation equations (2.9) to obtain a tentative estimate of the velocity field $\boldsymbol {u}^\star$.

5. (v) Solve the Poisson equation for pressure, (2.13), and project the tentative velocity field $\boldsymbol {u}^\star$ to ensure volume conservation using (2.14), obtaining a new estimate of the velocity field $\boldsymbol {u}$.

6. (vi) Calculate the intermediate volume fraction field using (2.3).

7. (vii) Calculate the liquid and vapour velocity fields $\boldsymbol {u}_l$ and $\boldsymbol {u}_v$ (individual phasic velocities) using our version of the extrapolation algorithm of Malan et al. (Reference Malan, Malan, Zaleski and Rousseau2021), as detailed in Bureš & Sato (Reference Bureš and Sato2021b).

8. (viii) Advect the intermediate volume fraction field using directional splitting, (2.4).

9. (ix) Ensure energy conservation by the simultaneous solution of (2.16), (2.17) and (2.19).

10. (x) Advance the time step.

11. (xi) Go back to step (i).

3.1. Discrepancy

Although the above assumption is generally acceptable for single-component multiphase systems (Ishii & Hibiki Reference Ishii and Hibiki2011), it is known that it can induce non-negligible error in the calculation of heat flow through the microregion and microlayer of an expanding bubble from a heated surface (Giustini et al. Reference Giustini, Jung, Kim and Walker2016). For illustration, consider a typical data point from the experiment of Bucci (Reference Bucci2020), with a measured wall superheat $\Delta T_{{wall}} \!=\! T_{{wall}} - T_{{sat}}$ in the range of 5–10 K. Though the microlayer thickness $d$ (m) was not explicitly measured in this experiment, during evaporation, values of $d = O(100~\text {nm})$ would be expected, as indicated by other recent experiments (e.g. Chen et al. Reference Chen, Hu, Hu, Utaka and Mori2020). Assuming steady-state 1-D conduction through the liquid microlayer beneath the expanding bubble (see § 4 for a discussion of the validity of this assumption), this gives for the heat flux $j_q$ (W m$^{-2}$):

(3.1)\begin{equation} j_q = \lambda_l\frac{\Delta T_{{wall}}}{d} = O(10~\text{MW}~\text{m}^{{-}2}). \end{equation}

Such values should thus be routinely observed during an experiment of this type, but the maximal values reported by Bucci (Reference Bucci2020) are in the region of 2 MW  m$^{-2}$ in the 5–10 K superheat range. Such a discrepancy is not exceptional, and has often been reported in the literature (e.g. Giustini et al. Reference Giustini, Jung, Kim and Walker2016), and appears to indicate that the simple assumption of thermostatic equilibrium does not capture the actual heat-transfer situation occurring in the microlayer.

3.2. Closure using statistical physics

An alternative closure model might be found in the framework of statistical physics. The Hertz–Knudsen relation expresses the interfacial heat flux in the following form (Knudsen Reference Knudsen1950):

(3.2)\begin{equation} j_{q,\gamma} = L\sqrt{\frac{M_v}{2{\rm \pi} R_g}}\left(\omega_e\frac{p_s(T_{\gamma,l})}{\sqrt{T_{\gamma,l}}} -\omega_c\frac{p_v}{\sqrt{T_{\gamma,v}}}\right). \end{equation}

Here, $M_v$ is the molar mass $({\rm kg}~{\rm mol}^{-1})$ of vapour, $R_g = 8.314~{\rm J}~({\rm mol}~{\rm K})^{-1}$ is the universal gas constant, $p_s(T_{\gamma,l})$ is the saturation pressure at the liquid interfacial temperature $T_{\gamma,l}$, while $p_v$ is the vapour pressure and $T_{\gamma,v}$ is the vapour interfacial temperature. In addition, $\omega _c$ and $\omega _e$ (–) represent the empirical accommodation coefficients for condensation and evaporation, introduced by Knudsen (Reference Knudsen1950) to account for the imperfect transport of molecules between the phases. Using statistical rate theory, Persad & Ward (Reference Persad and Ward2016) deduced closed-form expressions for $\omega _c$ and $\omega _e$ and demonstrated that

(3.3)\begin{equation} \omega_c \approx \omega_e \approx 1. \end{equation}

For small temperature variations in the domain, we can thus simplify the expression (3.2) accordingly:

(3.4)\begin{equation} j_{q,\gamma} \approx L\sqrt{\frac{M_v}{2{\rm \pi} R_g T_{{ref}}}} (p_s(T_{\gamma,l})-p_v), \end{equation}

with $T_{{ref}}$ being a representative temperature value, such as $T_{{sat}}$. This equation includes two independent parameters: the liquid interfacial temperature $T_{\gamma,l}$ and the vapour pressure $p_v$. Typically, it is assumed that the vapour is at saturation at the prevailing system pressure (Giustini et al. Reference Giustini, Jung, Kim and Walker2016), giving

(3.5)\begin{equation} j_{q,\gamma} \approx L\sqrt{\frac{M_v}{2{\rm \pi} R_g T_{{ref}}}} (p_s(T_{\gamma,l})-p_s(T_{{sat}})). \end{equation}

Using the linearised form of the Clausius–Clapeyron relation (Faghri & Zhang Reference Faghri and Zhang2006), we immediately deduce that

(3.6)\begin{equation} j_{q,\gamma} \approx \rho_v L^2 \sqrt{\frac{M_v}{2{\rm \pi} R_g T^3_{{sat}} }} \,(T_{\gamma,l}-T_{{sat}}) =\frac{T_{\gamma,l}-T_{{sat}}}{\mathcal{R}_\gamma}. \end{equation}

Here, we have introduced the factor $\mathcal {R}_\gamma$ $({\rm K~m}^2~{\rm W}^{-1})$ as the IHTR:

(3.7)\begin{equation} \mathcal{R}_\gamma = \frac{1}{\rho_v L^2} \sqrt{\frac{2{\rm \pi} R_g T_{{sat}}^3}{M_v}}. \end{equation}

However, this expression, when applied to water under atmospheric conditions, gives $\mathcal {R}_\gamma \simeq 1.3 \times 10^{-7}~{\rm K~m}^2~{\rm W}^{-1}$. This value is too small to change the results of the calculations of heat flux in the microlayer, such as those represented in (3.1). As discussed by Persad & Ward (Reference Persad and Ward2016), assumptions on vapour conditions can significantly impact the accuracy of any expression deduced from the Hertz–Knudsen relation (3.2).

To the best of our knowledge, a satisfactory closure model describing vapour conditions in the vicinity of the microlayer beneath a growing vapour bubble during boiling has not yet appeared in the literature, and thereby represents an unknown in the understanding of the physical processes taking place, which should be addressed in the future. For the present work, we are thus forced to resort to the classical approach, in which we assume that $p_v=p_s(T_{{sat}})$, and, although $\omega _c\approx \omega _e\equiv \omega$, the accommodation coefficient is in fact not taken to be unity, and is a priori unknown. This gives, following the same steps as above,

(3.8)\begin{equation} \mathcal{R}_\gamma = \frac{1}{\omega}\frac{1}{\rho_v L^2}\sqrt{\frac{2{\rm \pi} R_gT_{{sat}}^3}{M_v}}, \end{equation}

leaving the unknown $\omega$ yet to be determined.

3.3. Value of the accommodation coefficient

To evaluate $\omega$, we appeal to the available experimental data. The dataset of Bucci (Reference Bucci2020) includes measurements of both wall superheat and heat flux, with a spatial resolution of 30 $\mathrm {\mu }$m, and a temporal resolution (i.e. exposure time of the camera) of 200 $\mathrm {\mu }$s. Employing the same 1-D steady-state conduction model as used before (see figure 23 in § 4 and the corresponding discussion for a justification of the use of this model), we can collate the given measurements according to the following formula:

(3.9)\begin{equation} j_q = \frac{\Delta T_{{wall}}}{d/\lambda_l+\mathcal{R}_\gamma}. \end{equation}

If we now assume that, for each temporal snapshot, the highest locally measured value of $j_q$ corresponds to the situation for which $d\rightarrow 0$, we can deduce that, for each snapshot,

(3.10)\begin{equation} \mathcal{R}_\gamma \approx \frac{\Delta T_{{wall,loc}}}{j_{q,{max}}}, \end{equation}

where $\Delta T_{{wall,loc}}$ is the corresponding locally measured superheat value (see figure 4).

Figure 4.

Schematic representation of the $\mathcal {R}_\gamma$ evaluation procedure. Here, $j_{q,{max}}$ is the maximal heat flux for the given snapshot, and $T_{{wall,loc}}$ is the locally measured wall temperature.

Figure 5 shows the result of this processing of the available data. It can be observed that, despite deviations from the norm, a justifiable prediction of $\mathcal {R}_\gamma$ can be estimated, with a mean value of $\bar {\mathcal {R}}_\gamma \simeq 3.7 \times 10^{-6}~{\rm K~m}^2~{\rm W}^{-1}$. Note that this is at least one order of magnitude larger than the classical value ($\mathcal {R}_\gamma \simeq 1.3 \times 10^{-7}~{\rm K~m}^2~{\rm W}^{-1}$) based on the simplified Hertz–Knudsen relation (3.6). This new estimate corresponds to a value of the accommodation coefficient $\omega _A$ in (3.8) of $\sim$0.0345, a value consistent with the literature suggestions for water evaporation (Paul Reference Paul1962) and boiling (Giustini et al. Reference Giustini, Jung, Kim and Walker2016). Consequently, we have adopted this value for the remainder of the study presented here. Note that the equivalent liquid film thickness $d_{{eq}} = \lambda _l \bar {\mathcal {R}}_\gamma$ is approximately 2 $\mathrm {\mu }$m; i.e. the resistance of the interface has essentially the same importance as the conductive resistance of the microlayer itself, this being equal to $d/\lambda _l=O(1~\mathrm {\mu }\text {m})/\lambda _l = O(10^{-6}~\text {K~m}^2~\text {W}^{-1})$.

Figure 5.

The IHTR evaluated using (3.10), and based on the data of Bucci (Reference Bucci2020). The mean value has been obtained by an averaging procedure with respect to time, rather than one based on spatial position.

In the above analysis, it has been assumed that the experimental measurement technique can capture the location of the heat-flux extremum exactly. However, the IR camera used by Bucci (Reference Bucci2020) has finite spatial ($\Delta x = 30~\mathrm {\mu }$m) and temporal ($\Delta t = 200~\mathrm {\mu }$s) resolutions. Thus, $d$ in (3.9) should not be taken to be equal to zero, but rather has a value of

(3.11)\begin{equation} d = \frac{1}{\Delta t \Delta x} \int_{\Delta t}\int_{\Delta x} d(x',t') \,{\rm d}t' \,{{\rm d}\kern0.7pt x}', \end{equation}

the integrals representing averaging over a time step $\Delta t$ and a spatial segment with width $\Delta x$. Thus, we see that $\omega _A = 0.0345$ is a lower-bound estimate for the accommodation coefficient under the assumption $d\rightarrow 0$ in the above equation, and the actual value of $\mathcal {R}_\gamma$ should be given as

(3.12)\begin{equation} \mathcal{R}_\gamma = \frac{\Delta T_{{wall,loc}}}{j_{q,{max}}} - \frac{d}{\lambda_l} = \mathcal{R}_\gamma(d=0) - \frac{d}{\lambda_l} = \mathcal{R}_{\gamma,A} - \frac{d}{\lambda_l}. \end{equation}

To account for the uncertainty in $\mathcal {R}_\gamma$ due to the unknown value of $d$, we should also deduce an upper bound $\omega _B$. We have noted in our simulations in Bureš & Sato (Reference Bureš and Sato2021a,Reference Bureš and Satoc), as well as those presented in § 4, that the microlayer thickness does not gradually increase from zero near the contact line, but rather jumps abruptly to a non-zero value, $d_1$, as a result of the dynamics of dewetting. Thus, assuming the thickness of the microlayer to have a step-like profile in the vicinity of the contact line, the integral (3.11) reduces to $d = d_1$. Based on our preliminary simulations, we have chosen this value as $\sim$500 nm, a value we consider high enough to result in a reasonable upper-bound estimate $\omega _B$. Indeed, setting $d$ in (3.12) as $d_1 \simeq 500$ nm and evaluating $\mathcal {R}_\gamma$ corresponds to an increase of $\omega$ by approximately 33 %, giving $\omega _B \simeq 0.0460$. Although not derived in a fully rigorous manner, this upper bound should account for the uncertainty in the $\mathcal {R}_\gamma$ evaluation discussed here; results with both $\omega _A = 0.0345$ and $\omega _B = 0.0460$ are presented in § 4.

3.4. Full model implementation

To implement (liquid-side) IHTR into the PSI-BOIL code, the Dirichlet condition (2.18) must be replaced by a Robin-type (Papac, Gibou & Ratsch Reference Papac, Gibou and Ratsch2010) condition deduced from (3.6), with the interfacial heat flux evaluated using Fourier's law of heat conduction:

(3.13)\begin{gather} T_{\gamma,v} = T_{{sat}}\quad \text{at the vapour side of the interface}, \end{gather}

(3.14)\begin{gather}\frac{T_{\gamma,l}-T_{{sat}}}{\mathcal{R}_\gamma} = \boldsymbol{j}_{q}|_{\gamma,l}\boldsymbol{\cdot}\boldsymbol{n} = \lambda_l\boldsymbol{\nabla} T|_{\gamma,l}\boldsymbol{\cdot}\boldsymbol{n}\quad \text{at the liquid side of the interface}. \end{gather}

We further assume that the liquid-side heat flux in the direction tangential to the interface is negligible (which is typical both for boundary-layer situations at the interface within the bulk fluid and for 1-D conduction in the microlayer). Hence, $\,\boldsymbol {j}_{q}|_{\gamma,l} \approx j_{q}|_{\gamma,l}\,\boldsymbol {n}$, which is equivalent to the approximation introduced by Liu, Fedkiw & Kang (Reference Liu, Fedkiw and Kang2000) in the context of their solution approach for Poisson equations. This assumption immediately implies that

(3.15)\begin{equation} j_{q}|_{\gamma,l} = \frac{T_{\gamma,l}-T_{{sat}}}{\mathcal{R}_\gamma} \end{equation}

and

(3.16)\begin{equation} \boldsymbol{j}_{q}|_{\gamma,l} = \begin{pmatrix}j_{q,x}\\j_{q,z}\end{pmatrix}\bigg|_{\gamma,l} = \frac{T_{\gamma,l}-T_{{sat}}}{\mathcal{R}_\gamma}\begin{pmatrix} n_x\\n_z\end{pmatrix}. \end{equation}

Evidently, we can now decompose the Robin condition into two independent directional conditions as follows:

(3.17)\begin{gather} \frac{\partial T}{\partial x}\bigg|_{\gamma,l} = \frac{T_{\gamma,l}-T_{{sat}}}{\lambda_l\mathcal{R}_\gamma}n_x, \end{gather}

(3.18)\begin{gather}\frac{\partial T}{\partial z}\bigg|_{\gamma,l} = \frac{T_{\gamma,l}-T_{{sat}}}{\lambda_l\mathcal{R}_\gamma}n_z. \end{gather}

These conditions can be implemented into the implicit, sharp-interface discretisation scheme employed in PSI-BOIL to solve the energy conservation equation, by utilising conceptually similar methods to those described earlier in § 2.5.1, even though the resulting formulae are rather cumbersome.

Our preliminary results indicate that introducing the IHTR everywhere on the bubble surface results in an unrealistically small growth rate. This would indicate that IHTR only plays a significant role in the immediate vicinity of the heated surface. Although this might appear to be unphysical, i.e. we would expect that the interface has the same heat-transfer resistance everywhere, it is important to remember that the very concept of IHTR is speculative in itself. Indeed, referring back to (3.4), we should actually be concerned with the near-interface vapour conditions, which could vary around the bubble surface. Since we do not have any information on the vapour pressure distribution, we have decided to continue our numerical exercise using the approach we introduced in our conference paper (Bureš & Sato Reference Bureš and Sato2021a): i.e. we localise the effects of IHTR only in the microlayer, and then critically assess its success, or otherwise, against measured data. This approach is detailed below, while an example of results of a simulation with IHTR prescribed everywhere on the bubble surface is given in § 5, for comparison purposes, and where results of a simulation with $\omega =1$ are also given.

3.5. Near-wall implementation

Since any significant effect of IHTR is now assumed to be localised in the microlayer, for which conduction in the wall-normal direction represents the dominant heat-transfer mechanism, a simplified implementation of IHTR can be introduced directly at the solid wall boundary by means of a numerically equivalent, contact heat-transfer resistance factor $\mathcal {R}_c$ (see § 2.5.1). This approach automatically preserves the overall conduction resistance, since

(3.19)\begin{equation} j_q = \frac{T_{b,s}-T_{\gamma,l}}{d/\lambda_l} = \frac{T_{b,s}-T_{{sat}}}{d/\lambda_l + \mathcal{R}_\gamma} \equiv \frac{T_{b,s}-T_{{sat}}}{\mathcal{R}_c + d/\lambda_l}, \end{equation}

provided that $\mathcal {R}_c=\mathcal {R}_\gamma$. This simplification is illustrated schematically in figure 6. As can be seen, the temperature discontinuity due to heat-transfer resistance has been effectively moved from the liquid–vapour interface to the solid–liquid boundary; the overall temperature difference between the wall temperature $T_{b,s}$ and the vapour temperature $T_{{sat}}$ has been preserved in the process.

Figure 6.

Illustration of the concept of equivalent conductive resistance: (a) IHTR located at the liquid–vapour interface; and (b) IHTR located at the solid–liquid boundary. Note that $T_{\gamma,l}-T_{{sat}} = T_{b,s}-T_{b,l}$.

3.6. Apparent contact angle in the high-IHTR limit

In spite of the reported importance of the wetting conditions for microlayer formation (Urbano et al. Reference Urbano, Tanguy, Huber and Colin2018; Bureš & Sato Reference Bureš and Sato2021c), relevant experimental data remain sparse. Although static wetting conditions are usually reported – e.g. perfect wetting in the case of Bucci (Reference Bucci2020) – it is generally accepted (Raj et al. Reference Raj, Kunkelmann, Stephan, Plawsky and Kim2012; Janeček & Nikolayev Reference Janeček and Nikolayev2013; Fourgeaud et al. Reference Fourgeaud, Ercolani, Duplat, Gully and Nikolayev2016; Schweikert et al. Reference Schweikert, Sielaff and Stephan2019) that the presence of phase change, together with the motion of the contact line (see figure 1), significantly affect the wetting phenomenon. Note that in the case of DNS-BRM, the effects of the numerical slip length, as described by Afkhami et al. (Reference Afkhami, Buongiorno, Guion, Popinet, Saade, Scardovelli and Zaleski2018) and Bureš & Sato (Reference Bureš and Sato2021c), should also be taken into account, further complicating the issue.

The traditional approach to the problem of dynamic wetting of volatile liquids is to appeal to the so-called microregion model: a set of equations derived from 2-D lubrication theory extended to include heat transfer by including 1-D steady-state conduction through the film thickness. The classical result of Voinov (Reference Voinov1976) and Cox (Reference Cox1986) describes dynamic wetting under adiabatic conditions in the form of the Cox–Voinov law thus:

(3.20)\begin{equation} \theta(x)^3-\theta_0^3 = 9Ca_u \ln(x/l_s) = 9Ca_u \ln\chi . \end{equation}

Here $\theta (x)$ is the local interfacial slope, and $\theta _0$ is the microscopic contact angle, conventionally assumed to be constant (i.e. any hysteresis effects are neglected), and derived from force balance arguments at the triple line (Snoeijer et al. Reference Snoeijer, Andreotti, Delon and Fermigier2007). In addition, $l_s$ (m) is the slip length, $\chi =x/l_s$ is a dimensionless distance and $Ca_u$ is the standard capillary number,

(3.21)\begin{equation} Ca_u = \frac{\mu_l U_{{CL}}}{\sigma}, \end{equation}

with $U_{{CL}}$ being the contact-line velocity, taken positive for spreading of the liquid over the surface and negative for dewetting. In a like manner, Hocking (Reference Hocking1995) gives a closed-form solution for the apparent contact angle in the presence of phase change in the form

(3.22)\begin{equation} \theta(x)^4-\theta_0^4 = 12Ca_e\ln(K\chi), \end{equation}

where $Ca_e$ is the evaporation capillary number of Morris (Reference Morris2001), namely,

(3.23)\begin{equation} Ca_e = \frac{\mu_l\Delta T_{{wall}}}{\rho_l L\mathcal{R}_\gamma\sigma}. \end{equation}

Note that the parameter $K$ in the argument of the logarithm in (3.22) results from the special matching condition for droplets originally introduced by Hocking (Reference Hocking1995).

The asymptotic solution thus produced was shown to be valid in the high-IHTR limit by Janeček & Anderson (Reference Janeček and Anderson2016). In both these studies, the motion of the contact line was not taken into account. However, contact-line motion and evaporation were considered together in an earlier work by Janeček & Nikolayev (Reference Janeček and Nikolayev2013), who showed that the apparent contact-angle ($\theta _{mr}$ in figure 1) could be obtained from a combination of an inner numerical solution of the governing equations pertaining to the static microregion with phase change, and an outer, Cox–Voinov-type expression describing the dependence on velocity, overall giving

(3.24)\begin{equation} \theta_{outer}(\chi)^3-\theta_{inner}(\Delta T)^3 = 9Ca_u \ln\chi, \end{equation}

with $l_s$ from (3.20) replaced by a typical scale of the microregion, of the order of 10–100 nm (Janeček & Nikolayev Reference Janeček and Nikolayev2013).

The existence of analytical solutions at the high-IHTR limit (Hocking Reference Hocking1995; Morris Reference Morris2001), which is precisely the situation considered here, has motivated us to deduce an analytical solution for $\theta _{mr}(\Delta T,U_{{CL}})$ combining phase change and contact-line motion directly. Details are given in Appendix A, and the result may be expressed succinctly in the following form:

(3.25)\begin{equation} \theta_{mr} \simeq \frac{1}{|A|} = \frac{Ca_e}{|Ca_u|} = \frac{\Delta T_{{wall}}}{\rho_l|U_{{CL}}| L\mathcal{R}_\gamma}. \end{equation}

In general, the lower the value of the contact angle, the more is the contact-line mobility inhibited, and the role of hydrodynamic effects in microlayer destruction diminished: in other words, low values of the contact angle cause the microlayer to be destroyed predominantly by evaporation.

We have demonstrated this effect ourselves in our previous paper (Bureš & Sato Reference Bureš and Sato2021c). By employing (3.25), the resulting contact angle corresponds to a dynamic equilibrium situation (since $\theta _{mr}$ is influenced by $U_{{CL}}$, and vice versa), and the predicted values of the contact angle based on this equation have not exceeded $1^\circ$ in the simulations we have carried out, and will present in § 4. The resulting reduced contact-line mobility could be the cause of the discrepancy between experimental and calculated initial microlayer thickness (IMT) distributions, as discussed below, indicating that $\theta _{mr}$ is too small to represent the experimental conditions. Even though the exact value of the contact angle to be imposed in the simulations thus remains uncertain, it can be seen from (3.25) that a high value of $\mathcal {R}_\gamma$ leads to a decrease in the interfacial slope in the microregion, as measured by $\theta _{mr}$, since it is inversely proportional to $\mathcal {R}_\gamma$. Thus, it appears that the importance of the dewetting phenomena in the microlayer regime is somehow diminished during boiling. This is inconsistent with the results obtained from experiments on volatile liquid films, where dewetting, as reported by Fourgeaud et al. (Reference Fourgeaud, Ercolani, Duplat, Gully and Nikolayev2016), appears to play a decisive role in film destruction. It is possible that IHTR in volatile-film experiments might be less influential than during nucleate boiling, resulting in the increase in contact angle due to evaporation being much greater than in the situation considered here, and ultimately promoting the effects of dewetting, though this interpretation could be considered speculative.

4.1. Simulation set-up

Since the experimental observations report essentially perfect axial symmetry of the growing bubble, we have employed a cylindrical, axisymmetric domain for our simulations, incorporating three different grid resolutions (see table 1); the problem domain is illustrated schematically in figure 7. To reduce computational costs, both dimensions have been discretised uniformly only in the central region of the domain, indicated by darker colour in figure 7; a stretched grid has been employed for the remainder of the domain, with a gradual transition resulting in a maximum cell aspect ratio of 4.

Figure 7.

Schematic representation of the computational domain used for the validation simulations based on the experiment of Bucci (Reference Bucci2020), with $x$ representing the radial and $z$ the axial coordinate in an axisymmetric cylindrical system, respectively. The dimensions of the fluid domain are approximate and depend on the grid resolution used. Units in millimetres.

Table 1.

Domain characteristics for the grids considered: $\Delta x_{{min}}$ is the grid spacing in the uniformly discretised, fine-grid region of the domain (see figure 7) of dimensions $x_{uni}$ and $z_{uni}$. For the coarse grid, a typical CPU time is $\sim$20 000 core-hours; for the fine grid, it is $\sim$400 000 core-hours.

To facilitate faster convergence of the algebraic multigrid solver (Brandt & Livne Reference Brandt and Livne1984) used for the solution of the Poisson equation for pressure (2.13), and to economise on computational effort, the number of cells for each grid resolution was first chosen, and then the resultant domain size was computed. For this reason, the radial and axial extents of the fluid domain vary slightly for the different grid resolutions. Note that, for the coarse-mesh simulation, 70 cores have been used and, for the fine-mesh simulation, 336 cores have been used. Taking into account the increased number of time steps due to the reduction of step length, the performance of the code worsens by a factor of 6–7 between these two meshes, even though each core handles the same amount of information (i.e. $128\times 128$ computational cells). The main reason for this poor scaling is the worsening convergence of the pressure-Poisson equation multigrid solver. In future work, we will attempt to improve the computational efficiency of the code by introducing more advanced multigrid implementations. The medium-mesh simulation has been run on an older cluster and is, therefore, not considered in these comparisons.

With reference to figure 7, an axis-of-symmetry BC has been applied at $x = x_{min}$, and a free outflow BC has been imposed at $x = x_{max}$ and $z = z_{max}$, with a Neumann BC (zero heat flux) for temperature. A zero heat-flux BC has also been applied at the bottom of the substrate at $z = z_{min}$, where the material was in contact with air in the experimental set-up. The heater surface is treated as a no-slip wall. The titanium heater is included explicitly in the simulation, and the conjugate heat transfer between the fluid and solid phases is included in the calculation. A vapour seed bubble of radius $R_0=10~\mathrm {\mu }$m is placed at the origin of coordinates to initiate the bubble growth; see figure 7. The material properties used in the simulations are given in table 2.

Table 2.

Material properties used in our simulations up to three significant digits: $\beta$ is the coefficient of volumetric thermal expansion employed for the computation of the gravity force (Boussinesq approximation).

The initial temperature distribution was first obtained by solving the transient heating problem in an extended configuration of the domain, with uniform grid spacing $\Delta x=\Delta z = 1~\mathrm {\mu }$m: grid convergence was confirmed. The simulation was initiated under assumed saturation conditions, and continued until the nucleation temperature at the origin was attained. No discernible effects of natural convection could be seen at this stage. Although the experimental bubble-growth results were reported by Bucci (Reference Bucci2020) as fully axisymmetric, the heating power was not, due to the rectangular geometry of the heater. In addition, non-uniformity of the heating power was also reported. Thus, it may be seen that prescribing the heat flux uniformly within the heater at $425~{\rm kW}~{\rm m}^{-2}$ during the heating phase has resulted in the calculated surface temperature distribution not exactly matching that of the experiment, though the time to nucleation has been captured adequately.

For this reason, we have introduced an augmented heater power distribution, devised to reproduce both the time to nucleation and the transient temperature data reported by Bucci (Reference Bucci2020), with a special focus on the final surface temperature distribution at the onset of bubble growth. The applied power distribution is shown in figure 8, together with the step profile used originally. The two may be expressed as follows:

(4.1)\begin{gather} j_{step}\ ({\rm kW}~{\rm m}^{{-}2}) = \begin{cases} 425, & \text{if}\ x< x_{max}=\dfrac{1+\sqrt{2}}{2} \times 1.5 \text{mm}, \\ 0, & \text{otherwise}, \end{cases} \end{gather}

(4.2)\begin{gather}j_{lin}\ ({\rm kW}~{\rm m}^{{-}2}) = \max\left(\frac{4~\text{mm}-x~(\text{mm})}{4~\text{mm}}\times 481, 0\right). \end{gather}

The resulting temperature distributions at the onset of nucleation are shown in figure 9. Clearly, the linear-flux profile matches the experiment perfectly.

Figure 8.

Applied heat flux during the pre-nucleation transient heating problem, (4.1) and (4.2).

Figure 9.

Heater surface temperature at the onset of nucleation from the experiment and simulation with different applied heat fluxes. Times to nucleation are indicated in the legend.

In addition, figure 10 gives comparisons of the measured superheat over the heater surface against those derived numerically with the assumed linear heat-flux profile; again, the correspondence is excellent. The heater extent for the step distribution ($x_{max}$ in (4.1)) was chosen based on the experimental heater half-width of 1.5 mm. But, even if a different value had been selected, the qualitative shape of the temperature distribution measured by Bucci (Reference Bucci2020) could not be reproduced. Nevertheless, the ‘engineered’ approach, represented by (4.2), has allowed us to obtain an initial temperature distribution exactly matching that measured in the experiment, a prerequisite for the success of any subsequent validation exercise focused on bubble growth. Note that, during the simulation, a uniform heat flux has been applied; the short elapsed time of the simulation makes its impact on the numerical predictions essentially negligible. The temperature distribution in the domain at the onset of nucleation ($t=88.1$ ms) is displayed in figure 11. It should also be remarked that, due to the short heating time, the substrate is not uniformly superheated.

Figure 10.

Evolution of the heater surface temperature before the onset of nucleation at $t=88.1$ ms; simulated results obtained with the linear power distribution, (4.2).

Figure 11.

Shaded contours of the initial temperature distribution in the central region of the domain at the onset of nucleation. The white line indicates the solid–fluid boundary. Units in micrometres.

Throughout the simulation, a variable time step $\Delta t$ has been employed, with the upper limit imposed by the Courant number,

(4.3)\begin{equation} {CFL} = \frac{u_{{max}}\Delta t}{\Delta s} < 0.02, \end{equation}

where $u_{{max}}$ is the maximum directional velocity in the domain (i.e. in either the $x$- or $z$-direction), and $\Delta s$ is the corresponding local grid spacing in the given direction. A second upper limit is represented by the capillary-wave condition, as recommended by Brackbill et al. (Reference Brackbill, Kothe and Zemach1992) and Popinet (Reference Popinet2018):

(4.4)\begin{equation} \Delta t< C_\sigma \sqrt{\frac{(\rho_v+\rho_l) \Delta x_{{min}}^3}{\sigma}}, \end{equation}

with $C_\sigma =0.063$, a value suggested by Sussman & Ohta (Reference Sussman and Ohta2009). The absolute minimum of these two criteria is the time step actually used in the simulations. It transpires that, during the vigorous initial growth of the bubble, the ${CFL}$ condition dominates, while in the latter stages of the expansion, the capillary-wave condition plays the limiting role. It should be noted that both conditions are more stringent than usually employed in DNS computations. Indeed, we have typically used ${CFL}=0.05\text {--}0.25$ and $C_\sigma =0.10\text {--}0.25$ in our previous single-bubble, nucleate-boiling simulations (Sato & Niceno Reference Sato and Niceno2015; Bureš & Sato Reference Bureš and Sato2020). The need to limit the time step in a more stringent manner is a direct consequence of the fast propagation of the liquid–vapour interface through the fluid domain in the current simulation, and the very small grid spacing employed, both of these factors negatively impacting the stability of the simulation, and necessitating tighter time-step control.

Note that, in theory, the explicit treatment of the mass source (2.3) should result in another time-stepping restriction being introduced. However, implementation of $\dot {m}'''$ can be recast as motion of the interface with a velocity with magnitude $\dot {m}''/\rho _l$ (Malan et al. Reference Malan, Malan, Zaleski and Rousseau2021). This is in turn negligible in comparison with the Stefan velocity $\dot {m}''(1/\rho _v-1/\rho _l)$, which is guaranteed to occur somewhere in the computational domain during a bubble-growth simulation. As a result, the CFL condition should, in practice, be sufficient to guarantee stability of the explicit volume-of-fluid (VOF) scheme.

Preliminary simulations have indicated that a non-smooth distribution of the phase-change rate over the bubble surface, occasioned by the thinness of the interfacial temperature boundary layer in this application, has induced spurious capillary waves during the initial stages of the growth if the two finer grids in table 1 are employed. For this reason, we have introduced a global phase-change rate averaging procedure into the simulations, which remains active until the bubble radius exceeds 150 $\mathrm {\mu }$m, corresponding to 20–30 $\mathrm {\mu }$s after the onset of bubble growth, i.e. approximately 5 % of the total simulation time. Specifically, this entails the phase-change rate being calculated over the entire fluid domain, summed and uniformly redistributed along the interface. Note that, for water under atmospheric conditions and a superheat of 12.55 K, and based on conventional estimates of bubble growth predicted by analytical formulae for the heat-transfer-limited and inertial regimes (Faghri & Zhang Reference Faghri and Zhang2006), the bubble growth should be considered inertia-dominated for radii up to 200 $\mathrm {\mu }$m. Thus, it could be argued that stabilisation by global averaging employed in our simulations is only used in the period for which the results of our heat-transfer-limited, incompressible solver are anyway deviating from physical reality. A critical assessment of the effects of the averaging procedure (for the coarse-grid simulation) is presented in § 5.

The discussion above hints at limitations in the application of the sharp-interface GVOF technique when applied to volatile flows involving rapid propagation of the liquid–vapour interface. This is not surprising: indeed, the problem of bubble growth in a superheated, quiescent liquid as discussed by Scriven (Reference Scriven1959), which could be considered an antecedent of the problem of single-bubble nucleate boiling, has been successfully simulated using the sharp-interface GVOF method only very recently (Perez-Raya & Kandlikar Reference Perez-Raya and Kandlikar2018; Bureš & Sato Reference Bureš and Sato2021b; Malan et al. Reference Malan, Malan, Zaleski and Rousseau2021). Nevertheless, the GVOF method does offer several advantages over alternative approaches: principally, the ability to capture microlayer formation even with a coarse grid resolution ($\sim$1 $\mathrm {\mu }$m), as will be shown in § 4.2.

6.1. Modelling of the initial microlayer thickness

Typical analytical expressions for the IMT, i.e. $d_0$, distribution found in the literature are variations of the original square-root law proposed by Cooper & Lloyd (Reference Cooper and Lloyd1969):

(6.1)\begin{equation} d_0(t) \propto \sqrt{\frac{\mu_l}{\rho_l}t\,} = \sqrt{\nu_l t}, \end{equation}

with $\nu _l$ being the kinematic viscosity of the liquid $({\rm m}^2~{\rm s}^{-1})$; in addition to the referenced paper, see the works of van Ouwerkerk (Reference van Ouwerkerk1971) and van Stralen et al. (Reference van Stralen, Sohal, Cole and Sluyter1975) for more details. Furthermore, Koffman & Plesset (Reference Koffman and Plesset1983) and Yabuki & Nakabeppu (Reference Yabuki and Nakabeppu2017) have also correlated their experimental results with (6.1). Combining the above IMT expression with a Scriven-type dependence of bubble radius with respect to time, $R(t)\propto \sqrt {\alpha _l t}$, as discussed in our previous paper (Bureš & Sato Reference Bureš and Sato2021c), and originally derived by Olander & Watts (Reference Olander and Watts1969), the following dependence is suggested:

(6.2)\begin{equation} \frac{d_0(\xi)}{\xi} \propto \sqrt{\frac{\nu_l}{\alpha_l}} = \sqrt{\frac{c_{p,l}\mu_l}{\lambda_l}} = \sqrt{Pr_l}, \end{equation}

where $Pr_l$ is the liquid Prandtl number. Here, we are using $\xi (t)$ to denote the instantaneous outer edge of the microlayer. Thus, for two different fluids (under the restricted assumption of equal Jakob numbers):

(6.3)\begin{equation} \frac{d_{0,1}(\xi)}{d_{0,2}(\xi)} = \sqrt{\frac{Pr_{l,1}}{Pr_{l,2}}}. \end{equation}

For water and ethanol at atmospheric pressure, two working fluids commonly used for experimental studies of bubble growth, the ratio in (6.3) can be calculated using the respective properties at saturation to have the value 2.18. The original measurement of this ratio from the experiment of Koffman & Plesset (Reference Koffman and Plesset1983) was 1.6, though the more recent study of Utaka et al. (Reference Utaka, Kashiwabara and Ozaki2013) gives a value of 2.29.

We have decided to test the appropriateness of equation (6.3) by directly simulating bubble growth for both water and ethanol under atmospheric conditions, under the assumption of a common Jakob number of 30. Figure 34 gives a comparison of the resulting IMT distributions for the two fluids. The ratio between them, fitted using least-squares constant regression, is 1.91, a value reasonably close to the theoretical prediction of 2.18. Encouraged by this result, we have subsequently performed a set of simulations to evaluate the general efficacy of (6.2) using an artificial working fluid, also with ${Ja}=30$, the physical properties of which are listed in table 4; the molecular Prandtl number has been varied between $Pr_l=2$ and $Pr_l=5$. According to these simulations, essentially no change in the IMT distribution has occurred for this case. Thus, it appears that the model represented by (6.2) is not able to characterise the simulated IMT distribution.

Figure 34.

Initial microlayer thickness as a function of lateral position for water and ethanol with the common Jakob number ${Ja}=30$.

Table 4.

Material properties of the artificial working fluid used in our simulations.

A more detailed analysis of the IMT problem has been performed by Smirnov (Reference Smirnov1975), the result of which reads, in compact form,

(6.4)\begin{equation} d_{0,S}(t) = \sqrt{2\nu_l t\left/\left(K+\frac{4 n}{\textit{We}(t)}\right)\right.}, \end{equation}

with

(6.5)\begin{equation} K = 9(1-n) - 2\frac{(n-1)(n-2)}{n} + \frac{2n}{3} . \end{equation}

Here $n$ is the exponent in the formula describing the evolution of the bubble radius $R(t)= C_r t^n$, i.e. $n\approx 1/2$ for heat-transfer-limited bubble growth (giving $K\approx 11/6$), and ${We}(t)$ is the instantaneous microlayer formation Weber number, given as

(6.6)\begin{equation} {We}(t) = \frac{\rho_lR(t)}{\sigma} \left(\frac{\text{d}R}{\text{d}t}\right)^2 = \frac{\rho_l}{\sigma}n^2C_r^3t^{3n-2}. \end{equation}

Note that if we form a mean value of this Weber number over a given growth period, it becomes almost identical to the mean Bond number of Yabuki & Nakabeppu (Reference Yabuki and Nakabeppu2017):

(6.7)\begin{equation} \overline{{Bo}} = \frac{\rho_l \xi^2 \bar{U}}{\sigma t}=\frac{\rho_l \xi}{\sigma}\bar{U}^2=\frac{\rho_l \xi^3}{\sigma t^2}, \end{equation}

with $\xi$ being the position of the outer edge of the microlayer (i.e. the instantaneous location of microlayer formation) and $\bar {U}=\xi /t$. Note that, in Smirnov's analysis, it was assumed that $\xi \approx R$. Thus, the experimental result deduced by Yabuki & Nakabeppu (Reference Yabuki and Nakabeppu2017) as

(6.8)\begin{equation} d_{0,Y}(t) = C_Y\sqrt{\nu_l t} \end{equation}

with

(6.9)\begin{equation} C_Y = \min (0.34,0.13\overline{{Bo}}^{0.38}) \end{equation}

can be qualitatively explained by the fact that, for ${We}\gg 1$ (fast growth), no dependence of $d_0$ on surface tension should be observed, this being consistent with the classical film formation theory of Landau & Levich (Reference Landau and Levich1988).

We remark in passing that conceptually very similar results have been obtained by Moriyama & Inoue (Reference Moriyama and Inoue1996) for bubble growth in a constrained geometry with a Bond number exponent of 0.41, a value remarkably close to the 0.38 proposed by Yabuki & Nakabeppu (Reference Yabuki and Nakabeppu2017).

Smirnov's equation (6.4) has recently been re-examined by Jung & Kim (Reference Jung and Kim2018), who suggested the inclusion of an additional pre-multiplier to represent the deceleration of the flow in the microlayer:

(6.10)\begin{equation} d_{0,J}(t) = C_{{dec}} d_{0,S}(t), \end{equation}

for which $C_{{dec}}<1$. Additionally, these same authors reasoned that growing bubbles are not perfectly hemispherical, but rather of an oblate form, as illustrated in figure 35. Thus, in reality, the microlayer extent $\xi (t)$ will not be well approximated by the bubble lateral radius $R(t)$.

Figure 35.

Schematic representation of bubble geometry during growth: (a) the idealised configuration of Smirnov (Reference Smirnov1975), with $\xi \approx R$; and (b) the actual configuration as observed in experiments (e.g. Chen et al. Reference Chen, Haginiwa and Utaka2017).

Note that in the paper by Jung & Kim (Reference Jung and Kim2018), the Smirnov equation, and its subsequent modifications, are not reproduced correctly, possibly due to a typographic error: for example, equation (11) therein is inconsistent in terms of units, and consequently the Weber dimensionless group cannot be recovered. Also, in their analysis, the authors did not consider the fact that, by correctly assuming that $\xi (t)\not \approx R(t)$, one of the original assumptions of Smirnov (Reference Smirnov1975) is violated; indeed, Smirnov derived equation (6.4) from the relation

(6.11)\begin{equation} d_{0,S}(t) = \sqrt{2\nu_l \dot{\xi}\left/\left(-\frac{2}{\rho_l}\frac{\text{d}p}{\text{d}\xi} -\ddot{\xi}+\frac{2}{3}\frac{\dot{\xi}^2}{\xi}\right)\right.}, \end{equation}

the dot indicating a derivative with respect to time. It was assumed by Smirnov (Reference Smirnov1975) and Jung & Kim (Reference Jung and Kim2018) that the pressure derivative can be computed using the inviscid Rayleigh–Plesset equation (Faghri & Zhang Reference Faghri and Zhang2006), ultimately resulting in (6.4). However, as pointed out by Jung & Kim (Reference Jung and Kim2018), and also seen in our simulations, $\xi (t)$ can be only 30 % of $R(t)$, and the transition zone of Smirnov (Reference Smirnov1975), i.e. the region between $\xi$ and $R$ in figure 35, can be much larger than considered by the author. For this reason, we conclude that the Rayleigh–Plesset assumption is not valid under the present circumstances, i.e. $\xi \not \approx R$.

Alternatively, we could employ the standard, leading-order thin-film approximation (Faghri & Zhang Reference Faghri and Zhang2006) to obtain

(6.12)\begin{equation} \frac{1}{\rho_l}\frac{\text{d}p}{\text{d}\xi} ={-}\frac{2\sigma}{\rho_l}\frac{1}{R^2}\frac{\dot{R}}{\dot{\xi}}. \end{equation}

Note that in our analysis we have approximated the surface tension term in the interface-normal stress balance by assuming the bubble to have a uniform radius. Equation (6.12) can be compared with the original expression of Smirnov (Reference Smirnov1975):

(6.13)\begin{equation} \frac{1}{\rho_l}\frac{\text{d}p}{\text{d}\xi} \approx \frac{1}{\rho_l}\frac{\text{d}p}{\text{d}R}={-}\frac{2\sigma}{\rho_l}\frac{1}{R^2} + 4\ddot{R}+\frac{R\dddot{R}}{\dot{R}}. \end{equation}

This relation had been obtained from the inviscid Rayleigh–Plesset equation (Faghri & Zhang Reference Faghri and Zhang2006), and differs from the thin-film expression, (6.12), by the inclusion of terms arising both from inertial effects and from the $\xi (t)\approx R(t)$ approximation.

Finally, the modified Smirnov equation reads as

(6.14)\begin{equation} d_{0,{MS}}(t) = C_{{dec}} \sqrt{2\nu_l t\left/\left(\tilde{K}+\frac{4n}{\widetilde{{We}}(t)}\right)\right.}, \end{equation}

with

(6.15)\begin{gather} \tilde{K} = 1-m + \frac{2m}{3}, \end{gather}

(6.16)\begin{gather}\widetilde{{We}}(t)= \frac{\rho_lR(t)}{\sigma} \left(\frac{\text{d} \xi}{\text{d}t}\right)^2 = \frac{\rho_l}{\sigma}m^2C_rC_x^2t^{2m+n-2}, \end{gather}

in which

(6.17)\begin{equation} \xi(t)=C_xt^m \end{equation}

and

(6.18)\begin{equation} R(t)= C_r t^n .\end{equation}

It should be noted that the evaluation of the pressure gradient in (6.13) is questionable, as discussed in Appendix C. Nevertheless, the ability of the modified Smirnov equation to adequately reproduce the simulation results (detailed below), as well as the success of Jung & Kim (Reference Jung and Kim2018) in correlating existing measured data, has motivated the adoption of the equation in our analyses.

6.2. Parameter study for coefficient determination

In order to test the appropriateness of the Smirnov equation (6.14), we have run more than 30 simulations for different fluid properties and Jakob numbers. We have used water, ethanol and an artificial fluid (table 4) as the basis of the study, and examined the effects of the variation of parameters, both individually and in combination with each other. An illustrative sample from the overall test matrix is given in table 5. From this exercise, we have ascertained that $\rho _l$, $\mu _l$, $\sigma$, $\alpha _l$ and ${Ja}$ are indeed the key parameters affecting the IMT distribution; the last two parameters enter the problem through $C_r$ and $C_x$ (equations (6.17) and (6.18)), both of which are roughly proportional to $2\beta _g\sqrt {\alpha _l}$. Note that, for the cases considered, $\beta _g\approx Ja$. The parameters $C_r$, $C_x$, $n$ and $m$ have been estimated directly from the numerical predictions using least-squares regression. From this, we have deduced that

(6.19)\begin{align} n = 0.557 \pm 0.022, \end{align}

(6.20)\begin{align} m = 0.526 \pm 0.037. \end{align}

Evidently, the early growth of the bubble is somewhat faster than theory suggests, $n=0.5$, corresponding to $R(t)\propto \sqrt {t}$. At the same time, the small standard deviations indicate that the values of $n$ and $m$ might be considered universal.

Table 5.

Characteristics of the selected cases of the sensitivity study illustrated in figure 37. Here, $C_{{dec}}$ and $\tilde {C}_{{dec}}$ correspond to the fitted values of the pre-multiplier in (6.14), the latter being obtained if the Rayleigh–Plesset expression for pressure derivative is employed instead of (6.12).

Conversely, we have found that the ratio $\overline {\xi /R}$, averaged in time during the period of bubble growth, varies significantly between the individual cases:

(6.21)\begin{equation} \overline{\xi/R} = 0.59 \pm 0.13, \end{equation}

demonstrating strong dependence on the growth dynamics. In general, faster growth results in a higher value of $\overline {\xi /R}$ and vice versa, which could be expected, since slowly growing bubbles would have more time to relax their interface geometry to the equilibrium spherical shape under the influence of the surface tension. The mean value 0.59 derived in our numerical study is very close to the 0.6 suggested by Jung & Kim (Reference Jung and Kim2018), derived from their experimental data.

We assume that the bubble-shape relaxation to sphericity is principally governed by the interplay of the inertial and surface-tension forces. A simple and straightforward choice of the dimensionless group that could be used to devise a correlation for $\overline {\xi /R}$ is thus the mean Weber number derived for Scriven growth, i.e.

(6.22)\begin{equation} \overline{{We}}_s =\frac{2\rho_l \alpha_l^{1.5}}{\sigma\sqrt{t_g}}{Ja}^3 \approx\frac{2\rho_l \beta_g^3\alpha_l^{1.5}}{\sigma\sqrt{t_g}} = \frac{\rho_l R_{{scriven}}U_{{scriven}}^2}{\sigma}. \end{equation}

Figure 36 shows $\overline {\xi /R}$ plotted against this group for all the cases considered here, together with the following power-law fit:

(6.23)\begin{equation} \overline{\xi/R} = 0.52\overline{{We}}_s^{0.23}. \end{equation}

Note that this fit should not be used for values of $\overline {{We}}_s$ larger than those used in its derivation, the range being $0\lesssim \overline {{We}}_s\lesssim 10$, otherwise $\overline {\xi /R}\rightarrow \infty$ as $\overline {{We}}_s\rightarrow \infty$, even though the physical limits are $\overline {\xi /R}\rightarrow 1$ as $\overline {{We}}_s\rightarrow \infty$.

Figure 36.

Ratio $\overline {\xi /R}$ integrated over the growth period as a function of the mean Weber number for Scriven growth, (6.22), for all cases in the sensitivity study. The power-law fit is given by (6.23).

This result provides another perspective on microlayer formation: i.e. for vanishing Weber number, the bubbles retain their highly spherical initial shape during growth, and even though a microlayer should form at the bubble base, its lateral extent would remain essentially negligible until the point of bubble detachment. We have noted such behaviour in our conference paper (Bureš & Sato Reference Bureš and Sato2021a), in which we simulated a bubble growing in water for ${Ja}=3$, a value much smaller than that considered in the present study. Even though microlayer formation had been expected, according to the theory developed in Bureš & Sato (Reference Bureš and Sato2021c), $\beta _g\approx Ja$ being extremely small (and correspondingly the Weber number, (6.22), also) resulted in the bubble remaining spherical until detachment.

Figure 37 shows selected results (table 5 lists the characteristics of the cases displayed) obtained during the sensitivity study. Overall, very good agreement of the simulations with the modified Smirnov equation (6.14) has been achieved and a prediction of the deceleration coefficient with a rather low standard deviation has been deduced as

(6.24)\begin{equation} C_{{dec}} = 0.312 \pm 0.027. \end{equation}

In comparison, if the full Rayleigh–Plesset expression for pressure is used (i.e. one taking into account the inertial terms, similarly to (6.13)), the variation becomes much more uncertain:

(6.25)\begin{equation} \tilde{C}_{{dec}} = 0.63 \pm 0.19. \end{equation}

Figure 37.

Initial microlayer thickness as a function of time (solid lines) for selected cases in the sensitivity study performed, together with fits (dashed lines) derived from the modified Smirnov equation (6.14). The characteristics of the cases considered are listed in table 5.

For a comprehensive evaluation of the efficacy of (6.14) with the coefficient $C_{{dec}}$ taken as 0.312, figure 38 presents a histogram of the ratio of the microlayer thickness predicted by equation (6.14) against that obtained from the simulations. As can be seen, the majority of values fall within a $95\,\%$ confidence interval, given approximately as $\pm 17\,\%$.

Figure 38.

Histogram of the ratio of the microlayer thickness predicted by (6.14) and that obtained from the simulations, $d$((6.14))/$d_{sim}$. Dashed lines indicate the $95\,\%$ confidence level.