- Email: [email protected]
Bénard–Marangoni instability on evaporating menisci in capillary channels
International Journal of Heat and Mass Transfer 63 (2013) 239–248
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Bénard–Marangoni instability on evaporating menisci in capillary channels q Zhenhai Pan a, Hao Wang a,b,⇑ a b
Laboratory of Heat and Mass Transport at Micro-Nano Scale, College of Engineering, Peking University, Beijing 100871, China Beijing Key Laboratory for Solid Waste Utilization and Management, College of Engineering, Peking University, Beijing 100871, China
a r t i c l e
i n f o
Article history: Received 28 November 2012 Received in revised form 30 March 2013 Accepted 31 March 2013
Keywords: Meniscus Capillary channel Evaporation Bénard–Marangoni Instability
a b s t r a c t Despite numerous studies about Marangoni instabilities, there are very few of them conducted for the menisci in microscale channels or tubes. In this paper, the onset of Marangoni flow on an initially isothermal meniscus in a capillary channel is studied. The temperature gradient perpendicular to the meniscus is found arousing a Bénard–Marangoni instability when the channel size or the temperature gradient beyond certain values. The threshold Marangoni Number for the instability is obtained, which turns out to be dependent on the Biot Number. The evolution of the meniscus flow patterns with increasing Marangoni Number is studied. The influences of gravity are discussed. The present work is a preliminary step to a comprehensive understanding of Marangoni instability in capillary structures. The temperature gradient that is perpendicular to the meniscus should be paid attention to in the related applications. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Marangoni flow has been proved of importance in various applications like heat transfer [1–4], microfluidic control [5,6], particle deposition [7,8], etc. In Pearson’s pioneer analysis study [9] on Bénard’s experiments [10], when a liquid layer with an isothermal interface is heated from below, small perturbation gets enlarged and surface-tension-induced flow presents. Such a problem, where the Marangoni flow is induced by the temperature gradient perpendicular to the meniscus, is usually referred to as Bénard–Marangoni instability and has been studied extensively as a classical problem for decades [11–15]. Up to today, the study on Bénard–Marangoni instability is still going on for different situations and applications. For example, Doumenc et al. [16] studied the transition of Bénard– Marangoni flow on an infinite evaporating liquid layer with an adiabatic bottom. The numerical results agreed very well with previous Toussaint et al.’s experiments [17]. Touihri et al. [18] investigated the threshold and nonlinear variation of the Marangoni flow in a cylinder. Factors including Biot number, geometry aspect ratio and buoyancy force are considered. Benselama et al. [19] theoretically studied the Marangoni instability on an evaporating thin film
q Submitted for publication in International Journal of Heat and Mass Transfer, Nov. 2012. ⇑ Corresponding author at: Beijing Key Laboratory for Solid Waste Utilization and Management, College of Engineering, Peking University, Beijing 100871, China. Tel.: +86 10 8252 9010. E-mail address: [email protected] (H. Wang).
0017-9310/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.03.082
extending from a meniscus. The capillary pressure and disjoining pressure were found playing important roles in the hydrodynamic instabilities. With the fast development of microfluidics and MEMS applications, Marangoni flows in microscale/capillary structures have gained increasing attentions. Particularly, studies have been conducted for the evaporating menisci in capillary channels or tubes [20–27]. For most of the studies, the meniscus is initially non-isothermal, and thus Marangoni convection is directly induced by the temperature gradients along the meniscus. Usually, such Marangoni convection is also referred to as Thermalcapillary flow. For example, in Dhavaleswarapu et al.’s [24] and Chamarthy et al.’s [25] experiments on concave menisci, the evaporations were stronger at the meniscus edges. The meniscus edges were thus cooler than the centers, and the Marangoni convections were observed flowing from the meniscus centers to the edges. Similar results were found in Wang et al.’s [26] and Lan et al.’s [27] simulation studies. In Ward and Duan’s [28,29] and Buffone and Sefiane’s [30] experiments, the hot walls near the meniscus warmed the meniscus edges while the evaporation cooled the meniscus centers, inducing inward Marangoni convections. There are very few studies about the Marangoni instability in the capillary structures, even fewer for Bénard–Marangoni instability. Buffone et al. [20,31] noticed that the Marangoni flow was not stable in a horizontal tube and attributed it to the meniscus temperature difference near the contact line. Recently, Pan and Wang numerically [32] and experimentally [33] found that on a convex meniscus at the mouth of a microchannel, the initial
240
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
Nomenclature A Bo Bi C cp Cr D d hfg J k l Ma M mnet m00net ~ n p Pr r R Ra Sm Sh T T⁄ V q
area (m2) Bond number Biot number vapor molar concentration (mol/m3) heat capacity (J/kg K) Crispation number diffusion coefficient in air (m2/s) the channel width (m) latent heat of evaporation (J/kg) evaporating mole flux (mol/m2 s) thermal conductivity (W/ m K) the channel length (m) Marangoni number molecular weight (kg/kmol) mass flow rate (kg/s) mass flux (kg/m2 s) interface normal coordinate pressure (N/m2) Prandtl number radius (m) universal gas constant (J/mol K) Rayleigh number mass source (kg/m3 s) energy source (W/m3) temperature (K) non-dimensional time fluid velocity (m/s) heat flux (W/m2)
Greeks symbols a thermal diffusivity (m2/s) b thermal expansion coefficient (1/K) d thickness (m) h contact angle m kinematic viscosity (m2/s) l dynamic viscosity (N s/m2) q density (kg/m3) r surface tension coefficient (N/m) r interface accommodation coefficient s shear stress (N/m2) Subscripts am ambient Bous Boussinesq e evaporation Equ equilibrium g gas (vapor/air mixture) l liquid lv interface outlet outlet ref reference sat saturated v vapor x x component y y component z z component
symmetrical Marangoni flow is not stable and tends to evolve into an asymmetrical single vortex flow. The results also explained the abnormal flow pattern observed in Gazolla et al.’s experiment [34] where unexpected asymmetrical flow always presented under symmetrical boundary conditions. In the present paper we first study the onset of Marangoni flow on a flat, initially isothermal meniscus inside a capillary channel. It is found the perpendicular-to-meniscus temperature gradient (from the meniscus to the warm bulk) can induce Bénard–Marangoni instability as long as the channel size or the evaporation intensity beyond certain thresholds. The critical Ma values for the instability are obtained which turn out to be sensitive to the Biot number. Then the flow variations with Ma are presented and the influences of the gravity are discussed. It is suggested that the Bénard–Marangoni instability induced by the perpendicular temperature gradient
should be paid attention to in related microfluidic and MEMS designing. 2. Problem description A capillary channel with square cross section is shown in Fig. 1. Water at room temperature (298.15 K) is supplied from the channel inlet. A meniscus is formed in the channel. The channel outlet is open to the ambient air and water is volatizing i.e. evaporating on the meniscus. The meniscus temperature is lower than the bulk due to the latent heat dissipated by the evaporation. To study Bénard–Marangoni instability the meniscus has to be initially isothermal. In the present work to remove the geometry factors and catch the essentials, the meniscus shape is assumed to be flat. The evaporation on the meniscus is uniform, resulting
Meniscus
Outlet
Vapor transport
Inlet
Compensating flow
y
θ
z x
meniscus
Outlet Gas domain
lv
d
Inlet
Liquid domain
ll
Fig. 1. Schematic diagram of the open channel with an evaporating meniscus.
241
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
shear stress. Hence the Marangoni flow is produced. The shear stress at the meniscus is calculated as
Table 1 Fluid properties. Properties 3
Density (kg/m ) Thermal conductivity (W/m K) Thermal capacity (J/kg) Viscosity (kg/m s) Vapor molecular weight M (kg/mol) Thermal expansion coefficient (1/K) dr/dT (N/m K) Diffusion coefficient (m2/s) Accommodation coefficient for evaporation Saturated vapor pressure (Pa) Latent heat (J/kg)
Water
Air
997.05 at 298.15 K 0.6 4182 0.001003 18 103 0.000251 0.00015 2.55 105at 300 K 0.03
1.182 0.0242 1006.43 1.789 105 29 103 0.003355 / /
2334.6 at 293.15 K 2.26 106
/
sl jlv ¼
3.2.2. Meniscus deformation and profiles On a flat meniscus, the meniscus deformation is evaluated with Crispation number as suggested by Davis [35] and Doumenc et al. [16]. Crispation number is defined as
Cr ¼
in a uniform temperature distribution (details referring to the simulation results shown in Section 5). In practice, the flat meniscus could be true when the contact angle is nearly 90°, or the meniscus edge is pinned. The channel outlet is open to the ambient air and water is volatizing. The vapor transports through the channel outlet into the ambient as illustrated in Fig. 1. The vapor fraction on the channel outlet is assumed to be uniform and is determined by the relative humidity of the ambient. The channel width is d. The distances from the meniscus to the inlet and the outlet are ll and lv respectively. The ambient temperature Tam is 298.15 K. Fluid properties are listed in Table 1.
la rd
3.2.3. Evaporation Kinetic theory of gases has provided the molecular basis for understanding the evaporation for over a century. Hertz [36] and Knudsen [37] studied the evaporation of liquid mercury into vacuum and proposed Hertz–Knudsen theory. Schrage [38] further modified the theory by including the net macroscopic velocity of vapor and proposed Hertz–Knudsen–Schrage equation, that is, the evaporated mass flux across the interface at one moment is calculated as
m00net ¼
3.1. Liquid domain and gas domain
! qr V ¼ Sm !
!
ð1Þ !
q V r V ¼ rp þ r ðlr V Þ þ qBous~ g !
ð2Þ
qcp V rT ¼ kr2 T þ Sh
ð3Þ
qBous ¼ qref qref bðT T ref Þ
ð4Þ
3.2. Liquid–gas interface 3.2.1. Marangoni effect The temperature gradient along the meniscus would induce surface tension gradient. The tangential force due to the surface tension gradient on the meniscus should be balanced by the liquid
! 1=2 pv equ ðT lv Þ 2r M pv 2pR 2r T 1=2 T 1=2 v lv
ð7Þ
Assume that the interfacial vapor temperature Tv is equal to the interfacial liquid temperature Tlv, [23,32] leading to
m00net
Laminar flow of a Newtonian fluid is assumed in the liquid domain based on the small scale of the system. Boussinesq approximation is employed for modeling the natural convection such that the density changes with temperature only in the buoyancy term. The laminar and Newtonian fluid with Boussinesq approximation is also adapted for the air/vapor mixture since the pressure and temperature difference is small in the calculated gas domain. The continuity, momentum and energy equations and Boussinesq approach are given respectively as
ð6Þ
For a small Cr number, the surface deformation is negligible comparing with the channel width d. In the present study with temperature T 298 K, we have Cr < 2 105 with surface tension r 7.2 102 N/m, liquid viscosity l 1 103 kg/(m s), thermal diffusivity a 1.4 107 m2/s, channel width d > 1 104 m. Therefore, the deformation of the interface is neglected in the present simulations considering such small Cr values.
3. Mathematical modeling The heat and mass transport in the capillary channel includes the evaporation at the meniscus, vapor diffusion and convection in the gas domain, Marangoni and buoyancy flow in the liquid domain. A mathematical model is described in the following sections and all the factors will be solved together.
ð5Þ
where @T/@s is the temperature gradient along the meniscus.
/
/
dr @T dT @s
1=2 pv 2r M ¼ 2pR 2r
equ ðT lv Þ 1=2 T lv
pv
! ð8Þ
The capillary pressure caused by the meniscus curvature could affect the the equilibrium vapor pressure Pv_equ according to Kelvin equation, which is considerable when the curvature radius is smaller than several microns [46,47]. In our case, however, the curvature radius is 0.1 mm so the equilibrium vapor pressure is approximated as the corresponding saturation pressure.
pv
equ ðT lv Þ
psat ðT lv Þ ¼ psat;ref exp
Mhfg 1 1 T sat;ref T lv R
ð9Þ
^ is one source of uncertainty. As The accommodation coefficient r noted in [39], the evaporation and condensation coefficients of water that are obtained theoretically or experimentally scatter over a rather large range of more than two decades. Paul [40] estimated the value of the accommodation coefficient to be in the range of 0.02–0.04. Hickman [41] found that very high evaporation coefficients (close to unity) could only be obtained with moving liquids in a vacuum where the surface is continually refreshed by the moving stream and the evaporating molecules are removed immediately by the vacuum. Marek and Straub [39] commented that even small contaminations of the surface significantly reduce the interfacial mass transfer. The lowest value of the accommodation coefficient given in their review was 0.002 for water in a glass vessel. Badam et al. [42] estimated the value of the accommodation coefficient to be in the range of 0.028–0.1. In the present work, ^ is assumed to be a small value the accommodation coefficient r of 0.03 based on the previous works [32,40,42]. The possible
242
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
uncertainty induced by this coefficient value would be discussed in Section 4.
two layers of mesh cells adjacent to the interface as described specifically in [26]. The mass and energy source terms in the cells along the interface in both liquid and gas domains are calculated as
3.3. Coupling of vapor transport and evaporation
Sm;g ¼ The diffusion of vapor from the meniscus into the ambient air and its coupling to the evaporation has been described specifically in our previous work [23]. Briefly, the governing equation for the vapor transport in air is given as
@C v þ~ v rC v ¼ r ðDrC v Þ @t
ð10Þ
At the meniscus, the mass flux dissipated by the vapor transport can be obtained by
@C v m00v ¼ MJ ¼ M D þ v n C v @n lv
ð11Þ
The first term on the right side is the vapor transport due to mass diffusion and the second term is due to mass convection. This convection is due to the single direction of diffusion and known as Stefan flow [43]. After a series of derivations the evaporation at meniscus is coupled to the vapor transport in gas domain [23]:
1=2 MD @C v 2r M 1 ¼ ðpv 1=2 2pR 1 C v =C total @n lv 2 r T
equ ðT lv Þ
pv Þ
ð12Þ
lv
3.4. Boundary conditions The channel walls are treated as no-slip and adiabatic boundaries. A given stagnation-pressure inlet for the liquid and a given static pressure outlet [44] for the vapor are set respectively at the two ends of the tube. The temperature at both ends are assumed as a room temperature of 298.15 K. The vapor concentration is set based on the local vapor pressure.
p Cv ¼ v RT
ð13Þ
m00 Af ;g m00 Af ;l ; Sm;l ¼ V cell;g V cell;l
ð14Þ
Sh;l ¼ Sm;l hfg þ Sm;l cp;l ðT cell T ref Þ
ð15Þ
Sh;g ¼ Sm;g cp;g ðT cell T ref Þ
ð16Þ
Here, Af is the interface area of the cell. Vcell and Tcell are the volume and the temperature of the cell. The first term on the right side of Eq. (15) represents the latent heat consumed by the evaporation. The second term on the right side of Eq. (15) and the only term on the right side of Eq. (16) represent the sensible heat corresponding to the evaporative mass. Tref is the reference temperature used in ANSYS 12, which is 298.15 K. Similar mathematical modeling and numerical methodology have been validated in authors’ previous studies on evaporating menisci in microtube [26], V-groove [23], and at the mouth of capillary channel [32,33]. The obtained flow patterns and evaporation rates were all highly consistent with the corresponding experiments. For flat menisci, the mesh setup is illustrated in Fig. 1. Hexahedral elements are created for liquid and gas domains. The mesh independence test is performed. Four sets of meshes are checked under two different operating conditions as shown in Tables 2a and 2b. Mesh 3 is employed in the simulations. ^ is As mentioned in Section 3.2, the accommodation coefficient r a source of uncertainty. A series of tests are conducted to estimate the uncertainty. It is seen from Table 3 that the evaporation intensity is not sensitive to the choice of the coefficient values (the largest deviation is 5% when the coefficient value changes from 0.002 to 1). This result is reasonable because the primary resistance for water evaporation to open air relies on the vapor diffusion in the gas domain rather than the interfacial evaporation resistance [45,47], i.e. even when the accommodation coefficient is 0.002 the evaporation is still efficient enough compared to the vapor diffusion.
4. Numerical treatment 5. Results and discussions The numerical solution is obtained using the pressure based finite volume scheme. Software package ANSYS 12.0 (FLUENT SOLVER) enhanced with user-developed codes are employed. The SIMPLE algorithm is used for pressure–velocity coupling. To calculate the heat and mass exchange across the interface, it is assumed that the mass transport across the interface occurs only within the
5.1. Occurrence of instability when increasing channel size On the flat meniscus in the adiabatic channel as shown in Fig. 1, the evaporation distribution is uniform and consequently the temperature is uniform on the meniscus. However, Marangoni flow
Table 2a Mesh-independence studies under a large Ma value (d = 1000 lm, ll = 3 mm, lv = 1 mm, outlet humidity: 0%). Mesh setup (Nx Ny Nz)
Mesh 1 13 13 67
Mesh 2 19 19 95
Mesh 3 27 27 133
Mesh 4 40 40 140
Evaporation rate (kg/s) Deviation from mesh 4 Marangoni number Deviation from mesh 4
4.4424 1010 0.43% 1504.5 1.64%
4.4359 1010 0.28% 1513.0 1.09%
4.4289 1010 0.12% 1522.3 0.48%
4.4236 1010 0 1529.6 0
Table 2b Mesh-independence studies under a small Ma value (d = 600 lm, ll = 3 mm, lv = 1 mm, outlet humidity: 0%). Mesh setup (Nx Ny Nz)
Mesh 1 (13 13 67)
Mesh 2 (19 19 95)
Mesh 3 (27 27 133)
Mesh 4 (40 40 140)
Evaporation rate (kg/s) Deviation from mesh 4 Marangoni number Deviation from mesh 4
9.3867 1011 0.27% 356.43 1.02%
9.3742 1011 0.14% 358.16 0.53%
9.3661 1011 0.05% 359.29 0.22%
9.3610 1011 0 360.08 0
243
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248 Table 3 Tests of the accommodation coefficient (ll = 3 mm, lv = 1 mm, outlet humidity: 0%). Channel width (lm) 400 400 400 800 800 800 1000 1000 1000
Accommodation coefficient
Evaporation rate (kg/s) 11
0.002 0.03 1 0.002 0.03 1 0.002 0.03 1
6.62 10 6.91 1011 6.93 1011 2.69 1010 2.80 1010 2.81 1010 4.28 1010 4.43 1010 4.45 1010
occurs when the channel width d is greater than a critical value. Five channels with different d values, from 200 lm to 1000 lm, are simulated and the results of the temperature and flow fields are shown in Fig. 2. The ambient humidity at the channel outlet is fixed at 40%. The resulting evaporation flux is almost the same, 0.261 ± 0.002 g/(m2 s). Based on Fig. 2 when the channel width is below 600 lm, the meniscus is isothermal (the maximum temperature difference across the meniscus is about 0.0001 K). An extremely weak (velocity below 1 lm/s) flow is observed flowing from the meniscus center to the edges. This very weak flow arises due to the following reason: firstly, there is a compensating flow through the channel, from the inlet to the meniscus, to compensate the mass loss due to the evaporation; secondly, due to the friction of the channel walls, the compensating flow has higher velocity in the central than the edge of the channel, which brings more liquid to the meniscus center than to the meniscus edges.
295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272
Interfacial temperature (K)
Marangoni number
293.49 293.27 293.22 293.69 293.50 293.47 293.85 293.70 293.68
254.9 266.9 269.5 977.2 1018.8 1024.9 1471.0 1522.3 1529.3
However, when the channel width increases to 600 lm, a strong convection with a totally different pattern presents. The liquid flows from a corner to another on the meniscus with a max velocity up to 3900 lm/s, which is thousands times larger than the evaporation compensating flow. Corresondingly, the maximum temperature difference along the meniscus reaches 0.1 K, which is one thousand times larger than that of the 400 lm channel. For even bigger channels, i.e., 800 and 1000 lm in width, similar strong flows are obtained. 5.2. Occurrence of instability when increasing evaporation intensity If the channel width is fixed at 600 lm but the evaporation intensity is increased by reducing the ambient humidity at the channel outlet, the results are shown in Fig. 3. When the humidity is 80%, the evaporation on the meniscus is weak, with average evaporative mass flux of 8.8 105 kg/(m2 s). A uniform tempera-
0.001
0.000
0.000
d=200µm; Ma=39 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201
0.002
0.001
0.000
d=400µm; Ma=161 295.380 295.380 295.362 295.343 295.325 295.307 295.288 295.262 295.270 295.252 295.233 295.215 295.197 295.178 295.160 295.160
5.048
2.524
0.000
d=800 µm; Ma=622
295.300 295.300 295.287 295.273 295.260 295.247 295.233 295.220 295.220 295.207 295.193 295.180 295.167 295.153 295.140 295.140
3.907
1.963
0.019
d=600 µm; Ma=358 295.480 295.480 295.457 295.433 295.410 295.387 295.363 295.340 295.340 295.317 295.293 295.270 295.247 295.223 295.200 295.200
5.661
2.945
0.228
d=1000 µm; Ma=952
Fig. 2. Meniscus temperature (K) (left) and flow pattern (mm/s) (right) as channel width d increases (humidity = 40%, ll = 3 mm, ll/lv = 3).
244
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176
0.001
0.001
0.000
Humidity: 80%; Ma=120 296.210 296.210 296.205 296.200 296.195 296.190 296.185 296.180 296.180 296.175 296.170 296.165 296.160 296.155 296.150 296.150
295.300 295.300
1.628
295.287 295.273 295.260 295.247 295.233 295.220 295.220 295.207 295.193 295.180 295.167 295.153 295.140 295.140
0.815
0.000
Humidity: 60%; Ma=242 294.402 294.402 294.382 294.362 294.342 294.322 294.302 294.281 294.281 294.261 294.241 294.221 294.201 294.181 294.161 294.161
3.907
1.963
0.019
Humidity: 40%; Ma=358 293.485 293.485
5.273
293.462 293.438 293.414 293.391 293.368 293.344 293.344 293.320 293.297 293.274 293.250 293.226 293.203 293.203
2.666
0.038
6.215
3.142
0.069
Humidity: 0%; Ma=588
Humidity: 20%; Ma=473
Fig. 3. Meniscus temperature (K) (left) and flow pattern (mm/s) (right) as the humidity at the outlet decreases and the evaporation increases (ll = 3 mm, ll/lv = 3, d/ll = 0.2).
ture is obtained on the meniscus (the maximum temperature difference is about 0.0001 K). However, when the humidity decreases to 60% and the evaporation on the meniscus increases to 1.77 104 kg/(m2 s), the temperature difference jumps up to 0.1 K and a strong Marangoni convection occurs.
!
! V rT ¼ ar2 T
ð19Þ
At the meniscus where z = 0, assuming the latent heat of the evaporation is supplied from the liquid side, we have
k
!
rV ¼ 0
ð17Þ
!
ð18Þ
5.3. Bénard–Marangoni instability in capillary channels Above Marangoni convections occurring on the initially isothermal menisci could be categorized as one type of Bénard–Marangoni instability. The mechanism of the instability could be described by considering a perturbation occurs on an initially isothermal meniscus in a microchannel, as illustrated in Fig. 4a. As the meniscus keeps being cooled by the evaporation, the meniscus temperature T1 at the upriver of the perturbation flow must be higher than the downriver temperature T2. This temperature difference between T1 and T2 would then induce a Marangoni flow along the meniscus, whose flow direction is same to the perturbation. A positive feedback is then established. However, at the same time, the thermal diffusion tends to eliminate the temperature difference between T1 and T2 while the liquid viscosity tends to eliminate the induced Marangoni flow. If the former positive feedback effect is stronger than the thermal diffusion and liquid viscosity effect, i.e. Ma number is high enough, the Bénard–Marangoni instability occurs. As discussed in Section 5.1, the evaporation compensating flow is very weak compared to the Marangoni flow. By neglecting the evaporation compensating flow, the governing equations and the boundary conditions for the liquid domain are
!
q V r V ¼ rp þ lr2 V qgbðT T ref ÞÞe!g
l
@T ¼ hfg m00 ðTÞ @z
ð20Þ
@V x dr @T @V y dr @T ¼l ¼0 @z dT @x @z dT @y
ð21Þ
If choose characteristic length d, velocity a/d, pressure qa2/d2, and temperature DTd/ll in which DT = Tam Tlv,ave is the temperature drop from the inlet to the meniscus, we have dimensionless equations as
!
r V ¼ 0
ð22Þ
! ! ! ! V r V ¼ rp þ Prr2 V þ PrðRa1 T Ra2 Þeg
ð23Þ
! V rT ¼ r2 T
ð24Þ
in which the non-dimensional temperature T⁄ = (T Tlv,ave)ll/DTd. Tlv,ave is the average temperature on the meniscus. For the boundary conditions, at the meniscus, Eq. (20) can be linearized as 00
k
@T dm ¼ hfg m00 ðT lv ;av e Þ þ ðT lv ;av e ÞðT T lv ;av e Þ @z dT
ð25Þ
Meanwhile, according to the conduction through the liquid domain in channel
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
Evaporation cooling flux
296.4430 296.4430
0.007
296.4430 296.4430 296.4430 296.4430 296.4429 296.4429 296.4429 296.4429 296.4429 296.4429 296.4428 296.4428 296.4428 296.4428 296.4428
Meniscus
T1
T2 Channel wall
Channel wall Perturbation flow
Heat flux
245
0.004
0.000
Warm bulk
(a)
(b) d=1000 μm d=800 μm
210
ll/ lv=2, ll=2mm
d=600 μm
Critical Ma
ll/ lv=3, ll=3mm ll/ lv=6,ll=6mm
200
d=400 μm
ll/ lv=10, ll=10mm
d=200 μm
190
180 -6 10
d=100 μm
10-5
10-4
10-3
10-2
0.1
Bi number
(c) Fig. 4. The onset of the Marangoni instability: (a) the illustration of the mechanism; (b) an example showing the critical state (d = 600 lm; ll = 3 mm; Mac = 210); (c) the critical Ma values varying with Bi numbers.
hfg m00 ðT lv ;av e Þ k
DT ll
ð26Þ
with which Eq. (20) could finally be non-dimensionalized as
@T ¼ 1 þ BiT @z
ð27Þ
At the same time Eq. (21) is non dimensionalized as @V x @T @V y @T þ Ma ¼ þ Ma ¼ 0 @z @x @z @y
ð28Þ
In above dimensionless equations, Eqs. (22)–(24), (27) and (28), we have 4
gbd DT Ra1 ¼ ; amll 2
3
Ra2 ¼
dr d DT ; Ma ¼ dT lall
gbd DT
am
m ; Pr ¼ ; a 00
hfg d dm Bi ¼ ðT lv Þ k dT
ð29Þ
2
dr d hfg m00 dT la k
5.4. Influences of Biot number As discussed above, the critical Ma value for the onset of the Bénard–Marangoni instability should be dependent on the Bi number. Since the vapor transport is dominated by the vapor diffusion (as discussed in Section 4), the evaporating mass flux m00 could be approximated as
m00 MD
By combining with Eq. (26), Ma can also be written as
Ma ¼
distribution and interfacial flow pattern is significantly different from those before the Bénard–Marangoni instability occurs, the temperature gradient and the flow velocity are still very small (T 0.001 K; V 1 lm/s). This state is recogniazed as the critical state of the onset of the instability. One example is shown in Fig. 4b when d = 600 lm and Ma number comes to the threshold Mac = 210. Once the Ma number further increases and exceeds the critical value, the temperature and the velocity magnitude would increase sharply and establish strong flows as shown in Fig. 3.
ð30Þ
In order to obtain critical Ma values, the humidity at the outlet is lowered step by step (thus the evaporation intensity increases and Ma increases) for each case. As mentioned in Sections 5.1 and 5.2, the meniscus temperature is almost uniform before the Bénard–Marangoni instability occurs. Once the Ma value reaches the threshold, the temperature peak moves suddenly to one corner of the meniscus and a corresponding Marangoni flow is induced flowing from the location of the temperature peak to the other corners (Fig. 4b). At this critical state, even if the temperature
C lv C outlet lv
ð31Þ
Combining Eq. (13) with Eq. (31), we get
m00
MD pv ;lv pv ;outlet Rlv T lv T outlet
ð32Þ
Still because of the vapor transport is dominated by the vapor diffusion, the vapor pressure at the interface pv,lv is approximately equal to the saturation pressure at the interface psat(Tlv). Combining Eq. (32) with Eq. (9), we get
Mhfg @m00 MD 1 p ðT Þ 1 l v Rlv T 2lv sat @T lv RT lv
ð33Þ
246
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
Combining Eqs. (29), (31), and (33), we have Bi number
hfg d @m dhfg m Mhfg Bi ¼ p ðT Þ 1 lv sat 2 RT lv k @T kRT lv ðC sat ðT lv Þ C outlet Þ 00
00
ð34Þ
To understand the dependence of the critical Ma number to the Bi number described by Eq. (34), the critical Ma values are obtained for the channels with a wide range of geometry parameters, viz., different widths (d: 100 lm 1000 lm), aspect ratios (d/ll: 0.01 0.5) and length ratios (ll/lv: 2 10), as shown in Fig. 4c. For each case, the humidity at the outlet is lowered step by step. Once the instability occurs, the corresponding evaporation flux and interface temperature are recorded and the critical Ma value is obtained based on Eq. (29). Fig. 4c presents how critical Ma values vary with the Bi numbers. As Bi increases, the critical Ma increases. For small Bi numbers (Bi < 103), the critical Ma is almost a constant (180 < Ma < 186). 5.5. Flow patterns varying with Ma values
the Ma number is below the threshold. Once the Ma exceeds the threshold, the temperature peak on the meniscus moves to one of the four corners (Fig. 5a, point (a) to point (b)). The Marangoni convection is from the corner, radially outward. When Ma further increases, the flow pattern keeps varying. The temperature peak shifts along the diagonal, from the corner towards the center, point (b) toward (c). If the Ma number keeps increasing, a second instability would occur and the location of the temperature peak will leave the diagonal line and gets settled at a center line, point (c)–(d). The secondary critical Ma value varying with the Bi number is shown in Fig. 5b. The critical Ma number increases from 536 to 669 as the Bi number increases from 0.0003 to 0.1. 5.6. Influences of gravity The gravity does not significantly influence the instability onset or the Marangoni flow patterns if the capillary channel is vertical, 0.07
As presented in Figs. 2 and 3, there is a very weak evaporation compensating flow from the meniscus center to the edges when
0.04
b c
Gravity
d 0.00
(a)
a Diagonal line
28.0
Center
Secondary critical Ma
(a)
14.0
No Gravity
0.0
580
28.8 560
14.4 540
520
Gravity 0.0
10-4
10-2
10-3
0.1
(b)
Bi number
(b) Fig. 5. (a) The location of the temperature peak on the meniscus (the start point of the Marangoni flow) varies (from ‘a’ to ‘d’) as the Marangoni number increases; (b) The secondary critical Ma Number varies with the Bi number.
Fig. 6. The influences of the gravity (d = 800 lm; ll/lv = 3; d/ll = 0.27) (velocity is scaled by a/d). (a) Gravity brings forward the onset of the Marangoni instability (Ma = 105, Ra1 = 1.1, Bo1 = 0.010); (b) Gravity does not influence the flow patterns except the flow direction once the critical Ma value is reached (Ma = 622, Ra1 = 6.5, Bo1 = 0.011).
Table 4 Comparison between the Rayleigh number and Marangoni number. Channel width
Evaporation flux
Meniscus temperature
Ra1
Ra2
Ma
1/Bo1
1/Bo2
200 lm 200 lm 800 lm 800 lm
3.5 103 8.6 104 3.5 103 kg/ 8.6 104 kg/
21.1 °C 24.0 °C 21.2 °C 24.0 °C
0.035 0.009 8.693 2.285
0.52 0.13 32.6 8.57
53.4 13.3 831.3 218.5
1539 1539 95.6 95.6
102.6 102.6 25.5 25.5
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
no matter upward or downward. Since the gravity is perpendicular to the meniscus, the influence of the gravity manifests itself as the Rayleigh–Bénard instability. As shown in Table 4, both the Rayleigh number (Ra1 and Ra2) and the Bo number (Bo = Ra/Ma) are very small because of the small geometry of the capillary channels. Therefore, Rayleigh–Bénard instability does not occur in this case and the Marangoni effect is dominating. If the capillary channel is horizontal and the gravity direction is along one edge of the channel as marked in Fig. 6, the influence of the gravity is observable. The evaporation results in a cooling effect on the meniscus, which creates downward buoyancy force along the meniscus. As a result, the onset of Bénard–Marangoni instability is brought forward significantly. An example is shown in Fig. 6a, the Marangoni flow occurs when the Ma value is only about 105. More specific work is needed in the future to make a more comprehensive understanding of the threshold of Bénard–Marangoni instability in such conditions. The direction of the Marangoni flow pattern is also influenced by the gravity. As shown in Fig. 6b, the direction of the Marangoni flow is random when the gravity is zero, i.e., the temperature peak could locate at any corner of the meniscus. However, the Marangoni flow flows only downward on the meniscus if the gravity is imposed. 6. Conclusions Bénard–Marangoni instabilities on evaporating menisci in adiabatic capillary-channels are studied. On a flat, initially isothermal meniscus, even though the tangential temperature gradient on the meniscus is initially zero, the perpendicular temperature gradient can still arouse the instability leading to Marangoni flows when the evaporation intensity or the channel width increases to a certain critical value. The critical Ma values for the instability are obtained, ranging from 180 to 220 as the Bi number increases from 1 106 to 0.1. The flow pattern on the meniscus is varying as the Ma value increases. When Ma value reaches 600, a secondary instability occurs and the location of the temperatur peak moves from the diagonal lines to the center lines. The gravity does not have significant impact on the instability when its direction is perpendicular to the meniscus. However, the gravity along the meniscus tangential direction brings forward the threshold of the instability significantly, and influences the direction of the flow pattern. The present work is a preliminary step to a comprehensive understanding of Marangoni instability in capillary structures. The temperature gradient that is perpendicular to the meniscus should be paid attention to in the related applications. Acknowledgments This work was supported by the Natural Science Foundation of China (Grant No. 51276003) and Common Development Fund of Beijing. References [1] H.K. Dhavaleswarapu, J.Y. Murthy, S.V. Garimella, Numerical investigation of an evaporating meniscus in a channel, Int. J. Heat Mass Trans. 55 (2012) 915– 924. [2] Y.R. Li, H.R. Zhang, C.M. Wu, J.L. Xu, Effect of vertical heat transfer on thermocapillary convection in an open shallow rectangular cavity, Heat Mass Trans. 48 (2012) 241–251. [3] S.C. Zhao, Q.S. Liu, R. Liu, H. Nguyen-Thi, B. Billia, Thermal effects on Rayleigh– Marangoni–Bénard instability in a system of superposed fluid and porous layers, Int. J. Heat Mass Trans. 53 (2010) 2951–2954. [4] D. Melnikov, T. Takakusagi, V. Shevtsova, Enhancement of evaporation in presence of induced thermocapillary convection in a non-isothermal liquid bridge, Microgravity Sci. Tech. 2012 (online). [5] A.S. Basu, Y.B. Gianchandani, Virtual microfluidic traps, filters, channels and pumps using marangoni flows, J. Micromech. Microeng. 18 (2008) 115031.
247
[6] Z. Pan, H. Wang, Z. Yang, Marangoni bifurcation flow in a microchannel Tjunction and its micropumping effect: a computational study, Chin. Phys. Lett. 29 (2012) 074702. [7] R. Bhardwaj, X. Fang, D. Attinger, Pattern formation during the evaporation of a colloidal nanoliter drop: a numerical and experimental study, New J. Phys. 11 (2009) 075020. [8] H. Hu, R.G. Larson, Marangoni effect reverses coffee-ring depositions, J. Phys. Chem. B 110 (2005) 7090–7094. [9] J.R.A. Pearson, On convection cells induced by surface tension, J. Fluid Mech. 4 (1958) 489–500. [10] H. Bénard, Les tourbillons cellulaires dans une nappe liquide transportant de la chaleur par convection en regime permenant, Ann. Chim. Phys. 23 (1901) 62– 144. [11] J.A. Maroto, V. Perez-Munuzuri, M.S. Romero-Cano, Introductory analysis of Bénard–Marangoni convection, Euro. J. Phys. 28 (2007) 311. [12] Z.M. Tang, K. Li, W.R. Hu, Influence of free surface curvature of a liquid layer on the critical Marangoni convection, Int. J. Heat Mass Trans. 51 (2008) 5102– 5107. [13] H.S. Das, B.D. MacDonald, C.A. Ward, Stability of evaporating water when heated through the vapor and the liquid phases, Phys. Rev. E 81 (2010) 036318. [14] J.W. Peterson, G.F. Carey, Parallel adaptive solution of coupled Rayleigh– Bénard–Marangoni problems with the Navier-slip, Int. J. Numer. Methods Fluid 66 (2011) 428–451. [15] C. Ma, D. Bothe, Direct numerical simulation of thermocapillary flow based on the volume of fluid method, Int. J. Multiphase Flow 37 (2011) 1045–1058. [16] F. Doumenc, T. Boeck, B. Guerrier, M. Rossi, Transient Rayleigh–Bénard– Marangoni convection due to evaporation: a linear non-normal stability analysis, J. Fluid Mech. 648 (2010) 521–539. [17] G. Toussaint, H. Bodiguel, F. Doumenc, B. Guerrier, C. Allain, Experimental characterization of buoyancy- and surface tension-driven convection during the drying of a polymer solution, Int. J. Heat Mass Transfer 51 (2008) 4228– 4237. [18] R. Touihri, A. EI Gallaf, D. Henry, H. Ben Hadid, Instabilities in a cylindrical cavity heated from below with a free surface I effect of Biot and Marangoni numbers, Phys. Rev. E 84 (2011) 056302. [19] A.M. Benselama, S. Harmand, K. Sefiane, Thermocapillary effects on steadily evaporating contact line: a perturbative local analysis, Phys. Fluids 24 (2012) 072105. [20] C. Buffone, K. Sefiane, J.R.E. Christy, Experimental investigation of self-induced thermocapillary convection for an evaporating meniscus in capillary tubes using micro-particle image velocimetry, Phys. Fluids 17 (2005) 052104. [21] S.S. Panchamgam, A. Chatterjee, J.L. Plawsky, P.C. Wyner Jr., Comprehensive experimental and theoretical study of fluid flow and heat transfer in a microscopic evaporating meniscus in a miniature heat exchanger, Int. J. Heat Mass Trans. 51 (2008) 5368–5379. [22] X. Song, D.S. Nobes, Experimental investigation of evaporation-induced convection in water using laser based measurement techniques, Exp. Therm. Fluid Sci. 35 (2011) 910–919. [23] H. Wang, Z. Pan, S.V. Garimella, Numerical investigation of heat and mass transfer from an evaporating meniscus in a heated open groove, Int. J. Heat Mass Transfer 54 (2011) 3015–3023. [24] H.K. Dhavaleswarapu, P. Chamarthy, S.V. Garimella, J.R. Murthy, Experimental investigation of steady buoyant-thermocapillary convection near an evaporating meniscus, Phys. Fluids 19 (2007) 082103. [25] P. Chamarthy, H.K. Dhavaleswarapu, S.V. Garimella, J.Y. Murthy, S.T. Wereley, Visualization of convection patterns near an evaporating meniscus using lPIV, Exp. Fluids 44 (2008) 431–438. [26] H. Wang, J.Y. Murthy, S.V. Garimella, Transport from a volatile meniscus inside an open microtube, Int. J. Heat Mass Transfer 51 (2008) 3007–3017. [27] B. Lan, Y.R. Li, D.F. Ruan, Numerical simulation of thermocapillary flow induced by non-uniform evaporation on the meniscus in capillary tubes, Microgravity Sci. Technol. 23 (2011) 35–42. [28] C.A. Ward, F. Duan, Turbulent transition of thermocapillary flow induced by water evaporation, Phys. Rev. E 69 (2004) 056308. [29] F. Duan, C.A. Ward, Investigation of local evaporation flux and vapor-phase pressure at an evaporative droplet interface, Langmuir 25 (2009) 7424–7431. [30] C. Buffone, K. Sefiane, Controling evaporative thermocapillary convection using external heating: an experimental investigation, Exp. Therm. Fluid Sci. 32 (2008) 1287–1300. [31] C. Buffone, K. Sefiane, W. Easson, Marangoni-driven instabilities of an evaporation liquid-vapor interface, Phys. Rev. E 71 (2005) 056302. [32] Z. Pan, H. Wang, Symmetry-to-asymmetry transition of Marangoni flow at a convex volatizing meniscus, Microfluid. Nanofluid. 9 (2010) 657–669. [33] Z. Pan, F. Wang, H. Wang, Instability of marangoni toroidal convection in a micro channel and its relevance with the flowing direction, Microfluid. Nanofluid. 11 (2011) 327–338. [34] D. Gazzola, E. Franchi Scarselli, R. Guerrieri, 3D visualization of convection patterns in lab-on-chip with open microfluidic outlet, Microfluid. Nanofluid. 7 (2009) 659–668. [35] S.H. Davis, Thermocapillary instabilities, Annu. Rev. Fluid Mech. 19 (1987) 403–435. [36] H. Hertz, Ueber die verdungstung der flussigkeiten, inbesondere des quecksilbers, im luftleeren raume, Ann. Phys. Chem. 17 (1882) 177–200. [37] M. Knudsen, Die maximale verdampfungsgeschwindigkeit des quecksilbers, Ann. Phys. Chem. 47 (1915) 697–708.
248
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
[38] R.W. Schrage, A Theoretical Study of Interface Mass Transfer, Columbia University Press, New York, USA, 1953. [39] R. Marek, J. Straub, Analysis of the evaporation coefficient and the condensation coefficient of water, Int. J. Heat Mass Transfer 44 (2001) 39–53. [40] B. Paul, Complication of evaporation coefficients, ARS J. 32 (1962) 1321–1328. [41] K.C.D. Hickman, Maximum evaporation coefficient of water, Ind. Eng. Chem. Res. 46 (1954) 1442–1446. [42] V.K. Badam, V. Kumar, F. Durst, K. Danov, Experimental and theoretical investigations on interfacial temperature jumps during evaporation, Exp. Therm. Fluid Sci. 32 (2007) 276–292. [43] S.M. Yang, W.Q. Tao, Heat Transfer, Higher Education Press, Beijing, China, 2006.
[44] S.R. Mathur, J.Y. Murthy, Pressure boundary conditions for incompressible flow using unstructured meshes, Numer. Heat Transfer B 32 (1997) 283–298. [45] G.J. Dunn, S.K. Wilson, B.R. Duffy, S. David, K. Sefiane, The strong influence of substrate conductivity on droplet evaporation, J. Fluid Mech. 623 (2009) 329– 351. [46] Y.A. Cengel, Heat Transfer: A Practical Approach, McGraw-Hill Press, New York, USA, 2003. [47] S. Semenov, V.M. Starov, R.G. Rubio, M.G. Velarde, Computer simulations of evaporation of pinned sessile droplets: influence of kinetic effects, Langmuir 28 (2012) 15203–15211.
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Bénard–Marangoni instability on evaporating menisci in capillary channels q Zhenhai Pan a, Hao Wang a,b,⇑ a b
Laboratory of Heat and Mass Transport at Micro-Nano Scale, College of Engineering, Peking University, Beijing 100871, China Beijing Key Laboratory for Solid Waste Utilization and Management, College of Engineering, Peking University, Beijing 100871, China
a r t i c l e
i n f o
Article history: Received 28 November 2012 Received in revised form 30 March 2013 Accepted 31 March 2013
Keywords: Meniscus Capillary channel Evaporation Bénard–Marangoni Instability
a b s t r a c t Despite numerous studies about Marangoni instabilities, there are very few of them conducted for the menisci in microscale channels or tubes. In this paper, the onset of Marangoni flow on an initially isothermal meniscus in a capillary channel is studied. The temperature gradient perpendicular to the meniscus is found arousing a Bénard–Marangoni instability when the channel size or the temperature gradient beyond certain values. The threshold Marangoni Number for the instability is obtained, which turns out to be dependent on the Biot Number. The evolution of the meniscus flow patterns with increasing Marangoni Number is studied. The influences of gravity are discussed. The present work is a preliminary step to a comprehensive understanding of Marangoni instability in capillary structures. The temperature gradient that is perpendicular to the meniscus should be paid attention to in the related applications. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Marangoni flow has been proved of importance in various applications like heat transfer [1–4], microfluidic control [5,6], particle deposition [7,8], etc. In Pearson’s pioneer analysis study [9] on Bénard’s experiments [10], when a liquid layer with an isothermal interface is heated from below, small perturbation gets enlarged and surface-tension-induced flow presents. Such a problem, where the Marangoni flow is induced by the temperature gradient perpendicular to the meniscus, is usually referred to as Bénard–Marangoni instability and has been studied extensively as a classical problem for decades [11–15]. Up to today, the study on Bénard–Marangoni instability is still going on for different situations and applications. For example, Doumenc et al. [16] studied the transition of Bénard– Marangoni flow on an infinite evaporating liquid layer with an adiabatic bottom. The numerical results agreed very well with previous Toussaint et al.’s experiments [17]. Touihri et al. [18] investigated the threshold and nonlinear variation of the Marangoni flow in a cylinder. Factors including Biot number, geometry aspect ratio and buoyancy force are considered. Benselama et al. [19] theoretically studied the Marangoni instability on an evaporating thin film
q Submitted for publication in International Journal of Heat and Mass Transfer, Nov. 2012. ⇑ Corresponding author at: Beijing Key Laboratory for Solid Waste Utilization and Management, College of Engineering, Peking University, Beijing 100871, China. Tel.: +86 10 8252 9010. E-mail address: [email protected] (H. Wang).
0017-9310/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.03.082
extending from a meniscus. The capillary pressure and disjoining pressure were found playing important roles in the hydrodynamic instabilities. With the fast development of microfluidics and MEMS applications, Marangoni flows in microscale/capillary structures have gained increasing attentions. Particularly, studies have been conducted for the evaporating menisci in capillary channels or tubes [20–27]. For most of the studies, the meniscus is initially non-isothermal, and thus Marangoni convection is directly induced by the temperature gradients along the meniscus. Usually, such Marangoni convection is also referred to as Thermalcapillary flow. For example, in Dhavaleswarapu et al.’s [24] and Chamarthy et al.’s [25] experiments on concave menisci, the evaporations were stronger at the meniscus edges. The meniscus edges were thus cooler than the centers, and the Marangoni convections were observed flowing from the meniscus centers to the edges. Similar results were found in Wang et al.’s [26] and Lan et al.’s [27] simulation studies. In Ward and Duan’s [28,29] and Buffone and Sefiane’s [30] experiments, the hot walls near the meniscus warmed the meniscus edges while the evaporation cooled the meniscus centers, inducing inward Marangoni convections. There are very few studies about the Marangoni instability in the capillary structures, even fewer for Bénard–Marangoni instability. Buffone et al. [20,31] noticed that the Marangoni flow was not stable in a horizontal tube and attributed it to the meniscus temperature difference near the contact line. Recently, Pan and Wang numerically [32] and experimentally [33] found that on a convex meniscus at the mouth of a microchannel, the initial
240
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
Nomenclature A Bo Bi C cp Cr D d hfg J k l Ma M mnet m00net ~ n p Pr r R Ra Sm Sh T T⁄ V q
area (m2) Bond number Biot number vapor molar concentration (mol/m3) heat capacity (J/kg K) Crispation number diffusion coefficient in air (m2/s) the channel width (m) latent heat of evaporation (J/kg) evaporating mole flux (mol/m2 s) thermal conductivity (W/ m K) the channel length (m) Marangoni number molecular weight (kg/kmol) mass flow rate (kg/s) mass flux (kg/m2 s) interface normal coordinate pressure (N/m2) Prandtl number radius (m) universal gas constant (J/mol K) Rayleigh number mass source (kg/m3 s) energy source (W/m3) temperature (K) non-dimensional time fluid velocity (m/s) heat flux (W/m2)
Greeks symbols a thermal diffusivity (m2/s) b thermal expansion coefficient (1/K) d thickness (m) h contact angle m kinematic viscosity (m2/s) l dynamic viscosity (N s/m2) q density (kg/m3) r surface tension coefficient (N/m) r interface accommodation coefficient s shear stress (N/m2) Subscripts am ambient Bous Boussinesq e evaporation Equ equilibrium g gas (vapor/air mixture) l liquid lv interface outlet outlet ref reference sat saturated v vapor x x component y y component z z component
symmetrical Marangoni flow is not stable and tends to evolve into an asymmetrical single vortex flow. The results also explained the abnormal flow pattern observed in Gazolla et al.’s experiment [34] where unexpected asymmetrical flow always presented under symmetrical boundary conditions. In the present paper we first study the onset of Marangoni flow on a flat, initially isothermal meniscus inside a capillary channel. It is found the perpendicular-to-meniscus temperature gradient (from the meniscus to the warm bulk) can induce Bénard–Marangoni instability as long as the channel size or the evaporation intensity beyond certain thresholds. The critical Ma values for the instability are obtained which turn out to be sensitive to the Biot number. Then the flow variations with Ma are presented and the influences of the gravity are discussed. It is suggested that the Bénard–Marangoni instability induced by the perpendicular temperature gradient
should be paid attention to in related microfluidic and MEMS designing. 2. Problem description A capillary channel with square cross section is shown in Fig. 1. Water at room temperature (298.15 K) is supplied from the channel inlet. A meniscus is formed in the channel. The channel outlet is open to the ambient air and water is volatizing i.e. evaporating on the meniscus. The meniscus temperature is lower than the bulk due to the latent heat dissipated by the evaporation. To study Bénard–Marangoni instability the meniscus has to be initially isothermal. In the present work to remove the geometry factors and catch the essentials, the meniscus shape is assumed to be flat. The evaporation on the meniscus is uniform, resulting
Meniscus
Outlet
Vapor transport
Inlet
Compensating flow
y
θ
z x
meniscus
Outlet Gas domain
lv
d
Inlet
Liquid domain
ll
Fig. 1. Schematic diagram of the open channel with an evaporating meniscus.
241
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
shear stress. Hence the Marangoni flow is produced. The shear stress at the meniscus is calculated as
Table 1 Fluid properties. Properties 3
Density (kg/m ) Thermal conductivity (W/m K) Thermal capacity (J/kg) Viscosity (kg/m s) Vapor molecular weight M (kg/mol) Thermal expansion coefficient (1/K) dr/dT (N/m K) Diffusion coefficient (m2/s) Accommodation coefficient for evaporation Saturated vapor pressure (Pa) Latent heat (J/kg)
Water
Air
997.05 at 298.15 K 0.6 4182 0.001003 18 103 0.000251 0.00015 2.55 105at 300 K 0.03
1.182 0.0242 1006.43 1.789 105 29 103 0.003355 / /
2334.6 at 293.15 K 2.26 106
/
sl jlv ¼
3.2.2. Meniscus deformation and profiles On a flat meniscus, the meniscus deformation is evaluated with Crispation number as suggested by Davis [35] and Doumenc et al. [16]. Crispation number is defined as
Cr ¼
in a uniform temperature distribution (details referring to the simulation results shown in Section 5). In practice, the flat meniscus could be true when the contact angle is nearly 90°, or the meniscus edge is pinned. The channel outlet is open to the ambient air and water is volatizing. The vapor transports through the channel outlet into the ambient as illustrated in Fig. 1. The vapor fraction on the channel outlet is assumed to be uniform and is determined by the relative humidity of the ambient. The channel width is d. The distances from the meniscus to the inlet and the outlet are ll and lv respectively. The ambient temperature Tam is 298.15 K. Fluid properties are listed in Table 1.
la rd
3.2.3. Evaporation Kinetic theory of gases has provided the molecular basis for understanding the evaporation for over a century. Hertz [36] and Knudsen [37] studied the evaporation of liquid mercury into vacuum and proposed Hertz–Knudsen theory. Schrage [38] further modified the theory by including the net macroscopic velocity of vapor and proposed Hertz–Knudsen–Schrage equation, that is, the evaporated mass flux across the interface at one moment is calculated as
m00net ¼
3.1. Liquid domain and gas domain
! qr V ¼ Sm !
!
ð1Þ !
q V r V ¼ rp þ r ðlr V Þ þ qBous~ g !
ð2Þ
qcp V rT ¼ kr2 T þ Sh
ð3Þ
qBous ¼ qref qref bðT T ref Þ
ð4Þ
3.2. Liquid–gas interface 3.2.1. Marangoni effect The temperature gradient along the meniscus would induce surface tension gradient. The tangential force due to the surface tension gradient on the meniscus should be balanced by the liquid
! 1=2 pv equ ðT lv Þ 2r M pv 2pR 2r T 1=2 T 1=2 v lv
ð7Þ
Assume that the interfacial vapor temperature Tv is equal to the interfacial liquid temperature Tlv, [23,32] leading to
m00net
Laminar flow of a Newtonian fluid is assumed in the liquid domain based on the small scale of the system. Boussinesq approximation is employed for modeling the natural convection such that the density changes with temperature only in the buoyancy term. The laminar and Newtonian fluid with Boussinesq approximation is also adapted for the air/vapor mixture since the pressure and temperature difference is small in the calculated gas domain. The continuity, momentum and energy equations and Boussinesq approach are given respectively as
ð6Þ
For a small Cr number, the surface deformation is negligible comparing with the channel width d. In the present study with temperature T 298 K, we have Cr < 2 105 with surface tension r 7.2 102 N/m, liquid viscosity l 1 103 kg/(m s), thermal diffusivity a 1.4 107 m2/s, channel width d > 1 104 m. Therefore, the deformation of the interface is neglected in the present simulations considering such small Cr values.
3. Mathematical modeling The heat and mass transport in the capillary channel includes the evaporation at the meniscus, vapor diffusion and convection in the gas domain, Marangoni and buoyancy flow in the liquid domain. A mathematical model is described in the following sections and all the factors will be solved together.
ð5Þ
where @T/@s is the temperature gradient along the meniscus.
/
/
dr @T dT @s
1=2 pv 2r M ¼ 2pR 2r
equ ðT lv Þ 1=2 T lv
pv
! ð8Þ
The capillary pressure caused by the meniscus curvature could affect the the equilibrium vapor pressure Pv_equ according to Kelvin equation, which is considerable when the curvature radius is smaller than several microns [46,47]. In our case, however, the curvature radius is 0.1 mm so the equilibrium vapor pressure is approximated as the corresponding saturation pressure.
pv
equ ðT lv Þ
psat ðT lv Þ ¼ psat;ref exp
Mhfg 1 1 T sat;ref T lv R
ð9Þ
^ is one source of uncertainty. As The accommodation coefficient r noted in [39], the evaporation and condensation coefficients of water that are obtained theoretically or experimentally scatter over a rather large range of more than two decades. Paul [40] estimated the value of the accommodation coefficient to be in the range of 0.02–0.04. Hickman [41] found that very high evaporation coefficients (close to unity) could only be obtained with moving liquids in a vacuum where the surface is continually refreshed by the moving stream and the evaporating molecules are removed immediately by the vacuum. Marek and Straub [39] commented that even small contaminations of the surface significantly reduce the interfacial mass transfer. The lowest value of the accommodation coefficient given in their review was 0.002 for water in a glass vessel. Badam et al. [42] estimated the value of the accommodation coefficient to be in the range of 0.028–0.1. In the present work, ^ is assumed to be a small value the accommodation coefficient r of 0.03 based on the previous works [32,40,42]. The possible
242
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
uncertainty induced by this coefficient value would be discussed in Section 4.
two layers of mesh cells adjacent to the interface as described specifically in [26]. The mass and energy source terms in the cells along the interface in both liquid and gas domains are calculated as
3.3. Coupling of vapor transport and evaporation
Sm;g ¼ The diffusion of vapor from the meniscus into the ambient air and its coupling to the evaporation has been described specifically in our previous work [23]. Briefly, the governing equation for the vapor transport in air is given as
@C v þ~ v rC v ¼ r ðDrC v Þ @t
ð10Þ
At the meniscus, the mass flux dissipated by the vapor transport can be obtained by
@C v m00v ¼ MJ ¼ M D þ v n C v @n lv
ð11Þ
The first term on the right side is the vapor transport due to mass diffusion and the second term is due to mass convection. This convection is due to the single direction of diffusion and known as Stefan flow [43]. After a series of derivations the evaporation at meniscus is coupled to the vapor transport in gas domain [23]:
1=2 MD @C v 2r M 1 ¼ ðpv 1=2 2pR 1 C v =C total @n lv 2 r T
equ ðT lv Þ
pv Þ
ð12Þ
lv
3.4. Boundary conditions The channel walls are treated as no-slip and adiabatic boundaries. A given stagnation-pressure inlet for the liquid and a given static pressure outlet [44] for the vapor are set respectively at the two ends of the tube. The temperature at both ends are assumed as a room temperature of 298.15 K. The vapor concentration is set based on the local vapor pressure.
p Cv ¼ v RT
ð13Þ
m00 Af ;g m00 Af ;l ; Sm;l ¼ V cell;g V cell;l
ð14Þ
Sh;l ¼ Sm;l hfg þ Sm;l cp;l ðT cell T ref Þ
ð15Þ
Sh;g ¼ Sm;g cp;g ðT cell T ref Þ
ð16Þ
Here, Af is the interface area of the cell. Vcell and Tcell are the volume and the temperature of the cell. The first term on the right side of Eq. (15) represents the latent heat consumed by the evaporation. The second term on the right side of Eq. (15) and the only term on the right side of Eq. (16) represent the sensible heat corresponding to the evaporative mass. Tref is the reference temperature used in ANSYS 12, which is 298.15 K. Similar mathematical modeling and numerical methodology have been validated in authors’ previous studies on evaporating menisci in microtube [26], V-groove [23], and at the mouth of capillary channel [32,33]. The obtained flow patterns and evaporation rates were all highly consistent with the corresponding experiments. For flat menisci, the mesh setup is illustrated in Fig. 1. Hexahedral elements are created for liquid and gas domains. The mesh independence test is performed. Four sets of meshes are checked under two different operating conditions as shown in Tables 2a and 2b. Mesh 3 is employed in the simulations. ^ is As mentioned in Section 3.2, the accommodation coefficient r a source of uncertainty. A series of tests are conducted to estimate the uncertainty. It is seen from Table 3 that the evaporation intensity is not sensitive to the choice of the coefficient values (the largest deviation is 5% when the coefficient value changes from 0.002 to 1). This result is reasonable because the primary resistance for water evaporation to open air relies on the vapor diffusion in the gas domain rather than the interfacial evaporation resistance [45,47], i.e. even when the accommodation coefficient is 0.002 the evaporation is still efficient enough compared to the vapor diffusion.
4. Numerical treatment 5. Results and discussions The numerical solution is obtained using the pressure based finite volume scheme. Software package ANSYS 12.0 (FLUENT SOLVER) enhanced with user-developed codes are employed. The SIMPLE algorithm is used for pressure–velocity coupling. To calculate the heat and mass exchange across the interface, it is assumed that the mass transport across the interface occurs only within the
5.1. Occurrence of instability when increasing channel size On the flat meniscus in the adiabatic channel as shown in Fig. 1, the evaporation distribution is uniform and consequently the temperature is uniform on the meniscus. However, Marangoni flow
Table 2a Mesh-independence studies under a large Ma value (d = 1000 lm, ll = 3 mm, lv = 1 mm, outlet humidity: 0%). Mesh setup (Nx Ny Nz)
Mesh 1 13 13 67
Mesh 2 19 19 95
Mesh 3 27 27 133
Mesh 4 40 40 140
Evaporation rate (kg/s) Deviation from mesh 4 Marangoni number Deviation from mesh 4
4.4424 1010 0.43% 1504.5 1.64%
4.4359 1010 0.28% 1513.0 1.09%
4.4289 1010 0.12% 1522.3 0.48%
4.4236 1010 0 1529.6 0
Table 2b Mesh-independence studies under a small Ma value (d = 600 lm, ll = 3 mm, lv = 1 mm, outlet humidity: 0%). Mesh setup (Nx Ny Nz)
Mesh 1 (13 13 67)
Mesh 2 (19 19 95)
Mesh 3 (27 27 133)
Mesh 4 (40 40 140)
Evaporation rate (kg/s) Deviation from mesh 4 Marangoni number Deviation from mesh 4
9.3867 1011 0.27% 356.43 1.02%
9.3742 1011 0.14% 358.16 0.53%
9.3661 1011 0.05% 359.29 0.22%
9.3610 1011 0 360.08 0
243
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248 Table 3 Tests of the accommodation coefficient (ll = 3 mm, lv = 1 mm, outlet humidity: 0%). Channel width (lm) 400 400 400 800 800 800 1000 1000 1000
Accommodation coefficient
Evaporation rate (kg/s) 11
0.002 0.03 1 0.002 0.03 1 0.002 0.03 1
6.62 10 6.91 1011 6.93 1011 2.69 1010 2.80 1010 2.81 1010 4.28 1010 4.43 1010 4.45 1010
occurs when the channel width d is greater than a critical value. Five channels with different d values, from 200 lm to 1000 lm, are simulated and the results of the temperature and flow fields are shown in Fig. 2. The ambient humidity at the channel outlet is fixed at 40%. The resulting evaporation flux is almost the same, 0.261 ± 0.002 g/(m2 s). Based on Fig. 2 when the channel width is below 600 lm, the meniscus is isothermal (the maximum temperature difference across the meniscus is about 0.0001 K). An extremely weak (velocity below 1 lm/s) flow is observed flowing from the meniscus center to the edges. This very weak flow arises due to the following reason: firstly, there is a compensating flow through the channel, from the inlet to the meniscus, to compensate the mass loss due to the evaporation; secondly, due to the friction of the channel walls, the compensating flow has higher velocity in the central than the edge of the channel, which brings more liquid to the meniscus center than to the meniscus edges.
295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272 295.272
Interfacial temperature (K)
Marangoni number
293.49 293.27 293.22 293.69 293.50 293.47 293.85 293.70 293.68
254.9 266.9 269.5 977.2 1018.8 1024.9 1471.0 1522.3 1529.3
However, when the channel width increases to 600 lm, a strong convection with a totally different pattern presents. The liquid flows from a corner to another on the meniscus with a max velocity up to 3900 lm/s, which is thousands times larger than the evaporation compensating flow. Corresondingly, the maximum temperature difference along the meniscus reaches 0.1 K, which is one thousand times larger than that of the 400 lm channel. For even bigger channels, i.e., 800 and 1000 lm in width, similar strong flows are obtained. 5.2. Occurrence of instability when increasing evaporation intensity If the channel width is fixed at 600 lm but the evaporation intensity is increased by reducing the ambient humidity at the channel outlet, the results are shown in Fig. 3. When the humidity is 80%, the evaporation on the meniscus is weak, with average evaporative mass flux of 8.8 105 kg/(m2 s). A uniform tempera-
0.001
0.000
0.000
d=200µm; Ma=39 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201 295.201
0.002
0.001
0.000
d=400µm; Ma=161 295.380 295.380 295.362 295.343 295.325 295.307 295.288 295.262 295.270 295.252 295.233 295.215 295.197 295.178 295.160 295.160
5.048
2.524
0.000
d=800 µm; Ma=622
295.300 295.300 295.287 295.273 295.260 295.247 295.233 295.220 295.220 295.207 295.193 295.180 295.167 295.153 295.140 295.140
3.907
1.963
0.019
d=600 µm; Ma=358 295.480 295.480 295.457 295.433 295.410 295.387 295.363 295.340 295.340 295.317 295.293 295.270 295.247 295.223 295.200 295.200
5.661
2.945
0.228
d=1000 µm; Ma=952
Fig. 2. Meniscus temperature (K) (left) and flow pattern (mm/s) (right) as channel width d increases (humidity = 40%, ll = 3 mm, ll/lv = 3).
244
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176 297.176
0.001
0.001
0.000
Humidity: 80%; Ma=120 296.210 296.210 296.205 296.200 296.195 296.190 296.185 296.180 296.180 296.175 296.170 296.165 296.160 296.155 296.150 296.150
295.300 295.300
1.628
295.287 295.273 295.260 295.247 295.233 295.220 295.220 295.207 295.193 295.180 295.167 295.153 295.140 295.140
0.815
0.000
Humidity: 60%; Ma=242 294.402 294.402 294.382 294.362 294.342 294.322 294.302 294.281 294.281 294.261 294.241 294.221 294.201 294.181 294.161 294.161
3.907
1.963
0.019
Humidity: 40%; Ma=358 293.485 293.485
5.273
293.462 293.438 293.414 293.391 293.368 293.344 293.344 293.320 293.297 293.274 293.250 293.226 293.203 293.203
2.666
0.038
6.215
3.142
0.069
Humidity: 0%; Ma=588
Humidity: 20%; Ma=473
Fig. 3. Meniscus temperature (K) (left) and flow pattern (mm/s) (right) as the humidity at the outlet decreases and the evaporation increases (ll = 3 mm, ll/lv = 3, d/ll = 0.2).
ture is obtained on the meniscus (the maximum temperature difference is about 0.0001 K). However, when the humidity decreases to 60% and the evaporation on the meniscus increases to 1.77 104 kg/(m2 s), the temperature difference jumps up to 0.1 K and a strong Marangoni convection occurs.
!
! V rT ¼ ar2 T
ð19Þ
At the meniscus where z = 0, assuming the latent heat of the evaporation is supplied from the liquid side, we have
k
!
rV ¼ 0
ð17Þ
!
ð18Þ
5.3. Bénard–Marangoni instability in capillary channels Above Marangoni convections occurring on the initially isothermal menisci could be categorized as one type of Bénard–Marangoni instability. The mechanism of the instability could be described by considering a perturbation occurs on an initially isothermal meniscus in a microchannel, as illustrated in Fig. 4a. As the meniscus keeps being cooled by the evaporation, the meniscus temperature T1 at the upriver of the perturbation flow must be higher than the downriver temperature T2. This temperature difference between T1 and T2 would then induce a Marangoni flow along the meniscus, whose flow direction is same to the perturbation. A positive feedback is then established. However, at the same time, the thermal diffusion tends to eliminate the temperature difference between T1 and T2 while the liquid viscosity tends to eliminate the induced Marangoni flow. If the former positive feedback effect is stronger than the thermal diffusion and liquid viscosity effect, i.e. Ma number is high enough, the Bénard–Marangoni instability occurs. As discussed in Section 5.1, the evaporation compensating flow is very weak compared to the Marangoni flow. By neglecting the evaporation compensating flow, the governing equations and the boundary conditions for the liquid domain are
!
q V r V ¼ rp þ lr2 V qgbðT T ref ÞÞe!g
l
@T ¼ hfg m00 ðTÞ @z
ð20Þ
@V x dr @T @V y dr @T ¼l ¼0 @z dT @x @z dT @y
ð21Þ
If choose characteristic length d, velocity a/d, pressure qa2/d2, and temperature DTd/ll in which DT = Tam Tlv,ave is the temperature drop from the inlet to the meniscus, we have dimensionless equations as
!
r V ¼ 0
ð22Þ
! ! ! ! V r V ¼ rp þ Prr2 V þ PrðRa1 T Ra2 Þeg
ð23Þ
! V rT ¼ r2 T
ð24Þ
in which the non-dimensional temperature T⁄ = (T Tlv,ave)ll/DTd. Tlv,ave is the average temperature on the meniscus. For the boundary conditions, at the meniscus, Eq. (20) can be linearized as 00
k
@T dm ¼ hfg m00 ðT lv ;av e Þ þ ðT lv ;av e ÞðT T lv ;av e Þ @z dT
ð25Þ
Meanwhile, according to the conduction through the liquid domain in channel
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
Evaporation cooling flux
296.4430 296.4430
0.007
296.4430 296.4430 296.4430 296.4430 296.4429 296.4429 296.4429 296.4429 296.4429 296.4429 296.4428 296.4428 296.4428 296.4428 296.4428
Meniscus
T1
T2 Channel wall
Channel wall Perturbation flow
Heat flux
245
0.004
0.000
Warm bulk
(a)
(b) d=1000 μm d=800 μm
210
ll/ lv=2, ll=2mm
d=600 μm
Critical Ma
ll/ lv=3, ll=3mm ll/ lv=6,ll=6mm
200
d=400 μm
ll/ lv=10, ll=10mm
d=200 μm
190
180 -6 10
d=100 μm
10-5
10-4
10-3
10-2
0.1
Bi number
(c) Fig. 4. The onset of the Marangoni instability: (a) the illustration of the mechanism; (b) an example showing the critical state (d = 600 lm; ll = 3 mm; Mac = 210); (c) the critical Ma values varying with Bi numbers.
hfg m00 ðT lv ;av e Þ k
DT ll
ð26Þ
with which Eq. (20) could finally be non-dimensionalized as
@T ¼ 1 þ BiT @z
ð27Þ
At the same time Eq. (21) is non dimensionalized as @V x @T @V y @T þ Ma ¼ þ Ma ¼ 0 @z @x @z @y
ð28Þ
In above dimensionless equations, Eqs. (22)–(24), (27) and (28), we have 4
gbd DT Ra1 ¼ ; amll 2
3
Ra2 ¼
dr d DT ; Ma ¼ dT lall
gbd DT
am
m ; Pr ¼ ; a 00
hfg d dm Bi ¼ ðT lv Þ k dT
ð29Þ
2
dr d hfg m00 dT la k
5.4. Influences of Biot number As discussed above, the critical Ma value for the onset of the Bénard–Marangoni instability should be dependent on the Bi number. Since the vapor transport is dominated by the vapor diffusion (as discussed in Section 4), the evaporating mass flux m00 could be approximated as
m00 MD
By combining with Eq. (26), Ma can also be written as
Ma ¼
distribution and interfacial flow pattern is significantly different from those before the Bénard–Marangoni instability occurs, the temperature gradient and the flow velocity are still very small (T 0.001 K; V 1 lm/s). This state is recogniazed as the critical state of the onset of the instability. One example is shown in Fig. 4b when d = 600 lm and Ma number comes to the threshold Mac = 210. Once the Ma number further increases and exceeds the critical value, the temperature and the velocity magnitude would increase sharply and establish strong flows as shown in Fig. 3.
ð30Þ
In order to obtain critical Ma values, the humidity at the outlet is lowered step by step (thus the evaporation intensity increases and Ma increases) for each case. As mentioned in Sections 5.1 and 5.2, the meniscus temperature is almost uniform before the Bénard–Marangoni instability occurs. Once the Ma value reaches the threshold, the temperature peak moves suddenly to one corner of the meniscus and a corresponding Marangoni flow is induced flowing from the location of the temperature peak to the other corners (Fig. 4b). At this critical state, even if the temperature
C lv C outlet lv
ð31Þ
Combining Eq. (13) with Eq. (31), we get
m00
MD pv ;lv pv ;outlet Rlv T lv T outlet
ð32Þ
Still because of the vapor transport is dominated by the vapor diffusion, the vapor pressure at the interface pv,lv is approximately equal to the saturation pressure at the interface psat(Tlv). Combining Eq. (32) with Eq. (9), we get
Mhfg @m00 MD 1 p ðT Þ 1 l v Rlv T 2lv sat @T lv RT lv
ð33Þ
246
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
Combining Eqs. (29), (31), and (33), we have Bi number
hfg d @m dhfg m Mhfg Bi ¼ p ðT Þ 1 lv sat 2 RT lv k @T kRT lv ðC sat ðT lv Þ C outlet Þ 00
00
ð34Þ
To understand the dependence of the critical Ma number to the Bi number described by Eq. (34), the critical Ma values are obtained for the channels with a wide range of geometry parameters, viz., different widths (d: 100 lm 1000 lm), aspect ratios (d/ll: 0.01 0.5) and length ratios (ll/lv: 2 10), as shown in Fig. 4c. For each case, the humidity at the outlet is lowered step by step. Once the instability occurs, the corresponding evaporation flux and interface temperature are recorded and the critical Ma value is obtained based on Eq. (29). Fig. 4c presents how critical Ma values vary with the Bi numbers. As Bi increases, the critical Ma increases. For small Bi numbers (Bi < 103), the critical Ma is almost a constant (180 < Ma < 186). 5.5. Flow patterns varying with Ma values
the Ma number is below the threshold. Once the Ma exceeds the threshold, the temperature peak on the meniscus moves to one of the four corners (Fig. 5a, point (a) to point (b)). The Marangoni convection is from the corner, radially outward. When Ma further increases, the flow pattern keeps varying. The temperature peak shifts along the diagonal, from the corner towards the center, point (b) toward (c). If the Ma number keeps increasing, a second instability would occur and the location of the temperature peak will leave the diagonal line and gets settled at a center line, point (c)–(d). The secondary critical Ma value varying with the Bi number is shown in Fig. 5b. The critical Ma number increases from 536 to 669 as the Bi number increases from 0.0003 to 0.1. 5.6. Influences of gravity The gravity does not significantly influence the instability onset or the Marangoni flow patterns if the capillary channel is vertical, 0.07
As presented in Figs. 2 and 3, there is a very weak evaporation compensating flow from the meniscus center to the edges when
0.04
b c
Gravity
d 0.00
(a)
a Diagonal line
28.0
Center
Secondary critical Ma
(a)
14.0
No Gravity
0.0
580
28.8 560
14.4 540
520
Gravity 0.0
10-4
10-2
10-3
0.1
(b)
Bi number
(b) Fig. 5. (a) The location of the temperature peak on the meniscus (the start point of the Marangoni flow) varies (from ‘a’ to ‘d’) as the Marangoni number increases; (b) The secondary critical Ma Number varies with the Bi number.
Fig. 6. The influences of the gravity (d = 800 lm; ll/lv = 3; d/ll = 0.27) (velocity is scaled by a/d). (a) Gravity brings forward the onset of the Marangoni instability (Ma = 105, Ra1 = 1.1, Bo1 = 0.010); (b) Gravity does not influence the flow patterns except the flow direction once the critical Ma value is reached (Ma = 622, Ra1 = 6.5, Bo1 = 0.011).
Table 4 Comparison between the Rayleigh number and Marangoni number. Channel width
Evaporation flux
Meniscus temperature
Ra1
Ra2
Ma
1/Bo1
1/Bo2
200 lm 200 lm 800 lm 800 lm
3.5 103 8.6 104 3.5 103 kg/ 8.6 104 kg/
21.1 °C 24.0 °C 21.2 °C 24.0 °C
0.035 0.009 8.693 2.285
0.52 0.13 32.6 8.57
53.4 13.3 831.3 218.5
1539 1539 95.6 95.6
102.6 102.6 25.5 25.5
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
no matter upward or downward. Since the gravity is perpendicular to the meniscus, the influence of the gravity manifests itself as the Rayleigh–Bénard instability. As shown in Table 4, both the Rayleigh number (Ra1 and Ra2) and the Bo number (Bo = Ra/Ma) are very small because of the small geometry of the capillary channels. Therefore, Rayleigh–Bénard instability does not occur in this case and the Marangoni effect is dominating. If the capillary channel is horizontal and the gravity direction is along one edge of the channel as marked in Fig. 6, the influence of the gravity is observable. The evaporation results in a cooling effect on the meniscus, which creates downward buoyancy force along the meniscus. As a result, the onset of Bénard–Marangoni instability is brought forward significantly. An example is shown in Fig. 6a, the Marangoni flow occurs when the Ma value is only about 105. More specific work is needed in the future to make a more comprehensive understanding of the threshold of Bénard–Marangoni instability in such conditions. The direction of the Marangoni flow pattern is also influenced by the gravity. As shown in Fig. 6b, the direction of the Marangoni flow is random when the gravity is zero, i.e., the temperature peak could locate at any corner of the meniscus. However, the Marangoni flow flows only downward on the meniscus if the gravity is imposed. 6. Conclusions Bénard–Marangoni instabilities on evaporating menisci in adiabatic capillary-channels are studied. On a flat, initially isothermal meniscus, even though the tangential temperature gradient on the meniscus is initially zero, the perpendicular temperature gradient can still arouse the instability leading to Marangoni flows when the evaporation intensity or the channel width increases to a certain critical value. The critical Ma values for the instability are obtained, ranging from 180 to 220 as the Bi number increases from 1 106 to 0.1. The flow pattern on the meniscus is varying as the Ma value increases. When Ma value reaches 600, a secondary instability occurs and the location of the temperatur peak moves from the diagonal lines to the center lines. The gravity does not have significant impact on the instability when its direction is perpendicular to the meniscus. However, the gravity along the meniscus tangential direction brings forward the threshold of the instability significantly, and influences the direction of the flow pattern. The present work is a preliminary step to a comprehensive understanding of Marangoni instability in capillary structures. The temperature gradient that is perpendicular to the meniscus should be paid attention to in the related applications. Acknowledgments This work was supported by the Natural Science Foundation of China (Grant No. 51276003) and Common Development Fund of Beijing. References [1] H.K. Dhavaleswarapu, J.Y. Murthy, S.V. Garimella, Numerical investigation of an evaporating meniscus in a channel, Int. J. Heat Mass Trans. 55 (2012) 915– 924. [2] Y.R. Li, H.R. Zhang, C.M. Wu, J.L. Xu, Effect of vertical heat transfer on thermocapillary convection in an open shallow rectangular cavity, Heat Mass Trans. 48 (2012) 241–251. [3] S.C. Zhao, Q.S. Liu, R. Liu, H. Nguyen-Thi, B. Billia, Thermal effects on Rayleigh– Marangoni–Bénard instability in a system of superposed fluid and porous layers, Int. J. Heat Mass Trans. 53 (2010) 2951–2954. [4] D. Melnikov, T. Takakusagi, V. Shevtsova, Enhancement of evaporation in presence of induced thermocapillary convection in a non-isothermal liquid bridge, Microgravity Sci. Tech. 2012 (online). [5] A.S. Basu, Y.B. Gianchandani, Virtual microfluidic traps, filters, channels and pumps using marangoni flows, J. Micromech. Microeng. 18 (2008) 115031.
247
[6] Z. Pan, H. Wang, Z. Yang, Marangoni bifurcation flow in a microchannel Tjunction and its micropumping effect: a computational study, Chin. Phys. Lett. 29 (2012) 074702. [7] R. Bhardwaj, X. Fang, D. Attinger, Pattern formation during the evaporation of a colloidal nanoliter drop: a numerical and experimental study, New J. Phys. 11 (2009) 075020. [8] H. Hu, R.G. Larson, Marangoni effect reverses coffee-ring depositions, J. Phys. Chem. B 110 (2005) 7090–7094. [9] J.R.A. Pearson, On convection cells induced by surface tension, J. Fluid Mech. 4 (1958) 489–500. [10] H. Bénard, Les tourbillons cellulaires dans une nappe liquide transportant de la chaleur par convection en regime permenant, Ann. Chim. Phys. 23 (1901) 62– 144. [11] J.A. Maroto, V. Perez-Munuzuri, M.S. Romero-Cano, Introductory analysis of Bénard–Marangoni convection, Euro. J. Phys. 28 (2007) 311. [12] Z.M. Tang, K. Li, W.R. Hu, Influence of free surface curvature of a liquid layer on the critical Marangoni convection, Int. J. Heat Mass Trans. 51 (2008) 5102– 5107. [13] H.S. Das, B.D. MacDonald, C.A. Ward, Stability of evaporating water when heated through the vapor and the liquid phases, Phys. Rev. E 81 (2010) 036318. [14] J.W. Peterson, G.F. Carey, Parallel adaptive solution of coupled Rayleigh– Bénard–Marangoni problems with the Navier-slip, Int. J. Numer. Methods Fluid 66 (2011) 428–451. [15] C. Ma, D. Bothe, Direct numerical simulation of thermocapillary flow based on the volume of fluid method, Int. J. Multiphase Flow 37 (2011) 1045–1058. [16] F. Doumenc, T. Boeck, B. Guerrier, M. Rossi, Transient Rayleigh–Bénard– Marangoni convection due to evaporation: a linear non-normal stability analysis, J. Fluid Mech. 648 (2010) 521–539. [17] G. Toussaint, H. Bodiguel, F. Doumenc, B. Guerrier, C. Allain, Experimental characterization of buoyancy- and surface tension-driven convection during the drying of a polymer solution, Int. J. Heat Mass Transfer 51 (2008) 4228– 4237. [18] R. Touihri, A. EI Gallaf, D. Henry, H. Ben Hadid, Instabilities in a cylindrical cavity heated from below with a free surface I effect of Biot and Marangoni numbers, Phys. Rev. E 84 (2011) 056302. [19] A.M. Benselama, S. Harmand, K. Sefiane, Thermocapillary effects on steadily evaporating contact line: a perturbative local analysis, Phys. Fluids 24 (2012) 072105. [20] C. Buffone, K. Sefiane, J.R.E. Christy, Experimental investigation of self-induced thermocapillary convection for an evaporating meniscus in capillary tubes using micro-particle image velocimetry, Phys. Fluids 17 (2005) 052104. [21] S.S. Panchamgam, A. Chatterjee, J.L. Plawsky, P.C. Wyner Jr., Comprehensive experimental and theoretical study of fluid flow and heat transfer in a microscopic evaporating meniscus in a miniature heat exchanger, Int. J. Heat Mass Trans. 51 (2008) 5368–5379. [22] X. Song, D.S. Nobes, Experimental investigation of evaporation-induced convection in water using laser based measurement techniques, Exp. Therm. Fluid Sci. 35 (2011) 910–919. [23] H. Wang, Z. Pan, S.V. Garimella, Numerical investigation of heat and mass transfer from an evaporating meniscus in a heated open groove, Int. J. Heat Mass Transfer 54 (2011) 3015–3023. [24] H.K. Dhavaleswarapu, P. Chamarthy, S.V. Garimella, J.R. Murthy, Experimental investigation of steady buoyant-thermocapillary convection near an evaporating meniscus, Phys. Fluids 19 (2007) 082103. [25] P. Chamarthy, H.K. Dhavaleswarapu, S.V. Garimella, J.Y. Murthy, S.T. Wereley, Visualization of convection patterns near an evaporating meniscus using lPIV, Exp. Fluids 44 (2008) 431–438. [26] H. Wang, J.Y. Murthy, S.V. Garimella, Transport from a volatile meniscus inside an open microtube, Int. J. Heat Mass Transfer 51 (2008) 3007–3017. [27] B. Lan, Y.R. Li, D.F. Ruan, Numerical simulation of thermocapillary flow induced by non-uniform evaporation on the meniscus in capillary tubes, Microgravity Sci. Technol. 23 (2011) 35–42. [28] C.A. Ward, F. Duan, Turbulent transition of thermocapillary flow induced by water evaporation, Phys. Rev. E 69 (2004) 056308. [29] F. Duan, C.A. Ward, Investigation of local evaporation flux and vapor-phase pressure at an evaporative droplet interface, Langmuir 25 (2009) 7424–7431. [30] C. Buffone, K. Sefiane, Controling evaporative thermocapillary convection using external heating: an experimental investigation, Exp. Therm. Fluid Sci. 32 (2008) 1287–1300. [31] C. Buffone, K. Sefiane, W. Easson, Marangoni-driven instabilities of an evaporation liquid-vapor interface, Phys. Rev. E 71 (2005) 056302. [32] Z. Pan, H. Wang, Symmetry-to-asymmetry transition of Marangoni flow at a convex volatizing meniscus, Microfluid. Nanofluid. 9 (2010) 657–669. [33] Z. Pan, F. Wang, H. Wang, Instability of marangoni toroidal convection in a micro channel and its relevance with the flowing direction, Microfluid. Nanofluid. 11 (2011) 327–338. [34] D. Gazzola, E. Franchi Scarselli, R. Guerrieri, 3D visualization of convection patterns in lab-on-chip with open microfluidic outlet, Microfluid. Nanofluid. 7 (2009) 659–668. [35] S.H. Davis, Thermocapillary instabilities, Annu. Rev. Fluid Mech. 19 (1987) 403–435. [36] H. Hertz, Ueber die verdungstung der flussigkeiten, inbesondere des quecksilbers, im luftleeren raume, Ann. Phys. Chem. 17 (1882) 177–200. [37] M. Knudsen, Die maximale verdampfungsgeschwindigkeit des quecksilbers, Ann. Phys. Chem. 47 (1915) 697–708.
248
Z. Pan, H. Wang / International Journal of Heat and Mass Transfer 63 (2013) 239–248
[38] R.W. Schrage, A Theoretical Study of Interface Mass Transfer, Columbia University Press, New York, USA, 1953. [39] R. Marek, J. Straub, Analysis of the evaporation coefficient and the condensation coefficient of water, Int. J. Heat Mass Transfer 44 (2001) 39–53. [40] B. Paul, Complication of evaporation coefficients, ARS J. 32 (1962) 1321–1328. [41] K.C.D. Hickman, Maximum evaporation coefficient of water, Ind. Eng. Chem. Res. 46 (1954) 1442–1446. [42] V.K. Badam, V. Kumar, F. Durst, K. Danov, Experimental and theoretical investigations on interfacial temperature jumps during evaporation, Exp. Therm. Fluid Sci. 32 (2007) 276–292. [43] S.M. Yang, W.Q. Tao, Heat Transfer, Higher Education Press, Beijing, China, 2006.
[44] S.R. Mathur, J.Y. Murthy, Pressure boundary conditions for incompressible flow using unstructured meshes, Numer. Heat Transfer B 32 (1997) 283–298. [45] G.J. Dunn, S.K. Wilson, B.R. Duffy, S. David, K. Sefiane, The strong influence of substrate conductivity on droplet evaporation, J. Fluid Mech. 623 (2009) 329– 351. [46] Y.A. Cengel, Heat Transfer: A Practical Approach, McGraw-Hill Press, New York, USA, 2003. [47] S. Semenov, V.M. Starov, R.G. Rubio, M.G. Velarde, Computer simulations of evaporation of pinned sessile droplets: influence of kinetic effects, Langmuir 28 (2012) 15203–15211.