# Conjugate heat transfer effects on flow boiling in microchannels

F. Municchi1 I. El Mellas2 O. K. Matar3 M. Magnini2 1 Colorado School of
Mines, 1500 Illinois St., Golden, CO 80401
2 Department of Mechanical, Materials and Manufacturing Engineering,
University of Nottingham, Nottingham NG7 2RD, United Kingdom
3 Department of Chemical Engineering, Imperial College London, London SW7 2AZ,
United Kingdom

###### Abstract

This article presents a computational study of saturated flow boiling in non-
circular microchannels. The unit channel of a multi-microchannel evaporator,
consisting of the fluidic channel and surrounding evaporator walls, is
emulated and the conjugate heat transfer problem is solved. Simulations are
performed using OpenFOAM v2106 and the built-in geometric Volume Of Fluid
method, augmented with self-developed libraries to include liquid-vapour
phase-change and improve the surface tension force calculation. A systematic
study is conducted by employing water at atmospheric pressure, a channel
hydraulic diameter of $D_{h}=229\,\mathrm{\mu m}$, a uniform base heat flux of
$q_{b}=100\,\mathrm{kW/m^{2}}$, and by varying the channel width-to-height
aspect-ratio and channel fin thickness in the range $\epsilon=0.25-4$ and
$W_{f}=D_{h}/8-D_{h}$, respectively. The effects of conjugate heat transfer
and channel aspect-ratio on the bubble and evaporative film dynamics, heat
transfer, and evaporator temperature are investigated in detail. This study
reveals that, when the flow is single-phase, higher Nusselt numbers and lower
evaporator temperatures are achieved for $\epsilon<1$. In the two-phase flow
regime, the trends of the Nusselt number versus the aspect-ratio are mixed,
although for smaller channel fins an ascending trend of $\mathrm{Nu}$ for
increasing aspect-ratios is apparent. Nonetheless, due to conjugate heat
transfer, Nusselt numbers and evaporator base temperatures follow different
trends when varying the aspect-ratio, and channels with $\epsilon<1$ seem to
promote lower evaporator temperatures than higher aspect-ratio conduits,
despite exhibiting slightly worse two-phase convective heat transfer
performances.

###### keywords:

Conjugate heat transfer , Boiling , Microchannel , Two-Phase , Volume-Of-Fluid
, Bubbles

##  1 Introduction

Flow boiling in mini- and microchannels is recognised as one of the most
efficient cooling solutions for high-power-density applications. Recent
advances in manufacturing technology have enabled devices such as high-
performance computers, power electronics, lasers, avionics, electric vehicles,
batteries, photovoltaics, miniature fuel cells, energy conversion and storage
systems, evaporators, condensers and reactors, amongst other, to operate at
high power densities [1]. These applications involve heat fluxes of
$O(\mathrm{MW/m^{2}})$, while the heat removal capability of traditional
single-phase cooling is below $1\,\mathrm{MW/m^{2}}$ [2], which has resulted
in a dramatic and urgent demand for high-performance thermal management
systems that can transfer unprecedentedly high heat fluxes. Boiling two-phase
flows in microchannels [3]: (i) yield very high heat transfer coefficients and
maintain uniform surface temperatures, vital for the correct operation of
components; (ii) respond passively to alleviate localised ‘hot-spots’, as the
heat transfer coefficient increases with the heat flux without the need for
actively-controlled (higher) flow rates; (iii) offer a large surface-to-volume
ratio, benefitting the compactness of the system. As such, flow boiling in
microchannels has been studied extensively in recent years, with a focus on
flow pattern transitions, void fraction, pressure drop, heat transfer
coefficient, and critical heat flux [1, 4].

In order to increase the surface area, the two-phase flow is usually organised
into microevaporators, where multiple parallel microchannels are manufactured
into a thin metal die made of conductive material, which is placed in direct
contact with the surface to be refrigerated [5]. Heat is delivered to the
fluid via heat conduction through the solid walls of the evaporator and heat
convection at the contact surface between fluid and solid, thus forming a
conjugate heat transfer problem. Since heat is applied to one side of the
evaporator, the channels are subject to a nonuniform heating condition, as
heat is unequally distributed among the channel walls. The heat removal
capability of a microevaporator is characterised by the convective heat
transfer achieved by the two-phase flow, quantified by a boiling heat transfer
coefficient or Nusselt number, and by the value of the evaporator base
temperature, since the heat sink must limit the temperature of the device to
be refrigerated below a threshold value.

Within microchannels, after nucleation vapour bubbles grow quickly and occupy
the cross-section leading to slug or annular flow patterns, with bubbly flows
being suppressed already at very low values of vapour quality [6]. Surface
tension forces rearrange the liquid-vapour interface pattern into thick liquid
lobes at channel corners and thin liquid films at the centre [7], and the
thickness and morphology of this film strongly depends on the channel aspect-
ratio [8, 9]. The distribution of this liquid film has a direct impact on the
heat transfer coefficient, as both experimental [10, 11] and numerical [12,
13] studies have demonstrated that the local heat transfer coefficient $h$ is
inversely proportional to the film thickness $\delta$,
$h\approx\lambda_{l}/\delta$, with $\lambda_{l}$ being the liquid thermal
conductivity, with film dryout being highly detrimental to heat transfer
performance. As such, the channel aspect-ratio is expected to yield
significant impact on microchannel boiling heat transfer. The experimental
literature on the impact of the channel shape on boiling heat transfer is
rather vast, and comprehensive reviews have been carried out by Magnini and
Matar [13], Vontas et al. [14] and Al-Zaidi et al. [5]. It emerges that there
is still substantial disagreement on the effect of the channel aspect-ratio on
boiling heat transfer, with contrasting trends of heat transfer coefficient
and microchannel wall temperature versus aspect-ratio being reported in these
studies.

More recently, interface-resolving numerical methods have been employed to
investigate relevant fluid mechanics structures and heat transfer mechanisms
pertinent to flow boiling in noncircular microchannels. Magnini and Matar [13]
performed a systematic analysis of the impact of the channel aspect-ratio on
the bubble dynamics and heat transfer and concluded that square channels
performed better at lower flow rates while rectangular channels exhibited
larger Nusselt numbers at higher flow rates. However, their study did not
include conjugate heat transfer through the evaporator walls, and it
considered an idealised slug flow at low heat flux conditions ($q\sim
10\,\mathrm{kW/m^{2}}$), where liquid film dryout never occurred. Vontas et
al. [14] simulated flow boiling in single and multiple rectangular channels
for different channel hydraulic diameters also accounting for the evaporator
walls. They reported an ascending trend of the heat transfer coefficient when
reducing the channel size, but did not study the impact of the channel aspect-
ratio, which was maintained constant. Lin et al. [15] modeled flow boiling in
a single rectangular microchannel accounting for one evaporator wall. They
varied the material and thickness of the wall and observed that two-phase heat
transfer was enhanced by thicker walls, which exhibited higher temperature and
increased the bubble growth rate, and by highly-conductive materials. However,
their study considered only one microchannel wall and thus it is not
representative of an actual heat sink.

The literature review outlined above emphasises the fact that the impact of
the channel aspect-ratio and conjugate heat transfer on flow boiling in
microchannel evaporators is still unclear. Conjugate heat transfer is
particularly important in microgeometries, where solid walls are of thickness
comparable to the channel size, such that heat diffuses along all coordinate
directions [16] and the fluid is nonuniformly heated around the channel
perimeter. This article presents a computational study of the effect of the
geometrical features of channel and evaporator walls on heat transfer in both
single-phase and two-phase flow. The unit channel of a multi-microchannel
evaporator, consisting of channel and surrounding walls, is modelled and a
conjugate heat transfer problem is solved. Simulations are performed with
OpenFOAM v2106, using the built-in geometric Volume Of Fluid (VOF) solver
isoAdvector [17], augmented with self-developed functions implementing
thermally-driven liquid-vapour phase-change and improving the native surface
tension method. Simulations are run for a constant value of the channel
hydraulic diameter, $D_{h}=229\,\mathrm{\mu m}$, mass flux
($G=140\,\mathrm{kg/(m^{2}s)}$) and base heat flux
($q_{b}=100\,\mathrm{kW/m^{2}}$), using water at atmospheric pressure and
copper as working fluid and evaporator material, respectively. A constant heat
flux is provided to the outer wall of the evaporator base, and the resulting
single- and two-phase dynamics and heat transfer are investigated for a range
of aspect-ratios $\epsilon=0.25-4$ and of channel fin widths
$W_{f}=D_{h}/8-D_{h}$.

The rest of this article is organised as follows: the numerical framework is
described in Sec. 2 and results of validation benchmarks are illustrated in
Sec. 3; the results of the single-phase and flow boiling simulations are
presented in Sec. 4; Sec. 5 provides a discussion of the observed heat
transfer trends, and conclusions are summarised in Sec. 6.

##  2 Numerical framework

###  2.1 Governing equations

The numerical model is based on the solution of a conjugate heat transfer
problem in the fluid and solid regions of the domain. The fluid model solves
the Navier-Stokes and energy equations for the flow of two immiscible phases,
namely liquid and vapour, separated by an interface. The liquid and vapour
phases are both treated as incompressible, Newtonian fluids. A single-fluid
formulation is adopted and the two phases are treated as a single mixture
fluid with variable properties across the interface, such that a single field
of velocity, pressure, and temperature are sufficient to describe the flow,
and a single set of conservation equations holds throughout the domain [18].
Accordingly, the governing equations of mass, momentum, and energy, for a flow
with phase-change, are expressed as follows:

$\nabla\cdot\boldsymbol{u}=\frac{\dot{\rho}}{\rho}$ (1)
$\frac{\partial(\rho\boldsymbol{u})}{\partial
t}+\nabla\cdot(\rho\boldsymbol{u}\boldsymbol{u})=-\nabla
p+\nabla\cdot\mu\biggl{[}(\nabla\boldsymbol{u})+(\nabla\boldsymbol{u})^{T}\biggr{]}+\boldsymbol{F_{\sigma}}$
(2) $\frac{\partial(\rho c_{p}T)}{\partial t}+\nabla\cdot(\rho
c_{p}\boldsymbol{u}T)=\nabla\cdot(\lambda\nabla T)+\dot{h}$ (3)

where $\boldsymbol{u}$ indicates the fluid velocity, $\dot{\rho}$ the mass
flux due to phase-change, $\rho$ the mixture fluid density, $t$ the time, $p$
the pressure, $\mu$ the dynamic viscosity, $\boldsymbol{F_{\sigma}}$ the
surface tension force vector, $T$ the temperature, $c_{p}$ the constant
pressure specific heat, $\lambda$ the thermal conductivity, and $\dot{h}$ is
the enthalpy source due to phase-change. Details of surface tension and phase-
change models are provided in the subsections below. Gravitational effects are
neglected in this work.

By means of the VOF method, a volume fraction field $\alpha$ is defined to map
liquid and vapour phases throughout the flow domain. In each computational
cell of the domain, $\alpha$ identifies the fraction of the cell occupied by
the primary phase, which corresponds to liquid in the present case. Therefore,
the volume fraction takes values of 1 in the liquid, 0 in the vapour, and
$0<\alpha<1$ in cells that are cut by the interface. The volume fraction field
is evolved upon solution of the following transport equation:

$\frac{\partial\alpha}{\partial
t}+\nabla\cdot(\alpha\boldsymbol{u})=\frac{\dot{\rho}}{\rho}\alpha$ (4)

Using the volume fraction field, the properties of the mixture fluid can be
computed as an average over the two phases, e.g.
$\rho=\alpha\rho_{l}+(1-\alpha)\rho_{v}$, with the subscripts $v$ and $l$
denoting vapour- and liquid-specific properties, respectively. All the fluid-
specific properties (e.g. $\rho_{v}$, $\rho_{l}$, $\mu_{v}$, $\mu_{l}$, etc.)
are considered constant in this work.

A separate mesh is used to discretise the solid domain. The temperature field
$T_{s}$ in the solid domain is obtained by solving the following heat
conduction equation:

$\frac{\partial(\rho_{s}c_{p,s}T_{s})}{\partial
t}=\nabla\cdot(\lambda_{s}\nabla T_{s})$ (5)

where the subscript $s$ refers to solid. At the boundary between solid and
fluid regions, the fluid and solid temperature fields are coupled by imposing
continuity of temperature and heat flux, which is achieved by solving the
fluid and solid energy equations iteratively, and adjusting the boundary
conditions after every iteration.

###  2.2 Surface tension model

The surface tension force, $\boldsymbol{F_{\sigma}}$ in Eq. (2), is formulated
according to the Continuum Surface Force method [19] and computed as:

$\boldsymbol{F_{\sigma}}=\frac{2\rho}{\rho_{l}+\rho_{v}}\sigma\kappa|\nabla\alpha|$
(6)

where $\sigma$ is the surface tension coefficient (considered constant) and
$\kappa$ the local interface curvature; the term $2\rho/(\rho_{l}+\rho_{v})$
represents a density-correction factor that does not change the integral of
the surface tension across the interface, but redistributes the surface
tension towards the denser fluid to prevent unphysical accelerations in the
region occupied by the lighter fluid [19]; the impact of the density-
correction term on the simulation of flow boiling in microchannels will be
investigated in Sec. 3.2. The interface curvature is estimated by means of
derivatives of a smoothed volume fraction field $\tilde{\alpha}$,
$\kappa=\nabla\cdot(\nabla\tilde{\alpha}/|\nabla\tilde{\alpha}|)$ [20], where
the smoothed volume fraction field is obtained by interpolating $\alpha$ from
the computational cell centres to face centres and by averaging the resulting
$\alpha_{f}$ field back to cell centres according to:

$\tilde{\alpha}=\frac{\sum_{f}\alpha_{f}S_{f}}{\sum_{f}S_{f}}$ (7)

where the sum spans the $f$ faces (each of area $S_{f}$) of the control
volume. The use of a smoothed volume fraction field to calculate $\kappa$ was
already included in the original work of Brackbill et al. [19] and has the
effect of improving the accuracy of the calculation of $\kappa$, compared to
using the unsmoothed field $\alpha$, thus reducing the magnitude of the
spurious velocity caused by errors in the surface tension calculation. The
smoothing cycle can be repeated multiple times to further smoothen $\alpha$
before calculating $\kappa$, though its beneficial effect on mitigating the
spurious velocity saturates after a few cycles [20]; the optimal number of
smoothing cycles will be object of the analysis in Sec. 3.2.

A static contact angle is imposed at the wall boundaries where the liquid-
vapour interface comes in contact with solid walls. The contact angle is set
by adjusting the direction of the unit normal to the interface in boundary
cells that are cut by the contact line, using OpenFOAM’s built-in
implementation.

###  2.3 Phase-change model

The mass and enthalpy source terms due to evaporation, $\dot{\rho}$ in Eqs.
(1) and (4) and $\dot{h}$ in Eq. (3), are modelled according to the work of
Hardt and Wondra [21]. The evaporating mass flux at the liquid-vapour
interface, here denoted as $\dot{m}$, is calculated as a function of the local
interface superheat according to the Hertz-Knudsen-Schrage relationship [22],
and adopting the linearisation proposed by Tanasawa [23] for low values of the
superheat:

$\dot{m}=\frac{2\gamma}{2-\gamma}\biggl{(}\frac{M}{2\pi
R_{g}}\biggr{)}^{1/2}\frac{\rho_{v}h_{lv}(T_{lv}-T_{sat})}{T_{sat}^{3/2}}$ (8)

where $\gamma$ is the evaporation coefficient, $h_{lv}$ is the vaporisation
latent heat, $M$ is the molecular weight of the fluid, $R_{g}$ the universal
gas constant, $T_{lv}$ is the temperature at the liquid-vapour interface, and
$T_{sat}$ the saturation temperature of the fluid. In this work, the
evaporation coefficient is set to 1 according to preliminary test benchmarks
[24, 12]. An initial volumetric mass source $\dot{\rho}_{0}$ is estimated by
accounting for the evaporating mass flux calculated based on the temperature
on the liquid side of the interface:

$\dot{\rho}_{0}=N\alpha|\nabla\alpha|\dot{m}$ (9)

where $\alpha$ is the liquid volume fraction and $N$ is a normalisation factor
to ensure that the global evaporation rate is preserved [21]. A smoothed
evaporation mass source $\dot{\rho}_{1}$ is then obtained by solving a steady
diffusion equation as described in detail in [21, 12]; the smoothing of the
evaporation source term occurs over a few cells across the interface and
improves the numerical stability of the solver. The final volumetric source
$\dot{\rho}$ is obtained by redistributing $\dot{\rho}_{1}$ across the
interface on the vapor and liquid side according to:

$\dot{\rho}=\begin{cases}N_{v}(1-\alpha)\dot{\rho}_{1},\quad&\text{if}\,\alpha<\alpha_{cut}\\\
-N_{l}\alpha\dot{\rho}_{1},\quad&\text{if}\,\alpha>1-\alpha_{cut}\\\
0,\quad&\text{if}\,\alpha_{cut}<\alpha<1-\alpha_{cut}\end{cases}$ (10)

where $N_{l},N_{v}$ are normalisation factors ensuring that the masses of
liquid evaporated and vapour created are conserved by the redistribution step.
The threshold parameter $\alpha_{cut}$, here $\alpha_{cut}=10^{-3}$,
guarantees that the evaporation source term is nonzero only on vapour- or
liquid-full cells. The enthalpy source $\dot{h}$ accounts for the latent heat
dissipated by the evaporation process and is calculated from the initial
volumetric source term as $\dot{h}=-\dot{\rho}_{0}h_{lv}$.

###  2.4 Discretisation methods

The governing equations are solved using a customised implementation in
OpenFOAM, version v2106. The transport equations are discretised with a
finite-volume method on a collocated grid arrangement. OpenFOAM v2106 provides
both an algebraic (interFOAM; see [25, 26]) and a geometric (isoAdvector; see
[17, 27]) VOF solver. The numerical framework outlined in the previous
subsections has been implemented in both versions of VOF, which in our custom
solver are merged into a single solver, where the interface advection method
is specified by the user as a new solver option in fvSolution. This way, the
user can select whether to run the solver using the algebraic VOF (interFOAM
mode), or using the geometric VOF (isoAdvector mode) with any of the built-in
interface reconstruction methods (isoAlpha, isoRDF, plicRDF; see [27]). When
running the solver in interFOAM mode, the volume fraction Eq. (4) is
discretised with a first-order time-explicit scheme based on the MULES
(Multidimensional Universal Limiter for Explicit Solution) algorithm [25, 26],
set with nAlphaCorr 1, nAlphaSubCycles 3 and cAlpha 1. When running the solver
in isoAdvector mode, the advection method isoAlpha is selected, with
nAlphaSubCycles 3. All other equations are integrated in time with a first-
order implicit method. The divergence operators are discretised using second-
order TVD (Total Variation Diminishing) schemes [28], whereas Laplacian
operators are discretised with central finite-differences. The PISO (Pressure
Implicit Splitting of Operators) algorithm [29] is utilised to iteratively
update pressures and velocities (momentumPredictor yes; nCorrectors 3) within
each time-step. The residuals thresholds for the iterative solution of the
flow equations are set to $10^{-7}$ for the velocity, $10^{-9}$ for the
pressure, $10^{-8}$ for the volume fraction, and $10^{-10}$ for the
temperature and the evaporation source term smoothing equation. The time-step
of the simulation is variable and is calculated based on a maximum allowed
Courant number of $0.1$, unless otherwise specified.

##  3 Validation

###  3.1 Vapour bubble growth in superheated liquid

To assess the implementation of the evaporation model in both the algebraic
and geometric VOF methods, the first test case selected is the growth of a
spherical vapour bubble in an infinitely extended superheated liquid domain.
When a vapour bubble nucleates in a pool of superheated liquid
($T_{l}=T_{\infty}$) far from solid walls, there exists an initial growth
stage where $T_{v}=T_{\infty}$ and $p_{v}=p_{sat}(T_{\infty})$, the bubble
growth is limited only by the inertia of the surrounding liquid being
displaced and the bubble radius follows a linear trend with time. As the
bubble grows, the vapour pressure and temperature decrease until $p_{v}\approx
p_{\infty}$ and $T_{v}\approx T_{sat}(p_{\infty})$. This identifies a second
growth stage where the bubble growth rate is limited only by the rate at which
latent heat is supplied at the interface, and the bubble radius obeys a $\sim
t^{1/2}$ law [22].

Figure 1: Vapour bubble growth in superheated liquid. (a) Snapshot of the
bubble (grey surface), velocity field (top half) and temperature field (bottom
half) at the end of a simulation run using the geometric VOF and
$R_{b,0}/\Delta=40$. (b) Vapour bubble radius versus time for simulations run
with both the algebraic (MULES) and geometric (isoAdvector) VOF methods;
$R_{b,0}$ indicates the bubble radius at $t=0$ and $\Delta$ the cell size.

The latter heat-transfer-controlled growth stage is reproduced here as a
validation benchmark for the solver. A spherical steam bubble of initial
radius $R_{b,0}=50\,\mathrm{\mu m}$ and temperature $T_{v}=373.15\,\mathrm{K}$
is placed in a pool of superheated water at the system pressure
$p_{\infty}=1\,\mathrm{atm}$ and temperature $T_{\infty}=374.15\,\mathrm{K}$.
The fluid domain is a two-dimensional axisymmetric square box of side of
$1\,\mathrm{mm}$; using symmetry boundary conditions, only one quarter of the
domain is simulated. A thin thermal boundary layer surrounds the bubble during
its growth and the exact temperature profile within the liquid must be set at
$t=0$ for a coherent initial growth. Scriven [30] has derived an analytical
solution for the heat-transfer-controlled growth stage of the bubble,
providing expressions for the liquid temperature over time and distance from
the interface, and a law for the bubble radius over time
$R_{b}(t)=2\beta(at)^{1/2}$, with $\beta$ being a constant obtained from the
solution [30], and $a=\lambda_{l}/(\rho_{l}c_{p,l})$ being the liquid thermal
diffusivity. For the conditions presently simulated, the initial boundary
layer thickness (when, at $t=0$, $R_{b,0}=50\,\mathrm{\mu m}$) is about
$16\,\mathrm{\mu m}$ and the temperature profile resulting from the analytical
solution is set as an initial condition for the liquid temperature in the
simulation. The initial volume fraction field is set via OpenFOAM’s utility
setAlphaField, which enables accurate initialisation of volume fractions by
calculating the intersections between implicit functions defining the bubble
shape and the domain mesh. The simulations are run till
$t=0.0047\,\mathrm{s}$, where the bubble is about 4 times its initial size.
The maximum Courant number allowed for the simulation is set to 0.02.

The numerical solver is run both in algebraic and geometric VOF modes to
compare the two different interface advection methods. Surface tension is
disabled by setting $\sigma=0$, thus enabling us to test specifically the
performances of the phase-change model and interface advection methods,
without the influence of the surface tension algorithm. Four different meshes
are utilised, all structured with orthogonal and uniform hexahedrons, with
grid spacings of $\Delta=2.5\,\mathrm{\mu m}$ ($R_{b,0}/\Delta=20$),
$\Delta=1.25\,\mathrm{\mu m}$ (40), $\Delta=0.625\,\mathrm{\mu m}$ (80), and
$\Delta=0.3125\,\mathrm{\mu m}$ (160). A snapshot of the bubble growth
dynamics is shown in Fig. 1(a), where the velocity and temperature fields
surrounding the bubble at the last time instant of the simulation are
depicted. It can be seen that the bubble preserves the spherical shape during
its growth, a thin temperature boundary layer surrounds the bubble and the
liquid velocity is directed radially outward, with a maximum magnitude
identified in the proximity of the liquid-vapour interface. The comparison of
the bubble growth rate versus time for the different solver configurations
tested and the analytical solution are presented in Fig. 1(b). At low mesh
resolutions, all advection methods yield a faster growth rate during the
initial growth stage. This can be ascribed to the insufficient resolution of
the initial thermal boundary layer surrounding the bubble. At $t=0$, the
evaporation rate is calculated based on the temperatures at the centroids of
the first few liquid cells nearby the interface, and therefore when the mesh
is coarser these centroids are farther from the interface and experience
excessively high temperatures, resulting in higher evaporation rates. At low
mesh resolutions, the combination of phase-change model and algebraic VOF
method seems to perform better than the geometric VOF in terms of $R(t)$,
however the slopes of the $R(t)$ curve exhibited by the MULES method at later
growth stages are more far off the exact solution than the isoAdvector ones.
The deviations with the exact solution become of comparable magnitude between
the two methods when $R_{b,0}/\Delta=80$ and isoAdvector (interFOAM)
overpredicts (underpredicts) the exact growth rate by less than 2%, whereas
analytical and numerical solutions become almost indistinguishable when
$R_{b,0}/\Delta=160$. Therefore, both methods yield solutions that converge to
the analytical one as the mesh is refined.

In terms of computational overhead, the simulations performed with the
algebraic VOF were slightly faster than those run with isoAdvector, which can
be ascribed to the extra geometric interface reconstruction steps performed by
the latter. The simulations with $R_{b,0}/\Delta=80$ (160), featuring 640,000
(2,560,000) mesh cells, required 276 (3,100) core-hours for the algebraic VOF
and 320 (3,950) core-hours for the geometric VOF. These simulations were run
on UK’s Tier-2 supercomputer Sulis, which features Dell PowerEdge R6525
computing nodes each with two AMD EPYC 7742 (Rome) 2.25 GHz 64-core
processors, thus making 128 cores and 512 GB DDR4-3200 RAM per node.

###  3.2 Flow boiling in a microchannel

Figure 2: Vapour bubble growth in a microchannel, comparison with Mukherjee et
al. [31]. (a) Snapshot of the bubble, coloured by the evaporation rate, and
temperature field along a vertical centreline plane, at $t=1.7\,\mathrm{ms}$
for a simulation run using the geometric VOF, $n=4$ curvature smoothers and
density correction. Flow is from left to right. (b) Bubble equivalent diameter
versus time for simulations run with both the algebraic (MULES) and geometric
(isoAdv) VOF methods, different number of smoothers ($n$) in the curvature
calculation, with or without density-correction term ($\rho$ corr) for the
surface tension.

The second validation benchmark selected is the numerical simulation of the
growth of a vapour bubble at the wall of a heated microchannel, which emulates
the experiments and matching simulation performed by Mukherjee et al. [31].
Mukherjee et al. [31] performed an experiment in a microchannel of slightly
trapezoidal cross-section and hydraulic diameter $D_{h}=229\,\mathrm{\mu m}$,
where saturated water at atmospheric pressure was introduced with an average
speed of $0.146\,\mathrm{m/s}$. The microchannel was cut into a brass block
and heated from three sides, while a plexiglass cover allowed visualisation
from the fourth side. The temperature of the heated walls was reported as
$102.1\,\mathrm{{}^{\circ}C}$. They observed bubble nucleation and growth
along the microchannel and measured the bubble equivalent diameter versus
time, see data in Fig. 2(b). Additionally, they performed numerical
simulations of the flow by using a level-set method, modelling a perfectly
square microchannel cross-section, and initialising a steam bubble of diameter
of $50\,\mathrm{\mu m}$ over a heated wall with a contact angle of
$30^{\circ}$ (hydrophilic walls). Their numerical results are also included in
Fig. 2(b).

The experimental setup is here emulated as a square microchannel
($D_{h}=229\,\mathrm{\mu m}$) of $5D_{h}$ length, with the solid region being
disregarded. The channel is heated with a constant temperature of
$102.1\,\mathrm{{}^{\circ}C}$ imposed on three sides while the fourth wall is
set as adiabatic. Water at saturation temperature,
$100\,\mathrm{{}^{\circ}C}$, enters the channel with a uniform velocity of
$0.146\,\mathrm{m/s}$. No-slip and a static contact angle of $30^{\circ}$ are
set at all walls. A vapour bubble of diameter of $50\,\mathrm{\mu m}$ is
initialised at the heated wall, opposite the adiabatic one. Gravity is
neglected owing to the small spatial scales. To obtain realistic velocity and
temperature fields to be set as initial conditions at $t=0$, a preliminary
liquid-only simulation is run till steady-state is achieved. The computational
mesh is a structured orthogonal mesh made of uniform hexahedrons and grid
spacing of $2.29\,\mathrm{\mu m}$. A similar mesh was used by Mukherjee et al.
[31]; tests with a finer mesh did not show appreciable differences in the
results. The simulation is evolved in time until the bubble nose approaches
the outlet section of the microchannel.

The boiling solver is tested both in algebraic and geometric VOF modes. With
both methods, tests are performed to assess the impact of the surface tension
model with and without the density-correction term, and with different numbers
of smoothing cycles for the interface curvature calculation (see Section 2.2).
Figure 2(a) shows a snapshot of the flow dynamics little before the bubble
reaches the outlet section. A thermal boundary layer develops over the heated
walls as a result of the colder water coming through the inlet section. At the
time instant displayed, the bubble has grown sufficiently to become elongated
and dry patches form over the microchannel walls. The liquid evaporation rate
is maximum at the solid-liquid-vapour contact line, where temperature is the
highest. The bubble equivalent diameter versus time achieved with the
different simulation setups are displayed in Fig. 2(b). Without neither
density-correction nor smoothing cycles, both isoAdvector and MULES
(interFOAM) advection methods yield a significant overestimation of the bubble
growth. This can be ascribed to the presence of parasitic currents related to
errors in surface tension, which enhance convective flows near interfaces
[21]. isoAdvector seems to particularly suffer from parasitic currents, as it
was previously observed in benchmark tests by Gamet et al. [32] and Magnini et
al. [9], which is due to the sharper volume fraction changes at the interface
resulting from the geometric VOF advection, which are detrimental to the
interface curvature calculation when this is obtained via gradients of
$\alpha$. The situation improves when the density-correction term in Eq. (7)
is enabled, in particular from the instant when the bubble equivalent diameter
grows above $D_{h}$ and the bubble becomes elongated. Spurious currents
manifest in simulations of elongated bubbles in microchannels by generating
artificial vortices in the liquid ahead the bubble nose [33], which explains
the lower bubble growth rate trend observed in Fig. 2(b) as density-correction
is activated. The activation of smoothing cycles in the interface curvature
calculation does not yield significant further changes in the bubble growth
rate obtained with the MULES (though deviations appear at later stages), but
it makes a significant difference when using isoAdvector. Since isoAdvector
exhibits a sharper volume fraction field across the interface than MULES, a
test with 4 smoothing cycles was also performed, leading to a further
reduction of about $5\%$ of the bubble diameter achieved at the end of the
simulation.

Overall, the solver run with density-correction and isoAdvector mode with 4
smoothing cycles, or interFOAM mode with 2 smoothing cycles, yield results
that agree well with the reference data of Mukherjee et al. [31], with the
geometric VOF exhibiting a bubble growth rate closer to Mukherjee et al.
[31]’s simulation and the algebraic VOF following closely the experimental
data points. It is worth noting that, though the present simulations were set
similarly to the numerical setup of Mukherjee et al. [31], a few impactful
parameters from the experiment are unknown, e.g. the initial temperature field
in the liquid at the instant of bubble nucleation, frequency of bubble
nucleation and/or presence of other bubbles in the microchannel, actual
temperature distribution over the boundary walls; these have a significant
impact on the bubble growth rate, and thus deviations between experimental
data and simulation results are expected.

##  4 Results

A systematic analysis of the effect of the channel aspect-ratio and thickness
of the fins separating adjacent channels on the bubble dynamics and heat
transfer performance for flow boiling in microchannels was conducted and the
results are presented below, organised in subsections. First, the numerical
setup is described in Sec. 4.1; this is followed by an analysis of the grid
sensitivity of the numerical results in Sec. 4.2. Then, the results obtained
with liquid-only single-phase simulations are presented in Sec. 4.3, followed
by the final Sec. 4.4 where the flow boiling results are discussed.

###  4.1 Simulation setup

Figure 3: Schematic of the simulation setup and notation used in this work.
The microchannel has width $W_{ch}$ and height $H_{ch}$, with aspect-ratio
$\epsilon=W_{ch}/H_{ch}$. The microchannel base has width and height $W_{b}$
and $H_{b}$, respectively. The vertical walls (fins) have thickness $W_{f}$
but, as symmetry boundary conditions are used at the sides of the domain, only
half of the wall is modelled. At $t=0$, a vapour bubble is initialised at the
centre of the bottom wall, at a distance of $1D_{h}$ from the channel inlet;
the bubble has initial diameter of $0.2D_{h}$ and forms a contact angle of
$30^{\circ}$ with the wall. Heat is applied to the bottom boundary of the
base. The temperature profiles in both fluid and solid regions at $t=0$ are
obtained with a preliminary liquid-only simulation run till steady-state is
achieved. The image above refers to $\epsilon=2$ and $W_{f}=D_{h}/4$.

The geometrical configuration of multi-microchannel evaporators features
several parallel microchannels etched in a block of conductive material, which
are thus separated by solid walls (also called fins). The evaporator is
typically heated from below, by placing its base surface in contact with a
heater which applies a specified heat load. In experimental tests aimed at
characterising two-phase heat transfer [34], the top of the evaporator is
often covered by a poorly conductive transparent material (glass, pyrex) to
enable flow visualisation while limiting heat losses to the ambient.

In this numerical work, the geometrical configuration of a multi-microchannel
evaporator is emulated by considering a single channel unit composed of one
fluidic channel and the three connected solid wall regions at the two sides of
the channel and below it; a schematic of the flow configuration and notation
used in this work is provided in Fig. 3. Symmetry boundary conditions are
applied to the outer surface of the lateral walls to model the presence of
adjacent channels. A constant and uniform heat flux is applied at the
evaporator base.

We consider microchannels of square and rectangular cross-sections, of
constant hydraulic diameter $D_{h}=229\,\mathrm{\mu m}$ and length
$L=20D_{h}$. The microchannel has width denoted as $W_{ch}$ and height
$H_{ch}$, with the aspect-ratio defined as $\epsilon=W_{ch}/H_{ch}$; since the
hydraulic diameter is maintained constant, the channel width and height are
fully identified by $D_{h}$ and $\epsilon$: $W_{ch}=D_{h}(1+\epsilon)/2$,
$H_{ch}=D_{h}(1+\epsilon)/(2\epsilon)$. The microchannel base has width
$W_{b}$ and thickness $H_{b}$, so that $W_{b}-W_{ch}=W_{f}$ identifies the
thickness of the fins separating the channels. As boundary conditions, at the
fluid inlet $(x=0,y,z)$ a liquid-only fully-developed laminar velocity profile
of average speed $U_{l}$ is imposed, together with $T=T_{sat}$ and a zero-
gradient condition for the pressure; the solid boundary at $(x=0,y,z)$ is set
to adiabatic. At the fluid outlet $(x=L,y,z)$, zero-gradient conditions are
set for both velocity and temperature with a uniform pressure value, and the
solid boundary is adiabatic. On the bottom surface of the evaporator base,
$(x,y=-H_{b},z)$, a uniform heat flux is applied. The top boundary of the
domain $(x,y=H_{ch},z)$ is adiabatic for both fluid $(|z|\leq W_{ch}/2)$ and
solid regions $(|z|\geq W_{ch}/2)$. The three fluid-solid coupled boundaries
are identified as vertical walls $(x,y,|z|=W_{ch}/2)$ and horizontal walls
$(x,y=0,z)$; here, no-slip condition for the velocity, zero-gradient for the
pressure and a static contact angle of $30^{\circ}$ for the volume fraction
(hydrophilic walls) are imposed, as well as continuity of temperatures and
heat fluxes (see Sec. 2.1). Note that, since all the outer solid boundaries
except the base are adiabatic, all the heat load applied through the base
surface is dissipated via heat transfer to the fluid. The Bond number of the
flow is $\mathrm{Bo}=\rho_{l}gD_{h}^{2}/\sigma<0.01$, and thus gravitational
forces are neglected. To achieve fully-developed velocity and temperature
profiles as initial conditions for the two-phase simulation, a preliminary
single-phase case with only liquid is run till steady-state. Figure 3 displays
the solid walls of the channel, coloured with the steady temperature field,
for a representative case. The results of the single-phase runs will be
discussed in detail in Sec. 4.3, as they will be useful to rationalise the
two-phase results. At $t=0$, the two-phase simulation starts with steady-state
velocity and temperature fields, and a small spherical bubble of diameter
$D_{eq}=0.2D_{h}$ is initialised near the inlet, sitting along the centreline
over the bottom wall, with centroid of coodinates
$(x=D_{h},D_{eq}\cos(\pi/6)/2,0)$. Since the flow is expected to be symmetric
to the $z=0$ plane, only half of the domain is simulated and symmetry boundary
conditions are set at $(x,y,z=0)$. The two-phase simulation is run in time
till the bubble reaches the outlet section of the channel, which is on the
order of milliseconds.

To perform this study, the channel hydraulic diameter is maintained constant
and five different channel aspect-ratios are tested:
$\epsilon=0.25,0.5,1,2,4$. For each channel aspect-ratio, four different
values of the fin thickness are considered:
$W_{f}=D_{h},D_{h}/2,D_{h}/4,D_{h}/8$; these correspond to absolute fin
thicknesses ranging from $28.6\,\mathrm{\mu m}$ to $229\,\mathrm{\mu m}$. The
thickness of the base wall is maintained constant, $H_{b}=45.8\,\mathrm{\mu
m}$. The working fluid is water at $p=1\,\mathrm{atm}$ and
$T_{sat}=373.15\,\mathrm{K}$, with an inlet velocity of
$U_{l}=0.146\,\mathrm{m/s}$ (mass flux $G=140\,\mathrm{kg/(m^{2}s)}$); this
corresponds to a value of the capillary number of
$\mathrm{Ca}=\mu_{l}U_{l}/\sigma=0.0007$, as such very thin films and extended
dry vapour patches are expected to form between the elongated bubble and the
channel walls [9]. The Reynolds number of the flow is
$\mathrm{Re}=\rho_{l}U_{l}D_{h}/\mu_{l}=100$, therefore the flow regime is
laminar. A uniform heat flux of $q_{b}=100\,\mathrm{kW/m^{2}}$ is applied to
the outer wall of the microchannel base, corresponding to a base power (for
one single channel) ranging from $Q_{b}=0.08\,\mathrm{W}$ ($\epsilon=0.25$,
$W_{f}=D_{h}/8$) to $Q_{b}=0.37\,\mathrm{W}$ ($\epsilon=4$, $W_{f}=D_{h}$) and
a constant boiling number of $\mathrm{Bl}=q_{b}/(\rho_{l}U_{l}h_{lv})=0.0003$.
The material composing the solid regions of the evaporator is copper, taken
with constant properties $\rho_{s}=8940\,\mathrm{kg/m^{3}}$,
$c_{p,s}=385\,\mathrm{J/(kg\,K)}$ and $\lambda_{s}=380\,\mathrm{W/(m\,K)}$.

Following on from the results presented in Sec. 3, for this study the boiling
solver is set to run in isoAdvector mode (geometric VOF), with density-
correction enabled and 4 smoothing cycles for the interface curvature
calculation.

###  4.2 Mesh convergence analysis

The computational domain is meshed with structured orthogonal meshes made of
cubic hexahdrons. To identify the optimal mesh arrangement, a grid
independence analysis was performed for the case with $\epsilon=1$ and
$W_{f}=D_{h}/8$, with the domain length set to $L=10D_{h}$ to decrease the
computational cost. Three different meshes are tested, with grid spacing of
$\Delta=4.58\,\mathrm{\mu m},2.29\,\mathrm{\mu m},1.14\,\mathrm{\mu m}$,
corresponding to $50,100,200$ mesh elements per hydraulic diameter. Overall,
the three grids (with $L=10D_{h}$) have about 0.4, 3, and 24 million cells.

The grid independence analysis for the single-phase flow is first described.
The Nusselt number calculated over the three fluid-solid wall boundaries, at
steady-state, is considered. Local values of the Nusselt number are computed
from the simulation data as:

$\mathrm{Nu_{w}}(x,y,z)=\frac{q_{w}(x,y,z)}{[T_{w}(x,y,z)-T_{f}(x)]}\frac{D_{h}}{\lambda_{l}}$
(11)

where $q_{w}$ and $T_{w}$ are the heat flux and temperature on the horizontal
and vertical channel walls in contact with the solid region; see schematic in
Fig. 3. The fluid temperature along the channel is evaluated by integrating an
energy balance from $T=T_{sat}$ at the inlet:

$T_{f}(x)=T_{sat}+\frac{2}{\rho_{l}c_{p,l}U_{l}W_{ch}H_{ch}}\int\limits_{0}^{x}\left[\int\limits_{0}^{W_{ch}/2}q_{w}(x,y=0,z)\,dz+\int\limits_{0}^{H_{ch}}q_{w}(x,y,z=W_{ch}/2)\,dy\right]dx$
(12)

where the symmetry condition at $z=0$ has been used to simplify the integral.
From Eq. (11), a cross-sectional average Nusselt number,
$\overline{\mathrm{Nu_{w}}}(x)$, is obtained by averaging
$\mathrm{Nu_{w}}(x,y,z)$ around the cross-sectional perimeter at a specific
streamwise coordinate $x$:

$\overline{\mathrm{Nu_{w}}}(x)=\frac{1}{W_{ch}/2+H_{ch}}\left[\int\limits_{0}^{W_{ch}/2}\mathrm{Nu_{w}}(x,y=0,z)\,dz+\int\limits_{0}^{H_{ch}}\mathrm{Nu_{w}}(x,y,z=W_{ch}/2)\,dy\right]$
(13)

The profiles of the average Nusselt number along the channel for the three
meshes are plotted in Fig. 4(a). The three curves overlap and differences are
below $1\%$, thus suggesting that all meshes are sufficient to properly
resolve the single-phase flow.

Next, the performances of the grids for the two-phase flow are presented.
Owing to the bubble growth and elongation along the microchannel as time
elapses, wall temperatures and heat fluxes vary significantly during the
simulation and the Nusselt number is also a function of time,
$\mathrm{Nu_{w}}(x,y,z,t)$; note that the Nusselt number for the two-phase
flow is still calculated via Eq. (11), but replacing $T_{f}(x)$ with the
constant $T_{sat}$. To compare the grids, a spatially-averaged Nusselt number
is calculated during runtime, by averaging the Nusselt number over the three
fluid-solid walls at each time instant:

$\overline{\overline{\mathrm{Nu_{w}}}}(t)=\frac{1}{L(W_{ch}/2+H_{ch})}\int\limits_{0}^{L}\left[\int\limits_{0}^{W_{ch}/2}\mathrm{Nu_{w}}(x,y=0,z,t)\,dz+\int\limits_{0}^{H_{ch}}\mathrm{Nu_{w}}(x,y,z=W_{ch}/2,t)\,dy\right]dx$
(14)

where $L$ is the channel length. The average Nusselt number versus time for
the two-phase simulations and different meshes is plotted in Fig. 4(b). The
time in the abscissa is expressed in terms of the location of the most
downstream end of the bubble, or bubble nose, which is extracted at each time
instant and is denoted as $x_{N}$. The medium and most refined meshes yield
very similar bubble dynamics and the two curves of the Nusselt number versus
bubble nose location overlap throughout the bubble growth and elongation
process, whereas the coarsest mesh yields a $\sim 15\%$ smaller Nusselt
number, due to a larger extension of dry vapour patches which can be ascribed
to an insufficient resolution of the liquid film trapped between bubble and
walls. Therefore, the mesh with 100 cells per hydraulic diameter was selected
to run the analysis presented in the following sections.

Simulations were run on UK’s Tier-1 supercomputer ARCHER2, which features
computing nodes with two AMD EPYC 7742 2.25 GHz 64-core processors, thus
making 128 cores and 256 GB of RAM per node. Simulations were run using two
computing nodes and each two-phase run (with $L=20D_{h}$) required between
6,000 and 12,000 core-hours to complete.

Figure 4: Results of mesh convergence tests for (a) single-phase and (b) two-
phase simulations. (a) Single-phase Nusselt number at steady-state, averaged
around the cross-sectional perimeter, as a function of the streamwise
coordinate; see definition in Eq. (13). (b) Two-phase Nusselt number averaged
over all three microchannel walls, see definition in Eq. (14), plotted as a
function of the location of the bubble nose.

###  4.3 Single-phase results

In this section, the heat transfer results at steady-state obtained for
single-phase simulations run with only liquid, in the range of aspect-ratios
$\epsilon=0.25-4$ and fin thicknesses $W_{f}=D_{h}/8-D_{h}$, are presented.

####  4.3.1 Square channel and effect of fin width

Figure 5: Single-phase results for $\epsilon=1$ at steady-state. (a-f)
Contours of (a,b) temperature, (c,d) wall heat flux, and (e,f) Nusselt number,
over the (a,c,e) horizontal and (b,d,f) vertical walls of the channel, for
$W_{f}=D_{h}/8$; figures not to scale. Plots of (g) wall and base
temperatures, (h) wall heat flux, and (i) Nusselt number, averaged over the
cross-stream direction, as a function of the streamwise coordinate, for
different fin thicknesses. The legend in (i) applies also to (g) and (h).

Figure 5 illustrates contours of temperature, heat flux, and Nusselt number
over the horizontal and vertical channel walls, and averaged over the cross-
stream direction, for a representative case executed with $\epsilon=1$ and
$W_{f}=D_{h}/8$. The heat flux and Nusselt number contours emphasise the
development of thermal boundary layers over the coupled fluid-solid walls,
$|z|/W_{ch}=0.5$ in Fig. 5(c) and (e) and $y/H_{ch}=0$ in Fig. 5(d) and (f),
whereas smaller gradients appear nearby the channel top wall, $y/H_{ch}=1$ in
Fig. 5(d) and (f), owing to the adiabatic boundary condition. As a result of
the different boundary conditions, the temperature distribution over the
horizontal wall (Fig. 5(a)) is uniform in the cross-stream direction whereas
temperature decreases towards the top of the channel along the vertical wall
(Fig. 5(b)). Temperature, heat flux and Nusselt number are averaged around the
cross-sectional perimeter as indicated in Eq. (13) and the resulting
streamwise profiles are shown in Fig. 5(g)-(i), for different fin widths.
Figure 5(g) includes also the cross-stream average of the base temperature,
$T_{b}(x,z)=T(x,y=-H_{b},z)$, calculated as:

$\overline{\mathrm{T_{b}}}(x)=\frac{1}{W_{b}/2}\int\limits_{0}^{W_{b}/2}T_{b}(x,z)\,dz$
(15)

It can be seen that, though fluid enters the channel at $T_{sat}$, the average
temperature of the microchannel wall at the inlet ($x/D_{h}=0$) is a few
degrees higher than $T_{sat}$, due to the effect of axial heat conduction
which redistributes heat towards the inlet region where the fluid is colder,
and thus the temperature rise in the streamwise direction is milder than
expected in the case of negligible conjugate heat transfer effects. Since the
channel width $W_{ch}$ is maintained constant, thicker microchannel walls
($W_{f}$) yield larger wall and base temperatures as more heat is delivered to
the fluid. This is expected because the base heat flux $q_{b}$ is maintained
constant and a simple energy balance of the evaporator yields
$q_{b}W_{b}L=\overline{\overline{q_{w}}}(2H_{ch}+W_{ch})L$, such that
$\overline{\overline{q_{w}}}$ (average heat flux over all coupled fluid-solid
walls) must increase as $W_{b}$ increases, in order to deliver the increased
heat load to the fluid through the same exchange area $L(2H_{ch}+W_{ch})$. For
example, the energy balance suggests
$\overline{\overline{q_{w}}}=37.5\,\mathrm{kW/m^{2}}$ when $W_{f}=D_{h}/8$ and
$\overline{\overline{q_{w}}}=66.7\,\mathrm{kW/m^{2}}$ when $W_{f}=D_{h}$, as
confirmed by the wall heat flux magnitudes achieved in Fig. 5(h). Owing to the
high thermal conductivity of copper, wall and base temperatures are almost
overlapping in Fig. 5(g). A quick one-dimensional heat conduction calculation
across the evaporator base suggests that the base-wall temperature difference
should be $\approx q_{b}H_{b}/\lambda_{s}\approx 0.01\,\mathrm{K}$, thus
validating this observation. Nonetheless, the base-wall temperature difference
increases as the fin width is decreased, because the vertical wall temperature
decreases owing to the lower thermal resistance of the fin. Despite the
differences observed in boxes (g) and (h), the Nusselt number profiles in Fig.
5(i) exhibit little dependence on the fin width, which is expected as the
single-phase wall-fluid convective heat exchange is not directly impacted by
the wall thickness.

####  4.3.2 Effect of channel aspect-ratio

Figure 6: Single-phase results for $\epsilon=0.25$ and $4$, at steady-state.
The contours in (a,c,e,g) show (a,e) temperatures and (c,g) Nusselt numbers
over (a,c) horizontal and (e,g) vertical walls for $\epsilon=0.25$ and
$W_{f}=D_{h}/8$, while (b,d,f,h) show corresponding data for $\epsilon=4$ and
$W_{f}=D_{h}/8$; figures are not to scale. The plots in (i,j,k) present wall
and base temperatures, wall heat flux, and Nusselt number, averaged over the
cross-stream direction, versus the streamwise coordinate, for the same aspect
ratios and two selected values of the fin width. The legend in (k) applies
also to (i) and (j).

A comparison of the single-phase heat transfer performance for different
channel aspect-ratios and same fin width is provided with Fig. 6.
$H_{ch}>W_{ch}$ for $\epsilon=0.25$, and thus the shorter channel wall faces
the base, whereas $W_{ch}>H_{ch}$ for $\epsilon=4$, where the longer wall is
in contact with the evaporator base. The main differences between the two
aspect-ratios are related to the different thicknesses of the thermal boundary
layers developing over the channel walls. For $\epsilon=0.25$, the boundary
layer develops quickly over the shorter horizontal wall and the local Nusselt
number declines more rapidly, see Fig. 6(c), whereas heat convection over the
horizontal wall remains effective along the microchannel for $\epsilon=4$, see
Fig. 6(d), as the boundary layer over the wider wall developes more slowly.
The situation is reversed over the vertical channel wall, where heat transfer
is more effective for $\epsilon=0.25$, see Fig. 6(g) and (h), and the
heightwise temperature distribution becomes less uniform (Fig. 6(e)). When
comparing the contours of the Nusselt number over the wider wall, i.e. the
vertical wall for $\epsilon=0.25$ (Fig. 6(g)) and the horizontal wall for
$\epsilon=2$ (Fig. 6(d)), slightly better heat transfer is achieved for
$\epsilon=0.25$, owing to the adiabatic top boundary of the channel (at
$y/H_{ch}=1$) which allows for lower temperatures.

The streamwise profiles of cross-stream averaged temperature and heat flux in
Fig. 6(i) and (j) reveal that the wall and base temperatures are significantly
lower for the smaller aspect-ratio channel, mainly because the heat applied
through the evaporator base is smaller owing to the shorter channel width. For
example, the same energy balance for the solid region appied above yields
$\overline{\overline{q_{w}}}=13\,\mathrm{kW/m^{2}}$ when $\epsilon=0.25$ and
$\overline{\overline{q_{w}}}=70\,\mathrm{kW/m^{2}}$ when $\epsilon=4$, both
with $W_{f}=D_{h}/8$. Since $\overline{\overline{q_{w}}}$ increases with
increasing $\epsilon$, whereas the liquid mass flux is maintained constant,
higher aspect-ratio channels lead to higher fluid temperatures, thereby
resulting in larger evaporator temperatures.

The cross-stream average Nusselt number depicted in Fig. 6(k) confirms that
the fluid-solid heat transfer is relatively insensitive to the fin width, but
performances are substantially different for the two aspect-ratios despite the
fact that the hydraulic diameter and the ratio between the longer and shorter
channel sides are the same (i.e. four) for both configurations.
Counterintuitively, the configuration with $\epsilon=0.25$, where the shorter
wall is in contact with the hot evaporator base, yields a significantly higher
Nusselt number than that with the wider wall facing the heat source
($\epsilon=4$). This can be easily explained by inspection of the contours of
the Nusselt number in Fig. 6(d) and (g). The magnitudes of the Nusselt number
over the vertical wall for $\epsilon=0.25$ are comparable (or slightly higher)
to those detected over the horizontal wall for $\epsilon=4$. However, owing to
the three-side heating configuration, the heat transfer over the vertical wall
contributes twice to the cross-stream average of $\mathrm{Nu_{w}}$,
$\overline{\mathrm{Nu_{w}}}(x)=[W_{ch}\overline{\mathrm{Nu_{w,h}}}(x)+2H_{ch}\overline{\mathrm{Nu_{w,v}}}(x)]/(W_{ch}+2H_{ch})$,
and therefore the heat transfer performance over the vertical wall is twice as
important as that over a horizontal wall of same size.

This aspect is investigated in detail in Fig. 7, where the average Nusselt
number over the entire microchannel heated walls (and separately for vertical
and horizontal walls) is calculated using Eq. (14) and steady-state single-
phase data. Though the Nusselt number over the horizontal wall increases as
$\epsilon$ increases and that over the vertical wall decreases, and vice versa
as $\epsilon\rightarrow 0$, see Fig. 7(a), the overall Nusselt number is
larger at smaller aspect-ratios because, as mentioned above, the microchannel
wall with better Nusselt number (vertical wall) contributes twice to the
overall heat transfer. Shah and London [35] derived analytical solutions for
hydrodynamically- and thermally-developed laminar flows in rectangular
channels subject to nonuniform heating conditions and their results are
displayed in Fig. 7(b) together with the numerical data from the present work.
While the Shah and London [35] results emphasize that there would be no
difference in Nusselt numbers between $\epsilon=0.25$ and $4$ for uniform
heating (Case 1), the heat transfer performance increases monotonically with
$\epsilon$ for single-side heating (Case 4), whereas
$\overline{\overline{\mathrm{Nu_{w}}}}$ exhibits a non-monotonic trend for the
three-side heating identified as Case 2, with the Nusselt number being larger
at smaller aspect-ratios in the range investigated in this study. Our
numerical results follow very closely those for three-side heating, though the
magnitudes obtained with the simulations are slightly larger as the flow is
not thermally developed. The evident differences among the trends of
$\overline{\overline{\mathrm{Nu_{w}}}}$ versus $\epsilon$ for different
heating conditions observed in Fig. 7 emphasise the importance of including
conjugate heat transfer in the microchannel model.

Figure 7: Single-phase results for all aspect-ratios and fin widths, at
steady-state. (a) Nusselt number, averaged over the vertical and horizontal
walls, and averaged over all three microchannel walls (see Eq. (14)), for all
aspect-ratios and $W_{f}=D_{h}/8$. (b) Nusselt number, averaged over all three
microchannel walls, for different aspect-ratios and fin widths. Simulation
data are compared with Shah and London [35] predictions for fully-developed
laminar flow in a uniformly heated channel (Case 1), three-side heated channel
(Case 2) and single-side heated channel (Case 4).

####  4.3.3 Analysis of the entire single-phase database

In order to compare all results for different aspect-ratios and fin widths,
the whole numerical database of average base temperatures and wall heat fluxes
is compiled in Fig. 8(a) and (b). Base temperatures and heat fluxes follow
similar trends for the reasons explained above. Lower evaporator base
temperatures are achieved with smaller channel aspect-ratios and thinner fin
widths, because less heat per fluidic channel is applied to the evaporator.
The differences induced by the width of the fins are less apparent at smaller
aspect-ratios whereas they become more pronounced as $\epsilon$ increases. The
plots in Fig. 8 also include prediction curves obtained by means of a heat
transfer model for the evaporator and microchannel. Using a calculation
procedure usually applied to derive heat transfer coefficients from
temperature measurements in microchannel flow boiling experiments [5, 16], we
consider a model of the evaporator composed of two resistences in series, the
first related to heat conduction through the evaporator base,
$R_{b}=H_{b}/(\lambda_{s}W_{b}L)$, and the second related to heat convection
to the fluid. The latter is calculated as the result of two parallel thermal
resistances, one over the horizontal and the other over the vertical channel
walls. The resistance to convection at the horizontal wall is expressed as
$R_{w,h}=(h_{w,h}W_{ch}L)^{-1}$, with $h_{w,h}$ being the heat transfer
coefficient at the horizontal wall, that over the vertical wall is expressed
using the fin efficiency, $R_{w,v}=(\eta h_{w,v}2H_{ch}L)^{-1}$, with
$\eta=\text{tanh}(mH_{ch})/(mH_{ch})$ and
$m=\sqrt{2h_{w,v}/(\lambda_{s}W_{f})}$. The overall resistance to convection
is then calculated as $R_{w}=R_{w,h}R_{w,v}/(R_{w,h}+R_{w,v})$ and the total
resistance of the system is $R_{tot}=R_{b}+R_{w}$. Since the heat applied to
the system via the evaporator base is known, $Q_{b}=q_{b}W_{b}L$, the average
base temperature can be estimated as:

$\overline{\overline{T_{b}}}=\overline{\overline{T_{f}}}+Q_{b}R_{tot}$ (16)

where the average fluid temperature is calculated via an energy balance for
the flow:

$\overline{\overline{T_{f}}}=T_{sat}+\frac{1}{2}\frac{Q_{b}}{W_{ch}H_{ch}U_{l}\rho_{l}c_{p,l}}$
(17)

The predictions for the base temperature obtained with Eq. (16) are reported
in Fig. 8(a) as dashed lines. The only unknowns in the model are the heat
transfer coefficients $h_{w,h}$ and $h_{w,v}$, which are set by assuming a
constant Nusselt number value of $4$. Estimations of the average wall heat
flux can be obtained via the energy balance for the solid region used before:

$\overline{\overline{q_{w}}}=\frac{Q_{b}}{(2H_{ch}+W_{ch})L}$ (18)

which are plotted as dashed lines in Fig. 8(b). The heat transfer model
captures very well the base temperature and wall heat flux trends versus
aspect-ratio and fin width obtained with the numerical simulations, thus
confirming that lower base temperatures can be achieved by decreasing the
aspect-ratio of the microchannels. Figure 8(c) reports the ratio of heat
transferred through the horizontal or vertical wall and the base heat load,
when varying the aspect-ratio for $W_{f}=D_{h}/8$. The predictions displayed
as dashed lines are obtained as:

$Q_{w,h}=\frac{\overline{\overline{T_{w}}}-\overline{\overline{T_{f}}}}{R_{w,h}},\quad
Q_{w,v}=\frac{\overline{\overline{T_{w}}}-\overline{\overline{T_{f}}}}{R_{w,v}}$
(19)

where
$\overline{\overline{T_{w}}}=\overline{\overline{T_{b}}}-q_{b}H_{b}/\lambda_{s}$;
note that the model calculates an average wall temperature and does not
discriminate between horizontal and vertical wall temperatures. At smaller
aspect-ratios, most of the heat is transferred through the vertical wall
(above 90% for $\epsilon=0.25$). Even for the square channel ($\epsilon=1$),
the vertical wall contributes to about 70% of the heat dissipation to the
fluid, as the vertical wall counts twice on the overall heat transfer balance
owing to the three-side heating configuration. At larger aspect-ratios, heat
transfer through the horizontal wall takes over, though the vertical wall
still provides a considerable contribution, which explains the asymmetry on
the profiles of the Nusselt number around $\epsilon=1$ in Fig. 7.

Figure 8: Single-phase results for all aspect-ratios and fin widths, at
steady-state. (a) Base temperature, averaged over the entire base surface, (b)
wall heat flux averaged over all three microchannel walls, and (c) ratio of
heat transferred through the three microchannel walls and heat applied at the
base, for different aspect-ratio microchannels. The data in (c) refer to
$W_{f}=D_{h}/8$. The dashed lines in all figures depict the predictions
obtained with the heat transfer model; see Eqs. (16), (18) and (19).

###  4.4 Two-phase results

Figure 9: Snapshots of bubble growth and temperature field evolution over time
for a simulation run with $\epsilon=0.5$ and $W_{f}=D_{h}/2$. A shorter
channel ($L=10D_{h}$) was used for the sake of visualisation. Within each box,
the top image shows the bubble surface in grey and slight transparency, with
the contact lines over the horizontal and left vertical walls highlighted in
white. The right vertical wall of the channel is clipped to enable
visualisation. The middle image shows the temperature field on the inner
surface of the vertical wall, in contact with fluid. The bottom image shows
the temperature field on the outer surface, where symmetry boundary conditions
are applied.

This section presents the results obtained with two-phase flow simulations, in
the range of aspect-ratios $\epsilon=0.25-4$ and fin thicknesses
$W_{f}=D_{h}/8-D_{h}$. To illustrate the bubble dynamics and heat transfer for
a representative case, Fig. 9 shows snapshots of the bubble growth and
resulting microchannel temperatures as time elapses, for $\epsilon=0.5$ and
$W_{f}=D_{h}/2$. At $t=0$, the temperature field corresponds to the single-
phase steady-state solution. At the onset of the two-phase flow, the bubble
grows expanding the dry vapour region over the bottom surface of the channel
and, as its diameter reaches $W_{ch}$, a contact line is formed over the
vertical wall. The inner vertical wall responds instantaneously to the
presence of the contact line, see the central box in Fig. 9(c), as manifested
by the colder nearly-circular spot revealed by the temperature contours,
whereas this is not seen on the outer wall due to heat spreading. As the
bubble grows further, it elongates along the channel and the dry vapour
regions expand, leaving thick liquid lobes at the channel corners and thin
liquid layers between the bubble nose and the downstream ends of the contact
lines. As the bubble and contact lines progress along the channel, the wall
temperature reduces, exhibiting larger gradients in the proximity of the
contact lines. As the channel wall comes in contact with vapour in the dryout
regions, the wall temperature tends to become more uniform and to increase due
to the less efficient heat convection, though this happens with some delay due
to the thermal inertia of the solid.

####  4.4.1 Square channel and effect of fin width

Figure 10: Snapshots of bubble surface and temperature fields over the
evaporators walls as the bubble approaches the outlet section, for selected
channel aspect-ratios. One of the two vertical walls of the channel is clipped
to enable visualisation.

We begin the two-phase analysis with discussing the results obtained for a
reference case run with a square channel ($\epsilon=1$). Figure 10 presents
snapshots of the bubble surface and corresponding evaporator temperature
fields for $\epsilon=0.5,1,2$. For the square channel, the bubble front has a
cylindrical shape which flattens at the centre of the channel walls. Dry
vapour patches are formed near the tail of the bubble, at the centre of the
channel walls. The contours of temperature, wet fraction and Nusselt number
over the microchannel walls detected at $t=t_{end}$, when the bubble nose
reaches the outlet section, for $W_{f}=D_{h}/8$, are shown in Fig. 11(a)-(f).
At this stage, the bubble is elongated and dry vapour patches are clearly
visible on the channel walls in Fig. 11(c) and (d). Note that the wall dry
regions are identified as boundary faces where the face-interpolated liquid
volume fraction is $\alpha_{w}<0.5$, with $\alpha_{w}>0.5$ identifying wet
regions. Therefore, the boundaries between the blue and red regions in Fig.
11(c) and (d) do not have to be interpreted as sharp liquid-vapour boundaries.
As expected, the largest values of the Nusselt number coincide with contact
lines and very thin liquid films. The temperature field over the walls, Fig.
11(a) and (b), is very sensitive to the two-phase flow and the walls
progressively cool down as the contact line sweeps them. Figure 11(g)-(i)
reports the bubble equivalent diameter and the spatial average of wall wet
fraction and Nusselt number, both calculated as indicated in Eq. (14), as time
elapses for all the fin widths tested. Time in the abscissa is expressed in
terms of the location of the bubble nose $x_{N}$. The bubble growth rate
increases with the fin width, as a result of the larger amount of heat stored
by the solid regions during the previous single-phase stage, now dissipated in
the form of latent heat. As the bubble grows, dry vapour patches are formed
and expand over the walls, explaining the descending trends of
$\overline{\overline{\alpha_{w}}}$ with increasing time. For larger $W_{f}$,
the bubble nose propagates faster along the channel because the bubble grows
more rapidly, and thus thicker liquid films are left at the channel walls, as
expected from traditional lubrication theory [36]. Thicker films are less
prone to dewetting and dryout, thus mitigating the expansion of dry regions
and explaining the ascending trends of $\overline{\overline{\alpha_{w}}}$ with
$W_{f}$ observed in Fig. 11(h). The liquid film coverage of the wall explains
also the trends of the Nusselt number reported in Fig. 11(i).
$\overline{\overline{\mathrm{Nu_{w}}}}$ increases over time as the channel
wall becomes covered by a thin liquid film and both contact line and thin film
evaporation contribute to cool down the wall. However, the Nusselt number
increases at a larger rate for a thicker fin width, owing to the larger
fraction of channel wall that remains covered by a liquid film. Note that the
Nusselt number does not reach any steady-state or steady-periodic value,
because multiple nucleation cycles [37] and simulations over $100\times$
larger time-scales (order of $10^{-1}\,\mathrm{s}$) would be necessary for the
temperature field within the solid region to achieve a statistically steady
regime.

Figure 11: Two-phase results for $\epsilon=1$. (a-f) Contours of (a,b)
temperature, (c,d) liquid and vapour fraction, and (e,f) Nusselt number, over
the (a,c,e) horizontal and (b,d,f) vertical walls of the channel, for
$W_{f}=D_{h}/8$ at $t=t_{end}$; figures not to scale. The white lines in (c,d)
identify the bubble profiles extracted on the (c) $y=H_{ch}/2$ and (d) $z=0$
planes. (g) Bubble diameter versus time. (h) Wet area fraction and (i) Nusselt
number, averaged over all three microchannel walls, plotted as a function of
the location of the bubble nose along the channel. The legend in (g) applies
also to (h) and (i).

####  4.4.2 Effect of channel aspect-ratio

Figure 12: Two-phase results for $\epsilon=0.5$ and $2$, both with
$W_{f}=D_{h}/8$, at $t=t_{end}$. The contours in (a,c,e,g,i,k) show (a,g)
temperatures, (c,i) liquid and vapour fractions, and (e,k) Nusselt numbers
over (a,c,e) horizontal and (g,i,k) vertical walls for $\epsilon=0.5$, while
(b,d,f,h,j,l) illustrate the corresponding contours for $\epsilon=2$; figures
are not to scale. The white lines in (c,d,i,j) identify the bubble profiles
extracted on the (c,d) $y=H_{ch}/2$ and (i,j) $z=0$ planes.

Next, the results for two different aspect-ratios, $\epsilon=0.5$ and $2$, are
illustrated. Figure 10 shows bubble shapes and corresponding evaporator
temperature fields. Extended dry regions form over the wider channel walls,
whereas much smaller dry regions appear over the shorter walls; these features
are discussed below. Contours of wall temperatures, liquid and vapour
fraction, and Nusselt numbers for $\epsilon=0.5,2$ at same $W_{f}=D_{h}/8$ are
depicted in Fig. 12. The contours of the dry and wet area fraction, Fig.
12(c), (d), (i) and (j), confirm that larger dry regions develop over wider
walls, as observed in Fig. 10. This is a result of surface tension forces,
that arrange the bubble cross-section into circular arcs facing the shorter
walls, where a thicker liquid film develops, while the liquid-vapour interface
is rather flat along the wider walls, where a thinner liquid film is left [9,
13]. Since the capillary number of the flow is small, $\mathrm{Ca}=0.0007$,
the liquid film partially dries over both horizontal and vertical walls, but
the dry area fraction is smaller over the shorter wall owing to the thicker
film. It is interesting to note both in Fig. 12(d) and (i) that the dry region
at the upstream end of the contact line develops into two axial dry streaks at
its downstream end. This is due to the effect of capillary forces on the thin
film covering the larger wall, that create a saddle-like film profile on the
channel cross-section, which is thicker at the channel centre and exhibits a
dimple at the matching point with the static meniscus at the side [9, 38].
Dryout initiates in coincidence with this interfacial dimple where the film is
the thinnest [39, 40], as observed in Fig. 12(d) and (i) and, owing to the
hydrophilic walls, a narrow liquid ligament still exists along the wall
centreline between the two dry regions, although it eventually evaporates and
dries out, thus leaving an extended dry region. Dry patches are more extended
for $\epsilon=0.5$, because the bubble grows more slowly due to the lower wall
temperature. The white lines included in Fig. 12(c), (d), (i), and (j),
indicate the bubble profiles extracted along channel centreplanes, and reveal
that the liquid ligament (e.g., Fig. 12(c), $x/D_{h}=9-14$) is thicker than
the liquid film established in the downstream fully-wetted region
($x/D_{h}>14$). This is due to the fact that surface tension rearranges the
cross-sectional profile of the liquid-vapour interface from the fully-wetted
to the partially-dry region. In the downstream wet region, the film thickness
is determined by the meniscus at the channel corners [9], whereas the
thickness of the liquid ligament is determined by the contact angle.

Figure 13: Two-phase results for $\epsilon=0.5$ and $2$ and selected values of
fin widths. (a) Bubble diameter versus time, (b) wet area fraction and (c)
Nusselt number, averaged over all three microchannel walls, plotted as a
function of the location of the bubble nose along the channel. The legend in
(a) applies also to (b) and (c).

The temperature and Nusselt number contours in Fig. 12 are a result of the
liquid film and contact line distribution along shorter and larger channel
walls, as such the contours of $\overline{\overline{T_{w,h}}}$ (or
$\overline{\overline{T_{w,v}}}$) for $\epsilon=0.5$ are qualitatively similar
to those of $\overline{\overline{T_{w,v}}}$ (or
$\overline{\overline{T_{w,h}}}$) for $\epsilon=2$, with the Nusselt number
following analogous trends. The only apparent differences are on the contours
of temperature over the larger wall, Fig. 12(b) and (g), where the vertical
wall for $\epsilon=0.5$ restarts heating up in the dry vapour region.

A quantitative comparison of the results for $\epsilon=0.5$ and $2$ and
different fin widths is presented in Fig. 13. The vapour bubble grows more
rapidly for $\epsilon=2$ because the microchannel walls are warmer, as a
result of the previous single-phase steady-state temperature field, see Fig.
6. As a consequence, the larger aspect-ratio channel exhibits smaller dry
vapour patches and larger wet area fraction (Fig. 13(b)), although the curves
of $\overline{\overline{\alpha_{w}}}$ for the two aspect-ratios converge as
the fin width increases. Although the Nusselt numbers over the wider walls
(vertical wall for $\epsilon=0.5$, horizontal for $\epsilon=2$) are of
comparable magnitude, Fig. 13(c) shows that higher average Nusselt numbers are
achieved for a smaller aspect-ratio, since the vertical wall counts twice
towards $\overline{\overline{\mathrm{Nu_{w}}}}$. The situation is reversed for
$W_{f}=D_{h}/8$, see the green curves in Fig. 13(c), which is the case
illustrated in Fig. 12, because the bubble grows much faster for the larger
aspect-ratio channel and therefore the wet area fraction for $\epsilon=2$ is
considerably higher than that for $\epsilon=0.5$.

####  4.4.3 Conjugate heat transfer analysis

Figure 14: Two-phase temperature fields on the channel cross-section
($x=10D_{h}$) for different aspect-ratios and fin widths. Two different colour
scales are used for the temperature fields in the fluid and solid regions to
better reveal local gradients. The white lines in the fluid regions identify
the liquid-vapour interface on the channel cross-section.

To better reveal the nonuniform temperature distribution within the solid
regions, Fig. 14 shows cross-sectional temperature fields within both solid
and fluid, extracted half-way along the evaporator, for $\epsilon=0.5,1,2$,
and two selected values of fin thicknesses. The images also include the
profile of the liquid-vapour interface on the cross-section. When
$W_{f}=D_{h}/8$, Fig. 14(a), (c) and (e), the bubble grows at a slower rate
than $W_{f}=D_{h}$ and the liquid film over the wider walls is dry at the
displayed cross-section, except for $\epsilon=2$ where the bubble grows faster
owing to the warmer wall; thus, a narrow liquid ligament is left between the
dry regions at the sides, as observed in Fig. 12(d). Thicker liquid films
always cover the shorter walls at the selected cross-sections. When
$W_{f}=D_{h}$, Fig. 14(b), (d) and (f), the bubble grows more rapidly as the
channel walls are warmer and thin liquid films cover the channel walls for all
aspect-ratios; the saddle-like shape of the liquid-vapour interface over the
wider walls is clearly evident in Fig. 14(b) and (f), with minimum film
thickness regions appearing at the matching point between the static meniscus
at the corner and the thin film at the wall centre. By inspection of the
temperature field in the solid and fluid regions, it is evident that the
temperature is more uniform in the solid due to its much larger thermal
conductivity. Larger temperatures are measured in both the fluid and solid as
the aspect-ratio increases due to the increased heat load through the
evaporator base. The impact of contact lines and thin films on the solid wall
temperature are apparent. Temperature is the lowest near contact lines and
very thin films, e.g. in the correspondence of the interface dimples for
$\epsilon=0.5,2$, while temperatures are higher where thicker films cover the
walls. Wall temperatures remain relatively low in dry vapour regions, due to
the thermal inertia of the solid. Negligible temperature variation is observed
in the horizontal direction across the evaporator fins for thin walls
($W_{f}=D_{h}/8$), whereas horizontally-oriented gradients become evident for
thicker walls ($W_{f}=D_{h}$), depending on the local liquid film morphology.

####  4.4.4 Base temperatures and analysis of the entire two-phase database

The contours of temperature of the evaporator base for both single-phase
(steady-state) and two-phase ($t=t_{end}$) simulations, $\epsilon=0.5,1,2$ and
$W_{f}=D_{h}/8$ are presented in Fig. 15. As discussed previously,
temperatures increase when increasing $\epsilon$ and more heat must be
dissipated by the fluid. For single-phase flow, the spanwise distribution of
temperature is uniform. The formation of a contact line and an evaporating
film which are both advancing along the channel cool down the walls and this
effect is propagated to the evaporator base via heat conduction through the
solid. The cooling of the evaporator base is manifested by the isolines of
temperature shifting downstream in Fig. 15(b), (d) and (f). The spanwise
temperature profile for the two-phase flow stays rather uniform for
$\epsilon=0.5$ and 1 owing to the limited extension of the evaporating film,
whereas for $\epsilon=2$ the base centreline appears colder than the corners
due to the wide liquid film established over the horizontal wall; however,
part of this film eventually dries out and as time elapses the less effective
solid-vapour heat convection may lead to a local increase of temperatures.

Figure 15: Comparison of (a,c,e) single-phase and (b,d,f) two-phase base
temperatures, for $\epsilon=0.5,\,1,\,2$, and $W_{f}=D_{h}/8$; figures are not
to scale. The single-phase contours are exctracted at steady-state, whereas
the two-phase results are taken at $t=t_{end}$. The white profiles in the two-
phase temperature contours show the contact line formed by the bubble over the
horizontal wall at $t=t_{end}$.  Figure 16: Two-phase results for all aspect-
ratios and fin widths, extracted at $t=t_{end}$. Average (a) base and (b) wall
temperatures. (c,d,e) Heat fluxes, averaged over (c) the three microchannel
walls, (d) horizontal-only and (e) vertical-only walls. (f,g,h) Nusselt
numbers, averaged over (f) the three microchannel walls, (g) horizontal-only
and (h) vertical-only walls. (i,j,k) Wet area fractions, averaged over (i) the
three microchannel walls, (j) horizontal-only and (k) vertical-only walls. The
double-overline notation is here dropped for convenience.

To compare all the two-phase results at varying aspect-ratios and fin widths,
the entire two-phase database is compiled in the final Fig. 16. Base and wall
temperatures, Fig. 16(a) and (b), confirm that lower evaporator temperatures
are achieved by smaller channel aspect-ratios and fin widths, as less heat is
delievered to the evaporator through the narrower base surface, when the heat
flux is maintained constant. Therefore, the present results support the
conclusion that, for a given heat flux and evaporator width, accommodating
multiple low-aspect-ratio microchannels guarantees lower temperatures than a
few or one only high-aspect-ratio microchannel; this aspect is discussed
further in Sec. 5. This result is in agreement with the experiments of Al-
Zaidi et al. [5], who found that multi-microchannel evaporators with lower
aspect-ratio channels yield smaller wall superheat for the same base heat
flux. The heat flux at the wall-fluid boundary, Fig. 16(c)-(e), increases with
the fin width as more heat must be delivered to the fluid through the same
wall-fluid contact area. The heat flux through the horizontal wall, Fig.
16(d), increases monotonically with $\epsilon$ owing to the increasingly wider
evaporating liquid film. A similar trend is observed over the vertical wall as
the aspect-ratio is reduced from $\epsilon=4$ to $\epsilon=1$, see Fig. 16(e),
however a non-monotonic behavior of $\overline{\overline{q_{w,v}}}$ occurs as
$\epsilon\rightarrow 0$, because the bubble grows proressively more slowly due
to the lower wall temperatures, and thus extended dry vapour patches appear
over the vertical wall; see wet area fractions in Fig. 16(k). This effect is
amplified for thinner evaporator fins, as temperatures are lower thus
decreasing further the bubble growth rate. Note that the average two-phase
heat flux reaches values as high as $1\,\mathrm{MW/m^{2}}$, corresponding to
above $3\,\mathrm{W}$ dissipated by the fluid. This value is about 10 times
larger than the heat load applied to the evaporator base, which is possible
because the evaporator is not at steady-state, and liquid evaporation is fed
by the sensible heat stored within the evaporator walls. It is worth
inspecting the plots of the wet area fraction in Fig. 16(i)-(k) before
discussing those of the Nusselt number. Over the horizontal wall, dry patches
are small because when $\epsilon<1$ surface tension forces leave a thick film,
see for example Fig. 14(a) and (b), whereas when $\epsilon>1$ the bubble grows
faster thus depositing a thicker film that better resists dryout. Over the
vertical wall, channels with $\epsilon>1$ leave a thick film, whereas when
$\epsilon<1$ large dry patches appear (Fig. 14(a)) because the bubble
propagates more slowly and the thin film deposited by the bubble nose dewets
rapidly. The wet area fraction trends have a direct impact on the Nusselt
number plots in Fig. 16(f)-(h). Since dry patches over the horizontal wall are
small, $\overline{\overline{\mathrm{Nu_{w,h}}}}$ shows a monotonic ascending
trend when increasing $\epsilon$, because an increasingly larger fraction of
the channel perimeter is covered by a thin film which promotes heat transfer.
For the same reason, an analogous ascending trend occurs for
$\overline{\overline{\mathrm{Nu_{w,v}}}}$ when decreasing the aspect-ratio
from $\epsilon=4$ to about $\epsilon=0.5$. However, this trend changes when
$\epsilon<0.5$ due to a sudden drop of
$\overline{\overline{\mathrm{Nu_{w,v}}}}$. The latter happens due to the large
dry patches that form over the vertical walls when $\epsilon\rightarrow 0$, as
explained above. The fall of $\overline{\overline{\mathrm{Nu_{w,v}}}}$ at low
aspect-ratios is mitigated by larger evaporator fins, which maintain larger
temperatures thus promoting faster bubbles and thicker liquid films. The
resulting trends of the average Nusselt number versus $\epsilon$ reported in
Fig. 16(f) are mixed. For smaller evaporator fins, the overall Nusselt number
increases somewhat monotonically with the channel aspect-ratio, due to the
steep increase of the wall wet fraction. However, for thicker fins, the drop
of $\overline{\overline{\mathrm{Nu_{w,v}}}}$ as $\epsilon\rightarrow 0$ is
less severe and a nonmonotonic trend of
$\overline{\overline{\mathrm{Nu_{w}}}}$ for increasing aspect-ratios is
observed.

Al-Zaidi et al. [5] reported an increasing trend of the heat transfer
coefficient when increasing $\epsilon$ in the range $\epsilon=0.5-2$, using
the fluid HFE7100 in $D_{h}=0.46\,\mathrm{mm}$ channels with
$W_{f}=0.1\,\mathrm{mm}$. They calculated the heat transfer coefficient from
temperature and heat flux measurements via thermocouples installed below the
horizontal wall of the microchannel, and therefore their heat transfer trends
should be compared with our numerical trends for
$\overline{\overline{\mathrm{Nu_{w,h}}}}$, Fig. 16(g). The trends are in good
agreement and the numerical results confirm Al-Zaidi et al. [5] interpretation
that the higher heat transfer coefficient achieved by $\epsilon=2$ can be
explained with the presence of an extended evaporating film over the
horizontal wall.

##  5 Discussion

The results presented in the previous section emphasise that there are two
relevant parameters to describe the performance of an evaporator, the average
Nusselt number $\overline{\overline{\mathrm{Nu_{w}}}}$ and the base
temperature $\overline{\overline{T_{b}}}$, which do not necessarily follow the
same trends when varying the channel aspect-ratio $\epsilon$.

The single-phase results provided in Fig. 7 show that
$\overline{\overline{\mathrm{Nu_{w}}}}$ drops by about 33% when increasing the
aspect-ratio from $\epsilon=0.25$ to 4, to which it corresponds a
three/fourfold increase of the base superheat, see Fig. 8(a), due to the
combined effect of lower $\overline{\overline{\mathrm{Nu_{w}}}}$ and higher
fluid temperature. Nonetheless, the single-phase results exhibit a systematic
trend of increasing $\overline{\overline{\mathrm{Nu_{w}}}}$ and decreasing
$\overline{\overline{T_{b}}}$ when $\epsilon<1$, thus suggesting that
multichannel heat sinks with low aspect-ratio channels ensure better heat
transfer and lower base temperatures, when operating in single-phase regime.

The situation is less clear when the heat sink operates in two-phase flow
boiling conditions, because Fig. 16(f) indicates that the Nusselt number
exhibits non-monotonic trends versus $\epsilon$, and the base temperatures
reported in Fig. 16(a) are not at steady-state and thus some kind of model is
necessary to extrapolate their values to steady-state conditions. If the heat
transfer coefficients over the longer vertical walls for $\epsilon<1$ were of
similar magnitudes to those over corresponding longer horizontal walls for
$\epsilon>1$, it would be expected that
$\overline{\overline{\mathrm{Nu_{w}}}}$ increase monotonically for
$\epsilon\rightarrow 0$, as observed in the single-phase configuration.
However, the two-phase Nusselt numbers reported in Fig. 16(f) show mixed
trends when varying $\epsilon$ and the fin width $W_{f}$, as such it is not
possible to draw a firm conclusion about which aspect-ratio maximises two-
phase heat transfer. In microchannel two-phase flow, where slug flow and
annular flow are dominant flow patterns [41],
$\overline{\overline{\mathrm{Nu_{w}}}}$ is directly related to the thickness
and distribution of liquid films and dry vapour patches over the heated walls,
which depend on a number of hydrodynamic and thermodynamic parameters. Even
so, the beneficial effect of a higher Nusselt number for a specific value of
$\epsilon$ may be outweighed by the intrinsic thermal resistance of the whole
evaporator [5].

To further investigate this aspect, we extend here the heat transfer model of
the evaporator developed in Sec. 4.3.3 to account for a two-phase flow, where
we assume that $\overline{\overline{T_{f}}}=T_{sat}$, and consider an
evaporator of total base width $W_{tot}$ (base area $A_{b}=W_{tot}L$),
featuring $N_{ch}$ parallel microchannels. Therefore, Eq. (16) can be
rewritten as:

$\overline{\overline{T_{b}}}=T_{sat}+q_{b}A_{b}R_{tot}$ (20)

where $R_{tot}=R_{b}+R_{w}$ and these two thermal resistances are now
expressed for the entire evaporator width as:

$R_{b}=\frac{H_{b}}{\lambda_{s}A_{b}},\quad
R_{w}=\frac{R_{w,h}R_{w,v}}{R_{w,h}+R_{w,v}}\frac{1}{N_{ch}}$ (21)

where $N_{ch}=W_{tot}/(W_{ch}+W_{f})$, therefore the number of microchannels
that can be accommodated depends on their aspect-ratio and the fin width.
Equation (20) can therefore be rewritten as:

$\overline{\overline{T_{b}}}=T_{sat}+\frac{q_{b}H_{b}}{\lambda_{s}}+q_{b}A_{b}R_{w}$
(22)

where the heat convection resistance at the channel wall, $R_{w}$, is the only
parameter dependent on $\epsilon$. Using the same expressions for the heat
convection resistance over the vertical and horizontal walls of the channel
developed in Sec. 4.3.3, and taking $h_{w,h}=h_{w,v}=h_{w}$, $R_{w}$ can be
expressed as:

$R_{w}=\frac{\epsilon}{A_{b}}\frac{1+2W_{f}/[D_{h}(1+\epsilon)]}{h_{w}(\epsilon+2\eta)}$
(23)

which shows that $R_{w}\sim\epsilon$ such that the thermal resistance
decreases as the aspect-ratio is reduced, and so does the evaporator base
temperature by virtue of Eq. (20). It is possible to use Eq. (23) to estimate
the resistance to convective heat transfer per unit area, $R_{w}A_{b}$, for
the two-phase configuration studied in this paper. Using $W_{f}=D_{h}/8$ and
$h_{w}$ from the data in Fig. 16(f),
$R_{w}A_{b}=5.7\,\mathrm{[W/(mm^{2}\,K)]^{-1}}$ for $\epsilon=0.25$ and
$R_{w}A_{b}=8.7\,\mathrm{[W/(mm^{2}\,K)]^{-1}}$ for $\epsilon=4$. Therefore,
despite $\overline{\overline{\mathrm{Nu_{w}}}}\approx 27$ for $\epsilon=4$
while $\overline{\overline{\mathrm{Nu_{w}}}}\approx 12$ for $\epsilon=0.25$
(see data for $W_{f}=D_{h}/8$ in Fig. 16(f)), the thermal resistance to
convective heat transfer $R_{w}$ is still much lower for the lower aspect-
ratio channel, and Eq. (20) suggests that, at steady-state,
$\overline{\overline{T_{b}}}=373.7\,\mathrm{K}$ for $\epsilon=0.25$, versus
$\overline{\overline{T_{b}}}=374\,\mathrm{K}$ for $\epsilon=4$.

In summary, though our two-phase numerical results and steady-state model
outlined above seem to suggest that microchannels with $\epsilon<1$ promote
lower evaporator temperatures, no absolute answer emerges as to which channel
configuration yields the best heat transfer or lowest base temperature. This
is expected to depend on the interplay of heat transfer within the solid
regions of the evaporator, bubble dynamics and cross-sectional liquid film
morphology, which should all be incorporated into a novel three-zone model
[42, 43] to predict boiling heat transfer in noncircular channels, to provide
physics-based guidelines for the design of multichannel evaporators.

##  6 Conclusions

Flow boiling in a multi-microchannel evaporator was simulated by modelling one
single channel and the surrounding walls. The opensource software OpenFOAM
v2106 and the built-in geometric Volume Of Fluid method were employed, with
self-developed functions improving the estimation of the surface tension and
phase-change rate. Square and rectangular microchannels were considered, with
aspect-ratios varying in the range $\epsilon=0.25-4$ and widths
$W_{f}=D_{h}/8-D_{h}$ of the wall separating adjacent channels. The channel
hydraulic diameter, base heat flux, mass flux and fluid properties were
maintained constant throughout this work, and water at the saturation
temperature of $T_{sat}=100\,\mathrm{{}^{\circ}C}$ was the working fluid. The
channel walls were set as hydrophilic. The analysis of both single-phase and
two-phase fluid dynamics and heat transfer mechanisms led to the following
conclusions:

  * 1.

Heat is transferred to the fluid through the bottom wall and the two vertical
walls at the sides of the channel, thus the heating configuration corresponds
to a three-side heated channel. Conjugate heat transfer and channel shape have
a profound impact on heat transfer.

  * 2.

In the single-phase regime, the three-side configuration leads to increasing
Nusselt numbers as smaller aspect-ratio channels are considered, as the
vertical wall contributes twice to the average convective heat transfer
performance.

  * 3.

For the configuration studied, increasing heat loads are delivered to the
fluid when the channel aspect-ratio or fin width are increased. In single-
phase flow, this corresponds to increasingly higher evaporator temperatures,
and therefore the best heat removal performance are achieved when
$\epsilon<1$.

  * 4.

In the two-phase regime, the bubble quickly becomes elongated under the
conditions studied. Local heat transfer rates are the highest, and evaporator
temperatures the lowest, in coincidence with liquid-vapour-solid contact lines
and thin liquid films.

  * 5.

Extended thin liquid films form over the wider microchannel walls, while
thicker films are left over shorter walls. The thickness of the film over each
wall depends on the bubble speed. At high aspect-ratios, walls are warmer and
the bubble grows faster, leaving thicker liquid films that better resist
dryout. At low aspect-ratios, bubbles grow more slowly and the thin film over
the wider wall dewets leaving extended dry vapour patches which contribute
poorly to heat transfer.

  * 6.

The Nusselt number of the two-phase flow is directly related to the thickness
and morphology of liquid film and dry patches at the microchannel walls. The
trends of $\overline{\overline{\mathrm{Nu_{w}}}}$ versus $\epsilon$ are mixed,
although for smaller channel fins an ascending trend for increasing aspect-
ratios is apparent.

  * 7.

Due to the conjugate heat transfer, the heat transfer coefficient and
evaporator base temperature exhibit contrasting trends when varying channel
aspect-ratio. The present results and a steady-state heat transfer model for
the evaporator suggest that configurations with $\epsilon<1$ promote lower
evaporator temperatures even when their
$\overline{\overline{\mathrm{Nu_{w}}}}$ is below that achieved for
$\epsilon>1$, with the higher heat transfer coefficient of the latter being
outweighed by the larger overall thermal resistance of the system.

This work suggests that it is difficult to draw a general conclusion about
which channel aspect-ratio maximises boiling heat transfer, or minimises the
evaporator base temperature. This is the result of the interplay among
conjugate heat transfer in the solid, bubble, liquid film dynamics and two-
phase heat transfer, and only a prediction model incorporating all relevant
hydrodynamics and heat exchange processes can provide an optimised heat sink
configuration, which is expected to be case-dependent and vary according to
operating conditions and working fluid.

## Acknowledgements

This work is supported by the UK Engineering & Physical Sciences Research
Council (EPSRC), through the BONSAI (EP/T033398/1) grant. Calculations were
performed using ARCHER2 UK National Supercomputing Service (archer2.ac.uk),
and using the Sulis Tier-2 HPC platform hosted by the Scientific Computing
Research Technology Platform at the University of Warwick. Sulis is funded by
EPSRC Grant EP/T022108/1 and the HPC Midlands+ consortium.

## References

  * Karayiannis and Mahmoud [2017] T. G. Karayiannis, M. M. Mahmoud,  Flow boiling in microchannels: Fundamentals and applications,  Appl. Therm. Eng. 115 (2017) 1372 – 1397.
  * Tullius et al. [2011] J. F. Tullius, R. Vajtai, Y. Bayazitoglu,  A review of cooling in microchannels,  Heat Transf. Eng. 32 (2011) 527–541.
  * Agostini et al. [2007] B. Agostini, M. Fabbri, J. E. Park, L. Wojtan, J. R. Thome, B. Michel,  State of the art of high heat flux cooling technologies,  Heat Transf. Eng. 28 (2007) 258–281.
  * Cheng and Xia [2017] L. Cheng, G. Xia,  Fundamental issues, mechanisms and models of flow boiling heat transfer in microscale channels,  Int. J. Heat Mass Transf. 108 (2017) 97 – 127.
  * Al-Zaidi et al. [2021] A. H. Al-Zaidi, M. M. Mahmoud, T. G. Karayiannis,  Effect of aspect ratio on flow boiling characteristics in microchannels,  Int. J. Heat Mass Transf. 164 (2021) 120587.
  * Harirchian and Garimella [2009] T. Harirchian, S. V. Garimella,  The critical role of channel cross-sectional area in microchannel flow boiling heat transfer,  Int. J. Multiph. Flow 35 (2009) 904 – 913.
  * Wong et al. [1995] H. Wong, C. J. Radke, S. Morris,  The motion of long bubbles in polygonal capillaries. Part 1. Thin films,  J. Fluid Mech. 292 (1995) 71–94.
  * de Lózar et al. [2008] A. de Lózar, A. Juel, A. L. Hazel,  The steady propagation of an air finger into a rectangular tube,  J. Fluid Mech. 614 (2008) 173–195.
  * Magnini et al. [2022] M. Magnini, F. Municchi, I. El Mellas, M. Icardi,  Liquid film distribution around long gas bubbles propagating in rectangular capillaries,  Int. J. Multiph. Flow 148 (2022) 103939.
  * Han et al. [2012] Y. Han, N. Shikazono, N. Kasagi,  The effect of liquid film evaporation on flow boiling heat transfer in a microtube,  Int. J. Heat Mass Transf. 55 (2012) 547–555.
  * Rao and Peles [2015] S. R. Rao, Y. Peles,  Spatiotemporally resolved heat transfer measurements for flow boiling in microchannels,  Int. J. Heat Mass Transf. 89 (2015) 482–493.
  * Ferrari et al. [2018] A. Ferrari, M. Magnini, J. R. Thome,  Numerical analysis of slug flow boiling in square microchannels,  Int. J. Heat Mass Transf. 123 (2018) 928–944.
  * Magnini and Matar [2020] M. Magnini, O. K. Matar,  Numerical study of the impact of the channel shape on microchannel boiling heat transfer,  Int. J. Heat Mass Transf. 150 (2020) 119322.
  * Vontas et al. [2021] K. Vontas, M. Andredaki, A. Georgoulas, N. Miché, M. Marengo,  The effect of hydraulic diameter on flow boiling within single rectangular microchannels and comparison of heat sink configuration of a single and multiple microchannels,  Energies 14 (2021) 6641.
  * Lin et al. [2021] Y. Lin, J. Li, Y. Luo, W. Li, X. Luo, S. Kabelac, Y. Cao, W. J. Minkowycz,  Conjugate heat transfer analysis of bubble growth during flow boiling in a rectangular microchannel,  Int. J. Heat Mass Transf. 181 (2021) 121828.
  * Szczukiewicz et al. [2014] S. Szczukiewicz, M. Magnini, J. R. Thome,  Proposed models, ongoing experiments, and latest numerical simulations of microchannel two-phase flow boiling,  Int. J. Multiph. Flow 59 (2014) 84–101.
  * Roenby et al. [2016] J. Roenby, H. Bredmose, H. Jasak,  A computational method for sharp interface advection,  R. Soc. Open Sci. 3 (2016) 160405.
  * Tryggvason et al. [2011] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct numerical simulations of gas-liquid multiphase flows, Cambridge University Press, Cambridge, 2011.
  * Brackbill et al. [1992] J. U. Brackbill, D. B. Kothe, C. Zemach,  A continuum method for modeling surface tension,  J. Comput. Phys. 100 (1992) 335–354.
  * Hoang et al. [2013] D. A. Hoang, V. van Steijn, L. M. Portela, M. T. Kreutzer, C. R. Kleijn,  Benchmark numerical simulations of segmented two-phase flows in microchannels using the Volume of Fluid method,  Comput. Fluids 86 (2013) 28–36.
  * Hardt and Wondra [2008] S. Hardt, F. Wondra,  Evaporation model for interfacial flows based on a continuum-field representation of the source terms,  J. Comput. Phys. 227 (2008) 5871–5895.
  * Carey [1992] V. P. Carey, Liquid-vapor phase change phenomena, Taylor and Francis, 1992.
  * Tanasawa [1991] I. Tanasawa,  Advances in condensation heat transfer,  in: J. P. Hartnett, T. F. Irvine (Eds.), Advances in Heat Transfer, Academic Press, San Diego, 1991\.
  * Magnini et al. [2013] M. Magnini, B. Pulvirenti, J. R. Thome,  Numerical investigation of hydrodynamics and heat transfer of elongated bubbles during flow boiling in a microchannel,  Int. J. Heat Mass Transf. 59 (2013) 451–471.
  * Weller [2008] H. G. Weller,  A new approach to VOF-based interface capturing methods for incompressible and compressible flows,  OpenCFD Ltd. Report TR/HGW/04 (2008).
  * Deshpande et al. [2012] S. S. Deshpande, L. Anumolu, M. F. Trujillo,  Evaluating the performance of the two-phase flow solver interFoam,  Comput. Sci. Discov. 5 (2012) 1–36.
  * Scheufler and Roenby [2019] H. Scheufler, J. Roenby,  Accurate and efficient surface reconstruction from volume fraction data on general meshes,  J. Comput. Phys. 383 (2019) 1–23.
  * van Leer [1979] B. van Leer,  Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method,  J. Comput. Phys. 32 (1979) 101–136.
  * Issa [1985] R. I. Issa,  Solution of the implicitly discretized fluid flow equations by operator-splitting,  J. Comput. Phys. 62 (1985) 40–65.
  * Scriven [1959] L. E. Scriven,  On the dynamics of phase growth,  Chem. Eng. Sci. 10 (1959) 1–13.
  * Mukherjee et al. [2011] A. Mukherjee, S. G. Kandlikar, Z. J. Edel,  Numerical study of bubble growth and wall heat transfer during flow boiling in a microchannel,  Int. J. Heat Mass Transf. 54 (2011) 3702–3718.
  * Gamet et al. [2020] L. Gamet, M. Scala, J. Roenby, H. Scheufler, J.-L. Pierson,  Validation of volume-of-fluid OpenFOAM isoAdvector solvers using single bubble benchmarks,  Comput. Fluids 213 (2020) 104722.
  * Abadie et al. [2015] T. Abadie, J. Aubin, D. Legendre,  On the combined effects of surface tension force calculation and interface advection on spurious currents within Volume of Fluid and Level Set frameworks,  J. Comput. Phys. 297 (2015) 611–636.
  * Falsetti et al. [2018] C. Falsetti, M. Magnini, J. R. Thome,  Hydrodynamic and thermal analysis of a micro-pin fin evaporator for on-chip two-phase cooling of high density power micro-electronics,  Appl. Therm. Eng. 130 (2018) 1425–1439.
  * Shah and London [1971] R. K. Shah, A. L. London, Laminar Flow Forced Convection Heat Transfer and Flow Friction in Straight and Curved Ducts. A Summary of Analytical Solutions, Technical Report AD736260, Stanford University, 1971.
  * Bretherton [1961] F. P. Bretherton,  The motion of long bubbles in tubes,  J. Fluid Mech. 10 (1961) 166–188.
  * Magnini and Thome [2016] M. Magnini, J. R. Thome,  A CFD study of the parameters influencing heat transfer in microchannel slug flow boiling,  Int. J. Therm. Sci. 110 (2016) 119–136.
  * Hazel and Heil [2002] A. L. Hazel, M. Heil,  The steady propagation of a semi-infinite bubble into a tube of elliptical or rectangular cross-section,  J. Fluid Mech. 470 (2002) 91–114.
  * Khodaparast et al. [2017] S. Khodaparast, M. K. Kim, J. Silpe, H. A. Stone,  Bubble-driven detachment of bacteria from confined micro-geometries,  Environ. Sci. Technol. 51 (2017) 1340–1347.
  * Khodaparast et al. [2018] S. Khodaparast, O. Atasi, A. Deblais, B. Scheid, H. A. Stone,  Dewetting of thin liquid films surrounding long bubbles in microchannels,  Langmuir 34 (2018) 1363–1370.
  * Ong and Thome [2011] C. L. Ong, J. R. Thome,  Macro-to-microchannel transition in two-phase flow: Part 1 - two-phase flow patterns and film thickness measurements,  Exp. Therm. Fluid Sci. 35 (2011) 37–47.
  * Thome et al. [2004] J. R. Thome, V. Dupont, A. M. Jabobi,  Heat transfer model for evaporation in microchannels. Part I: Presentation of the model,  Int. J. Heat Mass Transf. 47 (2004) 3375–3385.
  * Magnini and Thome [2017] M. Magnini, J. R. Thome,  An updated three-zone heat transfer model for slug flow boiling in microchannels,  Int. J. Multiph. Flow 91 (2017) 296–314.