2.3.1. Defining nucleation-site locations and assigning activation temperatures

As a preliminary to the CFD simulation, we prescribe the locations of the nucleation sites on the heat transfer surface, together with a nucleation activation temperature $T_{act}$ . The methodology for defining the nucleation site locations, and the corresponding values for ${T}_{act}$ , is illustrated in flowchart form in figure 1(a). Initially, the counter for the nucleation site is set to 1. The activation temperature for Site #1 ( ${T}_{act,\#1}$ ) is obtained from the relationship between the number of active nucleation sites and the wall temperature, as read off using the trend line shown in figure 1(b), which is obtained from active nucleation site density (NSD) measured in an experiment. The abscissa represents the number of nucleation sites 1, 2, 3, $\ldots, m$ , spaced evenly, the final number $m$ being determined by the maximum activation temperature, $T_{max}$ , on the upper heater surface. Hence, the total number of nucleation sites ultimately depends on the magnitude of the applied heat flux and, as will be described later, this process is self-controlled, as implemented here. The nucleation site selection procedure begins with the first site, Site #1, its location being chosen using a (unbiased) random number generator, so could be positioned anywhere on the heat transfer surface. The associated activation temperature for this site is then read off from figure 1(b).

The second nucleation, Site #2, then needs to be defined, again via a random number generator. The counter for the number of nucleation sites is increased (i.e. $i = i+1$ ) and the associated activation temperature for this site is obtained from figure 1(b), as before. Note that ${T}_{act,\#2}$ is greater than ${T}_{act,\#1}$ . This process is continued until the activation temperature at Site $\#m$ (figure 1 b) exceeds ${T}_{max}$ , which is chosen higher than the maximum wall temperature anticipated in the simulation. It is to be recalled that the nucleation process is controlled by the local heat-transfer surface temperature, which is itself determined from the applied heat flux.

The active NSD is obtained from the measurements presented in figures 5–11 of Kossolapov (Reference Kossolapov2021), and is subsequently used in (3.2) and illustrated in figure 5. The maximum activation temperature ${T}_{max}$ was set to $T_{\textit{sat}} + 21.5\,\rm K$ , which corresponds to $NSD = {10}^{9}\,({\rm m}^{-2})$ . Here, 21.5 K is sufficiently higher than the 7 K wall superheat observed in the measured boiling curve. This concludes the pre-processing stage.

2.3.2. Activation criteria and seed-bubble initialisation

In the nucleation-site model, a small vapour seed bubble is introduced at a nucleation site when:

1. (i) the local temperature at the nucleation site reaches $T_{act}$ ; and

2. (ii) the fluid above the site is in the liquid phase.

The seed bubble is initially assumed to be spherical, with a diameter set to 2.4 times the grid cell width of the underlying grid $ ( =2.4\boldsymbol{\cdot }\max ( \Delta x,\Delta y,\Delta z ) )$ . Although the initial bubble shape is physically influenced by factors such as wettability, fluid properties, system pressure, applied heat flux and flow conditions, a spherical shape was chosen based on the observations of nucleate boiling of water at moderate pressure reported by Sakashita (Reference Sakashita2011) and Kossolapov et al. (Reference Kossolapov, Hughes, Phillips and Bucci2024). The dependence of the seed bubble diameter on the cell size inherently makes the model sensitive to the grid spacing chosen. Consequently, a grid refinement study is necessary to evaluate the impact of this assumption on the accuracy of the results.

Upon nucleation, the seed bubble is initialised as follows.

1. (i) The liquid volume fraction in the corresponding region is changed from liquid to vapour phase.

2. (ii) The velocity field in the seed bubble remains unchanged.

3. (iii) The temperature in the seed bubble is set to $T_{\textit{sat}}$ , assuming no superheated vapour inside an initial bubble.

4. (iv) The latent heat required for the seed bubble formation ( $\rho _v V_b L$ ) is subtracted from the solid cells beneath the bubble, where $V_b$ is the bubble volume.

We acknowledge that these procedures have room for improvement. Mass and energy conservations are violated, although a small seed bubble diameter (2.4 times the grid cell width) was chosen to minimise their impact. Additionally, the velocity field in the seed bubble could be approximated more physically.

In item (iv), performing latent heat subtraction in a single time step can lower the solid temperature below $T_{\textit{sat}}$ , which is considered unphysical and may lead to immediate condensation. To prevent this, latent heat subtraction is gradually applied over multiple time steps, introducing a heat sink in the energy equation, ensuring the solid temperature always remains above $T_{\textit{sat}}$ . Once the cumulative heat sink over the time integration reaches the latent heat ( $\rho _v V_b L$ ), then the seed bubble is planted as illustrated in the schematic in figure 2. The number of time steps required for the latent heat subtraction depends on the evolution of the wall temperature $T_w$ , which is influenced by the flow field, heater power, time increment $\Delta t$ and other simulation conditions. Thus, this number is not prescribed, but is dynamically determined during the simulation. In the large domain simulation (§ 4), the number of time steps required was 225 $\pm$ 40 (standard deviation).

Figure 2. Procedure of bubble nucleation.

2.3.3. Surface-condition effects

We conclude this section by detailing how our nucleation-site model accounts for surface conditions (roughness, oxidation) and for non-condensable gas dissolved in the liquid. The surface condition is known to significantly influence boiling heat transfer (Rohsenow Reference Rohsenow1971). Empirical correlations, such as those proposed by Jones, McHale & Garimella (Reference Jones, McHale and Garimella2009), indicate that roughened surfaces increase the number of nucleation sites and the size of surface cavities, thereby enhancing boiling heat transfer. In the present NSD model, this effect is accounted for indirectly by using the experimentally measured active NSD as input. By construction, the measured NSD already embodies the influence of surface conditions and dissolved non-condensable gas. In cases where experimental data are not available, correlations for the nucleation site density – for example, those proposed by Kocamustafaogullari & Ishii (Reference Kocamustafaogullari and Ishii1983) or Hibiki & Ishii (Reference Hibiki and Ishii2003) – may be used.

3.1.1. Computational domain

The computational domain size for the grid dependence study measures 1.25 $\times$ 0.63 $\times$ 0.39 mm $^3$ in the axial ( $x$ -), spanwise ( $y$ -) and wall-normal ( $z$ -) directions, respectively, as illustrated in figure 3. The origin of the coordinate system is located on the inlet plane, at the centre of the spanwise direction and on the wall surface. Gravity is applied in the negative $x$ -direction. This domain is smaller than the experimental test section due to the significant computational resources required for fine-grid simulations. Consequently, the computed results from the grid dependence study cannot be directly compared with experimental measurement. Instead, the results are compared across different grid sizes. Note that a simulation covering the full vertical length of the heated surface in the experiment, i.e. 10 mm, will be conducted in § 4.

Figure 3. Computational domain, boundary conditions and distribution of nucleation sites coloured with their activation temperature, $T_{act}$ , for the grid dependence study.

3.1.2. Boundary conditions

Inlet boundary. A velocity inlet boundary condition is imposed on the surface corresponding to the minimum $x$ -coordinate. Prior to the boiling simulations conducted with the PSI-BOIL code, a preliminary two-dimensional (2-D) steady-state, single-phase liquid flow simulation was performed using ANSYS Fluent 2022 (Ansys 2022). It is important to note that ANSYS Fluent was used exclusively to determine the inlet velocity; all subsequent simulations presented in this study were carried out using the PSI-BOIL code. The computational domain for the 2-D simulation measured 1178.0 $\times$ 5.89 ${\rm mm}^{2}$ in the streamwise and wall-normal directions, respectively, and was discretised using a structured mesh of 2000 $\times$ 200 rectangular cells. The streamwise domain length corresponds to 200 $h$ , where $h$ is the half-channel height, and was chosen to allow the turbulent flow to become fully developed before the outlet. At the inlet, a uniform velocity profile was prescribed, while a pressure outlet boundary condition was imposed at the outlet, where the pressure was fixed to zero. A no-slip wall boundary condition was applied at the wall and a symmetry condition was used on the opposite side. The SST $k$ – $\omega$ turbulence model was employed. The minimum mesh spacing near the wall was 5 $\unicode{x03BC}{\rm m}$ ; a non-dimensional wall distance of ${y}^{+}$ equals to 0.7. The computed friction Reynolds number at the outlet boundary was $Re_{tau} = 1840$ . From the profiles of velocity, turbulent viscosity, turbulent kinetic energy and turbulent dissipation rate at the outlet boundary computed with ANSYS Fluent, the transient inlet velocity field for the boiling simulations with PSI-BOIL was generated using the synthetic eddy method (Jarrin et al. Reference Jarrin, Benhamadouche, Laurence and Prosser2006). The implementation of the synthetic eddy method and its assessment for a single-phase liquid flow are provided in Appendix A.

The inlet temperature was set to a constant value of $(T_{\textit{sat}}-10)\,\rm K$ , with only the liquid phase entering the domain, i.e. no vapour phase was present at the inlet. Here, $T_{\textit{sat}}$ is the saturation temperature of water at 10.5 ${\rm bar}$ . The Neumann boundary condition is applied for pressure at the inlet.

Side and symmetry boundaries. A periodic boundary condition was applied to the surfaces at the minimum and maximum $y$ -directions, while a symmetry boundary condition was imposed on the surface at the maximum $z$ -direction, as illustrated in figure 3.

Outlet boundary. To ensure the stable release of bubbles without numerical instability from the computational domain, an outlet region was introduced in a manner similar to that described in § 6.1.4 of Sato & Niceno (Reference Sato and Niceno2013). The liquid phase is artificially transformed to vapour without volume change, thereby preventing the liquid–vapour interface from reaching the outlet boundary. Although the opposite treatment (vapour-to-liquid conversion) is also possible, we intentionally use the liquid-to-vapour conversion to avoid liquid-over-vapour stratification under gravity, which could promote Rayleigh–Taylor instability in the outlet region. The liquid volume fraction $\phi$ calculated with the Navier–Stokes solver $\phi _{\textit{cal}}$ is artificially transformed to $\phi _{\textit{trans}}$ using the next interpolation:

(3.1) \begin{equation} {{\phi }_{\textit{trans}}}=\left ( 1-c \right ){{\phi }_{\textit{cal}}}, \end{equation}

where the coefficient $c$ is defined as $c={ ( x-{{x}_{1}} )}/{ ( {{x}_{2}}-{{x}_{1}} )}$ , $x_1$ and $x_2$ being the $x$ -coordinates marking the start and end of the outlet region, respectively (see also figure 3). At the outlet boundary, corresponding to the maximum $x$ -surface, a Neumann condition is applied to the velocity, pressure and the temperature fields.

Boundaries for solid phase. The boundary conditions for the solid domain are as follows: an adiabatic condition is applied to all surfaces except for the wall surface adjacent to the fluid domain. Conjugate heat transfer is solved at the heat transfer surface by solving the energy balance (2.4). A heat source is assigned to the solid cells adjacent to the fluid domain, representing the ITO heater. The thickness of the ITO heater is 0.7 $\unicode{x03BC}{\rm m}$ , consistent with the experimental set-up. However, for simplicity, the physical properties of these cells are assumed to be those of sapphire.

Contact angle. The heat transfer surface is modelled as perfectly hydrophilic; that is, the contact angle in the CSF model (Brackbill et al. Reference Brackbill, Kothe and Zemach1992) is set to zero ( ${{\theta }_{eq}}=0$ ) for the following three reasons. First, vapour bubbles observed in experiments conducted under moderate pressure ( $\geqslant$ 10.5 ${\rm bar}$ ) exhibit very low apparent contact angles, as shown in figure 4. The bottom view (through the ITO and substrate) of vapour bubbles nucleated and growing on the ITO surface in our subcooled flow boiling experiment (Kossolapov Reference Kossolapov2021), shown in figure 4(a), features a small dry area compared with the bubble diameter, indicating a small apparent contact angle. It should be noted that the static contact angle of a water droplet on the ITO surface measured in our experiment exhibits slightly hydrophilic behaviour, i.e. the angle is slightly below $90^\circ$ . The apparent contact angles of the growing vapour bubbles on a nickel foil were directly captured by a high-speed camera (Sakashita Reference Sakashita2011), as shown in figure 4(b). Although this experiment was conducted under saturated pool boiling conditions at a system pressure of 27 ${\rm bar}$ – different from the present simulation conditions – it similarly exhibits a very small apparent contact angle.

Second, a sensitivity study of the contact angle conducted in our previous simulation (Bureš & Sato Reference Bureš and Sato2022) demonstrated that assuming a zero contact angle yielded results in close agreement with the experimental observations of Bucci (Reference Bucci2020). More precisely, the zero contact angle corresponds to a receding situation, since the simulation was performed only for the bubble growth phase. The corresponding experiment, conducted at MIT using an ITO heater (Bucci Reference Bucci2020), featured saturated pool boiling at atmospheric pressure. Although the numerical method employed in that study differed slightly from the present one – specifically, it used direct numerical simulation with the volume-of-fluid approach and a fine computational mesh of 0.5 $\unicode{x03BC}{\rm m}$ – the influence of the contact angle on the simulation results is considered consistent between the two studies.

Third, the validation case subsequently shown in figure 12 demonstrates that using the perfectly hydrophilic condition ( ${{\theta }_{eq}}=0$ ) overestimates the measured vapour area ratio on the wall. Increasing ${\theta }_{eq}$ would lead to even greater discrepancies. Therefore, the minimum value ${{\theta }_{eq}}=0$ is employed in this study.

Figure 4. Bubbles observed on the heat transfer surface under moderate pressure: (a) bottom view of subcooled flow boiling (Kossolapov Reference Kossolapov2021); and (b) side view of saturated pool boiling at 27 bar (Sakashita Reference Sakashita2011). Both experiments exhibited very small apparent contact angles.

Surface roughness. A no-slip wall boundary condition is applied to the heat-transfer surface without any slip length (Miksis & Davis Reference Miksis and Davis1994), implying that the effect of surface roughness on the velocity field is neglected. It should be noted, however, that the influence of surface roughness and oxidation on boiling heat transfer is indirectly incorporated into the NSD model, as discussed in § 2.3.

3.1.3. Initial conditions

The initial conditions were prescribed as follows: the temperature was initialised to $T_{\textit{sat}}-10\,\rm K$ in the fluid domain and $T_{\textit{sat}}$ in the solid domain. The streamwise velocity component was taken from the steady-state two-dimensional solution, and the remaining velocity components were set to zero everywhere.

3.1.4. Active nucleation site density and activation temperature $T_{act}$

The nucleation site density model described in § 2.3 requires the active nucleation site density to specify the nucleation activation temperature, as illustrated in figure 1. In the present study, this quantity was determined based on the experimental measurements reported by Kossolapov (Reference Kossolapov2021). The experimentally obtained active NSD ( ${m}^{-2}$ ) is expressed by the following linear relation:

(3.2) \begin{equation} \textit{NSD}=4.65\times {{10}^{7}}\boldsymbol{\cdot }\Delta {{T}_{w}}, \end{equation}

where $\Delta {{T}_{w}}$ is the wall superheat ( $=T_w - T_{\textit{sat}}$ ), $T_w$ being the wall surface temperature. The comparison between the measured NSD and the approximation in (3.2) is shown in figure 5, and the distribution of nucleation sites and their activation temperature are displayed in figure 3.

3.1.5. Number of grids and grid spacing

Three simulation cases of different grid resolution were calculated to evaluate the influence on the results. The grid parameters are listed in table 2, where ${\textit{Nx}}$ and ${\textit{Ny}}$ represent the number of cells in the $x$ - and $y$ -directions, respectively. Here, ${\textit{Nz}}_{\textit{fluid}}$ and ${\textit{Nz}}_{\textit{solid}}$ represent the number of cells in the $z$ -direction within the fluid and solid domains, respectively. The cell widths in the $x$ - and $y$ -directions ( $\Delta x$ and $\Delta y$ ) are uniform, whereas the grid size in the $z$ -direction ( $\Delta z$ ) is non-uniform. A small grid size for $\Delta z$ is used in the near-wall region to accurately resolve both the momentum and thermal boundary layers. The minimum and maximum cell widths in the $z$ -direction within the fluid domain are denoted as $\Delta z_{min}$ and $\Delta z_{max}$ , respectively. The cell thickness in the $z$ -direction within the solid domain adjacent to the fluid domain is set to 0.7 $\unicode{x03BC}{\rm m}$ in all cases, matching the thickness of the ITO heater used in the experiment. We introduce the Grid index defined as

(3.3) \begin{equation} \textit{Grid}\ \textit{index}={\Delta x}/{\Delta {{x}_{1}}}\;={\Delta y}/{\Delta {{y}_{1}}}\;={\Delta {{z}_{\textit{min}}}}/{\Delta {{z}_{\textit{min},1}},}\; \end{equation}

where the subscript 1 denotes the cell width of the finest grid. This gives grid index values of 1, 2 and 4 for the simulations conducted.

Table 2. Grid parameters for grid dependence study.

Figure 5. Measured active nucleation site density from the experiment (Kossolapov Reference Kossolapov2021) compared with the linear approximation given by (3.2).

4.3.1. Wall superheat, vapour coverage and bubble shapes

We begin by comparing the computed wall superheat with the experimental measurements. As shown in figure 11(a), the wall superheat approaches a value of 7.8 K, which represents the time-averaged value over the final 0.02 ${\rm s}$ of the 0.2 ${\rm s}$ transient simulation. This computed value falls within the uncertainty range of the measured value, 7.1 K. A comparison between the computed and measured wall superheat is shown in figure 11(b). According to the grid dependence study described in § 3.2, the wall superheat is expected to increase slightly with a finer grid, suggesting that the simulation may further overestimate the wall temperature.

Figure 12. Comparison of vapour area ratio on the heat transfer surface between the measurement and simulation.

The computed vapour area ratio on the wall is compared with the experimental measurement in figure 12. The simulation gives a value of approximately 3 %, obtained as a time-averaged result over the final 0.02 ${\rm s}$ of the simulation. The corresponding instantaneous distributions of vapour and liquid on the wall are shown later in figure 15(c).

In the experiment, the measured vapour area ratio is very small, of the order of $5\times 10^{-4}$ for a wall superheat of $\Delta T_w\approx 7\,\rm K$ . This is mainly because only a small number of nucleation sites are activated (approximately $4\times 10^8\,\rm m^{-2}$ ) under the present conditions. Although nucleation sites can be detected, the measured vapour area is strongly affected by the optical resolution, which corresponds to an uncertainty of roughly 5 % in vapour area ratio. This uncertainty is expected to decrease at higher heat fluxes, where individual vapour patches become larger.

Meanwhile, the simulated vapour area ratio also depends on the grid resolution, as shown in figure 7(b). In the present simulations, we use the mesh with $\textit{Grid index}=2$ . Extrapolating the results in figure 7(b) to $\textit{Grid index}\approx 0$ , corresponding to an infinitesimal grid size, increases the vapour area ratio by approximately 1 percentage point, indicating that the numerical uncertainty is of the order of 1 percentage point. Overall, the simulation value of 3 % is therefore consistent with the experimental measurement, given the uncertainties in both the experiment and the simulation.

Figure 13. Comparison of bubble shapes between (a) the experiment (10 $\times$ 10 ${\rm mm}^2$ ) and (b) simulation.

The computed bubble shapes at 0.2 ${\rm s}$ are compared with the experimental measurement in figure 13. In the lower region near the inlet, the smaller bubbles observed in the experiment are not captured in the simulation due to insufficient computational grid resolution. In the region 0.002 $\mathrm{m}\leqslant x \leqslant$ 0.005 $\mathrm{m}$ , the simulated bubble shapes closely resemble those observed in the experiment. However, in the downstream region beyond 0.005 $\mathrm{m}$ , the simulation may underestimate the bubble size. It should be emphasised that the comparison in figure 13 is qualitative, because both the experimental and simulated bubbles exhibit complex, non-spherical shapes and frequent coalescence, which makes an objective quantitative characterisation of individual bubbles difficult, particularly on the experimental side. A quantitative analysis of the bubble number density in size classes for the simulation is presented in § 4.3.10. Nonetheless, the overall bubble distribution and characteristic length scales are consistent between experiment and simulation. A more detailed quantitative shape analysis, for example, using image segmentation or volume-equivalent diameters, is left for future work.

4.3.2. Instantaneous overall flow field

The purpose of this section is to present a qualitative description of the instantaneous flow field to provide a physical context for the simulation results. A quantitative assessment follows in § 4.3.3.

Figure 14. Snapshot of the computed flow field at $t= 0.18$ ${\rm s}$ : (a) perspective view showing bubbles and temperature distributions in both the fluid and solid regions; (b) those in $x$ – $z$ view with the computational grid. The length dimensions are presented in metres ( $\mathrm{m}$ ).

A snapshot of the computed flow field at $t= 0.18$ ${\rm s}$ is shown in figure 14. The temperature distribution on the heat transfer surface, figure 14(a), exhibits small blue spots, associated with bubble nucleation, while the blue line at the inlet indicates the lower temperature region resulting from the subcooled inlet flow. As the flow progresses downstream, the bubbles grow larger due to bubble merging and vapourisation, though condensation also occurs simultaneously due to subcooling. The wall temperature also increases downstream, whereas the liquid temperature in the bulk region away from the wall remains relatively unchanged. A side view of figure 14(a) is shown in figure 14(b), where the height of the bubble crowd (dimension in the $z$ -direction) gradually increases downstream. The maximum bubble height reaches approximately 1.2 ${\rm mm}$ at $x = 0.01$ $\rm m$ , which is significantly lower than the domain height of 5.89 ${\rm mm}$ . This indicates that the influence of the domain height on the bubble growth simulation is minimal. To compare the grid and bubble sizes, the bubbles along with the grid are also shown in figure 14(b).

4.3.3. Instantaneous temperature and flow field near heated surface

To examine the computed wall temperature in more detail, the distribution of computed wall superheat on the heat transfer surface at $t = 0.18$ ${\rm s}$ is shown in figures 15(a) and 15(b), with and without bubbles, respectively. These figures apparently indicate that the temperature beneath the bubble is higher than that in the regions not covered by bubbles. To investigate the wall condition (i.e. wet or dry), the distribution of the liquid/vapour phase on the heat transfer surface is visualised in figure 15(c), along with the bubbles. Vapour occupies 3 % of the total area of the heat transfer surface, which is slightly larger than the measured value as shown in figure 12. In figure 15(c), dry patches are typically observed beneath the bubbles smaller than 0.0005 $\rm m$ (0.5 ${\rm mm}$ ) in size, but are rarely found under large bubbles ( $\gt$ 1 ${\rm mm}$ ). This indicates that the larger bubbles have already lifted off from the wall. To further investigate the liquid film beneath the bubbles, the distribution of liquid film thickness $\delta$ is presented in figure 15(d). In this figure, the red area represents a dry patch ( $\delta = 0$ ), while the white area indicates regions not covered by bubbles. From this, it can be deduced that most of the liquid film thickness beneath larger bubbles is less than 100 $\unicode{x03BC}{\rm m}$ (10 $^{-4}$ $\mathrm{m}$ ). Summarising the discussion so far, a liquid film of less than 100 $\unicode{x03BC}{\rm m}$ exists beneath large bubbles ( $\gt$ 1 ${\rm mm}$ in diameter) with no dry patch (see figure 15 c, d). However, the temperature beneath these large bubbles is noticeably higher than in other areas (see figure 15 a, b).

Figure 15. Snapshot of computed heat transfer surface at $t = 0.18$ ${\rm s}$ : (a, b) wall superheat distribution with and without bubble overlay; (c) liquid/vapour phase distribution with bubbles; (d) liquid film thickness $\delta$ , where the red area represents dry patches ( $\delta = 0$ ) and the white area indicates regions not covered by bubbles; (e, f) heat flux distribution with and without bubble overlay; and (g, h) heat flux distribution for regions with $\delta \,\leqslant \,100\,\unicode{x03BC}{\rm m}$ and $\delta \,\gt \,100\,\unicode{x03BC}{\rm m}$ , respectively. The length dimensions are measured in metres ( $\mathrm{m}$ ).

Figure 16. Magnified view of the flow field at t = 0.18 s. Sliced planes [A–A $'$ ], [B–B $'$ ] and [C–C $'$ ] are defined on (a) the liquid/vapour phase map and (b) the heat flux distribution on the wall. (c–h) Distribution of (c, e, g) velocity vectors and temperature, and (d, f, h) mass transfer rate on the sliced planes.

To evaluate the influence of liquid film thickness on heat transfer, the distribution of heat flux is presented in figure 15(e, f), with and without bubbles, respectively. While the overall distribution appears chaotic, the heat flux beneath large bubbles seems to be low and homogeneous. For better clarity, the heat flux distribution is further separated in figure 15(g, h): panel (g) shows the heat flux for liquid film thickness $\delta \lt 100\,\unicode{x03BC}{\rm m}$ , while panel (h) represents the heat flux for $\delta$ $\geqslant$ 100 $\unicode{x03BC}{\rm m}$ , including regions where the wall is not covered by vapour. The heat flux shown in figure 15(h) varies significantly over small spatial intervals. This variation may be attributed to the complex flow induced by bubble growth and shrinkage, as well as the presence of subcooled bulk liquid. Figure 15(g) also indicates that the heat flux beneath large bubbles (highlighted by purple circles) is noticeably lower than that beneath medium-sized bubbles (highlighted by brown circles). This difference is likely due to the liquid film thickness, as the film beneath large bubbles is thicker than that beneath medium-sized bubbles.

To investigate the interaction between the velocity, temperature and liquid film thickness, magnified views of sliced planes are presented in figure 16. We focus on three sliced planes – [A–A $'$ ], [B–B $'$ ] and [C–C $'$ ] – defined on the liquid/vapour phase map (figure 16 a) and the heat flux distribution on the wall, figure 16(b). The [A–A $'$ ] slice represents the region where the wall is in direct contact with both the bulk fluid and vapour (i.e. a dry patch). In this slice, superheated vapour is generated above the dry patch, while subcooled liquid flows towards the wall, reducing the temperature in the solid region, figure 16(c). Figure 16(d) shows that vapourisation occurs around the bubble surface in the vicinity of the wall, while condensation takes place in the bulk flow region. Once the bubble lifts off, as seen in the right bubble in figure 16(d), condensation occurs over almost the entire bubble surface. The [B–B $'$ ] slice (figure 16 e, f) visualises a medium-sized bubble, highlighted in figure 15(g). The medium-sized bubble detaches from the wall, and a thin liquid film with a thickness of 6–8 $\unicode{x03BC}{\rm m}$ forms beneath the entire bubble on the slice [B–B $'$ ], resulting in a high heat flux of approximately 9 $\times {10}^{5}$ $\rm W$ ${\rm m}^{-2}$ (figure 16 b). It is important to note that a high wall heat flux can be achieved in regions beneath a thin liquid film, as the liquid–vapour interface temperature is assumed to remain constant at the saturation temperature, resulting in a large temperature gradient between the wall and the interface. High intensity of vapourisation occurs at the interface of the thin liquid film (figure 16 f). The [C–C $'$ ] slice (figure 16 g, h) displays the flow around a large bubble, highlighted in figure 15(g). The bubble has lifted off and is sliding along the wall, with the liquid film exhibiting a non-uniform thickness, varying from 14 $\unicode{x03BC}{\rm m}$ to 70 $\unicode{x03BC}{\rm m}$ . Vapourisation occurs where the liquid film thickness is less than approximately 30 $\unicode{x03BC}{\rm m}$ , while condensation is observed in the thicker liquid film, as indicated in figure 16(h). Focusing on the liquid film indexed by $\delta =$ 14–30 $\unicode{x03BC}{\rm m}$ in figure 16(g), a thick superheated liquid film appears because the subcooled liquid cannot be provided from the bulk, which leads the low heat flux on [C–C $'$ ] in figure 16(b). The heat flux beneath the large bubble, where the liquid film thickness is approximately $\delta =$ 14–30 $\unicode{x03BC}{\rm m}$ , is estimated to be approximately 4 $\times$ ${10}^{5}$ $\rm W$ ${\rm m}^{-2}$ .

Additional discussion on thin liquid films. Referring to figure 16(e, f), the thin liquid film, with a thickness of 6–8 $\unicode{x03BC}{\rm m}$ , corresponding to 1–2 cells, is not considered to be fully resolved in terms of the velocity field. However, the flow inside such a thin film is expected to be nearly stagnant, except in cases where strong shear occurs at the liquid–vapour interface. In contrast, the temperature field and phase-change process are treated using an irregular star stencil that accounts for the subgrid position of the interface within the framework of the sharp-interface approach (Sato & Niceno Reference Sato and Niceno2013). This approach enables accurate evaluation of the temperature gradient to be obtained – and consequently the heat flux across the interface – even when the film is thinner than two grid cells.

Figure 17. Schematic illustration of the formation mechanisms of thick liquid films (a) in nucleate pool boiling, showing the so-called macrolayer, and (b) nucleate flow boiling.

Additional discussion on thick liquid films. A relatively thick liquid film was observed beneath large bubbles in the present simulation (see figure 16 g, h). A similar thick-film structure was also reported by Yajima et al. (Reference Yajima, Io, Miyazaki and Yabuki2022) and Io et al. (Reference Io, Tamura, Tanaka, Nakamura and Yabuki2023) for pool-boiling of water at atmospheric pressure. In their experiments, such a film – often referred to as the $\mathrm{macrolayer}$ – forms as the result of the merging of a $\mathrm{microlayer}$ and a remnant liquid clot trapped between coalescing bubbles and the heated wall, as illustrated in figure 17(a). Here, a $\mathrm{microlayer}$ can be formed for the case of rapid bubble growth with a smaller contact angle (i.e. high Jakob number and hydrophilic wall (Urbano et al. Reference Urbano, Tanguy, Huber and Colin2018)), and it does not develop if these conditions are not satisfied. The typical $\mathrm{macrolayer}$ thickness measured in the experiment of Yajima et al. (Reference Yajima, Io, Miyazaki and Yabuki2022) was approximately $30\,\unicode{x03BC}{\rm m}$ . In contrast, in flow boiling, the thick liquid film develops after bubble lift-off and is convected downstream parallel to the wall (see figure 17 b). Although the formation mechanisms differ between pool and flow boiling, the influence of the thick liquid film is similar in both cases: a reduction in the local heat transfer due to the presence of the confined liquid film between the bubble and the wall. The heat transfer within it is dominated by non-convective (conductive) heat flux across the film and, owing to its relatively large thickness, the heat flux is limited – unlike that in the much thinner $\mathrm{microlayer}$ . Yajima et al. (Reference Yajima, Io, Miyazaki and Yabuki2022) concluded that the heat-transfer coefficient associated with the macrolayer in the pool boiling was approximately $20\,{\rm kW}\,{\rm m}^{-2}\,\rm K^{-1}$ , which is lower than that in regions where liquid motion is possible. In the present simulation, the heat-transfer coefficient beneath a large bubble is approximately $47\,{\rm kW}\,{\rm m}^{-2}\,{\rm K}^{-1}$ , calculated from a local heat flux of $400\,{\rm kW\,m}^{-2}$ and a wall superheat of $8.5\,{\rm K}$ . The heat flux beneath the thick liquid film is lower than that in other regions, as shown in figure 15 (f, g).

4.3.4. Liquid film thickness and heat flux

To investigate the effect of liquid film thickness on heat flux quantitatively, the heat flux at each wall cell face is plotted as a function of its associated liquid film thickness in figure 18(a). The data used correspond to a single time step at $t = 0.18$ ${\rm s}$ . The number of the faces is 512 $\times$ 256, with each face measuring 19.5 $\times$ 19.5 $\unicode{x03BC}{\rm m}^{2}$ . Faces in dry patches, and those not covered by bubbles, are not shown in this figure. The minimum liquid film thickness detected in the simulation corresponds to half the height of the wall-adjacent cell, which is $\delta = 2.5$ $\unicode{x03BC}{\rm m}$ . In the thin liquid film region ( $\lt$ 10 $\unicode{x03BC}{\rm m}$ ) in figure 18(a), the heat flux decreases as the liquid film thickness increases, which is physically reasonable. However, in the region where $\delta$ $\gt$ 10 $\unicode{x03BC}{\rm m}$ , the heat flux does not follow the same trend, indicating the presence of convective heat transfer in addition to conduction.

Figure 18. (a) Heat flux at each wall cell face as a function of liquid film thickness. Each symbol ( $\circ$ ) represents an individual wall cell face (19.5 $\times$ 19.5 $\unicode{x03BC}{\rm m}^2$ ). (b) Histogram of the heat flux as a function of liquid film thickness, $\delta$ . The bars represent the mean heat flux for each $\delta$ range, with error bars indicating the standard deviation. (c, d) Histogram of area ratio and heat flux ratio, respectively, also as a function of liquid film thickness $\delta$ . All the data are taken from the result at $t = 0.18$ ${\rm s}$ .

Since the symbols in figure 18(a) overlap, histograms for heat flux, area ratio and heat flux ratio are presented in figure 18(b–d), respectively, as functions of liquid film thickness for better clarity. The data of dry patches ( $\delta = 0$ ) and regions not covered by bubbles ( $\delta$ $\gt$ 1024 $\unicode{x03BC}{\rm m}$ ) are included in these figures. The labels on the horizontal axis (e.g. 0–4) indicate that the liquid film thickness $\delta$ falls within the range 0 $\lt \delta \leqslant$ 4 $\unicode{x03BC}{\rm m}$ .

In figure 18(b), the bars represent the mean heat flux for each $\delta$ range, with error bars indicating the standard deviation. Although the applied heat flux, given as a boundary condition, is 1000 kW m^– $^2$ , the heat flux values displayed in the figure are lower than the applied heat flux, except for the range $0\lt \delta \leqslant 4\,\unicode{x03BC}{\rm m}$ . The wall heat flux is lower than the applied heat flux because $301\,\rm kW\,m^{-2}$ is extracted from the solid as a sink by the nucleation-site model (see figure 2). This sink supplies the latent heat required for nucleation of seed bubbles; only the remainder, ${\sim} 700\,{\rm kW\,m^{-2}}$ , crosses the heat-transfer surface. Further discussion of the heat sink will be provided in the section on heat flux partitioning (§ 4.3.6). The heat flux at dry patches is almost negligible, with a mean value of 200 W m ^– $^2$ . The maximum heat flux appears for the range 0 $\lt \delta \leqslant$ 4 $\unicode{x03BC}{\rm m}$ and decreases as the liquid film thickness increases. The mean heat flux remains relatively unchanged for $\delta$ $\gt$ 16 $\unicode{x03BC}{\rm m}$ .

Referring to figure 18(c), the area not covered by bubbles (indicated by $\delta \gt 1024\,\unicode{x03BC}{\rm m}$ ) is the highest histogram, accounting for 35 % of the total wall area. The area ratio for the range $0\lt \delta \leqslant 4\,\unicode{x03BC}{\rm m}$ remains low at 1.4 %, but those for the next four ranges (i.e. 4–8, 8–16, 16–32 and 64–128 $\unicode{x03BC}{\rm m}$ ) each exceed approximately 10 %, resulting in a total of more than 40 %. The heat flux ratio for each range, shown in figure 18(d), is calculated according to

(4.1) \begin{equation} {\textit{Heat}}\ {\textit{flux}}\ {\textit{ratio}}_i=\frac {\textit{Heat}\ {\textit{flux}}_i\times \textit{Area}\ {\textit{ratio}}_i}{\sum \limits _{i}{\left ( \textit{Heat}\ {\textit{flux}}_i\times \textit{Area}\ {\textit{ratio}}_i \right )}} \end{equation}

where the subscript $i$ indicates the range of the liquid film thickness (e.g. 0–4 or 4–8 $\unicode{x03BC}{\rm m}$ ). The highest histogram bar for the heat flux ratio appears in the range $\delta \,\gt$ 1024 $\unicode{x03BC}{\rm m}$ , accounting for 36 % of the total heat flux passing through the heat transfer surface. This total heat flux is not identical to the applied heat flux due to the nucleation site model, which directly subtracts the heat energy required to generate seed bubbles from the solid cells. In figure 18(c), each histogram bar for the thinner liquid film below 64 $\unicode{x03BC}{\rm m}$ remains under 15 %, though their sum (0–4, 4–8, 8–16, 16–32, and 32–64 $\unicode{x03BC}{\rm m}$ ) reaches 50 %.

4.3.5. Thermal-equivalent film thickness

To further investigate the effect of liquid film thickness on heat transfer, i.e. the results in figure 18(a), we introduce a thermal-equivalent film thickness, defined as

(4.2) \begin{equation} {{\delta }_{eq}}=\frac {{{\lambda }_{l}}\left ( {{T}_{w}}-{{T}_{\textit{sat}}} \right )}{{{q}''}}, \end{equation}

which represents the thickness of a liquid film that would yield the local heat flux $q''$ purely by heat conduction through the liquid film. Figure 19(a) is a scatter plot of the actual liquid film thickness $\delta$ versus the equivalent film thickness $\delta _{eq}$ at each wall cell face. Data points lying near the diagonal dashed line correspond to regions where heat transfer occurs primarily by conduction, whereas points on the right side of the dashed line indicate enhanced heat transfer, e.g. due to convective effects. Notably, points with a high wall superheat ( $\Delta T_{\textit{solid}}$ $\gt$ 10 K, shown in red) cluster around the diagonal line in the range 20 $\unicode{x03BC}{\rm m}\leqslant \delta \leqslant 50\,\unicode{x03BC}{\rm m}$ , as highlighted in figure 19(a). These cells, exhibiting high wall superheat with $\delta /\delta _{eq} \sim 1$ , are located beneath large bubbles, as highlighted by the circles in figure 19(b, c). These observations indicate that the heat transfer beneath large bubbles is governed primarily by conduction rather than convection. The presence of a thick liquid film under large bubbles leads to high wall superheat due to the low efficiency of heat conduction.

Figure 19. (a) Scatter plot of the equivalent liquid film thickness $\delta _{eq}$ as a function of the local liquid film thickness $\delta$ at each wall cell face. The colour indicates the local solid-wall superheat $\Delta T_{\textit{solid}}$ . (b) Spatial distribution of the wall temperature $\Delta T_{\textit{solid}}$ , with colour coding consistent with panel (a). (c) Distribution of the local film thickness ratio $\delta / \delta _{eq}$ , illustrating deviations from the thermal-equivalent film thickness.

4.3.6. Heat flux partitioning

To elucidate the mechanisms governing heat transfer in subcooled flow boiling, we perform a heat-flux partitioning analysis. The reported quantities are time-averaged over the final $0.02\,{\rm s}$ (0.18– $0.20\,{\rm s}$ ) of the $0.2\,{\rm s}$ transient simulation.

The total heat transfer rate through the wall $\dot {Q}_{w,\textit{total}}\,(W)$ is decomposed into four components: heat transfer to liquid-bulk ( $\dot {Q}_{w,\textit{liquid}\text{-}\textit{bulk}}$ ), liquid-film ( $\dot {Q}_{w,\textit{liquid}\text{-}\textit{film}}$ ), vapour ( $\dot {Q}_{w,vapour}$ ) and the contribution associated with initial bubble nucleation ( $\dot {Q}_{w,nucl}$ ), as illustrated in figure 20. The total heat transfer rate $\dot {Q}_{w,\textit{total}}\,(50\,\rm W)$ corresponds to the product of the applied heat flux ( $q^{\prime\prime}_w=1\,{\rm MW}\,{\rm m}^{-2}$ ) and the heater area ( $A_{w,\textit{total}}=10 \times 5\,{\rm mm}^2$ ). The terms $\dot {Q}_{w,\textit{liquid}\text{-}\textit{film}}$ and $\dot {Q}_{w,\textit{liquid}\text{-}\textit{bulk}}$ denote the heat transfer rates from the wall to the liquid film and to the liquid bulk, respectively. Here, the liquid film refers to the thin liquid layer beneath the bubble, while the liquid bulk corresponds to the surrounding liquid region outside this film. The term $\dot {Q}_{w,nucl}$ accounts for the latent heat required to generate an initial seed bubble with a diameter of 47 $\unicode{x03BC}{\rm m}$ ( $= 2.4 \times \Delta x$ ). In the nucleation-site model, $\dot {Q}_{w,nucl}$ for each bubble is subtracted from the ITO heater cells located directly beneath the activated bubble. It should be noted that if an infinitesimally small grid spacing were employed, $\dot {Q}_{w,nucl}$ would become negligible and the corresponding heat transfer rate would instead be attributed to $\dot {Q}_{w,\textit{liquid}\text{-}\textit{film}}$ or $\dot {Q}_{w,vapour}$ . The contributions of each component are $\dot {Q}_{w,\textit{liquid}\text{-}\textit{bulk}}=25.8\,\%$ $\dot {Q}_{w,\textit{liquid}\text{-}\textit{film}}=43.9\,\%$ , $\dot {Q}_{w,vapour}=0.1\,\%$ and $\dot {Q}_{w,nucl}=30.1\,\%$ .

Figure 20. Heat flux partitioning. The total heat transfer rate through the wall $\dot {Q}_{w,\textit{total}}\,(W)$ is divided into four parts: heat transfer to liquid-bulk ( $\dot {Q}_{w,\textit{liquid}\text{-}\textit{bulk}}$ ), liquid-film ( $\dot {Q}_{w,\textit{liquid}\text{-}\textit{film}}$ ), vapour ( $\dot {Q}_{w,vapour}$ ) and the contribution associated with bubble nucleation ( $\dot {Q}_{w,nucl}$ ). The heat transfer across the liquid–vapour interface is divided into bulk ( $\dot {Q}_{\textit{if,bulk}}$ ) and film ( $\dot {Q}_{\textit{if,film}}$ ) components. $A_{w,\textit{liquid}\text{-}\textit{bulk}}$ , $A_{w,\textit{liquid}\text{-}\textit{film}}$ and $A_{w,vapour}$ are the fractions of the heated surface that are in contact with the bulk liquid, covered by the wall-adjacent liquid film and vapour, respectively.

The liquid-film heat transfer rate ( $\dot {Q}_{w,\textit{liquid}\text{-}\textit{film}}$ ) is the dominant contribution, accounting for 43.9 % of the total wall heat transfer rate, while the area of the liquid-film region ( $A_{w,\textit{liquid}\text{-}\textit{film}}$ ) covers 61 % of the wall area ( $A_{w,\textit{total}}$ ). To quantify the heat transfer performance within this region, the heat transfer coefficient (HTC) is defined as

(4.3) \begin{equation} {h_{\textit{liquid}\text{-}\textit{film}}=\frac {\dot {Q}_{w,\textit{liquid}\text{-}\textit{film}}}{A_{w,\textit{liquid}\text{-}\textit{film}}(T_{w,\textit{liquid}\text{-}\textit{film}}- T_{\textit{sat}})}}, \end{equation}

where $T_{w,\textit{liquid}\text{-}\textit{film}}$ is the wall temperature averaged over the liquid-film region. In this equation, the wall superheat is used as the temperature difference, instead of $T_{w,\textit{liquid}\text{-}\textit{film}}- T_{\infty }$ ( $T_{\infty }$ is the bulk liquid temperature, which is identical to the inlet liquid temperature), because subcooled liquid rarely exists in the liquid film, as shown in figures 16(c), 16(e) and 16(g). Based on the wall superheat of $T_{w,\textit{liquid}\text{-}\textit{film}}-T_{\textit{sat}}=8.2\,\rm K$ , the HTC of the liquid-film region is $h_{liquid\text{-}film}= 88\,\rm kW\,{m^{-2}\,K^{-1}}$ .

The liquid-bulk heat transfer rate ( $\dot {Q}_{w,\textit{liquid}\text{-}\textit{bulk}}$ ) constitutes 25.8 % of the total wall heat transfer rate, with an area coverage of $A_{w,\textit{liquid}\text{-}\textit{bulk}} =$ 36 %. The HTC for the liquid bulk region is defined as

(4.4) \begin{equation} {h_{\textit{liquid}\text{-}\textit{bulk}}=\frac {\dot {Q}_{w,\textit{liquid}\text{-}\textit{bulk}}}{A_{w,\textit{liquid}\text{-}\textit{bulk}}(T_{w,\textit{liquid}\text{-}\textit{bulk}}- T_{\infty })}}, \end{equation}

where $T_{w,\textit{liquid}\text{-}\textit{bulk}}$ is the wall temperature averaged over the liquid-bulk region. Based on the computed wall superheat of $\Delta T_{w,\textit{liquid}\text{-}\textit{bulk}}=7.2\,\rm K$ and the degree of subcooling $\Delta T_{sub}=T_{\textit{sat}}-T_{\infty }=10.0\,\rm K$ , the HTC is calculated as $h_{\textit{liquid}\text{-}\textit{bulk}}= 42\,{\rm kW}\,{\rm m^{-2}\,K^{-1}}$ . Notably, this is significantly higher than that reported for single-phase liquid forced convection of water in the literature (Dittus & Boelter Reference Dittus and Boelter1930). A detailed discussion of the HTC is provided in the next section.

In addition to the wall heat transfer, the heat transfer across the liquid–vapour interface is also analysed. The interfacial heat transfer rate is partitioned into two components: $\dot {Q}_{\textit{if,film}}$ and $\dot {Q}_{\textit{if,bulk}}$ . The component $\dot {Q}_{\textit{if,film}}$ corresponds to heat transfer at the liquid–vapour interface between the vapour and liquid-film regions, as shown in figure 20, while $\dot {Q}_{\textit{if,bulk}}$ corresponds to the interface between the vapour and liquid-bulk region. The interfacial heat transfer rate $\dot {Q}_{\textit{if,film}}$ accounts for 14.2 % of the total wall heat transfer, directing from the liquid film to vapour phase, i.e. vapourisation of the liquid film. In contrast, $\dot {Q}_{\textit{if,bulk}}$ constitutes 20.7 %, corresponding to condensation of vapour into the subcooled liquid bulk.

Consequently, the total heat transfer rate at the wall ( $\dot {Q}_{w,\textit{total}}\,=50\,\rm W)$ is distributed as follows: 76.3 % ( $\approx$ $ 43.9 + 25.8 - 14.2 + 20.7$ ) is transferred to the liquid phase, while the remaining 23.7 % ( $\approx 0.1 + 30.1 + 14.2 - 20.7$ ) is transferred to the vapour phase, as shown in figure 20.

4.3.7. Discussion of the heat-flux partitioning and HTC

The heat-flux partitioning and heat-transfer coefficient (HTC) are analysed and compared with values reported in the literature. Following Rohsenow (Reference Rohsenow1953), the total wall heat flux in subcooled boiling is partitioned into single-phase convection and nucleate-boiling contributions:

(4.5) \begin{equation} q^{\prime\prime}_{\textit{total}}=q^{\prime\prime}_{\textit{spl}}+q^{\prime\prime}_{\textit{nb}}, \end{equation}

where $q^{\prime\prime}_{\textit{spl}}$ is the convective heat flux of the single-phase liquid and $q^{\prime\prime}_{\textit{nb}}$ is the nucleate-boiling contribution. The single-phase term $q^{\prime\prime}_{\textit{spl}}$ is calculated from

(4.6) \begin{equation} h_{\textit{spl}}=q^{\prime\prime}_{\textit{spl}}/(T_{w}- T_{\infty }), \end{equation}

where $h_{\textit{spl}}$ is the convective HTC. Assuming fully developed turbulent flow, it is estimated by the Dittus–Boelter correlation (Dittus & Boelter Reference Dittus and Boelter1930): $h_{\textit{spl}}=0.023(\lambda _{l}/D_{h})Re_l^{0.8}Pr_{l}^{0.4}$ , where $D_h$ is the hydraulic diameter, $Re_l$ the Reynolds number of the liquid phase based on $D_h$ and bulk velocity, and $Pr_l$ the Prandtl number of the liquid phase. Meanwhile, the $q^{\prime\prime}_{\textit{nb}}$ is modelled as

(4.7) \begin{equation} q^{\prime\prime}_{\textit{nb}}=\mu _l {\mathscr{h}}_{lv} \left [\frac {g(\rho _l-\rho _v)}{\sigma }\right ]^{1/2}Pr_l^{-3}\left [\frac {c_{pl}(T_w-T_{\textit{sat}})}{C_{sf}{\mathscr{h}}_{lv}}\right ]^{3}, \end{equation}

where $C_{sf}$ depends on the fluid-surface combination and the orientation of the axial flow. Since the water–ITO combination for vertical flow cannot be found in the database, we use $C_{sf}=0.013$ for water–copper in vertical flow here. We define the HTC for nucleate boiling term as

(4.8) \begin{equation} h_{\textit{nb}}=q^{\prime\prime}_{\textit{nb}}/(T_{w}- T_{\textit{sat}}). \end{equation}

Note that the temperature difference used in (4.8) is the wall superheat instead of $T_{w}- T_{\infty }$ , which is used in (4.6), because the nucleate boiling is more related to the wall superheat, as can be seen in the model equation for $q^{\prime\prime}_{\textit{nb}}$ .

The overall HTC can be defined as

(4.9) \begin{equation} h_{\textit{total}}=q^{\prime\prime}_{\textit{total}}/(T_{w}- T_{\infty }). \end{equation}

Empirical correlations for subcooled flow boiling have been extensively developed and many models have been proposed; a representative list is given by Devahdhanush & Mudawar (Reference Devahdhanush and Mudawar2022). Here, we select two recent models for additional comparison: Shah (Reference Shah2017) and Devahdhanush & Mudawar (Reference Devahdhanush and Mudawar2022). In both models, the Dittus–Boelter correlation is employed for $h_{\textit{spl}}$ , and $h_{\textit{total}}$ is obtained from the original empirical correlations of each model. Thus, following the concept of Rohsenow (Reference Rohsenow1953), we define $h_{\textit{nb}}$ as $h_{\textit{nb}}=(q^{\prime\prime}_{\textit{total}}-q^{\prime\prime}_{\textit{spl}})/(T_{w}- T_{\textit{sat}})$ .

To compare empirical heat-transfer coefficients (HTCs) with our simulation results, we define the HTCs in terms of the quantities introduced in § 4.3.6 as follows:

(4.10) \begin{align} h_{\textit{spl}}&=\frac {\dot {Q}_{w,\textit{liquid}\text{-}\textit{bulk}}}{A_{w,\textit{total}}(T_{w}- T_{\infty })}, \end{align}

(4.11) \begin{align} h_{\textit{nb}}&=\frac {\dot {Q}_{w,\textit{liquid}\text{-}\textit{film}}+\dot {Q}_{w,nucl}}{A_{w,\textit{total}}(T_{w}- T_{\textit{sat}})}, \end{align}

(4.12) \begin{align} h_{\textit{total}}&=\frac {\dot {Q}_{w,\textit{total}}}{A_{w,\textit{total}}(T_{w}- T_{\infty })}. \end{align}

In (4.11), $\dot {Q}_{w,\textit{liquid}\text{-}\textit{film}}$ and $\dot {Q}_{w,nucl}$ are included in $h_{\textit{nb}}$ because both heat-transfer rates are associated with nucleate boiling.

The comparisons of HTCs ( $h_{\textit{spl}}$ , $h_{\textit{nb}}$ and $h_{\textit{total}}$ ) between the empirical correlations and the present study are given in table 4. As mentioned earlier, the single-phase liquid heat-transfer coefficient $h_{\textit{spl}}$ predicted by Rohsenow (Reference Rohsenow1953), Shah (Reference Shah2017) and Devahdhanush & Mudawar (Reference Devahdhanush and Mudawar2022) is identical (all employ the Dittus–Boelter correlation) and is smaller than the value obtained in the present study. This difference arises from two effects: (i) the Dittus–Boelter correlation assumes a fully developed thermal boundary layer, whereas in the present configuration the thermal boundary layer starts from zero at the inlet and thickens in the streamwise direction; and (ii) the simulated $h_{\textit{spl}}$ is indirectly enhanced by the agitation of the liquid caused by bubble growth, shrinking due to condensation, and motion. The single-phase heat-transfer coefficient given by the Dittus–Boelter correlation, $h_{\textit{spl}} = 9.4\ {\rm kW}\,{\rm m}^{-2}\,{\rm K}^{-1}$ is much smaller than $h_{\textit{liquid}\text{-}\textit{bulk}} = 41.8\ {\rm kW}\,{\rm m}^{-2}\,{\rm K}^{-1}$ , as listed in table 4; the latter is obtained from the present simulations using (4.4). The coefficient $h_{\textit{liquid}\text{-}\textit{bulk}}$ is a more physically meaningful measure of single-phase convective heat-transfer efficiency in the context of multiphase flow, because (4.4) accounts for the area fraction of the liquid-bulk region rather than the total wall area.

Table 4. Comparison of heat transfer coefficients (HTCs) predicted by empirical correlations and the current study. All values are given in $\rm kW\,m^{-2}\,K^{-1}$ .

The nucleate-boiling HTC $h_{\textit{nb}}$ calculated from Rohsenow (Reference Rohsenow1953) is lower than both the values predicted by the other models and that obtained in the present simulations, presumably because Rohsenow’s correlation was originally developed for nucleate pool boiling. In the models of Shah and Devahdhanush & Mudawar, $h_{\textit{total}}$ is linearly proportional to $h_{\textit{spl}}$ , thus the underestimation of $h_{\textit{spl}}$ results in underestimation of both $h_{\textit{total}}$ and $h_{\textit{nb}}$ . Consequently, the computed values of $h_{\textit{nb}}$ and $h_{\textit{total}}$ are noticeably larger than those predicted by the empirical models. Nevertheless, the computed $h_{\textit{total}}$ in the present study is consistent with the measurements of Kossolapov (Reference Kossolapov, Phillips and Bucci2021), because the wall superheat shows good agreement, as illustrated in figure 11(b).

To give the computed heat transfer associated with nucleate boiling a more physically meaningful interpretation, we define $h_{liquid\text{-}film\,+\,nucl}$ as

(4.13) \begin{equation} h_{\textit{liquid}\text{-}\textit{film}\,+\,nucl}=\frac {\dot {Q}_{w,\textit{liquid}\text{-}\textit{film}}+\dot {Q}_{w,nucl}}{A_{w,\textit{liquid}\text{-}\textit{film}}(T_{w}- T_{\textit{sat}})}. \end{equation}

Here, $A_{w,\textit{liquid}\text{-}\textit{film}}$ is used because the phenomena related to nucleate boiling takes place mostly in the liquid-film region. As listed in table 4, this value reaches $148.2\ \rm kW\,m^{-2}\,K^{-1}$ , which indicates the high efficiency of nucleate boiling.

4.3.8. Time-averaged flow field

This section discusses time-averaged flow fields to give some general insights into the subcooled boiling phenomena. Figure 21 illustrates the distributions of the flow fields, averaged over time during the pseudo-steady state (0.18–0.2 ${\rm s}$ ) and across the $y$ -direction.

Figure 21. Distributions of fields averaged over time and across the $y$ -direction during the pseudo-steady state: (a) the liquid volume fraction $\phi$ ; (b) the void fraction $\alpha$ and vapour quality $\chi$ as a function of the local thermodynamic equilibrium quality $\chi _e$ ; (c) the mass transfer rate $\dot {m}$ ; (d) the vapourisation mass transfer rate ${\dot {m}}_{pos}$ ; (e) the condensation mass transfer rate ${\dot {m}}_{\textit{neg}}$ ; ( f) the degree of superheat/subcooling in the liquid phase $\Delta T_{liquid}$ ( $=T_{liquid}$ – $T_{\textit{sat}}$ ); (g) those in the vapour phase $\Delta T_{vapour}$ ( $=T_{vapour}$ – $T_{\textit{sat}}$ ); and (h) the wall superheat $\Delta T_{w}$ together with the local heat transfer coefficient at the heated surface.

Figure 21(a) shows the distribution of the liquid volume fraction $\phi$ . The 0.9 contour line gradually rises from the inlet to downstream, reaching $z \approx$ 0.001 $\rm m$ at $x = 0.01$ $\rm m$ with a slope of approximately 0.1. To quantify the vapour generation, the void fraction $\alpha$ and the vapour quality $\chi$ are evaluated as a function of the local thermodynamic equilibrium quality $\chi _e$ , as shown in figure 21(b). Here, $\alpha$ is the cross-sectional average of the vapour-phase volume fraction and $\chi$ is defined as $\chi =({\rho _v u_v A_v}/{\rho _v u_v A_v +\rho _l u_l A_l })$ , where $A_l$ and $A_v$ are the cross-sectional flow areas of the liquid and vapour phases, respectively. The local thermodynamic equilibrium quality is defined as $\chi _e= ({\mathscr{h}(x)-\mathscr{h}_{l,sat}}/{\mathscr{h}_{lv}})$ , where $\mathscr{h}(x)$ and $\mathscr{h}_{l,sat}$ denote the local enthalpy at position $x$ and the saturated-liquid enthalpy, respectively. The evolution of the void fraction $\alpha$ and the vapour quality $\chi$ in subcooled boiling flows in circular tubes has been extensively studied, as reviewed by Cai et al. (Reference Cai, Mudawar, Liu and Xi2021). However, because the present simulation covers only the entrance region of the heated section (10 ${\rm mm}$ ) in a square channel heated on one side only, no suitable correlation for direct comparison could be identified.

The time-averaged phase-change rate $\dot {m}$ (figure 21 c) exhibits both positive (vapourisation) and negative (condensation) values within the same region. Since averaging can obscure local contributions by cancellation, vapourisation and condensation rates are separated into their positive $\dot {m}_{pos}$ and negative $\dot {m}_{\textit{neg}}$ components in figure 21(d, e), respectively. Evapouration is confined to a very thin layer near the heated surface ( $z \lt$ 0.1 mm), where interfacial temperatures exceed the local saturation temperature. In contrast, condensation is distributed over a much larger region, extending into the channel core, driven by strong bulk subcooling. The highest condensation rates occur near the inlet, where liquid temperatures are lowest.

The temperature fields reinforce this interpretation. In the liquid phase (figure 21 f), a narrow superheated layer forms at the wall and tracks closely with the region of positive phase change, confirming that interfacial evapouration is thermally driven. Meanwhile, the vapour-phase temperature distribution (figure 21 g) remains close to $T_{\textit{sat}}$ , though slight superheating is observed near the wall where superheated vapour, generated by the wall, has insufficient time to equilibrate with the liquid.

The wall-surface superheat, $\Delta T_{w} = T_{w} - T_{{\textit{sat}}}$ , averaged over time and across the $y$ -direction, is shown in figure 21(h) together with the local HTC, which is evaluated from the temperature difference between the local wall superheat and the inlet subcooled-liquid temperature. The local HTC peaks at $x=0$ , whereas the vapourisation rate around $x=0$ is smaller than in the downstream region, as seen in figure 21(d). This indicates that the contribution of phase change to the HTC is lower in this region. The reason for the high HTC near $x=0$ is instead attributed to the subcooled liquid being in direct contact with the heated surface.

The distributions shown in figure 21(a,c–g) are intended only as qualitative visualisations using colour contours, whereas quantitative profiles of the liquid volume fraction and temperature in the wall-normal direction, together with the velocity field, will be presented and discussed in the next section.

4.3.9. Turbulence statistics

The HTC for the liquid bulk region attains a peak value $h_{\textit{liquid}\text{-}\textit{bulk}}=42\,{\rm kW}\,{\rm m}^{-2}\,{\rm K}^{-1}$ , as described in the previous sections, and the reasons were considered to be: (i) direct contact of subcooled liquid to the heat transfer surface especially near the inlet; and (ii) the agitation of liquid due to the motion of vapour bubbles being generated and condensed (Rohsenow & Clark Reference Rohsenow and Clark1951). To evaluate the intensity of the agitation and its influence on the heat transfer, turbulence statistics are presented in this section. Two simulation cases are presented here: single-phase liquid flow and two-phase boiling flow, both computed using LES with the PSI-BOIL code. The single-phase liquid-flow simulation was performed by deactivating the phase-change model and the nucleation-site model. All other conditions were identical to those of the boiling-flow simulation described in § 2 and further details are provided in Appendix A. Turbulence statistics for the boiling flow were computed over a period of 0.02 ${\rm s}$ using a sampling interval of $2.0 \times 10^{-5}\,{\rm s}$ during the pseudo-steady state (0.18–0.2 ${\rm s}$ ).

Figure 22. Velocity profiles averaged over time and in the $y$ -direction, $u^+_{\textit{mean}}$ versus $z^+$ , at selected streamwise locations ( $x=0,1,5,9\,{\rm mm}$ ): (a) single-phase liquid simulation (DNS shown for reference); (b) subcooled flow-boiling simulation (liquid phase). (c) Liquid volume-fraction profiles, $\phi (z)$ , at $x=1,5,9\,{\rm mm}$ .

Figure 22(a,b) shows the velocity profiles, averaged over time and in the $y$ -direction, for the single-phase case and for the liquid phase in the subcooled boiling-flow case, respectively. The velocity profiles are non-dimensionalised with respect to the friction velocity $u_{\tau }$ defined as $u_{\tau }=\sqrt {\tau _w/\rho }$ . The wall shear stress $\tau _w$ is evaluated from the near-wall velocity gradient using a one-sided finite difference based on the wall-adjacent cell value and the no-slip condition at the wall ((A10) and (A11)). For the boiling-flow case, the statistics are computed using only the liquid-phase velocity, i.e. values are sampled only in cells where the liquid volume fraction satisfies $\phi \gt 0.5$ . The computed velocity profiles for the single-phase case agree well with the reference DNS data of Lee & Moser (Reference Lee and Moser2015). In contrast, the velocity profiles for the boiling flow situation exhibit noticeably lower magnitudes (except at the inlet) owing to the increased wall friction induced by bubble growth and motion.

The existence of the vapour phase in the boiling simulation is confirmed in figure 22(c), which presents the liquid volume fraction $\phi (z)$ , the length scale being normalised by the half-channel height $h$ . The thickness of the vapour layer increases progressively in the downstream direction. The peak liquid volume fraction decreases with $x$ , reaching approximately $\phi =0.45$ at $x=9\,{\rm mm}$ , which corresponds to a vapour volume fraction of approximately 0.55.

Since the turbulence statistics, non-dimensionalised by the friction velocity $u_{\tau }$ , exhibit remarkably different characteristics from those of single-phase flow – making a direct comparison difficult – the subsequent analysis employs non-dimensionalisation based on the half-channel height $h$ and the bulk velocity $U_{bulk} = G/\rho _{l}$ , where $G$ is the mass flux.

Figure 23. Turbulence statistics of velocity at $x=1,5$ and $9\,{\rm mm}$ . Top row, single-phase liquid flow; middle row, liquid phase of the two-phase boiling flow; bottom row, vapour phase of the two-phase flow. From left to right, the mean streamwise velocity $u_{{mean}}/U_{\textit{bulk}}$ ; r.m.s. fluctuations $u^{\prime}_{{r{ms}}}/U_{\textit{bulk}}$ , $w_{{mean}}/U_{\textit{bulk}}$ , $w^{\prime}_{{r{ms}}}/U_{\textit{bulk}}$ ; Reynolds shear stress $\overline {u^{\prime}w^{\prime}}/U^2_{bulk}$ ; and turbulence intensity $TI=\sqrt {(\overline {{u^{\prime}}^2}+\overline {{v^{\prime}}^2}+\overline {{w^{\prime}}^2})/3}/U_{\textit{bulk}}$ .

Turbulent velocity statistics. Figure 23 shows the turbulence statistics of velocity at $x=1,5$ and $9\,{\rm mm}$ . The top row gives the results of single-phase flow, the middle and bottom rows those for the liquid and vapour phases of the two-phase boiling flow, respectively. In this figure, from left to right, the mean velocity in the $x$ -direction $u_{\textit{mean}}/U_{bulk}$ , the root-mean-square (r.m.s.) of the $u$ velocity fluctuation $u^{\prime}_{r{ms}}/U_{bulk}$ , the mean velocity in the $z$ -direction $w_{\textit{mean}}/U_{bulk}$ , the r.m.s. of the $w$ velocity fluctuation $w^{\prime}_{r{ms}}/U_{bulk}$ , the Reynolds stress $\overline {u'w'}/U^2_{bulk}$ , and the turbulence intensity $TI=\sqrt {(\overline {u\prime ^2}+\overline {v\prime ^2}+\overline {w\prime ^2})/3}/U_{\textit{bulk}}$ are displayed. Here, $\overline {(\boldsymbol{\cdot })}$ denotes averaging over time and the spanwise $y$ -direction at fixed $x$ and $z$ , and fluctuations are $u' = u-\overline {u}$ , etc. The statistics for the vapour phase show fluctuations because the number of bubbles appears in the region away from the wall is limited.

The mean streamwise velocity $u_{\textit{mean}}$ shows only a slight difference between the single-phase flow and the liquid phase of the boiling flow. In the vapour phase, $u_{v,\textit{mean}}$ exhibits significant spatial fluctuations; however, neglecting these variations, overall it is higher than the liquid-phase velocity $u_{l,mean}$ , owing to buoyancy-driven acceleration and the volume expansion of the vapour phase. That is, because the inlet velocity is fixed and the velocities at the heated wall and at the opposite side are constrained by symmetry boundary conditions, the vapour, which undergoes volumetric expansion, is accelerated in the streamwise direction.

The mean wall-normal velocity $w_{\textit{mean}}$ exhibits a clear difference between the single-phase flow and that of liquid phase of the boiling flow. In the single-phase flow case, $w_{\textit{mean}}$ remains very small, indicating that the flow is nearly parallel to the wall. In contrast, for the liquid phase of the boiling flow, $w_{l,mean}$ , exhibits a flow directed from the bulk towards the wall, driven by vapourisation of the liquid phase near wall. This current (i.e. negative $w_{l,mean}$ ) transports subcooled liquid from the bulk region to the heated surface, thereby enhancing the heat-transfer efficiency. Overall, the vapour phase velocity in the wall-normal direction $w_{v,\textit{mean}}$ is positive, corresponding to a flow toward the bulk region, which arises from bubble growth near the wall and the subsequent motion of detached bubbles.

The velocity fluctuation components in the liquid phase of the boiling flow, $u^{\prime}_{l,r{ms}}$ and $w^{\prime}_{l,r{ms}}$ , are approximately three and five times larger, respectively, than those for single-phase flow. This indicates that the bubble-induced agitation in the present subcooled boiling flow is approximately five times more effective than the turbulent velocity fluctuations in the wall-normal direction of the single-phase flow. The strong agitation in the boiling flow situation is also evidenced by the larger Reynolds stress term $\overline {u^{\prime}_lw^{\prime}_l}$ and $TI_l$ compared with those in the single-phase flow. The peak values of $\overline {u^{\prime}_lw^{\prime}_l}$ and $TI_l$ in the boiling flow appear at $z/h\approx 0.01$ , corresponding to $z\approx 60\,\unicode{x03BC}{\rm m}$ for all the streamwise positions.

Figure 24. Statistics related to temperature field for (a) the single-phase liquid flow and (b) the liquid phase of the two-phase subcooled-boiling flow. Left columns: Profiles of $(T_{{mean}}-T_\infty )/(T_w-T_\infty )$ , $T^{\prime}_{{r{ms}}}/(T_w-T_\infty )$ , and $\overline {w'T'}/(U_{\textit{bulk}}(T_w-T_\infty ))$ at $x=1,5,9\,{\rm mm}$ . Right columns: distribution of the liquid temperature $\Delta T_{liquid}=T_{liquid}-T_{\textit{sat}}\,(K)$ . All quantities are averaged over time and across the spanwise ( $y$ ) direction.

Temperature field and turbulent heat flux. Figure 24 presents the temperature and turbulent heat-flux distributions for the single-phase liquid flow (top row) and the liquid phase of the subcooled boiling flow (bottom row). From left to right, the profiles shown are the mean temperature, temperature fluctuation, and turbulent heat flux at $x=1,5$ and $9\,{\rm mm}$ . The temperature is non-dimensionalised by $T_w-T_{\infty }$ , where $T_w$ is the local wall temperature averaged over time and across the $y$ -direction. The rightmost panels display the distributions of the liquid temperature also averaged over the time and across the $y$ -direction. The wall superheat averaged over the entire heated wall in the single-phase liquid flow reaches approximately $62\,\rm K$ , whereas that in the boiling flow is $7.8\,\rm K$ .

As shown in figure 24 (rightmost panels), the thermal boundary layer in the single-phase liquid flow is noticeably thinner than that in the boiling flow. This difference is also evident in the temperature profiles shown in figure 24 (leftmost panels). The temperature fluctuation $T^{\prime}_{l,r{ms}}$ in the boiling flow exhibits a thicker distribution layer than that in the single-phase flow. This is because the bubble surface, i.e. the liquid–vapour interface, remains constant at the saturation temperature and bubbles with this temperature are convected into the subcooled liquid region. Presumably, as a result, the profiles of $T^{\prime}_{l,r{ms}}$ show a similar trend to the liquid–vapour volume-fraction distribution shown in figure 22(c).

The peak turbulent heat flux, $\overline {w^{\prime}_l T^{\prime}_l}$ , in the boiling flow at $x=5\,{\rm mm}$ and $9\,{\rm mm}$ is approximately eight times larger than in the single-phase flow, demonstrating the high efficiency of heat transfer in subcooled flow boiling. The peak occurs at $z/h\approx$ 0.01–0.015, corresponding to $z=60{-}90\,\unicode{x03BC}{\rm m}$ , for all streamwise positions.

4.3.10. Bubble statistics: size and spatial distribution

The flood-fill algorithm of Lafferty (Reference Lafferty2019), implemented in PSI-BOIL, identifies bubble sizes and centroids from the three-dimensional liquid volume-fraction field, $\phi$ . However, tracking temporal evolution is challenging: frequent coalescence and break-up require consistent bookkeeping of bubble IDs over time. Consequently, our analysis of bubble dynamics is restricted to instantaneous fields. Extending the method to track temporal evolution is left for future work.

Figure 25. Histogram of bubble number density, averaged over time and across the entire spatial domain.

Figure 25 presents the histogram of bubble number density, averaged over time from 0.18 ${\rm s}$ to 0.20 ${\rm s}$ and across the entire spatial domain, as a function of the equivalent bubble diameter, $D_e$ . The histogram peaks at an equivalent bubble diameter of $45\,\unicode{x03BC}{\rm m}\lt D_e\leqslant 68\,\unicode{x03BC}{\rm m}$ , which includes the seed bubble diameter of 47 $\unicode{x03BC}{\rm m}$ . Note that bubbles smaller than the seed bubble, but still larger than the grid size, may form due to condensation and break-up. Conversely, bubble growth from vapourisation and bubble coalescence lead to the formation of bubbles larger than the seed bubble.

The bubble-departure diameter (BDD) is widely used in mechanistic models of nucleate boiling. In pool boiling, the BDD can be defined unambiguously as the bubble diameter at the instant it departs from the wall, provided that coalescence with neighbouring bubbles does not occur beforehand. In flow boiling, however, as observed experimentally by Kossolapov et al. (Reference Kossolapov, Hughes, Phillips and Bucci2024) and also in the present simulation, wall-nucleated bubbles start sliding immediately after they appear, which makes it difficult to define a unique departure diameter. To circumvent this difficulty, Kossolapov et al. (Reference Kossolapov, Hughes, Phillips and Bucci2024) defined the BDD as the equivalent diameter of a bubble at the instant when its optical footprint no longer covers its original nucleation site, thereby allowing a new bubble to nucleate. Using this definition, the measured BDD under the conditions of the present simulation was $60\pm 6\,\unicode{x03BC}{\rm m}$ . In contrast, the BDD in the simulation remains approximately equal to the seed-bubble diameter, $47\,\unicode{x03BC}{\rm m}$ , because the seed bubble is immediately convected downstream before it can grow appreciably. For comparison, BDD estimates from conventional empirical correlations are $201\,\unicode{x03BC}{\rm m}$ (Cole & Rohsenow Reference Cole and Rohsenow1969); $256\,\unicode{x03BC}{\rm m}$ (Kocamustafaogullari Reference Kocamustafaogullari1983), assuming the contact angle $\theta =45^\circ$ ; and $62\,\unicode{x03BC}{\rm m}$ (Brooks & Hibiki Reference Brooks and Hibiki2015). Among these, the recent correlation of Brooks & Hibiki (Reference Brooks and Hibiki2015) is in closest agreement with the experimental value. To predict the BDD and accurately capture bubble dynamics with our simulation method, we would need to use a smaller seed bubble and a finer grid, possibly focusing on a single nucleation site within a small computational domain.

Figure 26. Distribution of bubble number density ( $\rm m^{-3}$ ) for selected ranges of equivalent bubble diameter $D_e$ . Each panel corresponds to a specific diameter range, increasing from top-left to bottom-right. The spatial resolution of each panel is set to the upper bound of its corresponding diameter range; for example, the resolution for the range 513 $\unicode{x03BC}{\rm m}\,\lt D_e\leqslant$ 769 $\unicode{x03BC}{\rm m}$ is 769 $\unicode{x03BC}{\rm m}$ .

Figure 26 shows the distribution of bubble number density, which can be used to inform the development and validation of Eulerian population-balance models. The distribution is time-averaged and spatially averaged across the $y$ -direction. The spatial resolution is adapted to the corresponding bubble-size class; for example, for the figure $513\,\unicode{x03BC}{\rm m}\lt D_e\leqslant 769\,\unicode{x03BC}{\rm m}$ , the resolution is set to 769 $\unicode{x03BC}{\rm m}$ .

Bubbles smaller than 68 $\unicode{x03BC}{\rm m}$ exhibit peak bubble number density near the inlet, in the vicinity of the wall. For the range $20\,\unicode{x03BC}{\rm m}\lt D_e\leqslant 30\,\unicode{x03BC}{\rm m}$ , a high bubble number density is observed also in the bulk region slightly away from the wall, indicating lifted-off bubbles. The bubble lift-off is also seen for the bubbles with the diameter of 30–45 $\unicode{x03BC}{\rm m}$ .

Because bubble lift-off was not measured in the experiment of Kossolapov et al. (Reference Kossolapov, Hughes, Phillips and Bucci2024), we instead compare the simulation results with available correlations from the literature. Using the empirical correlation for subcooled flow boiling proposed by Basu, Warrier & Dhir (Reference Basu, Warrier and Dhir2005), the predicted lift-off diameter under the present conditions is $208\,\unicode{x03BC}{\rm m}$ , which is larger than the value obtained in the simulation. This discrepancy may be attributed to the fact that the correlation of Basu et al. (Reference Basu, Warrier and Dhir2005) was derived from experiments conducted at pressures in the range 1–3 ${\rm bar}$ . Several other empirical correlations for lift-off diameter have been reported, but they are likewise limited to the ranges 1–3 ${\rm bar}$ (Prodanovic, Fraser & Salcudean Reference Prodanovic, Fraser and Salcudean2002), atmospheric pressure (Situ et al. Reference Situ, Hibiki, Ishii and Mori2005) and 1.2–1.3 ${\rm bar}$ (Okawa, Kubota & Ishida Reference Okawa, Kubota and Ishida2007).

Ahmadi, Ueno & Okawa (Reference Ahmadi, Ueno and Okawa2012) proposed a criterion for the transition between sliding on wall and direct lift-off, based on their experiments at 1–9 ${\rm bar}$ . They concluded that direct bubble lift-off – without prior sliding – occurs when the Jakob number, calculated from the wall superheat, exceeds 35. In our simulation, the Jakob number is 3 and, thus, direct lift-off should not occur according to this criterion. A possible reason for this inconsistency is that condensation in subcooled region is underestimated in the simulation due to insufficient grid resolution, which results in underestimation of temperature gradient and heat flux associated with condensation. In the subcooled bulk region, smaller bubbles should condense much more rapidly, whereas in the simulation, they persist in the bulk, as seen in figure 26 for $D_e \leqslant 45\,\unicode{x03BC}{\rm m}$ .

The peak bubble number density for $68\,\unicode{x03BC}{\rm m}\lt D_e\leqslant 101\,\unicode{x03BC}{\rm m}$ and $101\,\unicode{x03BC}{\rm m}\lt D_e\leqslant 152\,\unicode{x03BC}{\rm m}$ occurs downstream, rather than near $x=0$ , because bubbles larger than the seed bubble result from growth and coalescence. For bubbles with equivalent diameters in the range $513\,\unicode{x03BC}{\rm m}\lt D_e\leqslant 769\,\unicode{x03BC}{\rm m}$ , the peak occurs at $x =$ 2–3 ${\rm mm}$ , which is in good visual agreement with the experimental photograph in figure 13(a).