FB Crisis DNS and Porous Metasurfaces

International Journal of Heat and Mass Transfer 194 (2022) 123051
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer

journal homepage: www.elsevier.com/locate/hmt
Direct simulation of flow-boiling crisis and its porous-metasurface

control for very large dryout limit
Júlio Ferreira, Massoud Kaviany∗
Department of Mechanical Engineering, University of Michigan, Ann Arbor, MI, 48109, United States
a r t i c l e i n f o a b s t r a c t
Article history: The Zuber pool-boiling hydrodynamic instability theory establishes that the critical heat flux qCHF,Z is
Received 28 January 2022 based on rising surface vapor columns with liquid counterflow in a unit cell the of size of the critical
Revised 9 May 2022
wavelength λc determined by the Kelvin-Helmholtz and Rayleigh-Taylor interfacial instabilities. In the
Accepted 17 May 2022
flow boiling, the forced velocity ul,o creates a leading-edge region complementing these phase-buoyant
flows and breaking the periodicity of the Zuber unit cell. Here, with the direct numerical simulations
Keywords: (DNS), i.e., CFD-VOF-LES, of the Zuber unit cell for saturated water (one atm), the effect of ul,o on the
Flow boiling critical heat flux hydrodynamic dryout qCHF,h is examined. The results show that the upstream liquid flow penetrates the
Porous metasurface boiling region, forming surface liquid tracks meandering between the vapor columns, and upon qCHF,h
Zuber pool boiling CHF unit cell
these tracks become unstable causing surface dryout. For ul,o as small as 5 cm/s, the axial liquid inertia
Critical wavelength
deflects the vapor track and establishes the surface liquid track, raising qCHF,h noticeably over qCHF,Z and
Wavelength modulation
Hydrodynamic instability this is in good agreement with the available experiments. From the available theories, DNS-results and
Kelvin-Helmholtz instability dimensional analysis, a relation is found between ul,o and λc establishing the wavelength modulation
Rayleigh-Taylor instability enhancement of the qCHF,h by increasing ul,o which decreases λc . In this new relation qCHF,h is proportional
Capillary-viscous limit to u1l,/o6 . The limit on this modulation enhancement is the capillarity limit qCHF,c = 3.3 MW/m2 compared
Superheat limit to qCHF,Z = 1.1 MW/m2
Direct numerical simulation of flow boiling
Large eddy simulation of two-phase flow Next, the DNS results show that anisotropic arrangement of the vapor sites stabilizes the surface liq-
uid track, thus increasing the qCHF,h . This geometric control of the vapor sites is achieved with a 3-D
perforated porous coating, e.g., the flow-boiling canopy wick (FBCW). This geometric-modulated CHF en-
hancement is much larger than the wavelength-modulation enhancement and its vapor shear instability
limit is overcome with the use of levees (geometric confinement) allowing for q > 10 MW/m2 . However,
the FBCW has an internal capillary-viscous limit qCHF,c−v (about 5.2 MW/m2 using sintered-powder wick)
which favors small λc and currently is the bottleneck for achieving extremely large qCHF,h (which favors
large λc ).
© 2022 Elsevier Ltd. All rights reserved.
1. Introduction [6–8] have allowed for the direct numerical simulations (DNS) and
closer examination of the assumptions applied in the theoretical
In boiling crises, the competition between the liquid supply to treatments.
the heated surface and the escaping generated vapor results in lo- The Zuber two-phase hydrodynamic instability theory of sat-
cal surface dryout and this two-phase hydrodynamics has been the urated pool boiling postulates that the CHF phenomenon is trig-
subject of theoretical and experimental studies and reviews [1–3]. gered by a combination of the Kelvin-Helmholtz and the Rayleigh-
In recent years, in particular as related to the thermal manage- Taylor instabilities [9] and this was later expanded in [10,11]. Re-
ment of high heat flux devices, the control and enhancement of cently, the pool boiling CHF mechanism was simulated numeri-
the dryout limit (critical heat flux, CHF) has brought renewed in- cally (Lattice-Boltzmann Method), providing physical insights on
terest in understanding the boiling crises. In addition, advances in the triggering mechanism for the boiling crisis [12] and the roles
the computational two-phase flows [4,5] and turbulence modeling of surface wettability in reducing the surface superheat [13] and
vapor recoiling in the formation of hot spots [14]. In saturated
flow boiling, the inertia of the forced flow increases the CHF
∗
Corresponding author. above the pool boiling limit and comprehensive reviews on the
E-mail address: [email protected] (M. Kaviany). existing mechanistic CHF models are given in [15,16]. Among the
https://doi.org/10.1016/j.ijheatmasstransfer.2022.123051
0017-9310/© 2022 Elsevier Ltd. All rights reserved.
J. Ferreira and M. Kaviany International Journal of Heat and Mass Transfer 194 (2022) 123051
Nomenclature f fluid
g gas/vapor
A area (m2 ) h hydrodynamic
(ρ −ρ )gL2 i index
BoL Bond number, l σg
C inertial coefficient j index
cp heat capacity (J/kg-K) k index
D, d diameter (m) K Kolmogorov
(ρl −ρg )gL2 KH Kelvin-Helmholtz
EoL Eötvos number, σ
1/2 l liquid
ρl u2l,o lg liquid-gas phase change
FrL Froude number, g(ρl −ρg )L
max maximum
G/A thermal conductance [K/(W/m2 )] n normal
g gravity (m/s2 ) o outlet
H, h height (m) p particle, post
hlg heat of evaporation (kJ/kg) per perforation
K permeability (m2 ) ps plain surface
k thermal conductivity (W/m-K) RT Rayleigh-Taylor
L, l length (m) s surface
ρg u2g sh superheat
M dynamic pressure ratio,
ρl u2l Z Zuber
v
Mag Mach number, ug,ao
Others
p pressure (Pa)
¯ temporal average
P first velocity gradient invariant (1/s)
˜ filtered
Q heat flow rate (W), second velocity gradient invari-
perturbation/fluctuation
ant (1/s2 )
spatial average
q heat flux (W/m2 ) ∗ dimensionless
R curvature radius (m), third velocity gradient invari-
ant (1/s3 )
ρf uf L
ReL Reynolds number, flow-boiling CHF models is the interfacial lift-off model, which
μf
combines a separated two-phase model with the Kelvin-Helmholtz
S symmetric part of strain rate tensor (1/s)
and the Rayleigh-Taylor interfacial instabilities triggering the crisis
t time (s)
[17–19] with good agreement between the predictions and experi-
T temperature (K)
ments for rather narrow channels.
u velocity vector (m/s)
Although the bubble nucleation sites and bubble dynamics are
u axial velocity (m/s)
stochastic, the Zuber theory of boiling crises and its extensions are
ua speed of sound (m/s)
based on two-phase unit cells, with periodic arrangement of the
v lateral velocity (m/s)
vapor columns. This periodicity (multiple symmetry) is suitable to
V volume (m3 )
the DNS of the boiling crises. High-fidelity simulations have been
w perpendicular velocity (m/s)
applied to different aspects of the pool and flow boiling phenom-
W, w width (m)
ρl u2 L ena (e.g., cavity activation, bubble growth cycle) [5,20], however,
WeL Weber number, σl,o there are few studies of the flow boiling CHF. In [21,22], a model
x axial position (m) is proposed identifying CHF based on the local liquid film thick-
y lateral position (m) ness compared with the amplitude of the interfacial waves in mi-
z perpendicular position (m) crochannels and using a separated flow model to resolve the flow
Greek symbols field. In [23], a full 3-D CFD solver resolves the flow including mass
α void fraction transfer mechanisms, and the evaporation/condensation for differ-
δ liquid thickness (m), interfacial perturbation ent microchannel geometries, although before the CHF.
porosity Using the unit-cell based DNS, in addition to verifying the theo-
κ wavenumber (1/m) retical predictions of the flow boiling crises, allows for modulation
λ length scale (m) of the vapor column sites and explores the enhancement of the hy-
μ dynamic viscosity (Pa-s) drodynamic CHF through this control. These modulated vapor sites
ν kinematic viscosity (m2 /s) are created using a porous metasurface [24], and this is explored
ρ density (kg/m3 ) with 3-D capillary structures such as the flow-boiling canopy wick
σ surface tension (N/m) [25–27]. However, the porous metasurface has an internal capillary
ω vorticity (1/s) viscous limit which may not allow reaching the higher modulated
antisymmetric part of strain rate tensor (1/s) hydrodynamic limit.
Here, we review the unit-cell based, hydrodynamic instability
Subscripts theories of the pool and flow boiling (and their relations) and use
c capillary, critical them as the bases for relating the heat flux dryout limit to the
ca canopy vapor-site configuration. We use the DNS results to verify the theo-
ch channel, choking retical predictions and to explore the hydrodynamic CHF enhance-
CHF critical heat flux ments by vapor-site modulation. Snapshots of phase distribution
c −v capillary-viscous are used to monitor time variations of the liquid track width in
e evaporator the leading-edge close to the heated surface and identify its even-
tual dryout. Three hydrodynamic-limit regimes are defined, the
2
wavelength, the geometric, and the geometric-confined modula-

tion regimes, and their boundaries are explored analytically using
first principles. The related porous metasurfaces allowing for the
modulations of the vapor sites and liquid supply, and the capillary-
viscous limit are also treated. The quantitative results are for satu-
rated water at one atm.
2. Boiling hydrodynamic instability theories
2.1. Zuber pool boiling CHF theory and unit cell
Initially separated two-phase flow has interfacial instabilities

associated with it when there is relative motion between the
phases [28]. The instabilities depend on the flow condition and can
heavily impact the heat and mass transfer [29]. The linear stability
analysis considers the interface between liquid and vapor (gas) to
be planar and horizontal with the phases moving with axial veloc-
ities ug and ul . The derivation is summarized in Appendix A with
Fig. 1. Spatial distribution of vapor columns, with counter flowing liquid, on a hor-
details available in [29] and [30]. When the less dense phase (gas) izontal flat plate in the Zuber critical heat flux model. This is the Zuber unit cell.
is on top, the system is initially stable. The Kelvin-Helmholtz (K-H)
instability arises when the relative velocity is larger than a critical
value |ug − ul | = uc . The corresponding K-H critical wavelength is aλc , where a = 1/2. The square arrangement proposed by Zuber
1 / 2 was later confirmed by Sernas et al. [32]. The CHF is then calcu-
σ
λKH = λc = 2π . (1) lated as the heat removed from this rising vapor columns
g ( ρl − ρg )
Ag
The viscous effects are considered in [31]. qCHF = ρg vg,Z hlg , (6)
A
When the more dense phase is on top, the Rayleigh-Taylor (R-T) Z
instability arises with the critical wavelength (the most dangerous and the vapor velocity is found from K-H wavelength
wavelength)
1 / 2 1/2
2π σ
σ vg,Z = , (7)
λRT,d = 31/2 2π = 31/2 λc . (2) λKH ρg
g ( ρl − ρg )
So, the K-H and R-T critical wavelengths are related. The most dan- with the area ratio
gerous R-T 3-D waveforms are arranged in a square lattice with

Ag π a2 λ2c
each unit-cell side equal to the 2-D wavelength of Eq. (2) [29]. = . (8)
A 4λ2c
Equations (1) and (2) can be expressed with the Bond (or Eötvos) Z
number (would become equal to 2π and 31/2 2π ), which is the ra- Substituting this velocity and area ratio into Eq. (6) results in
tio of buoyancy and surface tension forces the pool-boiling CHF
g(ρl − ρg )λ2c Weλ 1 / 2
Boλ = = , (3) 2 π 3 a4 σ
σ 2
Frλ qCHF = ρg hlg , (9)
λKH ρg
where the Froude and Weber numbers are a ratio of the inertial
and gravity forces and a ratio of the inertial and surface tension with λKH = λRT,d = π dg . Substituting this into Eq. (6) we find the
forces, respectively Zuber CHF
1 / 2 1 / 4
ρl u2l,o ρl u2l,o λc π σ (ρl − ρg )g
Frλ = ; Weλ = . (4) qCHF,Z = ρg hlg . (10)
g(ρl − ρg )λc σ 24 ρg2
The capillary length lc is defined by the Bond number of unity, The wavelength adopted by Zuber corresponds to a Bond num-
Boλ =1, and signifies requiring an extra external force to keep the ber Boλ = 81 and for saturated water at one atm, qCHF,Z = 1.1
liquid and vapor phases separate MW/m2 .
1 / 2 When it is possible to modulate the critical wavelength, as in
σ
lc = . (5) [33] using a porous metasurface, the pool-boiling hydrodynamic
(ρl − ρg )g CHF can be enhanced according to
For saturated water it is lc = 2.5 mm. The K-H and R-T wavelengths π σρ 1/2
g
are proportional to the capillary length. qCHF,h = hlg , (11)
8 λc
Different pool-boiling critical heat flux (CHF) mechanisms are
summarized in [29]. The Helmholtz instability mechanism states where λc ≤ λRT,d is the modulated wavelength. The smaller wave-
that the generated bubbles create unstable columnar vertical length corresponds to smaller vapor columns, allowing for a larger
liquid-vapor interface and the CHF occurs when these columns im- vapor velocity when the K-H instability is triggered. The qCHF,h en-
pede the counterflow of liquid, ceasing the liquid supply. Zuber hancement achieved by the wavelength modulation [discussed be-
[9] adopted this mechanism as a basis for his theory for a flat low and for pool boiling reported in [33]] is similar to the pool-
plate. boiling heater-size effect. Measured qCHF,h over the qCHF,Z have
As shown in Fig. 1, the vapor columns are arranged in a square been reported using progressively smaller heaters, and the smallest
array with side λc = λRT,d from Eq. (2), and column diameter dg = heater size is considered to be the capillary length lc [10,34].
3
Fig. 2. Schematic representation of the liquid- and vapor-phase flow distribution in the transition from (i) pool boiling to (ii) flow boiling. (a) The 2-D interfacial lift-off
model [17,18] (narrow channel). (b) 3-D DNS leading-edge surface liquid track model (wide channel). The models are based on the Zuber unit cell, and the dominant forces
governing the CHF are the buoyancy, surface tension, inertial, and viscous.
2.2. Flow-boiling CHF theories terfacial pressure force, lifting the liquid away from the surface. In
order to perform this balance, it employs three submodels. (i) A
Differing from the pool boiling where the pressure gradient and separated-flow model to determine the velocities and thicknesses
the liquid momentum are in the gravity direction (z direction), the of liquid and vapor in the channel. (ii) A hydrodynamic stability
forced, axial inertial component of flow boiling breaks the symme- model which employs the results from the previous model to pre-
try of the counter, columnar two-phase flow of the Zuber unit cell dict the interfacial shape and the critical wavelength of the vapor
shown in Fig. 1. The smallest axial flow fundamentally alters the layer. (iii) A CHF trigger mechanism model, determined consider-
boundary conditions. ing both the critical wavelength of the vapor and the Helmholtz
Figure 2 illustrates this transition from pool to flow boiling, us- wavelength in the wetting front [19]. This is identical to the Zu-
ing the Zuber unit cell as the basis. Predicting the flow-boiling CHF ber formulation, where the surface is modelled as square unit
has traditionally been done empirically. Although the correlations cells containing vapor jets surrounded by liquid, with the jet di-
can provide useful results, they cannot be used outside the lim- ameter half the critical vapor wavelength from the second sub-
its of their experiments, so we consider the first-principles, phys- model, dg = λc /2. The implementation of this model is discussed
ical mechanisms triggering the CHF. The liquid supply mechanism in [26,27].
in flow boiling is essential in defining the mechanisms responsible The interfacial lift-off model is based on experimental measure-
for the boiling crises. Four different mechanisms have been pro- ments in a narrow channel (W < λc ) with a single row of vapor
posed and are widely accepted [15], namely, boundary-layer sep- jets. Figure 2(a) schematically illustrates the liquid and vapor flows
aration model [35], bubble-crowding model [36], sublayer dryout in a narrow channel for pool and flow boiling conditions. The va-
model [37], and interfacial lift-off model [17,18]. The last one states por expansion occupies the entire bottom portion of the channel,
that the liquid is supplied from the top in the Zuber unit cell, how- diverting the liquid upward and preventing the formation of a sur-
ever in a moving frame of reference [19]. This does not address the face liquid track. Liquid supply is limited to the wetting fronts. The
establishment of a leading-edge liquid track on the surface. This critical vapor wavelength λc corresponds to the vapor momentum
axial, forced surface liquid flow meanders around the vapor gener- lifting the interface, interrupting the liquid supply.
ated and escaping the surface. This fundamental change from the The axial velocities u f ( f = g or l for vapor and liquid) are found
pool to the flow boiling leads to a jump in the critical heat flux from the 2-D separated-flow submodel and employed in the fol-
with the smallest forced flow (verified with the DNS). lowing submodel (hydrodynamic stability theory) to determine λc .
The interfacial lift-off model and the leading-edge, liquid track The ascending vapor velocity and the liquid make-up velocity are
shear instability model, shown in Fig. 2, are discussed below. found from the third submodel. In [19], the effect of the axial flow
The interfacial lift-off model [Fig. 2(a)(ii)] postulates that CHF on the wetting fronts mimics the behavior of a moving wall gen-
occurs when the momentum of the generated vapor produced in erating columnar vapor jets, similar to the hydrodynamic instabil-
small portions of the surface called the wetting fronts, responsi- ity present in Zuber analysis (the only difference being the moving
ble for providing liquid to the heated surface, overcomes the in- frame of reference).
4
Similar to the pool boiling qCHF,h of Eq. (11), the interfacial lift-
off theory [17–19] of flow-boiling dryout results in
⎡ 1/2
⎤
ρl + ρg
πσ 1/2
⎢ ρl ρg ⎥
qCHF,h = ρg hlg ⎣ ρ ⎦, (12)
8 λc 1 + ρg (16π−π )
l
where λc is found from the 2-D separated flow model and interfa-
cial instability submodels [19]
⎡ ⎧ 2
ρl ρg (ug − ul )2 ⎨ ρl ρg (ug − ul )2
λc = 2π ⎣ +
2σ ρl + ρg ⎩ 2σ ρl + ρg
1/2 −1
(ρl − ρg )g cos(φ )
+ , (13)
σ
Fig. 3. Wavelength modulation of pool and flow-boiling dryout heat flux based on
where ρ f = ρ f tanh (κ H f ) is the modified phase density, H f is
−1
the Zuber unit-cell model. The pool-boiling results [33] are given by Eq. (11) and the
the height occupied by the phase. The liquid and vapor phase ve- flow-boiling interfacial lift-off theory results [17,18] are given by Eq. (12), and they
are identical. The water-adjusted pool-boiling experimental results of [33] using a
locities are calculated as modulated porous coating are also shown.
qx
ug = , (14)
ρg δ c p,l Tsub + hlg
this transition of the two-phase flow and distribution from the
ul,o H qx pool (i) to flow boiling (ii), under a wide-channel (Wch > λc ) treat-
ul = − , (15)
H−δ ρl (H − δ )(c p,l Tsub + hlg ) ment is rendered and fundamental forces outlined. The Zuber va-
por columns of pool boiling are deflected in the flow direction,
where δ is the thickness of the vapor layer, from the separated while a surface liquid track is formed in the leading-edge region
flow model. The last term on the right of Eq. (12) is from the as- and meanders around the vapor columns. This surface liquid track
cending vapor velocity and is related to the Helmholtz instability provides the axial liquid supply to the heated surface. This con-
2π σ 1/2 trasts with the liquid supply in the interfacial lift-off theory with
1/2
ρl + ρ g the submodel based on a narrow channel, where the liquid supply
vg,H − vl,H = . (16)
λH ρl ρg is from the top at the wetting fronts, triggering a Rayleigh-Taylor
instability (liquid on top of vapor).
The continuity equation gives
The two-phase flow rendered in Fig. 2(b)(ii) is hydrodynami-
ρg π cally stable, with the surface liquid track formed in the leading-
vl,H = v . (17)
ρl 16 − π g,H edge region beneath the vapor domain. As the vapor accelerates
downstream, however, the interface relative velocity reaches a crit-
The relationship between qCHF,h and the critical interfacial
ical value, triggering the Kelvin-Helmholtz instability. This mecha-
wavelength given in Eq. (12) for flow boiling in the interfacial lift-
nism is dubbed the leading-edge surface liquid track instability and
off model is functionally identical to Eq. (11), derived in [33] for
proposes the dryout is reached when the vapor shear causes the
pool boiling. Both express how qCHF,h can be modulated by λc
surface liquid track to shear break with the K-H instability down-
and we will later discuss the control of λc with metasurfaces.
stream. The wide-channel condition Wch > λc allows for the lateral
Equation (11) was derived from the momentum balnce of vapor
periodic treatment and formation of the liquid track, fundamen-
emerging from the surface with a modified wavelength λc .
tally altering the dryout mechanism. This is guided by direct nu-
This wavelength modulation, achieved experimentally through
merical simulation of flow boiling using the Zuber unit cell as the
porous metasurface with a periodic strucutral modulation of length
basis, and existing experimental results.
λc supersedes the dependence only on the thermophysical proper-
ties given by the Rayleigh-Taylor wavelength, Eq. (3). This allows
for the CHF enhancement beyond the Zuber limit for plain surface 3. Direct numerical simulations (DNS)
[33].
The interfacial lift-off theory identifies a similar dependency of Figure 2 (b)(ii) renders the computational domain used in the
the qCHF,h on the interfacial wavelength, obtained from the stabil- direct numerical simulations (DNS) of flow-boiling and its crisis
ity analysis of the 2-D two-phase flow using the separated flow (surface dryout). The inlet liquid velocity ul,o is prescribed up-
model. The forced flow ul,o alters the characteristic wavelength λc , stream, on the heated metasurface (z = 0) the vapor and liquid
enhancing the qCHF,h similar to that achieved by λc modulation in flows are prescribed for a given heat flux, lateral (x − z plane) pe-
pool boiling [33]. riodicity is used, shear-free top surface (z = Hch ), and continuous
Equations (11) and (12) have identical λc dependence as shown flow condition are used at the outlet downstream. Saturated wa-
in Fig. 3. In [33] a non-uniform porous coating was used with pen- ter (one atm) is assumed. The related two-phase fluid mechanics
tane as the fluid to modulate the λc and their experimental re- (transient, 3-D Navier-Stokes equations) are solved numerically us-
sults (adjusted for water properties) are shown in Fig. 3, with good ing ANSYS Fluent [38]. The treatment of turbulence is with Large
agreement to Eq. (11). Eddy Simulations. The unit-cell geometric parameters are related
While Fig. 3 shows the flow boiling qCHF,h is equivalent to to heat flux according to Eq. (11), which is for pool boiling, how-
the surface wavelength modulation of pool boiling, we expect the ever this critical heat flux depends on forced flow velocity ul,o .
anisotropic liquid inertia ul,o to alter the dryout crisis differently, The computational domain contains two axial and two lateral unit
and here we suggest a leading-edge surface liquid track insta- cells, with lateral periodicity. A leading and a trailing-edge re-
bility guided by direct numerical simulation results. In Fig. 2(b), gion are added for more realistic upstream and downstream flow
5
Table 1
CFD condition and dimensionless numbers for different velocity regimes in the wavelength modulation regime.
Velocity regime q, MW/m2 ul,o , m/s λc , mm Rel Reg Frλ Weλ Boλ Mag M
Low 1.35 0.05–0.25 17.7 3 × 103 –15 × 103 5.7 × 103 0.1–0.6 0.7-18 50 0.01 0.45
Moderate 1.75 0.25–1 9.4 7.9 × 103 –32 × 103 3.9 × 103 0.8–3.3 9.6-153 14 0.02 0.04-0.75
High 3 1–3 2.98 10 × 103 –30 × 103 2.1 × 103 5.8–17.5 49–439 1.4 0.03 0.015-0.14
conditions. The vapor-column metasurface velocity and countering combination of both [39]. In the VOF treatment, a single continu-
liquid flow are related to heat flux through ity and momentum conservation equations are used with an addi-
qA qA tional equation for the phase identifier, i.e., the volume fraction α
vg,o = ; vl,o = , (18)
ρg hlg Ag ρg hlg (A − Ag ) ∂
(αρg ) + ∇ · (αρg u˜ g ) = 0. (20)
where A is the unit-cell heated surface area and Ag the vapor col- ∂t
umn area, shown in Fig. 2(b). As an example, for saturated wa-
The phasic properties are used for the VOF properties, for example
ter at one atm and q = 1.75 MW/m2 , vg,o = 8.84 m/s and vl,o =
0.95 mm/s. As expected, the vapor velocity is about four orders-of-
ρ f = αρg + (1 − α )ρl , (21)
magnitude larger than the liquid.
An important distinction is that, in the channel, it is assumed
that the two-phase flow is under saturated liquid-vapor condi- μ f = αμg + (1 − α )μl . (22)
tion (no evaporation or condensation), and the predicted results
The simulation time step restriction is based on the convection
are compared with the saturation boiling experiments. No phase
term, namely the Courant-Friedrichs-Lewy (CFL) condition [40],
change is included in the two-phase flow, therefore empirical mod-
which states that the Courant number should be less than 0.25
els for evaporation and condensation are not used. For the case
of a plain surface, discussed in Section 4, the evaporation and u˜ f v˜ f ˜f
w
Co = t + + ≤ Comax . (23)
the return of liquid to the surface are based on the Zuber unit x y z
cell model/theory as surface boundary conditions under saturated
liquid-vapor condition. For the porous metasurfaces, discussed in
3.3. Large eddy simulation (LES)
Section 6, the vapor generation occurs in the evaporator layer wick
of the FBCW, where the liquid is supplied to the evaporator layer
For a more accurate calculation of the two-phase flow turbu-
through the porous posts which in turn are irrigated by the porous
lence, especially the flow adjacent to the heated metasurface and
and perforated canopy. The liquid-film evaporation in the evapo-
the state of the surface liquid track, the CFD simulations were con-
rator layer is modeled with one-dimensional heat conduction and
ducted using the LES [41]. In the LES, the large turbulent eddies are
the liquid film surface is assumed to be at the saturation tem-
resolved directly, while the small (i.e., subgrid) eddies are modeled.
perature. So, the superheat is within the monolayer wick which
The momentum and mass are transported mostly by the large ed-
is under local thermal equilibrium between the solid and liquid
dies which are significantly affected by the geometry and bound-
phase and the vapor is also at the saturated temperature (not su-
ary conditions, while the small eddies tend to be more isotropic,
perheated).
and the conventional turbulent models are more suited for their
3.1. Dimensionless numbers prediction.
A filter is applied to the Navier-Stokes equations, effectively cut-
The liquid and vapor inertial forces are scaled with their respec- ting off eddies smaller than the grid size. This filtering is analogous
tive viscosities in the Reynolds number for the liquid and vapor to the finite volume discretization.
phases. As mentioned earlier, the role of the liquid-vapor interfa- These filtered Navier-Stokes equations are
cial tension is presented by the Weber number using the liquid in- ∂ρ f ∂
ertia. For the role of gravity (buoyancy), the Froude number using + ρ f u˜ f,i = 0, (24)
∂t ∂ xi
the liquid inertia is used in Eq. (4). The compressibility limit of the
vapor is addressed with the Mach number using the vapor speed ∂ ∂ ∂σ ∂ p˜ ∂τi j
ρ u˜ + ρ u˜ u˜ = i j − − , (25)
of sound (for Mach number smaller than 0.3, the compressibility ∂ t f f,i ∂ x j f f,i f, j ∂xj ∂xj ∂xj
effects are neglected). The Zuber unit-cell size presented by wave-
length λc is used for the length scale. In addition to Eq. (4), the where σi j is the viscosity stress tensor component and τi j is
Reynolds numbers, Mach number [4], and vapor to liquid inertial the subgrid-scale stress component, described in Appendix B. The
ratio are defined as subgrid-scale resulting from the filtering operation requires mod-
ρl ul,o λc ρg vg,oλc eling. As in Reynolds-averaged Navier-Stokes (RANS) equations, the
Rel = , Reg = , Boussinesq hypothesis is employed and the stresses are [42]
μl μg
1
vg,per ρg u2g τi j − τkk δi j = −2μ f,t S̄i j . (26)
Mag,per = ,M= . (19) 3
ua ρl u2l
Alongside the void-fraction, the Q-criterion, a vortex identifi-
The range of these dimensionless numbers are summarized in cation method, also shows that vortices are limited to the vapor
Table 1 for the range of liquid velocity ul,o , typical heat flux, and phase. The vorticity is
the Zuber unit-cell presented by wavelength λc .
ω = ∇ × uf , (27)
3.2. Treatment of liquid-vapor interface: Volume of fluid (VOF) and its presentation requires the velocity gradient tensor invari-
method ants, P , Q, and R. These invariants are the coefficients of the cubic
characteristic polynomial det (∇ u f − λI ), i.e.,
Due to its mass conservation characteristics, the VOF method is
preferred over the level-set method, but authors have employed a λ3 + Pλ2 + Q λ + R = 0, (28)
6
Fig. 4. (a) Predicted distribution of the instantaneous, constant interfacial shear-flow turbulent kinetic energy Ē f,t contours in the axial midplane for plain surface for q =
1.75 MW/m2 , ul,o = 0.5 m/s. (b) Contour of constant time-averaged turbulence kinetic energy Ē f,t . Instantaneous contours of constant (c) invariant Q = 1 1/s2 , and (d) eddy
viscosity μ f,t . The results are for saturated water at one atm, and y = 0.
P 2 −Si j S ji +i j ji main within the vapor phase. Vortical structures dissipate turbu-
where P = −Sii , Q= 2 , and R=
lent kinetic energy [45]. Similarly, the eddy viscosity distribution is
−P 3 +3PQ−S i j S jk Ski −3i j jk ki
2 . For incompressible flows, the first shown in Fig. 4(d). It is calculated using Eq. (26), is proportional to
invariant is P =0 (since ∇ · u f =0), so the so-called Q-criterion is the Reynolds stresses normalized by the mean strain rate. Its peak
derived based on the second invariant [43], and is is observed just beneath the interface and is compatible with the
1 distribution of the vortices (and the large gradient components)
Q= 2 − S2 , (29) occurring mainly within the vapor phase. These analyses of the
2
turbulent quantities were repeated for all CFD runs with similar
where S is the symmetric part of the strain rate tensor and is
trends found regarding the turbulent vortical structure and distri-
its antisymmetric part (vorticity tensor). The Q-criterion states that
bution.
Q >0 represents a vortex. There are other methods for vortex iden-
tification available [44], however we use Eq. (29). Further discus-
sions of the turbulence treatment is given in Appendix B. 3.4. Dryout limit from the DNS
Figure 4 (a) to (d) show the distributions in the axial mid-
plane (y = 0) of the turbulence quantities for q = 1.75 MW/m2 and The direct simulations are used to verify the CHF wavelength-
ul,o = 0.5 m/s. The black line denotes the liquid-vapor interface. modulation regime for the plain surface, while varying the liquid
Figure 4(a) and (b) show that the turbulence is generated at the velocity ul,o . We begin by choosing a CHF and finding λc according
liquid-vapor interface and it is noticeably more intense in the va- to Eq. (11). Several iterations are required to determine the rela-
por phase, this is in agreement with recent studies of two-phase tion between ul,o and qCHF,h . In these iterations, either ul,o is varied
LES that show the occurrence of interface-generated turbulent ed- while q is kept constant or vice-versa (and consequentially unit-
dies [45–47]. The peak in the turbulent kinetic energy Ē f,t is in the cell is changed). The outcome is the ul,o corresponding to qCHF,h
vapor phase in the leading-edge region, where the incoming liq- calculated using Eq. (11), which depends on the thermophysical
uid flow is severely deflected and broken by the vapor columns. In properties and the Zuber unit cell, presented by the wavelength
this region, the vapor inertia overcomes the surface tension to en- λc . So, the inlet velocity affects the λc , thus altering qCHF,h .
ter the liquid channel flow. Near the liquid inlet, the relative liquid- So, for a given q, the λc is determined from Eq. (11) and for a
vapor velocity is the largest, because the liquid has not yet acceler- plain surface geometry, simulations are conducted with this wave-
ated and the vapor velocity has not changed from its large magni- length and progressively smaller ul,o to find the ul,o corresponding
tude prescribed by the heat flux boundary condition, four-order-of to the identified surface dryout (within the sensitivity of the prac-
magnitude larger than the liquid. tical ul,o or q steps in the simulation).
Figure 4 (c) shows isosurface of the Q-criterion for Q = 1 1/s2 To determine dryout, the width of liquid track wl in the sim-
superposed on the void fraction distribution in the axial mid- ulated unit-cell based domain of the heated surface is used. In-
plane. All the vortical structures identified by the Q-criterion re- stantaneous as well as running-time average values are used. The
7
Fig. 5. (a) Snapshots showing the liquid track width wl at the location marked by the red strip. Time variations of the normalized liquid track width for (b) q = 1.35
MW/m2 and ul,o = 0.05, 0.07 and 0.09 m/s, (c) q = 1.75 MW/m2 and ul,o = 0.25, 0.5 and 1.0 m/s, and (d) q = 3 MW/m2 and ul,o = 1, 1.5 and 2.5 m/s, respectively. The results
are for saturated water at one atm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
phase distribution in the x − y plane at elevation z = 0.25 mm is lack of it. Figure 5(b) to (d) show the dimensionless-time variations
used to detect the liquid phase and a phase marker is used to of the dimensionless liquid track width w∗l , for different heat
measure this width. The liquid-phase volume with α < 0.5 has a fluxes with progressively smaller inlet liquid velocities. In Fig. 5(b),
marker value of 1, while the vapor-phase is marked with a value the heat flux is q = 1.35 MW/m2 , corresponding to λc = 17.7 mm
of 0. The liquid track width is determined as the summation of [Eq. (11)], and ul,o = 0.05, 0.07 and 0.09 m/s. The running-time av-
the marked pixels, and presence of vapor reduces wl . The sum- erage from Eq. (32) is shown with broken lines. The surface dryout
mation is done over a rectangular domain (xy) and is then criterion of w∗l ≤ 0.1, marked by the red horizontal line, appears
normalized by the integral length of the liquid track which is to present the breakup/interruption of the liquid track. Based on
y = λc (31/2 − 1/2 ) (x = y/10), as shown in Fig. 5(a), i.e., this, the dryout is marked for ul,o = 0.05 m/s, since the running
mean is consistently below this dryout criterion. Figure 5(c) is for
x y
1 1 f = 0, α ≥ 0.5 q = 1.75 MW/m2 , corresponding to λc = 9.4 mm, and ul,o = 0.25,
w∗l = f , . (30)
x λc (3 − 1/2 )
1 / 2 f = 1, α < 0.5 0.5 and 1 m/s. Again, dryout is observed for the lowest velocity,
i i
i.e., 0.25 m/s. Figure 5(d) shows the liquid track width variations
Figure 5(b) to (d) show the time variations of the instantaneous for q = 3 MW/m2 , corresponding to λc = 2.98 mm, and ul,o = 1,
and running-time averaged dimensionless liquid track width. In or- 1.5 and 2.5 m/s. Dryout is identified to occur when the liquid ve-
der to account for the different velocities and wavelengths, a di- locity is ul,o = 1 m/s.
mensionless time, normalized by the fluid transit time is used, i.e., The DNS allows for using Zuber pool-boiling surface unit-cell
assumption and then adding the forced flow ul,o . The case of ul,o =
tul,o
t∗ = . (31) 0 is not realized, but liquid velocities as small as 5 cm/s can be
λc simulated without major computational time challenges. We will
The dashed lines denote the running-time average which is aver- further discuss this in the next section.
aging over elapsed time t of n steps Figure 6 shows the snapshots of the phase distributions for the
nine flow-boiling conditions shown in Fig. 5(b) to (d). The liquid-
1 ∗
w¯ ∗l = w l . (32) vapor interface is green and the surface liquid track is blue.
t n
When the liquid supply was interrupted for a time period long 3.5. Temporal and spatial resolution convergence
enough to trigger the rapid rise in the surface superheat observed
in the boiling crisis, dryout occurs. Here, this criterion for the dry- The finest grid size investigated was xi = 0.15 mm. Increas-
out is when w¯∗l < 0.1. These are marked in Fig. 5(b) to (d). ingly coarser grid sizes were evaluated and the value of xi =
Figure 5 (a) defines the downstream location in the computa- 0.25 mm was selected since it presented an error of 1.9% at a rea-
tional domain where the liquid track is evaluated for dryout or sonable computational cost. Although at a reduced computational
8
Fig. 6. Snapshots of phase distribution and the liquid track for the nine DNS conditions of Fig. 5, namely: q = 1.35 MW/m2 and ul,o = 0.05, 0.07 and 0.09 m/s, q = 1.75
MW/m2 and ul,o = 0.25, 0.5 and 1.0 m/s, and q = 3 MW/m2 and ul,o = 1, 1.5 and 2.5 m/s, respectively. Velocities marked in red indicate cases for which dryout was observed.
The results are for saturated water at one atm. The full video is available in the supplementary materials. (For interpretation of the references to colour in this figure legend,
the reader is referred to the web version of this article.)
Table 2 examine this as used in the interfacial lift-off submodel as it is

Variation of running-time averaged dimension-
found by the DNS.
less liquid track (Eq. (32)) for q = 1.75 MW/m2
and ul,o = 1 m/s. In the lift-off submodel, Eq. (13) gives the dependency of the
critical wavelength λc on the velocity of the liquid-phase ul , given
|w̄∗l −w̄∗l (xi =0.15 )|
xi (mm) w̄∗l w̄∗l (xi =0.15 ) (%) by Eq. (15). Equation (13) indicates that buoyancy and surface ten-
0.15 0.303 0.0 sion forces are balanced to determine the critical wavelength. At
0.25 0.309 1.9 ρl ρg (ug −ul )2
high velocities, when the gravity is negligible, i.e., >>
0.40 0.343 13.2 2σ ρ +ρg
0.50 0.376 24.1 1/2 l
(ρl −ρg )g cos(φ ) = 400 1/m, or when the relative velocity ug −
σ
ul ≥ 6.37 m/s, this becomes
cost, coarser grid sizes (xi = 0.40 and 0.50 mm) resulted in inac-
2π σ (ρl + ρg ) 2π σ
curate measurements (error too high). As an example, Table 2 sum- λc = ≈ , (33)
marizes the variations of the running-time averaged liquid track ρl ρg (ug − ul )2 ρg (ug − ul )2
w̄∗l , defined by Eq. (32), with respect to progressively coarser grid suggesting a Weber number relationship, or dominance of the sur-
size, for q = 1.75 MW/m2 and ul,o = 1 m/s for saturated water at face tension component.
one atm. In the interfacial lift-off submodel, the vapor momentum is de-
The time step was selected to satisfy the CFL condition which, termined by the balance between buoyancy and surface tension,
combined with the grid size, should result in Co ≤ 0.25. The Fig. 2(a)(ii), i.e., the Weber number and Froude number squared
adopted time step is t = 0.25 μs. (which is the Bond number). This treatment, similar to that used
by Zuber for pool boiling, does not include the liquid inertia, which
4. Wavelength-modulation regime qCHF,h based on surface further prevents the vapor from displacing the liquid. The complete
liquid track derivation is made available in Appendix A as presented in [29,30].
For flow boiling, the contribution of the forced axial flow cannot
In Fig. 3, it was suggested that using the Zuber model, the be neglected, as in Fig. 2(b)(ii) the surface liquid supply depends
wavelength variation of qCHF,h , i.e., its enhancement, is similarly on the axial flow. The forced liquid flow penetrates downstream
achieved in pool and flow boiling. In pool boiling this has been and through the vapor columns, based on the modulated Zuber
achieved by non-uniform (3-D modulated with pitch λc ) porous unit cell, ensuring the irrigation of the heated surface. This mod-
coating and in plain-surface flow boiling by the liquid velocity ul,o . ulation is responsible for the qCHF,h enhancement over the qCHF,Z
This is based on Zuber unit cell and surface vapor and liquid flow and is attributed to ul,o , as shown in Fig. 3.
prescriptions. The support for this assumption/treatment is the ex- The surface vapor columns disturb the surface liquid track and
pected tendency that under ul,o → 0, we should recover pool boil- this disturbance depends on the inertia of the forced liquid flow
ing. However, the slightest ul,o would break the streamwise sym- and the inertia of the escaping vapor columns, adjacent to the
metries and conditions. So, the relationship between λc and ul,o surface. Table 1 lists these liquid and vapor Reynolds number. So,
would address how the anisotropy is created by ul,o and here we when including these inertiae, in particular the effect of liquid
9
viscosity, which are not included in the K-H and R-T stability anal-
ysis, needs to be considered in the surface liquid track fluid me-
chanics. As outlined in Fig. 2(b)(ii) and with the snapshots shown
in Fig. 6, it is the continuous surface liquid track that irrigates the
heated surface, and its interruption leads to the flow-boiling dry-
out.
This effect can be included by multiplying the Bond number by
the Reynolds number of the liquid phase. So here we proposed in-
clusion of the liquid inertia, which is also present in the DNS, and
assume this product is a constant,
Weλ g(ρl − ρg )λ2c ρl ul,o
Reλ = Boλ Reλ = = C2 . (34)
Frλ 2 σ μl
The constant C 2 is found from Table 1, and for high velocities, C ≈
200. This relation between the critical wavelength and the inlet
velocity is of the form λc ≈ u−1
l,o
/3
. Rewriting Eq. (34) and replacing
the wavelength from Eq. (11), we have
1 / 6
π g(ρl − ρg )ρl ul,o
qCHF,h = hlg (σ ρg )1/2 . (35)
8 C 2 σ μl
Using the capillary length lc , the dimensionless qCHF,h becomes
qCHF,h π
1/6
= Re1λ/6 . (36)
ρ g hlg [σ g(ρl − ρg ) ]1/4 8C 1/3
This is comparable to the dimensionless Zuber CHF, Eq. (9), i.e.,
1/4
1/12 π ρ hlg σ ( ρl − ρg ) g
π u1l,/o6,min σ 24 g ρg2
= , (37)
8C 1 / 3
ν 1/6 g ( ρl − ρg ) ρg1/6 hlg [σ g(ρl − ρg )]1/4
l
where C = 200.
Solving for the liquid inlet velocity, the minimum ul,o giving
ν
to qCHF,Z , we find ul,o,min = C 2 729l l = 6 mm/s. Similarly, comparing
c
Eq. (36) to the CHF for flow boiling given by Eq. (11), for λc = lc ,
i.e., qCHF,c we have
1/12
1/12 π h ( σ ρg ) 1 / 2 g ( ρl − ρg )
π u1l,/o6,max σ 8 lg σ
= ,
8C 1/3 νl1/6 g ( ρl − ρg ) ρg1/6 hlg [σ g(ρl − ρg )]1/4
(38)
solving once again for liquid the inlet velocity, corresponding to
ν
qCHF,c is ul,o,max = C 2 l l = 5 m/s.
c
Figure 7 (a) shows the variations of the qCHF,h and critical wave-
length with respect to the inlet velocity for the interfacial lift-off
and the surface liquid-track dryout models. The available experi-
mental data for saturated water flow-boiling on plain surface CHF
and wider channels (Wch > λc ) [48–51] are also shown. Pool boil- Fig. 7. Variations of hydrodynamic CHF with the inlet velocity for saturated water
ing CHF data (corresponding to ul,o = 0 m/s and λc = λc,Z ) [52– (1 atm) in (a) linear and (b) log-log scales. (c) Variations of the critical wavelength
54] are marked on the y axis. The error bars denote the CHF mea- with the inlet velocity. The results are predictions form the interfacial lift-off and
surface liquid track models as well as the experimental and DNS results.
surement uncertainty and the blue band depicts a 25% uncertainty
in the value of the constant C in Eq. (34). The DNS results obtained
according to the iterative process described in Section 3.4 are also 5. Vapor-site geometric-modulation regime qCHF,h
shown. The three data points correspond to the lowest velocities
(ul,o = 0.05, 0.25 and 1 m/s) from Fig. 5(b) to (d), for qCHF,h = In the hexagonal surface vapor-site arrangement, the liquid
1.35, 1.75 and 3 MW/m2, respectively. The u1l,/o6 relation correlates track meanders and this disturbs and can break its streamwise
well the available data for the entire velocity range, presenting a continuous surface coverage. In flow boiling the liquid supply may
smaller maximum absolute error than the interfacial lift-off results, be exclusively from this axial flow, so the alignment of vapor sites
denoted by the red line. Figure 7(b) uses log-log scale to show the may matter.
low velocity results more clearly. Figure 7(c) show the variations Figure 8 (a) renders these vapor escape sites (also referred to
of the critical wavelength with the liquid velocity from the lift-off as surface perforations), following the Zuber unit-cell model, ar-
and surface liquid-track, DNS results. ranged in hexagonal (staggered) and square packing (inline, aligned
Both the high-velocity capillary limit, believed to be the up- with the flow direction). Figure 8(b) shows the time variations
per limit of wavelength modulation enhancement, and the zero- of the surface liquid track width wl for q = 1.75 MW/m2 and
velocity (pool-boiling) limit are also marked. Experiments with ul,o = 0.25 m/s. Dryout is clearly observed for the staggered ar-
flow velocities smaller than 1 cm/s are rather scarce in the liter- rangement. The inline arrangement presents an increase over the
ature, so the minimum ul,o from Eq. (37) cannot be verified. staggered, but less so than the inline rectangular vapor inlets.
10
is highlighted at two elapsed times, 75 and 150 ms. After enough

time, the liquid track in the latter portion of the channel starts to
dry and no wetting is observed from the middle of the channel
onwards. This marks the geometry modulation limitation.
Since the less dense fluid is on top, this condition is stable, but
a critical relative velocity can be found for this flow configuration
1/2 1/2
2 σ g(ρl∗ − ρg∗ )
|ug a − ul a | = , (39)
ρg∗ ρl∗ ρl∗ + ρg∗
where ρ ∗f = ρ f /(ρl + ρg ). From Eq. (39), we can infer the acceler-

ation of vapor will eventually trigger the K-H instability. The axial
vapor wavelength corresponding to the critical relative velocity is
found as
π ( ug − ul )2 ρl∗ ρg∗
λx = . (40)
g (ρl∗ − ρg∗ )
For liquid-vapor saturated water (1 atm) interfaces, λx =
15.7 mm. The relative velocity depends on the heat flux, so a func-
tional relationship of type qCHF,h ≈ λ1x /2 is found as
Fig. 8. (a) Schematic representation of vapor sites and surface liquid track snap- 1/2
shots in hexagonal arrangement, aligned with the inlet flow direction, and with λx (ρl∗ − ρg∗ ) ρg hlg 2
rectangular perforations. Red label indicates dryout. The full video is available in qCHF,KH = g + ul λ. (41)
the supplementary materials. (b) Time variations of liquid track width at different
π ρg ρl ∗ ∗ Aevap c
axial locations for q = 1.75 MW/m2 , ul,o = 0.25 m/s for saturated water (1 atm).
For the given λx = 15.7 mm, the maximum heat flux for marginal
stability is q = 5 MW/m2 . The relationship from Eq. (41) is graphi-
Although the inline circular vapor site arrangement does initially cally shown in Fig. 9(c) for constant λc = 4.5 mm.
establish a clear path for the liquid flow, this path is taken over
by the lateral expansion of vapor, thus nullifying any possible CHF 5.2. Geometric-confined liquid track
enhancements. The rectangular slot geometry, proposed originally
in [26,27,48], with the high aspect ratio L per : Wper helps direct the At high heat fluxes, the rectangular slot shape and the liquid
vapor towards the outlet. inertia might not be enough to prevent the vapor lateral expan-
Neighboring vapor slots do not coalesce immediately close to sion from destabilizing the liquid track. One natural idea is to
the surface (z = 0) due to a lateral vapor expansion prevention, further enhance the geometric confinement effects described in
imposed by the slot geometry, and axial liquid diversion, thus Section 5.1 by further delaying the merging of neighboring vapor
addressing the shortcomings of the square-aligned circular vapor slots.
slots, allowing the formation of a continuous liquid track in the The addition of levees is schematically shown in Fig. 10(a). Two
region between perforations. snapshots at locations from the surface, z = 0.25 and 2 mm, show
the vapor penetration is considerably reduced in the intralevee
5.1. Vapor shearing of liquid track spacings (aqueducts), it was first proposed in [26]. The liquid track
is shielded from the expanding vapor since a significantly smaller
Since the surface liquid track is preserved, no dryout is ob- fraction of the vapor diverts towards the liquid track to disturb it.
served. Further downstream, the surface liquid track continues un- A secondary effect of the addition of levees is the liquid anchor-
derneath the flowing vapor. This flow configuration can be ap- ing, due to the no-slip effect, the levees anchor the liquid, further
proximated by the case of two uniform fluids in relative horizon- stabilizing it and requiring additional shear from the vapor flow
tal motion separated by a horizontal boundary [28], as illustrates to disrupt it. CFD results show the addition of levees is capable of
Fig. 9(a). suppressing the onset of K-H instabilities since the vapor has to
Figure 9 (b) shows snapshots of phase distribution in the flow penetrate further into the intralevee subchannel, pushing the hy-
direction at different axial locations. The evolution of liquid track drodynamic CHF to 15 MW/m2 [26].
Fig. 9. (a) Rectangular vapor slots aligned with the flow direction. (b) Snapshots of the axial distribution of the liquid track at two distinct elapsed times (q = 5 MW/m2 ,
ul,o = 0.5 m/s). (c) Variations of interfacial vapor wavelength with Kelvin-Helmholtz CHF for λc = 4.5 mm. The axial K-H marginal stability is highlighted.
11
No dryout is observed even though the heat flux was increased

to 15 MW/m2 at the same inlet velocity of ul,o = 0.5 m/s. The
higher intensity turbulence is still stemming from the interface,
but further away from the liquid track, at the top portion of the
channel, where the high-velocity vapor deflects the liquid. As men-
tioned earlier, inside the intralevee subchannel, the vapor is at a
much lower velocity, so the flow is less turbulent, and the K-H sta-
bility is preserved.
6. Porous metasurface for geometric modulation:

Capillary-viscous bottleneck
Metasurfaces are unit-cell based, generally 3-D surface struc-

tures designed and fabricated for a targeted function. For example,
in electromagnetic and sound wave applications they have sub-
wavelength features capable of modulating the surface transmis-
sion, absorption and refection from phenomenon of interest [55–
57].
As discussed above, in boiling crises, the surface vapor gen-
eration triggers interfacial instabilities, and related metasurfaces,
including sub-critical wavelength surface modulations have been
used [58]. The evolution of boiling metasurfaces for enhancing the
CHF is shown in Fig. 11, employing capillary evaporator wicks. It
Fig. 10. (a) Snapshot of C-GMR video and top-view of phase distribution in the began with a uniform porous coating to reduce the surface super-
axial planes at elevations z = 0.25 and 2 mm for metasurface with levees. The full heat [59] and later to also increase the CHF over the plain sur-
video is available in the supplementary materials. (b) Time variation of liquid phase
face [60]. Figure 11 shows the schematic of metasurfaces and their
yz area fraction at x = 16 mm. q = 5, 10 and 15 MW/m2 , ul,o = 0.5 m/s for saturated
water (1 atm).
physical images, as well as the CHF enhancement compared to the
Zuber qCHF,Z . The shaded areas at the bottom shows the measured
enhancement achieved so far in the wavelength modulation, geo-
Figure 10 (b) shows the time variation of the dimensionless metric modulation, and geometric-confined modulation regimes.
area of liquid in the intralevee channel (spanning -1.25 mm≤ y ≤ The first porous metasurface sub-wavelength modulation of the
1.25 mm by 0 mm≤ z ≤ 3 mm) at axial location x = 16 mm for interfacial pool boiling wave was by Liter and Kaviany [33], reduc-
different heat fluxes (q = 5, 10 and 15 MW/m2 ). The dimension- ing the superheat and enhancing the CHF using sintered copper
less area is calculated according to particles. The 3-D capillary structures for phase separation using
x columnar and lateral arteries have been designed for irrigation and
1 1 wl Hl
f = 0, α ≥ 0.5 evaporative cooling of concentrated heat source [61–64]. In flow
A∗l = f dzdy, . (42)
x wl Hl 0 0 f = 1, α < 0.5 boiling, the porous metasurface provides separation for the liquid
i
supply and vapor escape paths, preventing the phase competition
Fig. 11. An evolution of porous metasurfaces for enhanced CHF [59], updated with the flow-boiling canopy wick metasurface [26,48]. The porous coatings are made of
sintered copper particles. The corresponding qCHF in the three enhancement regimes, wavelength-modulation, geometric modulation, and confined-geometric modulation,
are also shown.
12
Fig. 13. Variations of the capillary-viscous CHF with the total pressure drop, for
four different evaporator wick designs (Table 3) which affect the maximum capillary
pressure. Results for saturated water at one atm.
The maximum total pressure drop allowed corresponds to the

capillary-viscous heat flux limit (CHF controlled by the evapora-
tor maximum capillary pressure). Equation (43) can be rewritten
to relate the heat flux and the maximum capillary pressure, i.e.,

μl Lp2 (Lp − D p )/3 4 Hp
+ q
ρl hlg Ke δl (Lp + D p )/2 Kp D2p CHF,c-v
Fig. 12. Schematic illustration of the flow-boiling canopy wick. Adapted from Kim
et al. [48]. ( A∗ )2 C 2
+ q = pc,max , (44)
2ρg2 hlg CHF,c-v
close to the surface [25]. The wick superheat, capillary-vicous and (Wper +D p )(λc +D p )
where A∗ = and Lp =
(Wper +D p )(λc +D p −6π D2p /4 )
hydrodynamic qCHF are shown. !4 "1/2
Also shown in Fig. 11, is the recent work on raising qCHF,h π (Wper + D p )( p + D p ) . The derivation of this quadratic
by preventing the onset of Kelvin-Helmholtz instabilities with the equation is discussed in detail in [27].
flow-boiling canopy wick (FBCW), a 3-D (and multiscale) porous Figure 13 shows the variations of the capillary-viscous CHF with
metasurface that modulates vapor-escape sites (vapor perforations) the overall pressure drop for different wick designs. Results for
[25,26], Fig. 12. The FBCW separate of liquid and vapor phases, di- the blue curve have been experimentally verified in [48] for a bi-
rects the liquid from the liquid channel to the thin evaporator wick modal bilayer wick with particles sized de = 50 and 100 μm with a
with high capillary pressure capability, and allows for vapor space pc,max = 2.32 kPa and the current capillary-viscous CHF is qCHF,c-v =
over the evaporator for achieving large thermal conductance and 4.6 MW/m2 . The orange and green lines are design variations cor-
vapor escape sites (perforations) directing the vapor into the liq- responding to de = 30 and 60μm and de = 100 and 200 μm. The
uid channel. Its three main components are the evaporator, posts, smaller particles wick ( pc,max = 3.86 kPa) results in a CHF reduc-
and perforated canopy. All three porous components coordinately tion to qCHF,c-v = 2.1 MW/m2 while the large particles wick pro-
direct the liquid towards the heated surface. Liquid is assumed to vide an enhancement despite the pc,max reduction, qCHF,c-v = 5.5
be saturated in the channel and canopy and post wicks. The thin MW/m2 ( pc,max = 1.15 kPa). A large enhancement is obtained us-
wick evaporator is responsible for spreading and evaporating the ing a thicker evaporator layer (3 layers): qCHF,c-v = 6.3 MW/m2 . So,
liquid with low resistance (thermal and hydraulic), generating sat- search for raising qCHF,c-v , which is the bottleneck to higher qCHF,h ,
urated vapor. The posts separate the evaporator and the canopy continues. The corresponding pressure drops are listed in Table 3.
and form the vapor space. The perforated canopy is where the va- Further details are available in [48].
por wavelength modulation occurs, since it separates of the phases Figure 14 (a) and (b) show the variations of the hydrodynamic
and controls of vapor sites entering the liquid channel at regularly CHF with the critical wavelength λc and ul,o . The three regimes,
spaced perforations. The addition of levees further stabilizes the namely wavelength, geometry, and geometry-confined modulation
liquid track formed on the perforated canopy. The high velocity of with levees, are marked, as well as their transitions (boundaries).
the liquid being wicked across the canopy and posts leads to high The wavelength-modulation CHF is given by Eqs. (38) and
Péclet number, resulting in negligible upstream thermal conduc- (39) and is primarily controlled by ul,o and presented by the
tion [33]. At the thin wick evaporator, the sintered particles con- dashed gray line. The inlet velocity reduces the critical wavelength
tribute to an increased effective thermal conductivity, reducing the causing higher vapor velocity and increasing the heat flux required
thermal resistance [27]. to trigger the instability. This enhancement is limited by the small-
The liquid flow path from the channel to the heated surface est physically-limited wavelength, i.e., the capillary length, marked
goes through the canopy, posts, and evaporator layer, and each in Fig. 14(a). The 1/6 power-law dependency of the CHF on the in-
component has its associated pressure drop. Furthermore, the va- let velocity is shown in Fig. 14(b) and is in general agreement with
por escape through the perforations causes an additional pressure the low and high-velocity experimental data.
drop that must also be accounted for by the capillary wick. The The geometry modulation, decoupling the wavelength from the
combined pressure drops cannot exceed the maximum capillary inlet velocity, allows for a larger enhancement, as shown by the
pressure pc,max , i.e., strip regime in Fig. 14(a). The CFD predictions [26] and experi-
mental results [48–51] are also shown. This enhancement is con-
pi = pe + p p + pca + p per,g ≤ pc,max , (43) trolled by the vapor confinement provided by the rectangular per-
13
Table 3
Pressure drop components for the different evaporator layer designs. The corresponding qCHF,c-v are also listed.
#
d e , μm qCHF,c-v , MW/m2 G/A, MW/m2 -K pe , kPa p p , kPa pca , kPa p per,g , kPa pi , kPa pc,max , kPa
50,100 5.0 0.25 1.93 0.15 0.04 0.20 2.32 2.32

3060 2.1 0.42 3.67 0.08 0.03 0.08 3.86 3.86
100,200 5.5 0.15 0.58 0.22 0.05 0.28 1.14 1.14
50(3),100 6.3 0.08 1.63 0.29 0.06 0.34 2.32 2.32
Fig. 14. Regime diagram showing the variations of flow-boiling CHF with (a) wavelength and (b) inlet velocity. The three modulation regimes are highlighted: wavelength
modulation, geometric modulation, and geometric-confined modulation, as well as their boundaries. Experimental data and CFD data are also shown. Orange circles denote
experimental data for geometric modulation including length scale, inlet velocity and global thermal conductance.
forations, restricting the vapor lateral expansion and preserving boiling, hydrodynamic instability unit-cell theoretical construct, al-
the liquid tracks. The orange region in Fig. 14(b) is bound by the lows for circular vapor-column site wavelength modulation by the
qCHF − ul,o relation Eq. (34) at the bottom and the K-H instabil- liquid-flow velocity ul,o . In the DNS of flow boiling, the vapor and
ity Eq. (38) on the top. The capillary-viscous metasurface-imposed liquid velocities on the heated surface are prescribed, according to
limit is not shown since it depends on the characteristic wave- the ul,o = 0, i.e., the Zuber unit-cell construct. The results for ul,o >
length of the metasurface, which is decoupled from the inlet ve- 0 show a leading-edge, surface liquid track formed on the heated
locity. surface. This track meanders between the vapor columns and once
The current results show the limits for the K-H dominated the CHF is reached, it becomes unstable and dryout occurs within
geometry modulation qCHF,h and the porous metasurface qCHF,c-v the first (leading edge) unit cell.
are rather close, as confirmed by experiments [65], so no fur- The DNS results show with ul,o as small as 5 cm/s, the inertia is
ther enhancements can be obtained before further improvement to large enough to deflect the Zuber vertical vapor columns and no-
the porous metasurface, i.e., increasing qCHF,c-v . After that, to con- ticeably increase the critical heat flux (CHF) over the Zuber limit
tinue increasing the qCHF,h , the use of levees is required, i.e., the qCHF,Z . This is called the wavelength-modulation CHF enhance-
geometric-confined modulation regime, represented by the blue ment regime and ranges from this very-low ul,o limit qCHF,Z =
color in Fig. 14(a), creating the intralevee spacing which effectively 1.1 MW/m2 to the capillarity limit of about qCHF,c = 3.3 MW/m2 ,
raised qCHF,h > 10 MW/m2 . The wick superheat limit qCHF,sh which where the wavelength reaches the capillarity length based on the
marks vapor forming within the evaporation wick is reached there, critical Bond number. The predicted range is in good agreement
estimated as 13.3 MW/m2 for the current FBCW [48]. The mea- with experiments and a correlation is proposed based on the prod-
sured FBCW thermal conductance G/A which is rather very high is uct of Bond and liquid Reynolds number relating the wavelength
shown in Fig. 14 and listed in Table 3. and ul,o . In this new relation qCHF,h is proportional to u1l,/o6 .
The DNS results also show that to reach beyond the
7. Conclusions wavelength-modulation CHF enhancement regime, the vapor site
geometry should be anisotropic, with rectangular vapor sites
Direct numerical simulation (DNS), including turbulence, of sat- aligned with the flow direction to allow stable leading-edge liq-
urated water (one atm) flow boiling, based on the Zuber pool- uid track. This is called the geometric-modulation regime and upon
14
optimization of the site geometry can reach the next limit hydro- The velocities and pressure are decomposed into base flow and

dynamic limit qCHF,h and is over 10 MW/m2 . This is reached due perturbed components u f = ū f + uf , v f = v̄ f + vf , p = p̄ + p . Sub-
to vapor shearing (Kelvin-Helmholtz instability) of the liquid track stituting the decomposed flow quantities into Eqs. (A.1)–(A.3) and
between the vapor sites. To prevent this the liquid tracks are ge- performing the required simplifications, we get to
ometrically confined by levees placed around the vapor sites. This
geometric-confined CHF enhancement limit is rather very large. ∂ uf ∂vf
+ = 0, (A.4)
The geometric modulation, including liquid track confinement, ∂x ∂z
are achieved by porous metasurfaces, which are 3-D porous and
perforated wicks. One example is the flow-boiling canopy wick
∂ uf ∂ uf ∂ p
(FBCW) which delivers liquid from the above-mentioned liquid ρ + ū f =− , (A.5)
tracks formed on a perforated, porous canopy to porous posts and ∂t ∂x ∂x
to a thin evaporator wick. The vapor formed over the evaporator

flows between posts and escapes through the perforated canopy ∂vf ∂vf ∂ p
and mixes with the flowing liquid. While porous metasurfaces al- ρ + ū f =− . (A.6)
low for realization of the prescribed the surface liquid and vapor
∂t ∂x ∂z
velocities in the DNS, they introduce the internal hydrodynamic
Differentiating Eqs. (A.5) and (A.6) with respect to x and z, respec-
limit within the porous structure. This is the capillary-viscous limit
tively, and adding them together yields the Laplace equation
qCHF,c−v is controlled by the maximum capillary pressure and cur-
rently limits the saturated (one atm) water to about 5.2 MW/m2 . ∂ 2 p ∂ 2 p
The summary of the predicted three CHF-enhancement regimes, + = 0, (A.7)
∂ x2 ∂ z2
i.e., the wavelength, geometric, and geometric-confined modula-
tions, and their boundaries, along with the existing experimental and the base pressure field for each phase is simply hydrostatic
results, are shown in Fig. 14.
p̄i = po − ρi gz. (A.8)
Author Statement The initial interfacial perturbation is considered as a Fourier com-

ponent of the interfacial wavelength
Julio Ferreira: modeling, formal analysis, visualization, writing
- original draft; Massoud Kaviany: conceptualization, theory, su- δ (x, t = 0 ) = δoeiκ x , (A.9)
pervision, writing - review & editing administration. where δ is the perturbation and κ = 2π /λ the wave number. Sim-
ilarly,
Declaration of Competing Interest
δ = δoeiκ x+ωt ,
The authors declare that they have no known competing finan- vf = vˆ f (z )eiκ x+ωt ,
cial interests or personal relationships that could have appeared to
influence the work reported in this paper. p = pˆ (z )eiκ x+ωt , (A.10)
where ω = 2π / f is the angular frequency. Using this pressure in

Acknowledgements Eq. (A.7), we have
d2 pˆ
This work is supported by NSF USA (Thermal Transport and Pro- = κ 2 pˆ. (A.11)
cesses, CBET16235720). JF is thankful, acknowledging his scholar- dz2
ship supported by CAPES-Brazilian Federal Agency for Support and The solution is
Evaluation of Graduate Education, within the Ministry of Education
(88881.170629/2018-01). pˆ i = pi,oeκ z , (A.12)
Substituting these in Eq. (A.6), we obtain

Appendix A. Boiling hydrodynamic instability theory

! "−1 d pˆ
The results in Section 2 are based on the one treatments in vˆ f (z ) = − ρ ω + iκ ū f . (A.13)
dz
[29] and [30]. An initially stratified two-phase flow may experience
interfacial instabilities under relative phasic motion [28], depend- Evaluating the pressure gradient from Eq. (A.12), the velocities are
ing on the flow condition, this impacts the heat and mass transfer found and when using the functional form for the decomposed
[29]. The linear stability theory considers the liquid-vapor interface perturbed flow results in a differential equation for the perturbed
as planar and horizontal with the phases moving with axial veloc- axial velocity. This expression is integrated using the condition
ities ug and ul . The disturbance grows axially and with time, and uf → 0 far from the interface and gives
incompressible, inviscid, 2-D base flow is considered. The Navier-

Stokes equations are [29] i dvˆ f
uf = eiκ x+ωt . (A.14)
∂ u f ∂v f κ dz
+ = 0, (A.1)
∂x ∂z With this solution for the perturbed flow, the boundary conditions
are applied at the interface starting with
∂uf ∂uf ∂uf ∂p
ρ + uf + vf =− , (A.2) ∂δ ∂δ
∂t ∂x ∂z ∂x vf = + ū f , (A.15)
∂t ∂x

∂v f ∂v f ∂v f ∂p and imposing this condition and substituting Eqs. (A.10) and
ρ + uf + vf =− − ρ g. (A.3)
∂t ∂x ∂z ∂z (A.13) into Eq. (A.15) gives
15
Fig. A1. Evolution of interfacial instability on a horizontal liquid-vapor interface

and the Kelvin-Helmholtz instability and wavelength.
κ pg,o[ρg (ω + iκ ūg )]−1 = δoω + iδoκ ūg ,

− κ pl,o[ρl (ω + iκ ūl )]−1 = δoω + iδoκ ūl . (A.16)
Solving these equations for p f,o gives
ρg Fig. A2. Evolution of the interfacial instability on horizontal liquid-vapor interface
pg,o = (ω + iκ ūg )2 δo,
κ and the formation of the Rayleigh-Taylor instability and wavelength.
ρ
pl,o = − l (ω + iκ ūl )2 δo. (A.17)
κ with ρl + ρg ≈ ρl , we have
The momentum balance normal to the interface gives [30] 1 / 2
1 (ρl − ρg )g2π /λ σ (2π /λ )3
1 ∂vl ∂vg ω=± − . (A.24)
pl − pg = σ + + 2 μl − . (A.18) ρl ρl
R1 R2 ∂z ∂z
The maximum is found and the wavelength corresponding to
In the absence of phase change, the second term on the RHS is this maximum is the most dangerous wavelength λRT,d , given by
null and it reduces to the Young-Laplace equation, where 1/R1 is Eq. (2).
the interfacial Euler curvature The above analysis and results, although simplified, are experi-
d2 δ mentally verified. A similar decomposition including the thermo-
1 dx2
=
R1 dδ 2 3/2 , (A.19) mechanical wave propagation is presented in [30]. In the flow-
1+ dx
boiling DNS, we use 3-D fluid flow to track the phase and define
surface liquid track dryout.
and here (2-D), the other principal curvature is 1/R2 = 0. Substi-
tuting the curvature into Eq. (A.18), recalling that p f = p̄ f + pf , we Appendix B. LES turbulent treatment
obtain
! " In ANSYS Fluent [38], the finite-volume discretization provides
pl,o − pg,o = δo (ρl − ρg )g + σ κ 2 . (A.20)
the filtering operation for the LES as
Substituting Eq. (A.17) into Eq. (A.20) the angular frequency is 1
$ ! %1 / 2 " φ˜ (x ) = φ (x )dx , (B.1)
v
κ 2 ρl ρg (ūl − ūg ) − σ κ 3 + (ρl − ρg )gκ (ρl + ρg )
2 V
ω=± where V is the computational cell volume. The filter function is
ρl + ρg
then
iκ (ρl ūl + ρg ūg )
− . (A.21) 1
, x ∈ v
ρl + ρg G ( x, x ) = V . (B.2)
0, otherwise
Both the surface tension and gravitational acceleration are stabi-
lizing forces. The instability will occur when Eq. (A.21) has a real The filtered Navier-Stokes equations are given in Eqs. (24) and (25).
positive solution, i.e., The viscosity stress and the subgrid-scale stress tensors are given
⎧ ⎫1/2 as
⎨ 2πλσ + (ρl −ρg )gλ
( ρl + ρg ) ⎬
2π ∂ u˜ f,i ∂ u˜ f, j 2 ∂ u˜ f,l
ūl − ūg > . (A.22) σi j = μ f + − μf δ , (B.3)
⎩ ρl ρg ⎭ ∂xj ∂ xi 3 ∂ xl i j
A critical wavelength can be derived from Eq. (A.22) by setting it τi j = ρ f u f,i¯u f, j − ρ f ū f,i ū f, j . (B.4)
to zero. This type of instability is referred to as Kelvin-Helmholtz
To generate a time-independent inlet condition, a random 2-D vor-
instability, shown in Fig. A.1. Its wavelength is given in Eq. (1).
tex is considered. With this approach, a perturbation is added on a
When the velocities ul = ug = 0, with the change in sign of
specified mean velocity profile via a fluctuating vorticity field (i.e.,
gravity representing the dense fluid on top, we have the Rayleigh-
two-dimensional in the plane normal to the streamwise direction).
Taylor instability, as illustrates Fig. A.2. The critical wavelength for
The subgrid-scale resulting from the filtering operation requires
the Rayleigh-Taylor instability is identical to Eq. (2), but the angu-
modeling. As in RANS, the Boussinesq hypothesis is employed and
lar frequency is found by using these velocities and gravity term,
the stresses are calculated from Eq. (26).
and Eq. (A.22) becomes
The strain tensor is defined as
1 / 2
(ρl − ρg )g2π /λ − σ (2π /λ )3 1 ∂ ū f,i ∂ ū f, j
ω=± , (A.23) S̄i j = + . (B.5)
ρl + ρg 2 ∂xj ∂ xi
16
The interfacial force term in the momentum conservation equa-

tion is given by the continuum surface force (CSF) model as
ρ f κg ∇α
fs = σ , (B.9)
1
2 ( ρl + ρg )
where κg = ∇ · (∇α /α ) is the interface normal defined as the gra-
dient of the volume fraction.
In order to compute turbulent quantities, a decomposition is re-
quired and the velocity resolved by the grid can be decomposed
into mean and fluctuation quantities [42]
ū f = u˜ f − uf , (B.10)
where u˜ f is the filtered velocity field, ū f is the mean velocity field
averaged over a period of time t and uf is the fluctuation ve-
locity field. The turbulence kinetic energy resolved by the grid is
calculated as
Fig. B1. Turbulence kinetic energy spectra in the axial midplane for (a) plain sur-
1 ¯2
face q = 1.75 MW/m2 and ul,o = 0.5 m/s. The Kolmogorov inertial subrange -5/3 Ē f,t = u + v¯f2 + w¯f2 . (B.11)
[42] and the empirical dissipation range -3 [45,67] scales are shown. The subgrid 2 f
turbulence kinetic energy is highlighted in gray shade.
The unresolved portion of the turbulent kinetic energy is the sub-
grid component, in the WALE model [36], this portion of the tur-
In the Wall-Adapting Local Eddy-Viscosity (WALE) [66] the eddy bulent kinetic energy is calculated as
viscosity is modeled by 2
μ f,t
3/2 Ē f,t,sgs =
Cw ρ f Ls
. (B.12)
Sidj Sidj
μ f,t = ρ 2
f Ls 5/2 5/4 , (B.6)
Si j Si j + Sidj Sidj The total kinetic energy is, then
Ē f,t,total = Ē f,t + Ē f,t,sgs . (B.13)
where Ls is the mixing length for the subgrid scale
Ē f,t
Ls = min κ z, CwV 1/3 , (B.7) When ≥ 0.8, the mesh is well suited. For Fig. B.1 and similar
Ē f,t,total
results, the ratio is 0.96.
where κ is the von Karman constant, Cw = 0.325, V 1/3 is the local Figure B.1 shows the dimensionless turbulence spectra for
grid scale and q =1.75 MW/m2 on plain surface. The Kolmogorov -5/3 power is
2 2 2 shown as a blue dashed line in the inertial subrange. The subgrid
1 ∂ ū f,i ∂ ū f, j 1 ∂ ū f,k
Sidj = + − δi j . (B.8) turbulence kinetic energy decay is highlighted in gray. In the sub-
2 ∂xj ∂ xi 3 ∂ xk grid modeling range, for wavenumbers larger than the filter length,
Fig. B2. The distributions of the Reynolds stress tensor components in the axial midplane for q = 1.75 MW/m2 and ul,o = 0.5 m/s for saturated water (1 atm). The interface
is shown with the curved solid line.
17
i.e., in the dissipation subrange, the spectrum follows the -3 power [14] S. Gong, L. Zhang, P. Cheng, E.N. Wang, Understanding triggering mechanisms
of the dissipation scale observed experimentally for bubbly flows for critical heat flux in pool boiling based on direct numerical simulations, Int.
J. Heat Mass Transf. 163 (2020) 120546, doi:10.1016/j.ijheatmasstransfer.2020.
[45,67]. This is reasonable since the overall observation is that tur- 120546.
bulence originates from eddies formed at the liquid-vapor interface [15] C. Konishi, I. Mudawar, Review of flow boiling and critical heat flux in
and the smallest eddies would be encountered in the small vapor microgravity, Int. J. Heat Mass Transf. 80 (2015) 469–493, doi:10.1016/j.
ijheatmasstransfer.2014.09.017.
inclusions that detach from the escaping vapor stream. The highest [16] M. Bruder, G. Bloch, T. Sattelmayer, Critical heat flux in flow boiling-review of
wavenumber which can be resolved by the simulation corresponds the current understanding and experimental approaches, Heat Transf. Eng. 38
to twice the grid size and is denoted as (3) (2017) 347–360, doi:10.1080/01457632.2016.1189274.
[17] J. Galloway, I. Mudawar, CHF mechanism in flow boiling from a short heated
2π wall-i. examination of near-wall conditions with the aid of photomicrography
κ f,x = . (B.14) and high-speed video imaging, Int. J. Heat Mass Transf. 36 (10) (1993) 2511–
2Ls 2526, doi:10.1016/S0017-9310(05)80190-5.
The maximum wavenumber is related to the Kolmogorov length [18] J. Sturgis, I. Mudawar, Critical heat flux in a long, rectangular channel sub-
jected to one-sided heating-II. Analysis of critical heat flux data, Int. J. Heat
scale, given as [42] Mass Transf. 42 (10) (1999) 1849–1862, doi:10.1016/S0017- 9310(98)00275- 0.
1/4 [19] C.-N. Huang, C.R. Kharangate, A new mechanistic model for predicting flow
νg3 boiling critical heat flux based on hydrodynamic instabilities, Int. J. Heat Mass
λ f,K = . (B.15)
f,t Transf. 138 (2019) 1295–1309, doi:10.1016/j.ijheatmasstransfer.2019.04.103.
[20] K. Ling, S. Zhang, W. Liu, X. Sui, W. Tao, Interface tracking simulation for sub-
cooled flow boiling using VOSET method, Front. Energy Res. 8 (2021) 325,
These are also marked in Fig. B.1. These and the power-law com- doi:10.3389/fenrg.2020.526035.
parisons in Fig. B.1 support that the numerically predicted turbu- [21] R. Revellin, J.R. Thome, A theoretical model for the prediction of the critical
lent quantities have the expected theoretical behavior. heat flux in heated microchannels, Int. J. Heat Mass Transf. 51 (5) (2008) 1216–
1225, doi:10.1016/j.ijheatmasstransfer.20 07.03.0 02.
A similar trend is encountered with the components of the
[22] T.J. Zhang, S. Chen, E.N. Wang, A separated-flow model for predicting flow boil-
Reynolds stress tensor, as shows Fig. B.2. Results were averaged ing critical heat flux and pressure drop characteristics in microchannels, in:
over a similar time frame. Turbulence intensity is expected to be ASME 2012 10th International Conference on Nanochannels, Microchannels,
and Minichannels, in: International Conference on Nanochannels, Microchan-
more pronounced near the channel walls and the liquid-vapor in-
nels, and Minichannels, 2012, pp. 39–48, doi:10.1115/ICNMM2012-73046.
terface. Figure B.2 shows that the most intense stresses occur in [23] J. Broughton, Y.K. Joshi, Flow boiling in geometrically modified microchannels,
the vapor phase, indicating the production of turbulence is more Phys. Fluids 33 (10) (2021) 103308, doi:10.1063/5.0062585.
pronounced in the gas phase due to the large inertial force and [24] G. Liang, I. Mudawar, Review of channel flow boiling enhancement by surface
modification, and instability suppression schemes, Int. J. Heat Mass Transf. 146
interfacial shearing. (2020) 118864, doi:10.1016/j.ijheatmasstransfer.2019.118864.
[25] M. Kim, M. Kaviany, Flow-boiling canopy wick for extreme heat transfer, Int. J.
Supplementary material Heat Mass Transf. 117 (2018) 1158–1168, doi:10.1016/j.ijheatmasstransfer.2017.
10.079.
[26] J. Ferreira, M. Kaviany, Geometric-confinement suppression of flow-boiling in-
Supplementary material associated with this article can be stability using perforated wick: Part I CHF and conductance enhancement, Int.
found, in the online version, at 10.1016/j.ijheatmasstransfer.2022. J. Heat Mass Transf. 159 (2020) 120080, doi:10.1016/j.ijheatmasstransfer.2020.
120080.
123051. [27] J. Ferreira, M. Kaviany, Geometric-confinement suppression of flow-boiling in-
stability using perforated wick: Part II CHF limits and wick properties, Int.
References J. Heat Mass Transf. 159 (2020) 120079, doi:10.1016/j.ijheatmasstransfer.2020.
120079.
[1] D.E. Kim, D.I. Yu, D.W. Jerng, M.H. Kim, H.S. Ahn, Review of boiling heat trans- [28] S. Chandrasekhar, The stability of superposed fluids: the Kelvin-Helmholtz in-
fer enhancement on micro/nanostructured surfaces, Exp. Therm. Fluid Sci. 66 stability, in: Hydrodynamic and Hydromagnetic Stability, in: International Se-
(2015) 173–196, doi:10.1016/j.expthermflusci.2015.03.023. ries of Monographs on Physics, Clarendon Press, 1961, pp. 481–514.
[2] M. Dadhich, O.S. Prajapati, V. Sharma, A systematic review on the heat transfer [29] V.P. Carey, Convective Boiling in Tubes and Channels, in: Liquid-Vapor
investigation of the flow boiling process, Chem. Eng. Process. Process Intensif. Phase-Change Phenomena, Taylor & Francis, 1992, pp. 483–562.
165 (2021) 108425, doi:10.1016/j.cep.2021.108425. [30] M. Kaviany, Carrier Energy Transport and Transformation Theories, in: Heat
[3] K. Leong, J. Ho, K. Wong, A critical review of pool and flow boiling heat trans- Transfer Physics, 2nd ed., Cambridge University Press, 2014, pp. 119–172. https:
fer of dielectric fluids on enhanced surfaces, Appl. Therm. Eng. 112 (2017) 999– //doi.org/10.1017/CBO9781107300828.
1019, doi:10.1016/j.applthermaleng.2016.10.138. [31] H. Kim, J.C. Padrino, D.D. Joseph, Viscous effects on Kelvin-Helmholtz instabil-
[4] R. Balasubramaniam, E. Rame, J. Kizito, M. Kassemi, Two Phase Flow Modeling: ity in a channel, J. Fluid Mech. 680 (2011) 398–416, doi:10.1017/jfm.2011.206.
Summary of Flow Regimes and Pressure Drop Correlations in Reduced and Par- [32] V. Sernas, J. Lienhard, V. Dhir, The Taylor wave configuration during boiling
tial Gravity, Technical Report, National Center for Space Exploration Research, from a flat plate, Int. J. Heat Mass Transf. 16 (9) (1973) 1820–1821, doi:10.1016/
2006 https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20 060 0 08906.pdf. 0017-9310(73)90175-0.
[5] C.R. Kharangate, I. Mudawar, Review of computational studies on boiling and [33] S.G. Liter, M. Kaviany, Pool-boiling CHF enhancement by modulated porous-
condensation, Int. J. Heat Mass Transf. 108 (2017) 1164–1196, doi:10.1016/j. layer coating: theory and experiment, Int. J. Heat Mass Transf. 44 (22) (2001)
ijheatmasstransfer.2016.12.065. 4287–4311, doi:10.1016/S0017-9310(01)0 0 084-9.
[6] N.I. Kolev, Multiphase Flow Dynamics 4, Springer, 2012 https://link.springer. [34] A. Bar-Cohen, A. McNeil, Parametric effects on pool boiling critical heat flux in
com/book/10.1007/978- 3- 642- 20749- 5. dielectric liquids, in: Pool and External Flow Boiling, 1992, pp. 171–175.
[7] B. Magolan, E. Baglietto, C. Brown, I.A. Bolotnov, G. Tryggvason, J. Lu, Multi- [35] S.S. Kutateladze, A.I. Leont’ev, Some applications of the asymptotic theory of
phase turbulence mechanisms identification from consistent analysis of direct the turbulent boundary layer, in: International Heat Transfer Conference, vol. 3,
numerical simulation data, Nucl. Eng. Technol. 49 (6) (2017) 1318–1325, doi:10. 1966, pp. 1–6.
1016/j.net.2017.08.001. Special Issue on International Conference on Mathe- [36] J. Weisman, B. Pei, Prediction of critical heat flux in flow boiling at low
matics and Computational Methods Applied to Nuclear Science and Engineer- qualities, Int. J. Heat Mass Transf. 26 (10) (1983) 1463–1477, doi:10.1016/
ing 2017 (M&C 2017) S0 017-9310(83)80 047-7.
[8] T.Q.D. Pham, J. Jeon, D. Jo, S. Choi, Two-phase flow simulations using 1D [37] C. Lee, I. Mudawwar, A mechanistic critical heat flux model for subcooled flow
centerline-based C- and U-shaped pipe meshes, Appl. Sci. 11 (5) (2021), doi:10. boiling based on local bulk flow conditions, Int. J. Multiphase Flow 14 (6)
3390/app11052020. (1988) 711–728, doi:10.1016/0301- 9322(88)90070- 5.
[9] N. Zuber, Hydrodynamic aspects of boiling heat transfer, University of Califor- [38] Fluent Reference Manual, ANSYS, 2009.
nia, Los Angeles, 1959 Ph.D. thesishttps://www.osti.gov/biblio/4175511. [39] G. Tomar, G. Biswas, A. Sharma, A. Agrawal, Numerical simulation of bubble
[10] J.H. Lienhard, V.K. Dhir, Hydrodynamic prediction of peak pool-boiling heat growth in film boiling using a coupled level-set and volume-of-fluid method,
fluxes from finite bodies, J. Heat Transf. 95 (2) (1973) 152–158, doi:10.1115/ Phys. Fluids 17 (11) (2005) 112103, doi:10.1063/1.2136357.
1.3450013. [40] R. Courant, F. K, H. Lewy, On the partial difference equations of mathematical
[11] P.J. Berenson, Film-Boiling heat transfer from a horizontal surface, J. Heat physics, Math. Ann. 100 (1) (1928) 32–74, doi:10.1063/1.2136357.
Transf. 83 (3) (1961) 351–356, doi:10.1115/1.3682280. [41] P. Sagaut, Application to Navier—Stokes Equations, in: Large Eddy Simulation
[12] L. Zhang, S. Gong, Z. Lu, P. Cheng, E.N. Wang, Boiling crisis due to bub- for Incompressible Flows: An Introduction, Springer Berlin Heidelberg, Berlin,
ble interactions, Int. J. Heat Mass Transf. 182 (2022) 121904, doi:10.1016/j. Heidelberg, 2006, pp. 45–81. https://doi.org/10.1007/3- 540- 26403- 5_3.
ijheatmasstransfer.2021.121904. [42] S.B. Pope, Turbulent Flows, Cambridge University Press, 20 0 0, doi:10.1017/
[13] Q.æ. Li, Y.ä. Yu, Z.X.æ. Wen, How does boiling occur in lattice Boltzmann sim- CBO9780511840531.
ulations? Phys. Fluids 32 (9) (2020) 093306, doi:10.1063/5.0015491.
18
[43] J. Jeong, F. Hussain, On the identification of a vortex, J. Fluid Mech. 285 (1995) [55] S.S. Bukhari, J.Y. Vardaxoglou, W. Whittow, A metasurfaces review: definitions
69–94, doi:10.1017/S0 0221120950 0 0462. and applications, Appl. Sci. 9 (13) (2019), doi:10.3390/app9132727.
[44] J.-m. Zhan, Y.-t. Li, W.-h.O. Wai, W.-q. Hu, Comparison between the Q criterion [56] T. Gric, O. Hess, Chapter 5 - metasurfaces, in: T. Gric, O. Hess (Eds.), Phe-
and Rortex in the application of an in-stream structure, Phys. Fluids 31 (12) nomena of Optical Metamaterials, Micro and Nano Technologies, Elsevier, 2019,
(2019) 121701, doi:10.1063/1.5124245. pp. 131–154, doi:10.1016/B978- 0- 12- 813896- 0.0 0 0 05-5.
[45] S. Long, J. Yang, X. Huang, G. Li, W. Shi, M. Sommerfeld, X. Yang, Large- [57] A. Lininger, A.Y. Zhu, J.-S. Park, G. Palermo, S. Chatterjee, J. Boyd, F. Ca-
eddy simulation of gas-liquid two-phase flow in a bubble column reactor passo, G. Strangi, Optical properties of metasurfaces infiltrated with liquid
using a modified sub-grid scale model with the consideration of bubble- crystals, Proc. Natl. Acad. Sci. 117 (34) (2020) 20390–20396, doi:10.1073/pnas.
eddy interaction, Int. J. Heat Mass Transf. 161 (2020) 120240, doi:10.1016/j. 2006336117.
ijheatmasstransfer.2020.120240. [58] M. Yang, P. Sheng, Sound absorption structures: from porous media to acous-
[46] M. Saeedipour, S. Puttinger, N. Doppelhammer, S. Pirker, Investigation on tur- tic metamaterials, Annu. Rev. Mater. Res. 47 (1) (2017) 83–114, doi:10.1146/
bulence in the vicinity of liquid-liquid interfaces - large eddy simulation and annurev- matsci- 070616- 124032.
PIV experiment, Chem. Eng. Sci. 198 (2019) 98–107, doi:10.1016/j.ces.2018.12. [59] G. Hwang, M. Kaviany, High-heat-flux distributed capillary evaporators, in:
040. K. Vafai (Ed.), Handbook of Porous Media, Taylor and Francis & CRC Press, 2015,
[47] M. Saeedipour, S. Schneiderbauer, Favre-filtered LES-VOF of two-phase flows pp. 535–555, doi:10.1201/b18614.
with eddy viscosity-based subgrid closure models: an a-posteriori analysis, [60] N.H. Afgan, L.A. Jovic, S.A. Kovalev, V.A. Lenykov, Boiling heat transfer from
Int. J. Multiphase Flow 144 (2021) 103780, doi:10.1016/j.ijmultiphaseflow.2021. surfaces with porous layers, Int. J. Heat Mass Transf. 28 (2) (1985) 415–422,
103780. doi:10.1016/0017- 9310(85)90074- 2.
[48] T.K. Kim, J. Ferreira, H. Jo, M. Kaviany, Flow-boiling canopy wick capillary- [61] G.-S. Hwang, M. Kaviany, Critical heat flux in thin, uniform particle coat-
viscous limit, Int. J. Heat Mass Transf. 181 (2021) 121999, doi:10.1016/j. ings, Int. J. Heat Mass Transf. 49 (5) (2006) 844–849, doi:10.1016/j.
ijheatmasstransfer.2021.121999. ijheatmasstransfer.2005.09.020.
[49] H.S. Ahn, H. Kim, H. Jo, S. Kang, W. Chang, M.H. Kim, Experimental study [62] G. Hwang, Y. Nam, E. Fleming, P. Dussinger, Y. Ju, M. Kaviany, Multi-artery heat
of critical heat flux enhancement during forced convective flow boiling of pipe spreader: experiment, Int. J. Heat Mass Transf. 53 (13) (2010) 2662–2669,
nanofluid on a short heated surface, Int. J. Multiphase Flow 36 (5) (2010) 375– doi:10.1016/j.ijheatmasstransfer.2010.02.046.
384, doi:10.1016/j.ijmultiphaseflow.2010.01.004. [63] D. Min, G. Hwang, Y. Usta, O. Cora, M. Koc, M. Kaviany, 2-D and 3-D modulated
[50] T. Lee, D.H. Kam, J.H. Lee, Y.H. Jeong, Effects of two-phase flow conditions on porous coatings for enhanced pool boiling, Int. J. Heat Mass Transf. 52 (11)
flow boiling CHF enhancement of magnetite-water nanofluids, Int. J. Heat Mass (2009) 2607–2613, doi:10.1016/j.ijheatmasstransfer.2008.12.018.
Transf. 74 (2014) 278–284, doi:10.1016/j.ijheatmasstransfer.2014.03.028. [64] G. Hwang, E. Fleming, B. Carne, S. Sharratt, Y. Nam, P. Dussinger, Y. Ju, M. Ka-
[51] P. Weber, K. Johanssen, Study of the critical heat flux condition at convective viany, Multi-artery heat-pipe spreader: lateral liquid supply, Int. J. Heat Mass
boiling of water: temperature and power-controlled experiments, in: 9th In- Transf. 54 (11) (2011) 2334–2340, doi:10.1016/j.ijheatmasstransfer.2011.02.029.
ternational Heat Transfer Conference, New York, 1990. [65] T.K. Kim, Enhancement of critical heat flux in a flow channel with capillary
[52] D. Cooke, S.G. Kandlikar, Effect of open microchannel geometry on pool boiling wicking structure, Pohang University of Science and Technology, Pohang, Re-
enhancement, Int. J. Heat Mass Transf. 55 (4) (2012) 1004–1013, doi:10.1016/j. public of Korea, 2022 Ph.D. thesis.
ijheatmasstransfer.2011.10.010. [66] F. Nicoud, F. Ducros, Subgrid-Scale stress modelling based on the square of
[53] A. Jaikumar, S.G. Kandlikar, Enhanced pool boiling heat transfer mechanisms the velocity gradient tensor, Flow Turbul. Combust. 62 (1999) 183–200, doi:10.
for selectively sintered open microchannels, Int. J. Heat Mass Transf. 88 (2015) 1023/A:10 099954260 01.
652–661, doi:10.1016/j.ijheatmasstransfer.2015.04.100. [67] Y. Murai, X. qun Song, T. Takagi, M. aki Ishikawa, F. Yamamoto, J. Ohta, Inverse
[54] P.A. Raghupathi, S.G. Kandlikar, Effect of thermophysical properties of the energy cascade structure of turbulence in a bubbly flow: PIV measurement and
heater substrate on critical heat flux in pool boiling, J. Heat Transf. 139 (11) results, JSME Int. J. Ser. B 43 (2) (20 0 0) 188–196, doi:10.1299/jsmeb.43.188.
(2017) 111502, doi:10.1115/1.4036653.
19

FB Crisis DNS and Porous Metasurfaces

Uploaded by

Copyright:

Available Formats

FB Crisis DNS and Porous Metasurfaces

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

FB Crisis DNS and Porous Metasurfaces

Uploaded by

Copyright:

Available Formats

International Journal of Heat and Mass Transfer 194 (2022) 123051

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer

Direct simulation of ﬂow-boiling crisis and its porous-metasurface

wavelength, the geometric, and the geometric-conﬁned modula-

2. Boiling hydrodynamic instability theories

2.1. Zuber pool boiling CHF theory and unit cell

Initially separated two-phase ﬂow has interfacial instabilities

Table 2 examine this as used in the interfacial lift-off submodel as it is

is highlighted at two elapsed times, 75 and 150 ms. After enough

where ρ ∗f = ρ f /(ρl + ρg ). From Eq. (39), we can infer the acceler-

No dryout is observed even though the heat ﬂux was increased

6. Porous metasurface for geometric modulation:

Metasurfaces are unit-cell based, generally 3-D surface struc-

The maximum total pressure drop allowed corresponds to the

50,100 5.0 0.25 1.93 0.15 0.04 0.20 2.32 2.32

Author Statement The initial interfacial perturbation is considered as a Fourier com-

where ω = 2π / f is the angular frequency. Using this pressure in

Substituting these in Eq. (A.6), we obtain

Fig. A1. Evolution of interfacial instability on a horizontal liquid-vapor interface

κ pg,o[ρg (ω + iκ ūg )]−1 = δoω + iδoκ ūg ,

The interfacial force term in the momentum conservation equa-

You might also like