# Analytical Solutions of Heat Transfer and Film Thickness with Slip Condition Effect in Thin-Film Evaporation for Two-Phase Flow in Microchannel

Authors: Ahmed Jassim Shkarah; Mohd Yusoff Bin Sulaiman; Md. Razali bin Hj Ayob
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Physical and mathematical model has been developed to predict the two-phase flow and heat transfer in a microchannel with evaporative heat transfer. Sample solutions to the model were obtained for both analytical analysis and numerical analysis. It is assumed that the capillary pressure is neglected (Morris, 2003). Results are provided for liquid film thickness, total heat flux, and evaporating heat flux distribution. In addition to the sample calculations that were used to illustrate the transport characteristics, computations based on the current model were performed to generate results for comparisons with the analytical results of Wang et al. (2008) and Wayner Jr. et al. (1976). The calculated results from the current model match closely with those of analytical results of Wang et al. (2008) and Wayner Jr. et al. (1976). This work will lead to a better understanding of heat transfer and fluid flow occurring in the evaporating film region and develop an analytical equation for evaporating liquid film thickness.


---

## Body

## 1. Introduction

Over the last decade, micromachining technology has been increasingly used to develop highly efficient heat sink cooling devices due to advantages such as lower coolant demands and smaller machinable dimensions. One of the most important micromachining technologies is the ability to fabricate microchannels. Hence, the studies of fluid flow and heat transfer in microchannels which are two essential parts of such devices have attracted attention due to their broad potential for solving both engineering and medical problems. Heat sinks are classified as either single-phase or two-phase according to whether liquid boiling occurs inside the microchannels. The primary parameters that determine the single-phase and two-phase operating regimes are the heat flux through the channel wall and the coolant flow rate. For a fixed heat flux (heat load), the coolant may maintain its liquid state throughout the microchannels. For a lower flow rate, the liquid coolant flowing inside the microchannel may reach its boiling point, causing flow boiling to occur, resulting in a two-phase heat sink [1–3]. Since the initial conception of microheat pipes in 1984 by Cotter [4], a number of analytical and experimental investigations have been conducted. Most of these investigations have concentrated on the capillary heat transport capability. With rapid advancements in microelectronic devices, their required total power and power density are significantly increasing. Thin-film evaporation plays an important role in these modern highly efficient heat transfer devices. When thin-film evaporation occurs in the thin-film region, most of the heat transfers through a narrow area between the nonevaporation region and intrinsic meniscus region as shown in Figure 1. The flow resistance of the vapor phase during thin-film evaporation is very small compared with the vapor flow in the liquid phase in a typical nucleate boiling heat transfer configuration. In addition, the superheat needed for the phase change in the thin-film region is much smaller than that for a bubble growth in a typical nucleate boiling, in particular, at the initial stage of the bubble growth. The heat transfer efficiency of thin-film evaporation is much higher than the nucleate boiling heat transfer. Because thin-film evaporation occurs in a small region, increasing thin-film regions and maintaining its stability are very important for heat transfer enhancement. It is known that the thermodynamic properties of the liquid thin-film region are very different from those of the macroregion. The effect of the intermolecular forces between the liquid thin film and wall can be characterized by the disjoining pressure, which controls the wettability and stability of liquid thin film formed on the wall. A better understanding of the evaporation mechanisms governed by the disjoining pressure in the thin-film region, especially for high heat flux, is very important in the development of highly efficient heat transfer devices. As early as 1972, Potash Jr. and Wayner Jr. [5] expanded the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory [6] to describe evaporation and fluid flow from an extended meniscus. Following this work, extensive investigations have been conducted to further understand mechanisms of fluid flow coupled with evaporating heat transfer in thin-film region. Stephan and Busse [7] developed a mathematical model based on the theoretical analysis presented by Wayner Jr. et al. [5, 8] to investigate the heat transfer coefficient occurring in small triangular grooves and found that the interface temperature variation plays an important role in the thin-film evaporation. Schonberg and Wayner Jr. [9] developed an analytical model by ignoring capillary pressure and found an analytical solution for the maximum heat evaporation from the thin film. Ma and Peterson [10] studied the thin-film profile, heat transfer coefficient, and temperature variation along the axial direction of a triangular groove. Hanlon and Ma [11] found that when particles become smaller, thin-film region can be significantly increased. More recently, Shao and Zhang [12] considered the effect of thin-film evaporation on the heat transfer performance in an oscillating heat pipe. Wang et al. [13, 14] established a simplified model based on the Young-Laplace equation and obtained an analytical solution for the total heat transfer in the thin-film region. In the current investigation, a mathematical model is established and its analytical solutions are obtained to evaluate the heat flux, total heat transport per unit length along the thin-film profile, thin-film thickness, and location for the maximum heat evaporation in the thin-film region and maximum heat transfer rate per unit length by thin-film evaporation.Figure 1
Schematic of an evaporating thin film.


## 2. Theoretical Analysis

Figure1 illustrates a schematic of an evaporating thin film formed on a wall. For the current investigation, it is assumed that fluid flow in the thin-film region is two-dimensional and pressure in the liquid film is a function of the x-coordinate only. The wall temperature, T
w, is greater than the vapor temperature, T
v. The momentum equation governing the fluid flow in the thin film can be found by
(1)
d
p
l
d
x
-
μ
l
∂
2
u
∂
y
2
=
0
,
where
(2)
u
=
u
wall
=
β
d
u
d
x
wall
y
=
0
,
δ
u
δ
y
=
0
y
=
δ
.
Solving (1) with boundary condition (2), we get
(3)
u
=
1
μ
l
d
p
l
d
x
y
2
2
-
δ
μ
l
d
p
l
d
x
1
y
+
u
wall
,
u
wall
=
β
d
u
d
x
wall
,
where β is the slip coefficient. If β is equal to zero, then no slip boundary condition is obtained.So we can find(4)
u
=
1
μ
l
d
p
l
d
x
y
2
2
-
δ
y
-
β
δ
,
where
(5)
β
=
β
0
1
1
-
γ
/
γ
c
,
where β
0 is the limiting slip length and γ
c represents the critical value of the shear rate. The slip coefficient under the condition of this study turns out to be approximately 1
*
10
-
9 m. From (4), the mass flow rate at a given location x can be found as
(6)
m
˙
x
=
∫
0
δ
ρ
l
u
l
d
y
=
ρ
l
μ
l
d
p
l
d
x
-
2
δ
3
-
6
β
δ
2
6
.
Taking a derivative of (6),
(7)
d
m
˙
x
d
x
=
d
d
x
ρ
l
μ
l
d
p
l
d
x
-
2
δ
3
-
6
β
δ
2
6
.
But
(8)
m
˙
e
=
-
d
m
˙
x
d
x
so the net evaporative mass transfer can be obtained as
(9)
m
˙
e
=
-
d
m
˙
x
d
x
=
d
d
x
ρ
l
μ
l
d
p
l
d
x
2
δ
3
+
6
β
δ
2
6
.
The heat transfer rate by evaporation occurring at the liquid-vapor interface in the thin-film region can be determined by
(10)
q
′′
=
m
˙
e
h
f
g
=
h
f
g
d
d
x
ρ
l
μ
l
d
p
l
d
x
2
δ
3
+
6
β
δ
2
6
.
But at the same time q
′
′ in (10) is equal to the heat transfer rate through the liquid thin film; that is,
(11)
q
′′
=
k
l
T
wall
-
T
δ
δ
.
From expanding the Clausius-Clapeyron equation, we can find
(12)
d
p
d
T
sat
⁡
=
h
f
g
T
v
1
ρ
v
-
1
ρ
l
,
(13)
T
δ
=
T
v
1
+
Δ
p
ρ
v
h
f
g
.
The pressure difference between vapor and liquid, Δ
p, at the liquid-vapor interface is due to both the capillary pressure and disjoining pressure and is expressed using the augmented Young-Laplace equation:
(14)
Δ
p
=
p
v
-
p
l
=
p
c
+
p
d
.
The disjoining pressure for a nonpolar liquid is expressed as
(15)
p
d
=
A
δ
3
,
where A is the dispersion constant and δ is the film thickness. The capillary pressure is the product of interfacial curvature K and surface tension coefficient σ:
(16)
p
c
=
σ
K
,
(17)
K
=
d
2
δ
d
x
2
1
+
d
δ
d
x
2
3
/
2
.
Substituting (13) into (11) the heat flux can be rewritten as
(18)
q
′′
=
T
wall
-
T
v
1
+
Δ
p
/
ρ
v
h
f
g
δ
/
k
l
.
Substituting (14) into (16) the heat flux can be rewritten as
(19)
q
′′
=
T
wall
-
T
v
1
+
p
c
+
p
d
/
ρ
v
h
f
g
δ
/
k
l
.
Substituting (19) into (10) gives
(20)
T
wall
-
T
v
1
+
p
c
+
p
d
/
ρ
v
h
f
g
δ
/
k
l
=
h
f
g
d
d
x
ρ
l
μ
l
d
p
l
d
x
2
δ
3
+
6
β
δ
2
6
.
We can find d
p
l
/
d
x by differentiating (14) with respect to x, and assume uniform vapor pressure, p
v, along the meniscus:
(21)
d
p
l
d
x
=
-
d
p
c
+
p
d
d
x
.
Substituting (21) into (20) yields
(22)
T
wall
-
T
v
1
+
p
c
+
p
d
/
ρ
v
h
f
g
δ
/
k
l
=
-
h
f
g
d
d
x
ρ
l
μ
l
d
p
c
+
p
d
d
x
2
δ
3
+
6
β
δ
2
6
.
For the evaporating thin-film region, the disjoining pressure is one dominant parameter, which governs the fluid flow in the evaporating thin-film region. And in the evaporating thin film region, the absolute disjoining pressure is much larger than the capillary pressure especially when the curvature variation along the meniscus is very small. In order to find the primary factor affecting the thin-film evaporation in the evaporating thin-film region, it is assumed that the capillary pressure, p
c, is neglected.Then, (22) becomes
(23)
T
wall
-
T
v
1
+
p
d
/
ρ
v
h
f
g
δ
/
k
l
=
-
h
f
g
d
d
x
ρ
l
μ
l
d
p
d
d
x
2
δ
3
+
6
β
δ
2
6
.
From (15) we find p
d and d
p
d
/
d
x;
(24)
d
p
d
d
x
=
-
3
A
δ
4
d
δ
d
x
.
Substituting (15) and (24) into (23) yields
(25)
T
wall
-
T
v
1
+
A
/
δ
3
/
ρ
v
h
f
g
δ
/
k
l
=
-
h
f
g
d
d
x
ρ
l
μ
l
-
3
A
δ
4
d
δ
d
x
2
δ
3
+
6
β
δ
2
6
.
We can rewrite the heat flux as
(26)
q
′
′
=
h
f
g
d
d
x
ρ
l
A
μ
l
d
δ
d
x
1
δ
+
3
β
δ
2
.
Simplifying (25) gives
(27)
1
h
f
g
k
l
T
wall
δ
-
k
l
T
v
δ
-
k
l
T
v
A
ρ
v
h
f
g
δ
4
=
d
d
x
ρ
l
A
μ
l
d
δ
d
x
1
δ
+
3
β
δ
2
.
We need to solve (27) to find d
δ
/
d
x.By multiplying 2 sides by2
*
d
δ
/
d
x
(
1
/
δ
+
3
β
/
δ
2
) and integrating two sides with respect to x we get
(28)
d
δ
d
x
1
δ
+
3
β
δ
2
2
+
D
=
2
μ
l
h
f
g
A
ρ
l
k
l
T
v
A
4
ρ
v
h
f
g
δ
4
+
3
β
k
l
T
v
A
5
ρ
v
h
f
g
δ
5
k
l
δ
T
v
-
T
wall
+
3
k
l
β
2
δ
2
T
v
-
T
wall
vvvvvvvvvvv
+
k
l
T
v
A
4
ρ
v
h
f
g
δ
4
+
3
β
k
l
T
v
A
5
ρ
v
h
f
g
δ
5
+
C
.
D, C are integration constants.So(29)
d
δ
d
x
2
=
1
1
/
δ
+
3
β
/
δ
2
2
×
2
μ
l
h
f
g
A
ρ
l
3
β
k
l
T
v
A
5
ρ
v
h
f
g
δ
5
k
l
δ
T
v
-
T
wall
vvvvvvvvvvvv
+
3
k
l
β
2
δ
2
T
v
-
T
wall
vvvvvvvvvvvv
+
k
l
T
v
A
4
ρ
v
h
f
g
δ
4
+
3
β
k
l
T
v
A
5
ρ
v
h
f
g
δ
5
+
C
-
D
.
Let
(30)
B
=
C
-
D
,
(31)
d
δ
d
x
=
1
1
/
δ
+
3
β
/
δ
2
2
vvvv
×
2
μ
l
h
f
g
A
ρ
l
+
k
l
T
v
A
4
ρ
v
h
f
g
δ
4
+
3
β
k
l
T
v
A
5
ρ
v
h
f
g
δ
5
vvvvvvv
×
3
β
k
l
T
v
A
5
ρ
v
h
f
g
δ
5
k
l
δ
T
v
-
T
wall
+
3
k
l
β
2
δ
2
T
v
-
T
wall
vvvvvvvvvv
+
k
l
T
v
A
4
ρ
v
h
f
g
δ
4
+
3
β
k
l
T
v
A
5
ρ
v
h
f
g
δ
5
+
B

1
1
/
δ
+
3
β
/
δ
2
2
1
/
2
.
To find constant B we apply the boundary condition
(32)
d
δ
d
x
=
0
δ
=
δ
o
.
So we get
(33)
B
=
-
2
μ
l
h
f
g
A
ρ
l
×
k
l
δ
o
T
v
-
T
wall
+
3
k
l
β
2
δ
o
2
T
v
-
T
wall
vvvvv
+
k
l
T
v
A
4
ρ
v
h
f
g
δ
o
4
+
3
β
k
l
T
v
A
5
ρ
v
h
f
g
δ
o
5
.
From (26) and (27),
(34)
q
′′
=
h
f
g
d
d
x
ρ
l
A
μ
l
d
δ
d
x
1
δ
+
3
β
δ
2
=
k
l
T
wall
δ
-
k
l
T
v
δ
-
k
l
T
v
A
ρ
v
h
f
g
δ
4
.
So
(35)
μ
l
k
l
/
h
f
g
ρ
l
A
T
w
-
T
v
-
μ
l
k
l
T
v
/
h
f
g
2
ρ
l
ρ
v
A
1
/
δ
3
δ
=
d
d
x
d
δ
d
x
1
δ
+
3
β
δ
2
.
Then
(36)
d
d
x
d
δ
d
x
(
1
δ
+
3
β
δ
2
)
=
k
l
υ
l
/
h
f
g
A
T
w
-
T
v
-
k
l
υ
l
T
v
/
h
f
g
2
ρ
v
δ
-
3
δ
.
Note that
(37)
υ
l
=
μ
l
ρ
l
.
The optimum thickness of the evaporating thin film, δ
optimum, where the heat flux reaches its maximum, can be found by taking a derivative of x for heat flux equation (34):
(38)
d
q
′′
d
x
=
0
=
d
k
l
T
wall
/
δ
-
k
l
T
v
/
δ
-
k
l
T
v
A
/
ρ
v
h
f
g
δ
4
δ
=
δ
optimum
d
x
.
So
(39)
-
T
wall
δ
optimum
2
+
T
v
δ
optimum
2
+
4
T
v
A
ρ
v
h
f
g
δ
optimum
5
=
0
.
Then
(40)
δ
optimum
=
4
A
T
v
ρ
v
h
f
g
T
w
-
T
v
3
.
For the nonevaporating film region, the heat flux is zero. Clearly, the interface temperature is equal to the wall temperature. The equilibrium thickness δ
o can be readily found by
(41)
q
′′
=
k
l
T
wall
δ
o
-
k
l
T
v
δ
o
-
k
l
T
v
A
ρ
v
h
f
g
δ
o
4
=
0
.
So
(42)
δ
0
≥
A
T
v
ρ
v
h
f
g
T
w
-
T
v
3
,
δ
0
min
⁡
=
A
T
v
ρ
v
h
f
g
T
w
-
T
v
3
.
From (31), (33), and (34) we get
(43)
δ
optimum
=
4
3
δ
o
.
The total heat transfer rate per unit width along the meniscus, q
t, can be calculated by
(44)
q
t
=
∫
0
x
q
′′
d
x
.
From (34)(45)
q
t
=
∫
0
x
h
f
g
d
d
x
ρ
l
A
μ
l
d
δ
d
x
(
1
δ
+
3
β
δ
2
)
d
x
.
So
(46)
q
t
=
h
f
g
ρ
l
A
μ
l
d
δ
d
x
1
δ
+
3
β
δ
2
.
Dimensionless thickness of the evaporating region can be defined as
(47)
δ
^
=
δ
δ
o
,
where δ
o is the thickness of the nonevaporating region.A dimensionless position also can be defined as(48)
ψ
=
x
δ
o
.
For dimensionless heat flux,
(49)
ϕ
=
q
′′
q
o
′′
,
where
(50)
q
o
′′
=
k
l
(
T
w
-
T
v
)
δ
o
is the heat flux at the interface temperature equal to the vapor temperature.Substituting (36), (34), and (50) in (49) we found
(51)
ϕ
=
h
f
g
ρ
l
A
μ
l
×
k
l
υ
l
/
h
f
g
A
T
w
-
T
v
-
k
l
υ
l
T
v
/
h
f
g
2
ρ
v
δ
-
3
/
δ
k
l
(
T
w
-
T
v
)
/
δ
o
=
δ
o
δ
-
δ
o
4
δ
4
,
or
(52)
ϕ
=
1
δ
^
1
-
1
δ
^
3
.
By considering d
ϕ
/
d
δ
^
=
0, we can determine the maximum dimensionless heat flux ϕ
max
⁡. The local heat flux through the evaporating thin film reaches its maximum when δ
^
=
4
1
/
3. Letting δ
^
=
4
1
/
3 in (52), the maximum dimensionless heat flux, δ
^
max
⁡, can be found as
(53)
δ
^
max
⁡
=
3
4
4
/
3
≈
0.473
.
Equation (53) indicates that the maximum heat flux occurring in the evaporating thin-film region is not greater than 0.473 times of the characteristic flux heat.Consider the heat transfer coefficient of(54)
q
′
′
=
h
T
s
-
T
v
.
To solve (31) we have two methods: (1)
exact solution or analytical solution,
(2)
numerical solution.We are going to solve (31) by using both methods


### 2.1. Analytical Solution

We can rewrite (31) in the form
(55)
d
δ
d
x
=
Σ
σ
δ
1
-
δ
0
δ
+
Ψ
σ
δ
1
-
δ
0
2
δ
2
vvvv
+
Ω
σ
δ
δ
0
4
δ
4
-
1
+
Π
σ
δ
δ
0
5
δ
5
-
1
1
/
2
,
where
(56)
σ
δ
=
δ
0
2
δ
0
/
δ
+
3
β
/
δ
0
δ
0
2
/
δ
2
2
,
Σ
=
2
μ
l
k
l
T
w
-
T
v
h
f
g
A
ρ
l
1
δ
0
,
Ψ
=
3
μ
l
β
T
w
-
T
v
h
f
g
A
ρ
l
1
δ
0
2
,
Ω
=
μ
l
k
l
T
v
2
h
f
g
2
ρ
v
1
δ
0
4
,
Π
=
6
μ
l
k
l
β
T
v
5
h
f
g
2
ρ
v
1
δ
0
5
,
Λ
=
3
β
δ
0
.
Since Σ
≫
Ψ
,
Π
,
Ω so we can rewrite (55) as
(57)
d
δ
d
x
=
Σ
δ
0
2
δ
0
/
δ
+
3
β
δ
0
2
/
δ
0
δ
2
vvvvvvv
×
1
-
δ
0
δ
+
Ψ
Σ
1
-
δ
0
2
δ
2
vvvvvvvvvv
+
Ω
Σ
δ
0
4
δ
4
-
1
+
Π
Σ
δ
0
5
δ
5
-
1
δ
0
2
δ
0
/
δ
+
3
β
δ
0
2
/
δ
0
δ
2
1
/
2
.
We introduce δ
^
=
δ
/
δ
0 and x
^
=
x
Σ to get
(58)
d
δ
^
d
x
^
=
δ
^
2
1
+
Λ
1
/
δ
^
2
vvvv
×
1
-
1
δ
^
+
Ψ
Σ
1
-
1
δ
^
2
vvvvvvvv
+
Ω
Σ
1
δ
^
4
-
1
+
Π
Σ
1
δ
^
5
-
1
δ
^
2
1
+
Λ
1
/
δ
^
2
1
/
2
.
Since 3
β
/
δ
0
≪
1, we can make the following approximation:
(59)
1
1
+
Λ
1
/
δ
^
2
=
1
-
2
Λ
δ
^
+
Λ
δ
^
2
.
To get
(60)
d
δ
^
d
x
^
=
δ
^
2
1
-
2
Λ
δ
^
+
3
Λ
δ
^
2
vvvv
×
1
-
1
δ
^
+
Ψ
Σ
1
-
1
δ
^
2
vvvvvvv
+
Ω
Σ
1
δ
^
4
-
1
+
Π
Σ
1
δ
^
5
-
1
1
-
2
Λ
δ
^
+
3
Λ
δ
^
2
1
/
2
,
d
δ
^
d
x
^
=
1
-
2
Λ
δ
^
+
3
Λ
δ
^
2
vvvv
×
1
δ
^
2
-
δ
^
2
+
Π
Σ
1
δ
^
3
-
δ
^
2
δ
^
2
-
δ
^
+
Ψ
Σ
δ
^
2
-
1
vvvvvvv
+
Ω
Σ
1
δ
^
2
-
δ
^
2
+
Π
Σ
1
δ
^
3
-
δ
^
2
1
-
2
Λ
δ
^
+
3
Λ
δ
^
2
1
/
2
we have
(61)
d
δ
^
d
x
^
=

Ψ
Σ
-
2
Λ
+
⋯
1
δ
^
+
⋯
1
δ
^
2
+
⋯
δ
^
2
1
+
Ψ
Σ
-
Ω
Σ
-
Π
Σ
vvvv
-
δ
^
1
+
2
Λ
1
+
Ψ
Σ
-
Ω
Σ
-
Π
Σ
vvvv
-
Ψ
Σ
-
2
Λ
+
⋯
1
δ
^
+
⋯
1
δ
^
2
+
⋯
1
/
2
.
Let
(62)
Γ
=
1
+
Ψ
Σ
-
Ω
Σ
-
Π
Σ
.
Because δ
^
≥
1, we neglect 1
/
δ
^
,
1
/
δ
^
2
,
…
. Therefore, the solution is accurate for δ
^
≫
1 as
(63)
d
δ
^
d
x
^
=
δ
^
2
Γ
-
δ
^
1
+
2
Λ
Γ
-
Γ
-
1
-
2
Λ
Γ
.
The general solution of (63) by using the initial condition
(64)
δ
^
=
1
x
^
=
0
is
(65)
δ
^
=
1
4
×
Γ
/
1
+
2
Λ
Γ
×
1
2
×
Γ
/
1
+
2
Λ
Γ
-
1
×
e
-
x
^
Γ
1
+
2
Γ
1
+
2
Λ
Γ
-
1
e
-
x
^
Γ
2
4
Γ
1
+
2
Λ
Γ
Γ
1
+
2
Λ
Γ
-
1
vvvvvvvvv
+
1
+
2
Γ
1
+
2
Λ
Γ
-
1
e
-
x
^
Γ
2
.
And from (65) we can evaluate an analytical solution for q
t  δ
→
∞ as
(66)
q
t
=
h
f
g
ρ
l
A
μ
l
Σ
.


### 2.2. Numerical Solution

The solution can be readily obtained using the fourth order Runge-Kutta method for the evaporating thin film profile. The governing equation (31) is solved with the use of a Runge-Kutta (4) method; the solution procedure is iterative. As a first guess δ
0, the values from the previous step are used, and the calculated values are returned from the Runge-Kutta solver and compared to the guess values. A comparison of the guess and calculated values is performed and looped until reaching convergence criteria for both film thicknesses. The numerical solver is coded in MATLAB. With these initial conditions, we have
(67)
δ
=
δ
0
x
=
0
,
d
δ
d
x
=
0
x
=
0
.


## 2.1. Analytical Solution

We can rewrite (31) in the form
(55)
d
δ
d
x
=
Σ
σ
δ
1
-
δ
0
δ
+
Ψ
σ
δ
1
-
δ
0
2
δ
2
vvvv
+
Ω
σ
δ
δ
0
4
δ
4
-
1
+
Π
σ
δ
δ
0
5
δ
5
-
1
1
/
2
,
where
(56)
σ
δ
=
δ
0
2
δ
0
/
δ
+
3
β
/
δ
0
δ
0
2
/
δ
2
2
,
Σ
=
2
μ
l
k
l
T
w
-
T
v
h
f
g
A
ρ
l
1
δ
0
,
Ψ
=
3
μ
l
β
T
w
-
T
v
h
f
g
A
ρ
l
1
δ
0
2
,
Ω
=
μ
l
k
l
T
v
2
h
f
g
2
ρ
v
1
δ
0
4
,
Π
=
6
μ
l
k
l
β
T
v
5
h
f
g
2
ρ
v
1
δ
0
5
,
Λ
=
3
β
δ
0
.
Since Σ
≫
Ψ
,
Π
,
Ω so we can rewrite (55) as
(57)
d
δ
d
x
=
Σ
δ
0
2
δ
0
/
δ
+
3
β
δ
0
2
/
δ
0
δ
2
vvvvvvv
×
1
-
δ
0
δ
+
Ψ
Σ
1
-
δ
0
2
δ
2
vvvvvvvvvv
+
Ω
Σ
δ
0
4
δ
4
-
1
+
Π
Σ
δ
0
5
δ
5
-
1
δ
0
2
δ
0
/
δ
+
3
β
δ
0
2
/
δ
0
δ
2
1
/
2
.
We introduce δ
^
=
δ
/
δ
0 and x
^
=
x
Σ to get
(58)
d
δ
^
d
x
^
=
δ
^
2
1
+
Λ
1
/
δ
^
2
vvvv
×
1
-
1
δ
^
+
Ψ
Σ
1
-
1
δ
^
2
vvvvvvvv
+
Ω
Σ
1
δ
^
4
-
1
+
Π
Σ
1
δ
^
5
-
1
δ
^
2
1
+
Λ
1
/
δ
^
2
1
/
2
.
Since 3
β
/
δ
0
≪
1, we can make the following approximation:
(59)
1
1
+
Λ
1
/
δ
^
2
=
1
-
2
Λ
δ
^
+
Λ
δ
^
2
.
To get
(60)
d
δ
^
d
x
^
=
δ
^
2
1
-
2
Λ
δ
^
+
3
Λ
δ
^
2
vvvv
×
1
-
1
δ
^
+
Ψ
Σ
1
-
1
δ
^
2
vvvvvvv
+
Ω
Σ
1
δ
^
4
-
1
+
Π
Σ
1
δ
^
5
-
1
1
-
2
Λ
δ
^
+
3
Λ
δ
^
2
1
/
2
,
d
δ
^
d
x
^
=
1
-
2
Λ
δ
^
+
3
Λ
δ
^
2
vvvv
×
1
δ
^
2
-
δ
^
2
+
Π
Σ
1
δ
^
3
-
δ
^
2
δ
^
2
-
δ
^
+
Ψ
Σ
δ
^
2
-
1
vvvvvvv
+
Ω
Σ
1
δ
^
2
-
δ
^
2
+
Π
Σ
1
δ
^
3
-
δ
^
2
1
-
2
Λ
δ
^
+
3
Λ
δ
^
2
1
/
2
we have
(61)
d
δ
^
d
x
^
=

Ψ
Σ
-
2
Λ
+
⋯
1
δ
^
+
⋯
1
δ
^
2
+
⋯
δ
^
2
1
+
Ψ
Σ
-
Ω
Σ
-
Π
Σ
vvvv
-
δ
^
1
+
2
Λ
1
+
Ψ
Σ
-
Ω
Σ
-
Π
Σ
vvvv
-
Ψ
Σ
-
2
Λ
+
⋯
1
δ
^
+
⋯
1
δ
^
2
+
⋯
1
/
2
.
Let
(62)
Γ
=
1
+
Ψ
Σ
-
Ω
Σ
-
Π
Σ
.
Because δ
^
≥
1, we neglect 1
/
δ
^
,
1
/
δ
^
2
,
…
. Therefore, the solution is accurate for δ
^
≫
1 as
(63)
d
δ
^
d
x
^
=
δ
^
2
Γ
-
δ
^
1
+
2
Λ
Γ
-
Γ
-
1
-
2
Λ
Γ
.
The general solution of (63) by using the initial condition
(64)
δ
^
=
1
x
^
=
0
is
(65)
δ
^
=
1
4
×
Γ
/
1
+
2
Λ
Γ
×
1
2
×
Γ
/
1
+
2
Λ
Γ
-
1
×
e
-
x
^
Γ
1
+
2
Γ
1
+
2
Λ
Γ
-
1
e
-
x
^
Γ
2
4
Γ
1
+
2
Λ
Γ
Γ
1
+
2
Λ
Γ
-
1
vvvvvvvvv
+
1
+
2
Γ
1
+
2
Λ
Γ
-
1
e
-
x
^
Γ
2
.
And from (65) we can evaluate an analytical solution for q
t  δ
→
∞ as
(66)
q
t
=
h
f
g
ρ
l
A
μ
l
Σ
.


## 2.2. Numerical Solution

The solution can be readily obtained using the fourth order Runge-Kutta method for the evaporating thin film profile. The governing equation (31) is solved with the use of a Runge-Kutta (4) method; the solution procedure is iterative. As a first guess δ
0, the values from the previous step are used, and the calculated values are returned from the Runge-Kutta solver and compared to the guess values. A comparison of the guess and calculated values is performed and looped until reaching convergence criteria for both film thicknesses. The numerical solver is coded in MATLAB. With these initial conditions, we have
(67)
δ
=
δ
0
x
=
0
,
d
δ
d
x
=
0
x
=
0
.


## 3. Results and Discussion

### 3.1. Comparison of Analytical Solution and Full Model

As presented above, a mathematical model for predicting evaporation and fluid flow in thin-film region is developed. Utilizing dimensionless analysis, analytical and numerical solutions are obtained for the heat flux distribution, total heat transfer rate per unit length, location of the maximum heat flux, and ratio of the conduction to convection thermal resistance in the evaporating film region. In order to verify the analytical solution derived herein, results predicted by Wang et al. [13] and numerical solution by Wayner Jr. et al. [8] are used. Figure 2 shows the comparison of analytical and numerical results of the total heat transfer rate through thin-film region with results presented by Wang et al. [13] and Wayner Jr. et al. [8]. Total heat flux is presented in function of the superheat temperature. Our numerical and analytical solutions are compared to the ones of Wang and Schonberg and Wayner. It shows good agreement for moderate superheat temperatures; however, our analytical solution tends to underestimate total heat flux at large superheat temperatures. In addition, the model presented herein can be used to predict analytically and numerically all of the equilibrium film thickness, heat flux distribution, film thickness variation of evaporating film region, maximum total heat transfer rate through the evaporating film region, and ratio of the conduction to convection thermal resistance. The following calculations and predictions are based on the thermal properties and operating conditions shown in Table 1.Table 1
Liquid properties and operating conditions.
Liquid
Water
A
10−20 J
T
v
353 K°
ρ
v
0.083 kg/m3
υ
l
0.4996 * 10−6 m2/s
h
f
g
2382700 J/kg
k
l
0.65 w/m⋅KFigure 2
Comparison of the total heat transfer rate through thin film region with results presented by Wang et al. [13] and Wayner Jr. et al. [8].


### 3.2. Comparison of Analytic Equation forδ with Previous Studies

It was seen from Table2 that, in many studies on a wide range of time, we did not find any researcher who found the analytical equation for δ even only approximately, but we found that they have numerical studies. So according to this, we can say that (65) is the first analytical equation for δ or at least approximately.Table 2
Comparison of previous studies on evaporating extended meniscus.
Authors
Numerical solution
analytical solution
Finding analytical equation forδ
Slip condition
Potash and Wayner [5]
ο
x
x
x
Moosman and Homsy [15]
ο
ο
x
x
Schonberg and Wayner [9]
ο
ο
x
x
Stephan and Busse [7]
ο
x
x
x
Schonberg et al. [16]
ο
ο
x
x
Shikhmurzaev [17]
x
ο
x
x
Ma  and Peterson [10]
ο
x
x
x
Pismen and Pomeau [18]
x
ο
x
x
Catton and Stroes [19]
ο
ο
x
x
Choi et al. [20]
ο
x
x
ο
Qu and Ma [21]
ο
ο
x
x
Choi et al. [22]
ο
x
x
ο
Park and  Lee [23]
ο
x
x
x
By Morris [24]
x
ο
x
x
Demsky and Ma  [25]
ο
x
x
x
Jiao et al. [26]
x
ο
x
x
Na et al. [27]
ο
ο
x
x
Sultan et al. [28]
x
ο
x
x
Wang et al. [14]
ο
x
x
x
Ma et al. [29]
ο
x
x
x
Wang et al. [13]
ο
ο
x
x
Zhao et al. [30]
ο
x
x
ο
Zhao et al.  [31]
ο
x
x
ο
Benselama et al. [32]
x
ο
x
x
Biswal et al. [33]
ο
x
x
ο
Liu et al.  [34]
x
ο
x
x
Bai et al. [35]
ο
ο
x
x
Biswal et al. [36]
ο
x
x
x
Thokchom et al. [37]
ο
x
x
x
Yang et al. [38]
ο
x
x
x
Current study
ο
ο
ο
ο


### 3.3. Distribution of Evaporative Film Thickness

Figures3, 4, 5, 6, 7, and 8 compare our numerical and analytical predictions for the dimensionless film thickness as function of the dimensionless position for 0.1 K, 0.5 K, 1 K, 2 K, 5 K, and 10 K superheat temperatures. They clearly show that the error of the analytical approximation is decreasing with increasing superheat temperature. We can see that for large positions the solution is dominated by exponential growth.Figure 3
Dimensionless evaporative film thickness profile at a superheat of 0.1 K.Figure 4
Dimensionless evaporative film thickness profile at a superheat of 0.5 K.Figure 5
Dimensionless evaporative film thickness profile at a superheat of 1 K.Figure 6
Dimensionless evaporative film thickness profile at a superheat of 2 K.Figure 7
Dimensionless evaporative film thickness profile at a superheat of 5 K.Figure 8
Dimensionless evaporative film thickness profile at a superheat of 10 K.Figures9, 10, 11, 12, 13, and 14 compare our numerical and analytical predictions for the dimensionless heat flux as function of the dimensionless position for 0.3 K, 0.7 K, 1 K, 2 K, 5 K, and 10 K superheat temperatures. Similar property can be seen; that is, the error of the analytical approximation is decreasing with increasing superheat temperature. The maximum of the heat flux is predicted correctly for all values of the superheat temperature analytically; however, the location of the maximum is predicted with significant error at low superheat temperatures.Figure 9
Dimensionless heat flux profile at a superheat of 0.3 K.Figure 10
Dimensionless heat flux profile at a superheat of 0.7 K.Figure 11
Dimensionless heat flux profile at a superheat of 1 K.Figure 12
Dimensionless heat flux profile at a superheat of 2 k.Figure 13
Dimensionless heat flux profile at a superheat of 5 K.Figure 14
Dimensionless heat flux profile at a superheat of 10 K.Figures15, 16, 17, and 18 show that nondimensional heat flux is presented as function of the nondimensional position for different values of β: 0, 0.5, 1 and 3
×
10
-
9. As β gets larger, the error in the analytical heat flux approximation is getting larger.Figure 15
Dimensionless heat flux profile at aβ
=
0 nm.Figure 16
Dimensionless heat flux profile at aβ
=
0.5 nm.Figure 17
Dimensionless heat flux profile at aβ
=
1 nm.Figure 18
Dimensionless heat flux profile at aβ
=
3 nm.Figures19, 20, 21, and 22 compare our numerical and analytical predictions for the dimensionless film thickness as function of the dimensionless position for different values of the slip coefficient, β: 0, 0.5, 1 and 3
×
10
-
9. Comparing the figures indicates similar conclusions as previously mentioned: smaller slip coefficient yields smaller error in the analytical approximation.Figure 19
Dimensionless evaporative film thickness profile at aβ
=
0 nm.Figure 20
Dimensionless evaporative film thickness profile at aβ
=
0.5 nm.Figure 21
Dimensionless evaporative film thickness profile at aβ
=
1 nm.Figure 22
Dimensionless evaporative film thickness profile at aβ
=
3 nm.


## 3.1. Comparison of Analytical Solution and Full Model

As presented above, a mathematical model for predicting evaporation and fluid flow in thin-film region is developed. Utilizing dimensionless analysis, analytical and numerical solutions are obtained for the heat flux distribution, total heat transfer rate per unit length, location of the maximum heat flux, and ratio of the conduction to convection thermal resistance in the evaporating film region. In order to verify the analytical solution derived herein, results predicted by Wang et al. [13] and numerical solution by Wayner Jr. et al. [8] are used. Figure 2 shows the comparison of analytical and numerical results of the total heat transfer rate through thin-film region with results presented by Wang et al. [13] and Wayner Jr. et al. [8]. Total heat flux is presented in function of the superheat temperature. Our numerical and analytical solutions are compared to the ones of Wang and Schonberg and Wayner. It shows good agreement for moderate superheat temperatures; however, our analytical solution tends to underestimate total heat flux at large superheat temperatures. In addition, the model presented herein can be used to predict analytically and numerically all of the equilibrium film thickness, heat flux distribution, film thickness variation of evaporating film region, maximum total heat transfer rate through the evaporating film region, and ratio of the conduction to convection thermal resistance. The following calculations and predictions are based on the thermal properties and operating conditions shown in Table 1.Table 1
Liquid properties and operating conditions.
Liquid
Water
A
10−20 J
T
v
353 K°
ρ
v
0.083 kg/m3
υ
l
0.4996 * 10−6 m2/s
h
f
g
2382700 J/kg
k
l
0.65 w/m⋅KFigure 2
Comparison of the total heat transfer rate through thin film region with results presented by Wang et al. [13] and Wayner Jr. et al. [8].


## 3.2. Comparison of Analytic Equation forδ with Previous Studies

It was seen from Table2 that, in many studies on a wide range of time, we did not find any researcher who found the analytical equation for δ even only approximately, but we found that they have numerical studies. So according to this, we can say that (65) is the first analytical equation for δ or at least approximately.Table 2
Comparison of previous studies on evaporating extended meniscus.
Authors
Numerical solution
analytical solution
Finding analytical equation forδ
Slip condition
Potash and Wayner [5]
ο
x
x
x
Moosman and Homsy [15]
ο
ο
x
x
Schonberg and Wayner [9]
ο
ο
x
x
Stephan and Busse [7]
ο
x
x
x
Schonberg et al. [16]
ο
ο
x
x
Shikhmurzaev [17]
x
ο
x
x
Ma  and Peterson [10]
ο
x
x
x
Pismen and Pomeau [18]
x
ο
x
x
Catton and Stroes [19]
ο
ο
x
x
Choi et al. [20]
ο
x
x
ο
Qu and Ma [21]
ο
ο
x
x
Choi et al. [22]
ο
x
x
ο
Park and  Lee [23]
ο
x
x
x
By Morris [24]
x
ο
x
x
Demsky and Ma  [25]
ο
x
x
x
Jiao et al. [26]
x
ο
x
x
Na et al. [27]
ο
ο
x
x
Sultan et al. [28]
x
ο
x
x
Wang et al. [14]
ο
x
x
x
Ma et al. [29]
ο
x
x
x
Wang et al. [13]
ο
ο
x
x
Zhao et al. [30]
ο
x
x
ο
Zhao et al.  [31]
ο
x
x
ο
Benselama et al. [32]
x
ο
x
x
Biswal et al. [33]
ο
x
x
ο
Liu et al.  [34]
x
ο
x
x
Bai et al. [35]
ο
ο
x
x
Biswal et al. [36]
ο
x
x
x
Thokchom et al. [37]
ο
x
x
x
Yang et al. [38]
ο
x
x
x
Current study
ο
ο
ο
ο


## 3.3. Distribution of Evaporative Film Thickness

Figures3, 4, 5, 6, 7, and 8 compare our numerical and analytical predictions for the dimensionless film thickness as function of the dimensionless position for 0.1 K, 0.5 K, 1 K, 2 K, 5 K, and 10 K superheat temperatures. They clearly show that the error of the analytical approximation is decreasing with increasing superheat temperature. We can see that for large positions the solution is dominated by exponential growth.Figure 3
Dimensionless evaporative film thickness profile at a superheat of 0.1 K.Figure 4
Dimensionless evaporative film thickness profile at a superheat of 0.5 K.Figure 5
Dimensionless evaporative film thickness profile at a superheat of 1 K.Figure 6
Dimensionless evaporative film thickness profile at a superheat of 2 K.Figure 7
Dimensionless evaporative film thickness profile at a superheat of 5 K.Figure 8
Dimensionless evaporative film thickness profile at a superheat of 10 K.Figures9, 10, 11, 12, 13, and 14 compare our numerical and analytical predictions for the dimensionless heat flux as function of the dimensionless position for 0.3 K, 0.7 K, 1 K, 2 K, 5 K, and 10 K superheat temperatures. Similar property can be seen; that is, the error of the analytical approximation is decreasing with increasing superheat temperature. The maximum of the heat flux is predicted correctly for all values of the superheat temperature analytically; however, the location of the maximum is predicted with significant error at low superheat temperatures.Figure 9
Dimensionless heat flux profile at a superheat of 0.3 K.Figure 10
Dimensionless heat flux profile at a superheat of 0.7 K.Figure 11
Dimensionless heat flux profile at a superheat of 1 K.Figure 12
Dimensionless heat flux profile at a superheat of 2 k.Figure 13
Dimensionless heat flux profile at a superheat of 5 K.Figure 14
Dimensionless heat flux profile at a superheat of 10 K.Figures15, 16, 17, and 18 show that nondimensional heat flux is presented as function of the nondimensional position for different values of β: 0, 0.5, 1 and 3
×
10
-
9. As β gets larger, the error in the analytical heat flux approximation is getting larger.Figure 15
Dimensionless heat flux profile at aβ
=
0 nm.Figure 16
Dimensionless heat flux profile at aβ
=
0.5 nm.Figure 17
Dimensionless heat flux profile at aβ
=
1 nm.Figure 18
Dimensionless heat flux profile at aβ
=
3 nm.Figures19, 20, 21, and 22 compare our numerical and analytical predictions for the dimensionless film thickness as function of the dimensionless position for different values of the slip coefficient, β: 0, 0.5, 1 and 3
×
10
-
9. Comparing the figures indicates similar conclusions as previously mentioned: smaller slip coefficient yields smaller error in the analytical approximation.Figure 19
Dimensionless evaporative film thickness profile at aβ
=
0 nm.Figure 20
Dimensionless evaporative film thickness profile at aβ
=
0.5 nm.Figure 21
Dimensionless evaporative film thickness profile at aβ
=
1 nm.Figure 22
Dimensionless evaporative film thickness profile at aβ
=
3 nm.


## 4. Conclusions

This paper presents a mathematical model for predicting evaporation and fluid flow in thin-film region. The thin-film region of the extended meniscus is delineated. Utilizing analytical solutions were obtained for heat flux distribution, total heat transfer, and liquid film thickness in the evaporating film region. The mathematical model also develops an analytic equation forδ. So according to this, we can say that (65) is the first analytical equation for δ or at least approximately. Also we can conclude that there is small effect of slip condition β on the thin-film evaporation for two-phase flow in microchannel. The dimensionless heat flux through thin-film region is a function of dimensionless thickness. Also the results showed the assumption that neglecting the capillary pressure is acceptable.


## Highlights

(i)
Analytical two-phase flow for microchannel heat sink is studied.
(ii)
Numerical two-phase flow for microchannel heat sink is studied.
(iii)
New evaporating film thickness equation is developed.



---

*Source: 369581-2015-01-01.xml*