- Email: [email protected]
Stability of an evaporating meniscus: Part I – Theoretical analysis
International Journal of Thermal Sciences 107 (2016) 209e219
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Stability of an evaporating meniscus: Part I e Theoretical analysis John Polansky*, Tarik Kaya Carleton University, Department of Mechanical and Aerospace Engineering, Ottawa ON, Canada
a r t i c l e i n f o
a b s t r a c t
Article history: Received 19 November 2015 Received in revised form 29 March 2016 Accepted 30 March 2016 Available online 22 April 2016
The stability of an evaporating meniscus in a narrow channel is theoretically investigated. A non-linear differential equation is derived using the lubrication approximation and taking into account evaporation, long-range molecular forces, surface tension, vapor recoil, thermocapillarity and viscous forces. Using a linear stability analysis, the susceptibility of the film to instability is studied. Three different base state profiles are used: flat replenished thin film, modified Young-Laplace profile and curved thin film. The last base state most accurately reflects an evaporating meniscus; however, its implementation requires the numerical solution of the non-linear thin-film equations. A shooting technique together with a fourth order Runge-Kutta method is used to obtain the curved thin-film profiles. The results obtained for each base state are analyzed and limitations are discussed. The curved thin-film base profile predicts strong instability for higher wall superheats. The growth rates are found to be positive (destabilizing) in the thin-film region and become negative (stabilizing) towards the intrinsic meniscus. The results indicate that the instability most likely originates in the evaporating thin-film region. Part II of the presented paper is devoted to experimental investigation. © 2016 Elsevier Masson SAS. All rights reserved.
Keywords: Meniscus instability Evaporating meniscus Thin films
1. Introduction The stability of evaporating thin films are important for several technological applications such as the cooling of power plants and advanced electronics, microfluidic devices, and heat pipes. In these applications, the heat from a solid substrate is conducted through a thin liquid film to its free surface to take advantage of the high heat transfer rates associated with phase change at the liquid-vapor interface. However, if the interface equilibrium is disturbed, a potential runaway condition may lead to the eventual rupture or collapse of the thin film upon the solid wall. Instabilities leading to film rupture would constitute a failure of the device as damaging dry hot spots would result in the inability to transport the applied heat. Thus, to appropriately design and implement two-phase heat transfer devices it is important to characterize these instabilities. The stability of an evaporating liquid film is governed by a complex interaction of mechanisms spanning mass loss due to evaporation, long-range molecular forces, surface tension, vapor recoil, thermocapillarity and viscous forces. The combination of these factors is the basis from which instabilities may arise. In the case of an evaporating meniscus, the geometric constraints and associated thin-film geometry, see Fig. 1, can further complicate the
* Corresponding author. E-mail address: [email protected] (J. Polansky). http://dx.doi.org/10.1016/j.ijthermalsci.2016.03.021 1290-0729/© 2016 Elsevier Masson SAS. All rights reserved.
problem. Since the scales and associated flow rates are small, the gravitational effects and liquid inertia may be considered as negligible. To explain a potential instability mechanism of an evaporating meniscus, we will assume that a local small disturbance results in a surface depression of the liquid-vapor interface, causing a departure of the system from equilibrium. The depressed interfacial region closest to the solid wall becomes warmer and experiences a local increase in evaporation. The resulting increased vapor recoil then forces the liquid towards the wall, amplifying the disturbance. The spatial changes in interface thickness lead to spatial temperature gradients and in turn thermocapillarity by way of surface tension variation (Marangoni instability). The thermocapillarity serves to further enhance the runaway effect as warmer liquid from the depression is drawn toward the cooler regions of the interface. As the film gets thinner, the increasing disjoining pressure will also provide positive feedback. The principal stabilizing mechanism here appears to be the capillary force as it will increase with increasing deformation of the interface, thereby forcing the system back toward equilibrium. It is thus possible to reach a new equilibrium state by means of rolling cell structures with a wavy interface in a similar fashion observed nard or Be nard-Marangoni for fluid instabilities such as Rayleigh-Be cells. Alternatively, there is another potential mechanism promoting equilibrium: as the interface gets closer to the solid wall,
210
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
Nomenclature a A A ;B;C b cp C Ca D E h hfg hp H H i j J k L M Ma b n p p P Pe R Re Rem t b t T T u U v V V
Clapeyron coefficient [kg/m2s K] Dispersion constant [J] Coeff. of growth rate equation [e] Kelvin coefficient [kg/m2s Pa] Specific heat [J/kg K] Accommodation coefficient [e] Capillary number [e] Difference [e] Evaporation number [e] Film thickness [m] Latent heat [J/kg] Perturbation thickness [m] Base state film thickness [e] Dimensionless film thickness [e] Imaginary unit [e] Mass flux [kg/m2 s] Dimensionless mass flux [e] Thermal conductivity [W/m K] and disturbance wavenumber Thin-film length [m] Molar mass [kg/kmol] Marangoni number [e] Normal unit vector [e] Pressure [Pa] Pressure vector [Pa] Dimensionless pressure [e] Peclet number [e] Universal gas constant [J/mol K] Reynolds number [e] Modified Reynolds number [e] Time [s] Tangential unit vector [e] Temperature [K] Cauchy stress tensor [Pa] Velocity in x [m/s] Dimensionless velocity in x [e] Velocity in y [m/s] Velocity vector [m/s] Dimensionless velocity in y [e]
increased evaporation may cause local cooling of the depressed region thereby reducing the vapor recoil and thermocapillary effects. It is such complex interplay among these various mechanisms that needs to be studied together so as to accurately characterize the potential for meniscus instability. The problem of flat thin-film stability has been investigated through a series of theoretical studies aimed at describing the film and its propensity to destabilize. The earliest stability studies focused on long range molecular forces in non-evaporating films [1e3]. Utilizing a thermodynamic approach Sheludko [1] was able to determine that a perturbation upon a non-evaporating liquid thinfilm can destabilize due to long range molecular forces. The destabilization of the thin film was correlated with a characteristic film thickness. This finding was further supported by the findings of Vrij [2], who correlated the critical film thickness with a critical perturbation wavenumber. The characterization of critical wavenumbers was further supported by the work of Ruckenstein and Jain [3]. Following this, Williams and Davis [4] constructed a non-linear stability model for which they performed a numerical analysis for a finite perturbation. A comparison of their non-linear model with
Vl W x, y X, Y
Molar volume [m3/mol] Channel width [m] Cartesian coordinates [m] Dimensionless coordinates [e]
Greek symbols Tangent to Cartesian x angle [rad] Surface tension coefficient [N/m K] Dispersion number [e] Non-evaporating film thickness [m] Perturbation amplitude [e] Wall superheat [K] Characteristic ratio [e] Dimensionless temperature [e] L Liquid velocity coefficient [e] k Local curvature [m1] m Dynamic viscosity [Pa s] n Kinematic viscosity [m2/s] r Density [kg/m3] s Surface tension [N/m] s0 Reference surface tension [N/m] t Dimensionless time [e] f Coordinate transformation coeff. [e] F Dimensionless coord. trans. coeff. [e]
a g G d0 dH DT z q
Subscripts c Characteristic and cutoff cap Capillary d Disjoining lv Liquid-vapor interface max Maximum peak Peak value of thin film property ref Reference t Derivative in dimensional time t Derivative in dimensionless time v Vapor w Wall x Derivative in dimensional x X Derivative in dimensionless X y Derivative in dimensional y Y Derivative in dimensionless Y
that of a linear model yielded a reduction in the time to rupture in addition to shifting the required perturbation wavelength to that of smaller values. The findings of Williams and Davis [4] are further supported by Sharma and Ruckenstein [5] who also found the nonlinear effects to reduce the rupture time. Their work also included the addition of Marangoni effects by way of surfactants. The nonlinear effects attributed to Marangoni motion were found to enhance the film's stability. The consideration of Marangoni effects on thin-film stability was continued in Davis' [6] work in which thermocapillary effects were the principal focus. This analysis highlighted the thermal gradients as the source of energy for hydrodynamic instabilities, which manifest as two and three-dimensional surface waves. The influence of phase change on thin-film stability was presented by Bankhoff's [7] study of evaporating and condensing thin liquid films on an inclined plane. The results of the study found the effects of phase change to be destabilizing for evaporation and stabilizing for condensation. Combining the aforementioned influences on thin-film stability was the work of Burelbach, Davis, and Bankhoff [8]. Their mathematical formulation was constructed around a one-sided model,
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
211
stability analysis. Different thin-film base state profiles are presented for the analysis. The results and limitations of the different base states are then discussed. The more realistic curved thin-film base profile shows that the film becomes more susceptible to meniscus rupture with increasing wall superheats. The results also suggest that the instability most likely originates in the evaporating thin-film region. This research is performed in two parts, with the second part being devoted to experimental investigation, which is presented in a separate paper. 2. Mathematical formulation
Fig. 1. Schematic of an extended meniscus in a channel.
which simplified the free surface boundary condition, yielding a single non-linear differential equation. Using a perturbation analysis, they were able to determine the potential for thin-film instability. The numerical results of the analysis demonstrated the film was liable to destabilize faster for an evaporating film than for a condensing film. The conditions necessary for instability were related to the nature of the perturbation by way of a critical wavenumber as is consistent with previous works. Following this work, Shklyaev and Fried [9] developed a modified linear stability model, which removed the singularity condition found at the point of rupture as discussed by Burelbach, Davis, and Bankhoff [8]. Instead of rupture, the film was found to thin until such a point that a non-evaporating thin-film developed on the substrate. Benselama, Harmand and Sefiane [10], developed a mathematical model to determine the stability of an evaporating curved interface. The analysis included long range molecular forces, surface tension and thermocapillarity effects. Their analysis was limited to moderate evaporation rates so as to neglect vapor recoil effects. This work used a domain perturbation analysis technique, while retaining higher order terms pertaining to the lubrication problem. The results indicated the potential for instability around the non-evaporating thin-film region, dependent on the sign of the Marangoni number. Additional work by Avramenko et al. [11], considered curved thin-film stability with thermocapillarity and long range molecular forces. The analysis demonstrated changes in interfacial temperature gradients shift the location of the peak replenishing fluid velocity. The introduction of curved thin-film geometries has served to introduce a spatial component to thinfilm stability. It is this spatial aspect that may provide an indication as to the location of the instability for cases such as that of an evaporating meniscus. Comprehensive reviews on thin-film stability can be found in Oron, Davis, and Bankoff [12], as well as Craster and Matar [13]. The main objective of the present work is to develop a mathematical model to study the stability of an evaporating meniscus. For the evolution of the thin-film thickness, a non-linear differential equation is derived using the lubrication approximation. Evaporation, long-range molecular forces, surface tension, vapor recoil, thermocapillarity and viscous forces are accounted for in the derivation. The governing equations with the associated boundary conditions and scaling of these equations are explained thereafter. The susceptibility of the film to instability is studied using a linear
A liquid meniscus constrained within a rectangular channel with immiscible walls of width W and infinite depth is considered as shown in Fig. 1. Heat is applied to the sides of the rectangular channel walls. The symmetry condition at the center of the channel reduces the problem to that of a thin film of increasing thickness. A microscopic region with a flat non-evaporating thin-film near the wall merges into the evaporating thin-film region. At the far end, a macroscopic region hereafter referred to as the intrinsic meniscus is dominated by the curvature of the liquid-vapor interface. The extended meniscus profile can therefore be described in a Cartesian coordinate system with the origin positioned at the transition between planar non-evaporating thin-film region (for x < 0) and the evaporating thin-film region (for x 0). For this geb and b ometry, the normal and tangential unit vectors ( n t) can be written as:
b¼ n
1 hx f 1
1 1 b t¼ f hx
(1)
(2)
where,
f¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ h2x
(3)
2.1. The governing equations The working fluid is assumed to be incompressible, Newtonian and non-polar with its surface tension varying linearly with temperature. The channel geometry is limited to sub-millimeter sizes such that the Bond number is sufficiently small so as to ignore gravitational contributions. The fluid is assumed to be of sufficient thickness that continuum mechanics hold, while viscous dissipation effects are neglected. The variables associated with the vapor and interface are denoted with a subscript v and lv, respectively, and the subscript l for liquid phase being omitted for simplicity. Under these assumptions, the two-dimensional continuity, Navier-Stokes and the energy equations can be written as:
V$V ¼ 0
(4)
rðV t þ V$VVÞ ¼ Vp þ mV2 V
(5)
rcp ðTt þ V$VTÞ ¼ kV2 T
(6)
The boundary conditions at the wall (y ¼ 0) are imposed as a no slip, no penetration and constant temperature.
u¼0
(7)
212
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
v¼0
(8)
T ¼ Tw
(9)
The boundary conditions at the liquid-vapor interface, (y ¼ hðx; tÞ), address the jump conditions associated with evaporation. As the interface is curved, the interfacial jump conditions are b and b described in the n t coordinate system. The mass jump balance at the interface can be written as:
b ¼ rv ðV v V lv Þ$ n b j ¼ rðV V lv Þ$ n
Lastly, an evaporation model is needed to relate the evaporative mass flux to the interfacial conditions. The Kelvin-Clapeyron constitutive equation is used for this purpose: j ¼ aðTlv Tv Þ þ bðp pv Þ [15]. The second term in this equation is ignored as it was also done in [8,11] considering that this term will not affect the results significantly but introduce unnecessary complexity. The constitutive equation then becomes:
3 2sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 M Pv h M C fg 5ðTlv Tv Þ j¼4 2pRTlv RTlv Tv
(23)
(10)
Taking into account the velocity components in both x- and ydirections for a curved interface, the liquid and interface velocities can be written as:
V ¼ ½u v
(11)
V lv ¼ ½ulv vlv
(12)
Considering that the density, viscosity and thermal conductivity of the vapor phase are much smaller than those of the liquid phase, the vapor flow field is decoupled from the liquid (one sided model), as it was done in earlier related works such as Burelbach, Bankoff and Davis [8]. For the one sided model, the normal and shear stresses and energy jump conditions are then simplified to:
j2
rv
where,
T
b Þ $n b ðp pv Þ ¼ þ ðT$ n
shxx f3
þ
A h3
(24)
ulv ¼ sinðaÞht
(13)
b ÞT $b t ¼ Vs$b t ðT$ n
(25)
vlv ¼ cosðaÞht
(14)
b¼0 jhfg þ kVT$ n
(26)
a ¼ tan1 ðhx Þ
(15)
Substituting Eq. (11)e(15) into Eq. (10) and simplifying, we get the mass jump condition for a curved interface as:
rðfht þ hx u vÞ
j¼
(16)
f
The normal-stress jump condition at the interface is: T
b ðp pv Þ ¼ pcap þ pd b þ ððT T v Þ$ n b Þ $n jðV V v Þ$ n
(17)
The vapor field is assumed isotropic with a pressure equal to that of the saturated vapor pressure. The capillary pressure is proportional to the surface tension s and the local curvature k giving:
pcap ¼ sk ¼
shxx
(18)
f3
The surface tension is represented as a function of temperature as:
s ¼ s0 g Tlv Tref
(19)
The disjoining pressure is expressed as [14]:
pd ¼
2.2. Scaling of the equations The governing equations and boundary conditions are nondimensionalized using scales appropriate to the problem [8,10]. The characteristic scale in the x-direction is assigned as xc and comes as a result of subsequent scaling. The thickness of the nonevaporating thin-film layer d0 is used as the characteristic scale in the y-direction, where the non-evaporating thickness is obtained from the Kelvin-Clapeyron evaporation model as: d0 ¼ ðAbDT=aÞ1=3 . By using the continuity equation, Eq. (4), a relation can be found between the transversal vc and longitudinal uc characteristic velocities as vc =uc ¼ d0 =xc . For the transversal characteristic velocity, we use vc ¼ aDT=r where DT ¼ Tw Tv . Thus the characteristic velocity in x-direction then takes the form:
uc ¼
aDTxc
The characteristic pressure Pc can be deduced from the y-direction momentum equation as Pc ¼ muc xc =d20 . Knowing that the characteristic pressure Pc can also be obtained from Eq. (20), it is possible to deduce a relation for the longitudinal characteristic length as:
A h3
(20)
where A > 0 is the dispersion constant. A no slip shear-stress jump condition at the interface can be written as:
b Þ $b ððT T v Þ$ n t ¼ Vs$b t T
b¼0 jhfg þ ðkVT kv VTv Þ$ n
xc ¼
(21)
(22)
A naDT
1=2 (28)
The viscous scale is used for non-dimensionalizing time as follows:
tc ¼
The energy jump condition is simplified assuming kinetic energy and viscous dissipation terms to be negligible yielding:
(27)
rd0
xc d 0
n
(29)
Lastly, the temperature is scaled as the ratio of temperature differences as follows:
q¼
Tw T Tw Tv
(30)
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
Applying the characteristic scales to Eqs. (4)e(6), we obtain the non-dimensional forms of continuity, Navier-Stokes, and energy equations as:
UX þ VY ¼ 0
(31) 2
RezðUUX þ VUY Þ ¼ PX þ z UXX þ UYY 3
2
Rez ðUVX þ VVY Þ ¼ PY þ z
2
z VXX þ VYY
(33) (34)
where the dimensionless coordinates are defined as: X ¼ x=xc and Y ¼ y=d0 . The characteristic ratio (z), Reynolds (Re) and Peclet (Pe) numbers are defined by:
d0
(35)
xc
ruc d0 aDTxc ¼ m m
Re ¼
Ca ¼
G¼
Pez2 ðU qX þ V qY Þ ¼ z2 qXX þ qYY
z¼
Rem ¼
(32)
(36)
213
(47)
maDT rs0
(48)
rA m2 d0
Ma ¼
Е¼
d0 vc r n rv
(49)
gDTxc rcP mk
(50)
d0 ahfg
(51)
k
Here Rem , Ca, G, Ma, and Е are the dimensionless modified Reynolds, capillary, dispersion, Marangoni and evaporation numbers respectively. Knowing that the characteristic ratio is small, z≪1, we can reduce the governing equations, Eqs. (31)e(34), thereby obtaining the leading-order governing equations as:
UX þ VY ¼ 0
(52)
PX ¼ UYY
(53)
For the heated wall (Y ¼ 0) boundary conditions, Eqs. (7)e(9), the non-dimensional forms are then:
PY ¼ 0
(54)
U¼0
(38)
qYY ¼ 0
(55)
V ¼0
(39)
q¼0
(40)
The interfacial jump conditions and constitutive equation, Eqs. (41)e(45), are also reduced on the basis that characteristic ratio is small. The mass jump and constitutive equations, Eq. (41) and Eq. (45), remain unchanged. We selectively retain the vapor recoil, surface tension, long range molecular forces and thermocapillary contributions in the normal and shear stress boundary conditions so as to determine their effect on the stability of an evaporating thin film. Therefore, we obtain the leading-order interfacial stress and energy jump conditions as:
Pe ¼
rcp uc xc k
cp aDTx2c ¼ d0 k
(37)
The interfacial height and time are non-dimensionalized as H ¼ h=d0 and t ¼ t=tc , respectively. The jump conditions and constitutive equation, Eqs. (16), (23)e(26), at the interface, Y ¼ H ðX; tÞ, then come as:
J¼
ReðV H X UÞ FH t ReF
(41)
2
J 2 Rem z ¼
2z2 z2 H X ðVX H X UX Þ þ H X UY þ UX 2
F
ðP Pv Þ
z4 H XX Gz þ ReH 3 CaF3
J 2 Rem z2 ðP Pv Þ ¼
z4 H XX Gz þ ReH 3 CaF3
(56)
Pe F2 2 UY ¼ MazFðqX þ H X qY Þ
(57)
ЕJ F qY ¼ 0
(58)
(42) 2 PeUY F2 2 þ Pez 4H X UX þ VX F2 2
2.3. The evolution equation
¼ MazFðH X qY þ qX Þ
(43)
JЕF qY þ z H X qX ¼ 0
2
(44)
J ¼1q
(45)
where,
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F ¼ 1 þ ðzH X Þ2
The system of differential equations, Eqs. (52)e(55), together with boundary conditions described above are solved herein so as to obtain a single differential equation describing the evolution of an evaporating interface. We begin by integrating Eq. (55) twice, subject to the wall boundary and interfacial jump conditions, Eq. (40) and Eq. (58), respectively. The dimensionless temperature is then found to be:
q ¼ ЕJ FY (46)
(59)
To obtain a relation for the mass flux, we can substitute Eq. (45) into Eq. (59) and evaluate at the interface, Y ¼ H ðX; tÞ, thereby obtaining the dimensionless mass flux relation as:
214
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
2.4. Linear stability analysis
1 J¼ ЕFH þ 1
(60)
The substitution of Eq. (60) into Eq. (59) yields a new relation for the dimensionless temperature as:
q¼
Е FY Е FH þ 1
(61)
Next, we solve the Navier-Stokes equation in the Y-direction by integrating Eq. (54) and obtaining the constant of integration from the normal stress jump condition, Eq. (56), and substituting the previously obtained mass flux relation, Eq. (60), yielding:
P ¼ Pv þ
Rem z2 ðЕFH þ 1Þ
2
z4 H XX Gz ReH 3 CaF3
(62)
Given the pressure, we can now solve Navier-Stokes in the Xdirection by differentiating Eq. (62) with respect to X and substituting into Eq. (53). The resulting substitution yields a result such that the right hand side is only a function of X and t. Equation (53) is then integrated twice, subject to no slip at the wall, Eq. (38), and the shear stress jump condition, Eq. (57), at the interface. Finally, the X-component of the liquid velocity is obtained as:
2EMazFFX 1 ЕF2 H X L Y 2 þ EMazF2 H X ðЕ FH þ 1ÞY U¼ Pe F2 2 ðЕ FH þ 1Þ2 (63) where
L¼
z4 ð3FX H XX FH XXX Þ 2Rem z2 ðЕFH X þ ЕFX H Þ CaF4 ðЕFH þ 1Þ3 3GzH X þ
ReH 4
(64)
To analyze the stability of an evaporating film, we perturb the base state film thickness H by a small perturbation hp :
H ðX; tÞ ¼ HðX; tÞ þ hp ðX; tÞ
(66)
The perturbation is assumed to be of the form:
hp ðX; tÞ ¼ dHðtÞ expðikXÞ
(67)
Substituting Eq. (66) into Eq. (A1) and removing the base state by way of binomial series and simplifying, we obtain the resulting leading-order differential equation of the following form:
hpt ¼ A k4 þ B k2 þ C hp
(68)
Equation (68) describes the growth rate of a perturbation with the coefficients provided in Appendix B. By using Eq. (68), it is then possible to describe the range of disturbance wavenumbers for which a small disturbance will grow (instability) or decay (stability). 3. Analysis and discussion The results of the linear stability analysis for the thin film are presented in this section. All the results are obtained for n-pentane as the working fluid. The base profile is that of a steady-state evaporating curved thin film, which can be computed from the steady-state evolution equation, Eq. (A1). Given the non-linearity, no closed form solution is possible to be used as a base profile HðX; tÞ, making the linear stability analysis highly complex. Therefore, we first consider an evaporating flat thin film and then a modified Young-Laplace profile as the base states. Then, we consider the more realistic, though more complex, curved evaporating thin film for the base state. 3.1. Flat thin film
We can now solve the continuity equation, by differentiating Eq. (63) with respect to X and substituting into Eq. (52) and integrate once with respect to Y subject to the no penetration boundary condition, Eq. (39). Using the transverse velocity component V in conjunction with Eq. (60) and Eqs. (63) and (64), we can solve the mass jump condition, Eq. (41), at the interface, thereby obtaining the evolution equation. We note that the evolution equation contains F, which is a function of X and t, as denoted in Eq. (46). It is known that the thin film slope and characteristic ratio are sufficiently small such that, zH x ≪1 to make the following approximation:
A flat thin film base state was considered as an approximation to study the stability of an evaporating curved thin film. The quasiparallel theory was also used in previous works, for example in Unsal [16] for the stability of condensing flows and in Pratt, Brown and Hallinan [17] for the evaporating meniscus. For an evaporating thin film to be at steady state, a flow of fluid into the thin film is needed to replenish that which is lost to evaporation. If we consider an evaporating flat film of a steady height, the base state will be of the following form:
FðX; tÞy1
H¼1
(69)
HX ¼ HXX ¼ HXXX ¼ HXXXX ¼ 0
(70)
(65)
Therefore, we simplify the evolution equation, subject to Eq. (65) thereby yielding a thin-film evolution equation. Thus, we obtain a single non-linear differential equation which describes how the time rate of change of the film thickness is related to mass loss due to phase change, long-range molecular forces, surface tension, vapor recoil, thermocapillarity and viscous forces. The resulting evolution equation is sizable and presented in Appendix A. Because of the large size of the resulting equation, all parametric derivation is carried out in Maple and the resulting differential equation is solved using MATLAB.
Substituting the base state conditions, Eqs. (69) and (70), into the growth rate equation, Eq. (68), and after some manipulation, we obtain the coefficients as:
z4 Re
A ¼
6Ca
(71)
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
B ¼
zEð2PeReRem z 9PeG 3MaReÞ 6PeðE þ 1Þ3
C ¼
zE2 ð3PeGðE þ 3Þ þ MaReðE þ 4ÞÞ 6PeðE þ 1Þ3
3PeGz 6PeðE þ 1Þ3
ERe ðЕ þ 1Þ2
H¼ (72)
(73)
Here we note that the stability of the meniscus is not spatially dependent as all of the base state derivatives with respect to X are zero, thereby leaving the solution as a function of the perturbation (disturbance) wavenumber k, and the applied superheat DT. Fig. 2 shows a plot of the growth rate as a function of the disturbance wavenumber for a range of wall superheats.
Fig. 2. Growth rate as a function of disturbance wavenumber.
The negative growth rates indicate stability. As the superheat increases, the range of wavenumbers capable of destabilizing the flat film decreases. In the unstable region, there is a cross over after which the higher wall superheats lead to larger growth rates. For a given wall superheat, the growth rate continuously increases toward a maximum as k/0.
W 2
215
ffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 W þ d0 ðXxc Þ 2
d0
(74)
Substituting the base state profile, Eq. (74), into the growth rate equation, Eq. (68), we obtain the fourth-order polynomial growth rate equation for the modified Young-Laplace base case. As the modified Young-Laplace base state is a function of position and is continuous and differentiable, the stability analysis becomes a function of position, unlike that of a flat thin-film case. The resulting growth rate equation is notably more complex than that of the flat film case and not reproduced here as it is comprised of 875 operands. For the modified Young-Laplace base state we plot the growth rate for a typical case, DT ¼ 0:05, W ¼ 10 mm and k ¼ 0:5; as it is shown in Fig. 3. It should be noted that in the proximity of the symmetry point at x ¼ W=2, HX /∞, and the assumption zHX ≪1 breaks down. For the typical case presented in Fig. 3, the solutions should be sufficiently reliable for the meniscus positions less than 4.5 mm where zHX is of an order of magnitude of 102 or smaller. It can be seen from Fig. 3 that the film is increasingly stable as one moves away from the evaporating thin-film region towards the intrinsic meniscus. The positive growth rates observed in the evaporating thin-film region indicates instability as illustrated in Fig. 3. The growth rate rapidly changes sign as a function of the position in the evaporating thin-film region. As the wavenumber decreases asymptotically toward zero, k/0, the scatter in the thin-film region reduces. For a range of wavenumbers less than some critical value, take for example kc ¼ 0:5, the growth-rate curve becomes smooth without scatter and the thin-film remains stable for DT < 105 K in a channel of W ¼ 10 mm. For the same critical wavenumber, as the channel width increases, the interfacial temperature jump at which the instability is first observed also decreases suggesting that the large channels are more susceptible to instability. While the growth rates increase with increasing wavenumber, the intrinsic meniscus region remains generally unaffected. Thus, the growthrate analysis indicates an inherent stability in the intrinsic meniscus region while the evaporating thin-film region is susceptible to instability. Considering the local heat transfer and evaporation rates peak in the evaporating thin-film region [18], the finding that the thin-film region could be the origin of instability is not surprising. This point will be further discussed later.
3.2. Modified Young-Laplace profile
3.3. Curved thin film
A base state described by the Young-Laplace equation is of interest as it has a closed form and is easy to implement for the linear stability analysis. The modified Young-Laplace equation represents a more realistic base state than that of a flat thin film. The modified Young-Laplace base state is described by the segment of a circle, centered on the symmetry plane of the channel and reduced in width by twice the non-evaporating thinfilm thickness. While the non-evaporating thin-film thickness is a function of the wall superheat, its scale as compared to the channel width places emphasis on the geometry over the applied superheat. Thus, we define the modified Young-Laplace base state as:
As there is no closed form solution for the steady-state evaporating curved thin-film base state, the governing equations need to be solved numerically to obtain the thin-film thickness and the associated derivatives at discrete locations along the extended meniscus. Then, these values can be substituted into the growthrate equation, Eq. (68), at each corresponding location to obtain the growth rate along the entire meniscus. Note that the base state is timeeindependent as the evaporated fluid will be replenished by the flow in the capillary channel. To obtain the thin film profile, the steady-state form of the evolution equation, Eq. (A1), needs to be solved. Using a slightly modified version of Eq. (A1), we first solve the governing non-
216
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
Fig. 3. Growth rate as function of position for the modified Young-Laplace base state, DT ¼ 0:05 K, W ¼ 10 mm and k ¼ 0:5.
linear equations using a fourth order Runge-Kutta method. Solving the evaporating thin-film equation requires a set of initial conditions to achieve the target curvature of 2=W at the far field with prescribed constant curvature tolerance of ± 0:02=W. As the sensitivity of the non-linear equation lies with the first derivative, hx , a shooting technique is employed to find the required initial slope to obtain the desired far field boundary condition. The results are found to be extremely sensitive to the initial conditions. The details of the solution methodology were presented in Ball,
Polansky and Kaya [19]. As the thin-film profile is a function of wall superheat, the profiles must be obtained over a range of wall superheats for a given channel. To ensure the validity of the stability model in the thin-film region, the slopes of the thin-film profiles are checked against the assumption: zHX ≪1. For the cases tested, we found the assumption to hold. The obtained thin-film thickness and its derivatives are shown in Fig. 4. Given the strong dependence on superheat, the thin-film curvature profiles transition from an asymptotic rise to 2=W for small superheats to an overshoot before achieving the far field curvature as illustrated in Fig. 4(b). The location of the overshoot peak is noted to shift spatially toward the non-evaporating thinfilm region with increasing superheats. The higher order derivatives also change sign in response to the change in curvature as shown in Fig. 4(c) and (d). These changes in sign on the higher derivative values along the thin film directly affect the growth rate results through Eq. (68). The resulting disturbance cutoff wavenumbers kc obtained using the curved thin-film solutions are shown in Fig. 5 for a range of wall superheats and meniscus locations in a channel 100 mm wide. The cutoff wavenumber is calculated by assigning the growth rate to zero, hpt =hp ¼ 0. As the thin-film solution is dependent upon the fluid properties, superheat and position so too is the cutoff wavenumber. The broken appearance of the surface as the cutoff wavenumber approaches to zero, kc /0; is a consequence of the sharp gradients coupled with the numerical resolution of thin-film solutions in this region. The stable region corresponds to the negative growth values and occupies the top portion of the curved plane.
Fig. 4. (a) Thin film thickness (b) local curvature (c) third derivative (d) fourth derivative for wall superheats ranging from 0.01 to 0.20 K, W ¼ 100 mm.
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
217
Looking at the most critical characteristics, we consider the locations where the maximum evaporative mass flux, interface temperature gradient, liquid pressure gradient and thin-film curvature occur in the base case profile. Fig. 7 shows the difference between the locations where the growth rate peaks ðxpeak Þ and a thin-film characteristic reaches a maximum ðxmax Þ, normalized by the length of the extended thin-film ðL Þ, as D ¼ ðxpeak xmax Þ=L . The local curvature has the smallest degree of offset and the differences continuously decrease with increasing wall superheat. The local curvature appears to be the most influential parameter. However, the differences shown in Fig. 7 are insufficiently varied to draw a sound conclusion. It is clear that all these thin-film characteristics reach a maximum in the evaporating thin-film region. As it was stated previously, the evaporating thin-film region seems to be the most likely region of the onset of instabilities.
Fig. 5. The disturbance cutoff wavenumber as a function of wall superheat and meniscus location, W ¼ 100 mm.
4. Conclusions
The cutoff wavenumber surface illustrated in Fig. 5 suggests that there are two regions of an evaporating thin film with a positive growth rate. The first occurs nearest the non-evaporating thin-film region, and the second occurs for a small portion of the thin film before tending back to negative growth rates towards the intrinsic meniscus region. The shape of the curve suggests that there is a peak cutoff wavenumber above which meniscus is always stable. Fig. 6 shows the growth rates for a range of wall superheats (0.01e0.20 K) for a channel width of 100 mm and a disturbance wavenumber of 0.15. The growth rates for the thin-film base case are observed to be smooth and continuous as compared to the sporadic nature found in the modified Young-Laplace base case. As the applied superheat increases, so too does the peak value of the growth rate. This suggests that higher superheats are associated with strong instability and will drive the meniscus to rupture faster than lower superheats for a given perturbation. Similar to the modified YoungLaplace base case, the growth rate starts off positive (destabilizing) in the thin film region with a peak positive growth rate and transitions to a negative (stabilizing) further along the meniscus. The position of the maximum positive growth rate shifts towards the nonevaporating thin-film region with increasing superheats. The location where the growth rate peaks should drive the instability. Thus, it would be interesting to investigate what other thin-film characteristics reach a maximum around this location.
The stability of an evaporating meniscus is studied using a linear stability analysis. For the base state, three different profiles are considered; a flat replenished thin film, modified Young-Laplace profile and a curved thin film. The flat thin-film analysis shows that the range of destabilizing wavenumbers decreases for increasing wall superheats. For the modified Young-Laplace base profile, the instability is principally driven by the position along the meniscus and the perturbation wavenumber. Though the solution is subject to the assumption zHX ≪1, the validity extends well into the intrinsic meniscus where the growth rate is found to be negative thereby indicating strong stability. The growth rate analysis also indicates that the evaporating thin-film region is susceptible to instability, where the local heat transfer and evaporation rates peak. The first two base profiles have the advantage of leading to the closed form solutions; however, they are not a very accurate representation of an evaporating meniscus. For a more accurate analysis, the steady-state thin-film profiles are numerically solved by using a fourth order Runge-Kutta method. It is found that there exists a cutoff wavenumber above which the meniscus always remains stable. The results predict strong instability for higher superheats. The growth rates are found to be positive (destabilizing) in the thin film region where they reach a maximum and become negative (stabilizing) towards the intrinsic meniscus. As the superheat increases, the position of the maximum positive growth rate shifts towards the non-evaporating thin-film region. No direct
Fig. 6. Growth rates calculated as a function of meniscus location for a range of wall superheats; W ¼ 100 mm and k ¼ 0:15.
Fig. 7. Comparison of the locations where the growth rate and thin-film characteristics peak for a channel of W ¼ 100 mm.
218
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
correlation is found between the location of the peak growth rate and individual thin-film characteristics. The results obtained using modified Young-Laplace and curved thin-film base profiles both indicate that the evaporating thin-film region is the most likely origin of the instability.
Appendix B. The perturbation growth rate equation coefficients
4
4
Acknowledgements
The support from Natural Sciences and Engineering Research Council of Canada under Discovery Grants Program is greatly appreciated.
5E2 Rez4 H5 3CaðEH þ 1Þ
B ¼
Appendix A. The thin film evolution equation
19E2 MaRezFFX H 3 H X 4E3 MaRezF2 FX H 4 H X 2 2 6Pe F 2 ðEFH þ 1Þ 3Pe F2 2 ðEFH þ 1Þ3
4
þ
2E2 MaRezF3 FX H 3 H X Rez H 3 H XXXX 2 6CaF4 3Pe F2 2 ðEFH þ 1Þ2 E2 MaRezF2 H 3 H
2
XX 2
3Pe F 2 ðEFH þ 1Þ
2H
EMaRezFH XX 2Pe F2 2 ðEFH þ 1Þ
þ
EMaRezFH H 2XX E2 MaRezFFXX H 4 2 2 3Pe F 2 ðEFH þ 1Þ Pe F2 2 ðEFH þ 1Þ
þ
2E3 MaRezF3 H 3 H 2X 2EMaRezFF2X H 3 2 2 3 3Pe F2 2 ðEFH þ 1Þ 3Pe F 2 ðEFH þ 1Þ
2E2 MaRezF2 F2X H 4 2E3 MaRezFF2X H 5 2 2 3 3Pe F2 2 ðEFH þ 1Þ2 3Pe F 2 ðEFH þ 1Þ 2EMaRezFX H 2 H X EMaRezF2 FX H 2 H X þ 2 2 Pe F 2 ðEFH þ 1Þ Pe F2 2 ðEFH þ 1Þ
þ
2
þ
þ
2CaF5 E2 ReRem z F2X H 5
FðEFH þ 1Þ
4
3
3FðEFH þ 1Þ
þ
GzH XX GzH 2X 2FH 2FH 2
2
2
5EReRem z FX H 3 H X 3
3FðEFH þ 1Þ
2
þ þ þ
E2 ReRem z FH 3 H 2X ðEFH þ 1Þ4 Rez4 FX H 3 H XXX 5
CaF
Rez4 FXX H 3 H XX 2CaF
5
2H 2 X 3
EReRem z H Re E FH þ 1 ðEFH þ 1Þ
5
2ðEH þ 1Þ
PeðEH þ 1Þ
5E4 GzH3 2ðEH þ 1Þ
ðEH þ 1Þ
5
þ
þ
Rez4 H 2 H X H XXX
E3 MaRezH2 HX2
5
PeðEH þ 1Þ
2CaF
(A1)
5E2 MaRezH 3 3PeðEH þ 1Þ5
5ERez4 H2 HX HXXX CaðEH þ 1Þ
CaðEH þ 1Þ5
10E2 Rez H3 HX HXXX CaðEH þ 1Þ EMaRezHHXX
5
þ
5
PeðEH þ 1Þ
2PeðEH þ 1Þ5 E4 ReH3 ðEH þ 1Þ5
þ
þ
E5 GzH2 HX2 5 ðEH þ 1Þ
þ
13E3 MaRezH3 HXX 6PeðEH þ 1Þ5
6PeðEH þ 1Þ
þ
3E2 ReH ðEH þ 1Þ5
10E2 GzHX2 5
5E4 GzH2 HXX 5
CaðEH þ 1Þ
ERe ðEH þ 1Þ5
Rez H2 HXXXX 2CaðEH þ 1Þ5
10E3 GzHX2 ðEH þ 1Þ5
2CaðEH þ 1Þ5
5E3 Rez H5 HXXXX CaðEH þ 1Þ5 5E2 GzHXX 5
ðEH þ 1Þ
5EGzHXX 5
þ
5E4 Rez H6 HXXXX
4
ðEH þ 1Þ5
4
5
2CaðEH þ 1Þ
þ
EReRem z2 H2 HXX
ðEH þ 1Þ5 H2
Rez HHX HXXX
5
5EGzHX2
þ
ðEH þ 1Þ H
E5 Rez H7 HXXXX
6PeðEH þ 1Þ5
4
2ðEH þ 1Þ5 H 2
2CaðEH þ 1Þ5
E5 MaRezH 5 HXX
ðEH þ 1Þ5
4
þ
E2 ReRem z2 H 3 HXX
GzHXX
5ERez H3 HXXXX
2ðEH þ 1Þ
CaðEH þ 1Þ5
XX 5
ðEH þ 1Þ5
ðEH þ 1Þ5 H 3 ðEH þ 1Þ
þ
3E3 ReH 2
GzHX2
CaðEH þ 1Þ5
10E3 Rez H4 HX HXXX
5E4 MaRezH4 H
5E2 MaRezH 2 HXX
E5 Rez4 H6 HX HXXX
5
4
5E4 Rez H 5 HX HXXX
þ
4
ðEH þ 1Þ5
2
5E4 GzHHX2 5
3
5E2 GzH
2ðEH þ 1Þ5 H
4
EReRem z2 H 3 H XX
6PeðEH þ 1Þ5
Gz
4
CaF6
E5 MaRezH 6
5
EMaRezH2 2PeðEH þ 1Þ
3ðEH þ 1Þ5
5E3 GzH 2
5
4
2Rez4 F2X H 3 H XX
3ðEFH þ 1Þ
5
4
þ
EMaRezFXX H 3 E2 MaRezF2X H 4 þ 3Pe F2 2 ðEFH þ 1Þ Pe F2 2 ðEFH þ 1Þ2
þ
ðEFH þ 1Þ4
EReRem z FXX H 4
5
5EGz
2
2EReRem z HHX2 2E2 ReRem z H 2 HX2 EMaRezHX2 C ¼ þ þ 5 5 ðEH þ 1Þ ðEH þ 1Þ PeðEH þ 1Þ5
2E2 ReRem z FX H 4 H X 2
E5 GzH4
E4 MaRezH5
(B1)
(B2)
2
3Rez FX H 2 H X H XX
PeðEH þ 1Þ
5
6CaðEH þ 1Þ5
EReRem z H3
þ
5
2
3E2 MaRezF2 H 2 H 2X EMaRezF2X H 3 þ 2 3FPe F 2 ðEFH þ 1Þ 2Pe F2 2 ðEFH þ 1Þ2
4
3ðEH þ 1Þ
2E3 MaRezH4
2ðEH þ 1Þ
6CaðEH þ 1Þ
2E2 ReRem z H4
þ
5
Rez4 H 3
5
2
E3 ReRem z H 5 3ðEH þ 1Þ
5ERez4
5
2
Ht¼
4
E5 Rez H 8 5ReE4 z H7 10E3 Rez H6 A ¼ 5 5 6CaðEH þ 1Þ 6CaðEH þ 1Þ 6CaðEH þ 1Þ5
2ðEH þ 1Þ H
4
5E2 Rez H4 HXXXX CaðEH þ 1Þ5
E5 GzH3 HXX 2ðEH þ 1Þ5 5E3 GzHHXX ðEH þ 1Þ5 (B3)
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
References [1] Sheludko A. Thin liquid films. Adv Colloid Interface Sci 1967;1(4):391e464. [2] Vrij A. Possible mechanism for the spontaneous rupture of thin, free liquid films. Discuss Faraday Soc 1966;42:23e33. [3] Ruckenstein E, Jain R. Spontaneous rupture of thin liquid films. J Chem Soc Faraday Trans 2 Mol Chem Phys 1974;70:132e47. [4] Williams M, Davis S. Nonlinear theory of film rupture. J Colloid Interface Sci 1982;90(1):220e8. [5] Sharma A, Ruckenstein E. An analytical nonlinear theory of thin film rupture and its application to wetting films. J Colloid Interface Sci 1986;113(2):456e79. [6] Davis S. Thermocapillary instabilities. Annu Rev Fluid Mech 1987;19:403e35. [7] Bankoff S. Stability of liquid flow down a heated inclined plane. Int J Heat Mass Transf 1971;14:377e85. [8] Burelbach J, Bankoff S, Davis S. Nonlinear stability of evaporating/condensing liquid films. J Fluid Mech 1988;195:463e94. [9] Shklyaev O, Fried E. Stability of an evaporating thin liquid film. J Fluid Mech 2007;584:157e83. [10] Benselama A, Harmand S, Sefiane K. Thermocapillary effects on steadily evaporating contact line: a perturbative local analysis. Phys Fluids 2012;24(7).
219
pp. 072105e1-38. [11] Avramenko A, Shevchuk I, Harmand S, Tyrinov A. Thermocapillary instability in an evaporating two-dimensional thin layer film. Int J Heat Mass Transf 2015;91:77e88. [12] Oron A, Davis S, Bankoff S. Long-scale evolution of thin liquid films. Rev Mod Phys 1997;69(3):931e80. [13] Craster R, Matar O. Dynamics and stability of thin liquid films. Rev Mod Phys 2009;81(3):1131e98. [14] Wayner P. Intermolecular forces in phase-change heat transfer: 1998 Kern award review. Am Inst Chem Eng J 1999;45(10):2055e68. [15] Wanyer P. The effect of interfacial mass transport on flow in thin liquid films. Colloids Surf. 1991;52:71e84. [16] Unsal M. Non-linear stability of forced vapour flow condensation. Int J Heat Mass Transf 1988;31(8):1619e26. [17] Pratt D, Brown J, Hallinan K. Thermocapillary effects on the stability of a heated, curved meniscus. ASME J Heat Transf 1998;120(1):220e6. [18] Carey VP. Liquid-Vapor Phase-Change Phenomena. Taylor & Francis; 1992. [19] Ball G, Polansky J, Kaya T. Investigation of particular features of the numerical solution of an evaporating thin film in a channel. Front Heat Mass Transf 2013;4: 1e9.
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Stability of an evaporating meniscus: Part I e Theoretical analysis John Polansky*, Tarik Kaya Carleton University, Department of Mechanical and Aerospace Engineering, Ottawa ON, Canada
a r t i c l e i n f o
a b s t r a c t
Article history: Received 19 November 2015 Received in revised form 29 March 2016 Accepted 30 March 2016 Available online 22 April 2016
The stability of an evaporating meniscus in a narrow channel is theoretically investigated. A non-linear differential equation is derived using the lubrication approximation and taking into account evaporation, long-range molecular forces, surface tension, vapor recoil, thermocapillarity and viscous forces. Using a linear stability analysis, the susceptibility of the film to instability is studied. Three different base state profiles are used: flat replenished thin film, modified Young-Laplace profile and curved thin film. The last base state most accurately reflects an evaporating meniscus; however, its implementation requires the numerical solution of the non-linear thin-film equations. A shooting technique together with a fourth order Runge-Kutta method is used to obtain the curved thin-film profiles. The results obtained for each base state are analyzed and limitations are discussed. The curved thin-film base profile predicts strong instability for higher wall superheats. The growth rates are found to be positive (destabilizing) in the thin-film region and become negative (stabilizing) towards the intrinsic meniscus. The results indicate that the instability most likely originates in the evaporating thin-film region. Part II of the presented paper is devoted to experimental investigation. © 2016 Elsevier Masson SAS. All rights reserved.
Keywords: Meniscus instability Evaporating meniscus Thin films
1. Introduction The stability of evaporating thin films are important for several technological applications such as the cooling of power plants and advanced electronics, microfluidic devices, and heat pipes. In these applications, the heat from a solid substrate is conducted through a thin liquid film to its free surface to take advantage of the high heat transfer rates associated with phase change at the liquid-vapor interface. However, if the interface equilibrium is disturbed, a potential runaway condition may lead to the eventual rupture or collapse of the thin film upon the solid wall. Instabilities leading to film rupture would constitute a failure of the device as damaging dry hot spots would result in the inability to transport the applied heat. Thus, to appropriately design and implement two-phase heat transfer devices it is important to characterize these instabilities. The stability of an evaporating liquid film is governed by a complex interaction of mechanisms spanning mass loss due to evaporation, long-range molecular forces, surface tension, vapor recoil, thermocapillarity and viscous forces. The combination of these factors is the basis from which instabilities may arise. In the case of an evaporating meniscus, the geometric constraints and associated thin-film geometry, see Fig. 1, can further complicate the
* Corresponding author. E-mail address: [email protected] (J. Polansky). http://dx.doi.org/10.1016/j.ijthermalsci.2016.03.021 1290-0729/© 2016 Elsevier Masson SAS. All rights reserved.
problem. Since the scales and associated flow rates are small, the gravitational effects and liquid inertia may be considered as negligible. To explain a potential instability mechanism of an evaporating meniscus, we will assume that a local small disturbance results in a surface depression of the liquid-vapor interface, causing a departure of the system from equilibrium. The depressed interfacial region closest to the solid wall becomes warmer and experiences a local increase in evaporation. The resulting increased vapor recoil then forces the liquid towards the wall, amplifying the disturbance. The spatial changes in interface thickness lead to spatial temperature gradients and in turn thermocapillarity by way of surface tension variation (Marangoni instability). The thermocapillarity serves to further enhance the runaway effect as warmer liquid from the depression is drawn toward the cooler regions of the interface. As the film gets thinner, the increasing disjoining pressure will also provide positive feedback. The principal stabilizing mechanism here appears to be the capillary force as it will increase with increasing deformation of the interface, thereby forcing the system back toward equilibrium. It is thus possible to reach a new equilibrium state by means of rolling cell structures with a wavy interface in a similar fashion observed nard or Be nard-Marangoni for fluid instabilities such as Rayleigh-Be cells. Alternatively, there is another potential mechanism promoting equilibrium: as the interface gets closer to the solid wall,
210
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
Nomenclature a A A ;B;C b cp C Ca D E h hfg hp H H i j J k L M Ma b n p p P Pe R Re Rem t b t T T u U v V V
Clapeyron coefficient [kg/m2s K] Dispersion constant [J] Coeff. of growth rate equation [e] Kelvin coefficient [kg/m2s Pa] Specific heat [J/kg K] Accommodation coefficient [e] Capillary number [e] Difference [e] Evaporation number [e] Film thickness [m] Latent heat [J/kg] Perturbation thickness [m] Base state film thickness [e] Dimensionless film thickness [e] Imaginary unit [e] Mass flux [kg/m2 s] Dimensionless mass flux [e] Thermal conductivity [W/m K] and disturbance wavenumber Thin-film length [m] Molar mass [kg/kmol] Marangoni number [e] Normal unit vector [e] Pressure [Pa] Pressure vector [Pa] Dimensionless pressure [e] Peclet number [e] Universal gas constant [J/mol K] Reynolds number [e] Modified Reynolds number [e] Time [s] Tangential unit vector [e] Temperature [K] Cauchy stress tensor [Pa] Velocity in x [m/s] Dimensionless velocity in x [e] Velocity in y [m/s] Velocity vector [m/s] Dimensionless velocity in y [e]
increased evaporation may cause local cooling of the depressed region thereby reducing the vapor recoil and thermocapillary effects. It is such complex interplay among these various mechanisms that needs to be studied together so as to accurately characterize the potential for meniscus instability. The problem of flat thin-film stability has been investigated through a series of theoretical studies aimed at describing the film and its propensity to destabilize. The earliest stability studies focused on long range molecular forces in non-evaporating films [1e3]. Utilizing a thermodynamic approach Sheludko [1] was able to determine that a perturbation upon a non-evaporating liquid thinfilm can destabilize due to long range molecular forces. The destabilization of the thin film was correlated with a characteristic film thickness. This finding was further supported by the findings of Vrij [2], who correlated the critical film thickness with a critical perturbation wavenumber. The characterization of critical wavenumbers was further supported by the work of Ruckenstein and Jain [3]. Following this, Williams and Davis [4] constructed a non-linear stability model for which they performed a numerical analysis for a finite perturbation. A comparison of their non-linear model with
Vl W x, y X, Y
Molar volume [m3/mol] Channel width [m] Cartesian coordinates [m] Dimensionless coordinates [e]
Greek symbols Tangent to Cartesian x angle [rad] Surface tension coefficient [N/m K] Dispersion number [e] Non-evaporating film thickness [m] Perturbation amplitude [e] Wall superheat [K] Characteristic ratio [e] Dimensionless temperature [e] L Liquid velocity coefficient [e] k Local curvature [m1] m Dynamic viscosity [Pa s] n Kinematic viscosity [m2/s] r Density [kg/m3] s Surface tension [N/m] s0 Reference surface tension [N/m] t Dimensionless time [e] f Coordinate transformation coeff. [e] F Dimensionless coord. trans. coeff. [e]
a g G d0 dH DT z q
Subscripts c Characteristic and cutoff cap Capillary d Disjoining lv Liquid-vapor interface max Maximum peak Peak value of thin film property ref Reference t Derivative in dimensional time t Derivative in dimensionless time v Vapor w Wall x Derivative in dimensional x X Derivative in dimensionless X y Derivative in dimensional y Y Derivative in dimensionless Y
that of a linear model yielded a reduction in the time to rupture in addition to shifting the required perturbation wavelength to that of smaller values. The findings of Williams and Davis [4] are further supported by Sharma and Ruckenstein [5] who also found the nonlinear effects to reduce the rupture time. Their work also included the addition of Marangoni effects by way of surfactants. The nonlinear effects attributed to Marangoni motion were found to enhance the film's stability. The consideration of Marangoni effects on thin-film stability was continued in Davis' [6] work in which thermocapillary effects were the principal focus. This analysis highlighted the thermal gradients as the source of energy for hydrodynamic instabilities, which manifest as two and three-dimensional surface waves. The influence of phase change on thin-film stability was presented by Bankhoff's [7] study of evaporating and condensing thin liquid films on an inclined plane. The results of the study found the effects of phase change to be destabilizing for evaporation and stabilizing for condensation. Combining the aforementioned influences on thin-film stability was the work of Burelbach, Davis, and Bankhoff [8]. Their mathematical formulation was constructed around a one-sided model,
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
211
stability analysis. Different thin-film base state profiles are presented for the analysis. The results and limitations of the different base states are then discussed. The more realistic curved thin-film base profile shows that the film becomes more susceptible to meniscus rupture with increasing wall superheats. The results also suggest that the instability most likely originates in the evaporating thin-film region. This research is performed in two parts, with the second part being devoted to experimental investigation, which is presented in a separate paper. 2. Mathematical formulation
Fig. 1. Schematic of an extended meniscus in a channel.
which simplified the free surface boundary condition, yielding a single non-linear differential equation. Using a perturbation analysis, they were able to determine the potential for thin-film instability. The numerical results of the analysis demonstrated the film was liable to destabilize faster for an evaporating film than for a condensing film. The conditions necessary for instability were related to the nature of the perturbation by way of a critical wavenumber as is consistent with previous works. Following this work, Shklyaev and Fried [9] developed a modified linear stability model, which removed the singularity condition found at the point of rupture as discussed by Burelbach, Davis, and Bankhoff [8]. Instead of rupture, the film was found to thin until such a point that a non-evaporating thin-film developed on the substrate. Benselama, Harmand and Sefiane [10], developed a mathematical model to determine the stability of an evaporating curved interface. The analysis included long range molecular forces, surface tension and thermocapillarity effects. Their analysis was limited to moderate evaporation rates so as to neglect vapor recoil effects. This work used a domain perturbation analysis technique, while retaining higher order terms pertaining to the lubrication problem. The results indicated the potential for instability around the non-evaporating thin-film region, dependent on the sign of the Marangoni number. Additional work by Avramenko et al. [11], considered curved thin-film stability with thermocapillarity and long range molecular forces. The analysis demonstrated changes in interfacial temperature gradients shift the location of the peak replenishing fluid velocity. The introduction of curved thin-film geometries has served to introduce a spatial component to thinfilm stability. It is this spatial aspect that may provide an indication as to the location of the instability for cases such as that of an evaporating meniscus. Comprehensive reviews on thin-film stability can be found in Oron, Davis, and Bankoff [12], as well as Craster and Matar [13]. The main objective of the present work is to develop a mathematical model to study the stability of an evaporating meniscus. For the evolution of the thin-film thickness, a non-linear differential equation is derived using the lubrication approximation. Evaporation, long-range molecular forces, surface tension, vapor recoil, thermocapillarity and viscous forces are accounted for in the derivation. The governing equations with the associated boundary conditions and scaling of these equations are explained thereafter. The susceptibility of the film to instability is studied using a linear
A liquid meniscus constrained within a rectangular channel with immiscible walls of width W and infinite depth is considered as shown in Fig. 1. Heat is applied to the sides of the rectangular channel walls. The symmetry condition at the center of the channel reduces the problem to that of a thin film of increasing thickness. A microscopic region with a flat non-evaporating thin-film near the wall merges into the evaporating thin-film region. At the far end, a macroscopic region hereafter referred to as the intrinsic meniscus is dominated by the curvature of the liquid-vapor interface. The extended meniscus profile can therefore be described in a Cartesian coordinate system with the origin positioned at the transition between planar non-evaporating thin-film region (for x < 0) and the evaporating thin-film region (for x 0). For this geb and b ometry, the normal and tangential unit vectors ( n t) can be written as:
b¼ n
1 hx f 1
1 1 b t¼ f hx
(1)
(2)
where,
f¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ h2x
(3)
2.1. The governing equations The working fluid is assumed to be incompressible, Newtonian and non-polar with its surface tension varying linearly with temperature. The channel geometry is limited to sub-millimeter sizes such that the Bond number is sufficiently small so as to ignore gravitational contributions. The fluid is assumed to be of sufficient thickness that continuum mechanics hold, while viscous dissipation effects are neglected. The variables associated with the vapor and interface are denoted with a subscript v and lv, respectively, and the subscript l for liquid phase being omitted for simplicity. Under these assumptions, the two-dimensional continuity, Navier-Stokes and the energy equations can be written as:
V$V ¼ 0
(4)
rðV t þ V$VVÞ ¼ Vp þ mV2 V
(5)
rcp ðTt þ V$VTÞ ¼ kV2 T
(6)
The boundary conditions at the wall (y ¼ 0) are imposed as a no slip, no penetration and constant temperature.
u¼0
(7)
212
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
v¼0
(8)
T ¼ Tw
(9)
The boundary conditions at the liquid-vapor interface, (y ¼ hðx; tÞ), address the jump conditions associated with evaporation. As the interface is curved, the interfacial jump conditions are b and b described in the n t coordinate system. The mass jump balance at the interface can be written as:
b ¼ rv ðV v V lv Þ$ n b j ¼ rðV V lv Þ$ n
Lastly, an evaporation model is needed to relate the evaporative mass flux to the interfacial conditions. The Kelvin-Clapeyron constitutive equation is used for this purpose: j ¼ aðTlv Tv Þ þ bðp pv Þ [15]. The second term in this equation is ignored as it was also done in [8,11] considering that this term will not affect the results significantly but introduce unnecessary complexity. The constitutive equation then becomes:
3 2sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 M Pv h M C fg 5ðTlv Tv Þ j¼4 2pRTlv RTlv Tv
(23)
(10)
Taking into account the velocity components in both x- and ydirections for a curved interface, the liquid and interface velocities can be written as:
V ¼ ½u v
(11)
V lv ¼ ½ulv vlv
(12)
Considering that the density, viscosity and thermal conductivity of the vapor phase are much smaller than those of the liquid phase, the vapor flow field is decoupled from the liquid (one sided model), as it was done in earlier related works such as Burelbach, Bankoff and Davis [8]. For the one sided model, the normal and shear stresses and energy jump conditions are then simplified to:
j2
rv
where,
T
b Þ $n b ðp pv Þ ¼ þ ðT$ n
shxx f3
þ
A h3
(24)
ulv ¼ sinðaÞht
(13)
b ÞT $b t ¼ Vs$b t ðT$ n
(25)
vlv ¼ cosðaÞht
(14)
b¼0 jhfg þ kVT$ n
(26)
a ¼ tan1 ðhx Þ
(15)
Substituting Eq. (11)e(15) into Eq. (10) and simplifying, we get the mass jump condition for a curved interface as:
rðfht þ hx u vÞ
j¼
(16)
f
The normal-stress jump condition at the interface is: T
b ðp pv Þ ¼ pcap þ pd b þ ððT T v Þ$ n b Þ $n jðV V v Þ$ n
(17)
The vapor field is assumed isotropic with a pressure equal to that of the saturated vapor pressure. The capillary pressure is proportional to the surface tension s and the local curvature k giving:
pcap ¼ sk ¼
shxx
(18)
f3
The surface tension is represented as a function of temperature as:
s ¼ s0 g Tlv Tref
(19)
The disjoining pressure is expressed as [14]:
pd ¼
2.2. Scaling of the equations The governing equations and boundary conditions are nondimensionalized using scales appropriate to the problem [8,10]. The characteristic scale in the x-direction is assigned as xc and comes as a result of subsequent scaling. The thickness of the nonevaporating thin-film layer d0 is used as the characteristic scale in the y-direction, where the non-evaporating thickness is obtained from the Kelvin-Clapeyron evaporation model as: d0 ¼ ðAbDT=aÞ1=3 . By using the continuity equation, Eq. (4), a relation can be found between the transversal vc and longitudinal uc characteristic velocities as vc =uc ¼ d0 =xc . For the transversal characteristic velocity, we use vc ¼ aDT=r where DT ¼ Tw Tv . Thus the characteristic velocity in x-direction then takes the form:
uc ¼
aDTxc
The characteristic pressure Pc can be deduced from the y-direction momentum equation as Pc ¼ muc xc =d20 . Knowing that the characteristic pressure Pc can also be obtained from Eq. (20), it is possible to deduce a relation for the longitudinal characteristic length as:
A h3
(20)
where A > 0 is the dispersion constant. A no slip shear-stress jump condition at the interface can be written as:
b Þ $b ððT T v Þ$ n t ¼ Vs$b t T
b¼0 jhfg þ ðkVT kv VTv Þ$ n
xc ¼
(21)
(22)
A naDT
1=2 (28)
The viscous scale is used for non-dimensionalizing time as follows:
tc ¼
The energy jump condition is simplified assuming kinetic energy and viscous dissipation terms to be negligible yielding:
(27)
rd0
xc d 0
n
(29)
Lastly, the temperature is scaled as the ratio of temperature differences as follows:
q¼
Tw T Tw Tv
(30)
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
Applying the characteristic scales to Eqs. (4)e(6), we obtain the non-dimensional forms of continuity, Navier-Stokes, and energy equations as:
UX þ VY ¼ 0
(31) 2
RezðUUX þ VUY Þ ¼ PX þ z UXX þ UYY 3
2
Rez ðUVX þ VVY Þ ¼ PY þ z
2
z VXX þ VYY
(33) (34)
where the dimensionless coordinates are defined as: X ¼ x=xc and Y ¼ y=d0 . The characteristic ratio (z), Reynolds (Re) and Peclet (Pe) numbers are defined by:
d0
(35)
xc
ruc d0 aDTxc ¼ m m
Re ¼
Ca ¼
G¼
Pez2 ðU qX þ V qY Þ ¼ z2 qXX þ qYY
z¼
Rem ¼
(32)
(36)
213
(47)
maDT rs0
(48)
rA m2 d0
Ma ¼
Е¼
d0 vc r n rv
(49)
gDTxc rcP mk
(50)
d0 ahfg
(51)
k
Here Rem , Ca, G, Ma, and Е are the dimensionless modified Reynolds, capillary, dispersion, Marangoni and evaporation numbers respectively. Knowing that the characteristic ratio is small, z≪1, we can reduce the governing equations, Eqs. (31)e(34), thereby obtaining the leading-order governing equations as:
UX þ VY ¼ 0
(52)
PX ¼ UYY
(53)
For the heated wall (Y ¼ 0) boundary conditions, Eqs. (7)e(9), the non-dimensional forms are then:
PY ¼ 0
(54)
U¼0
(38)
qYY ¼ 0
(55)
V ¼0
(39)
q¼0
(40)
The interfacial jump conditions and constitutive equation, Eqs. (41)e(45), are also reduced on the basis that characteristic ratio is small. The mass jump and constitutive equations, Eq. (41) and Eq. (45), remain unchanged. We selectively retain the vapor recoil, surface tension, long range molecular forces and thermocapillary contributions in the normal and shear stress boundary conditions so as to determine their effect on the stability of an evaporating thin film. Therefore, we obtain the leading-order interfacial stress and energy jump conditions as:
Pe ¼
rcp uc xc k
cp aDTx2c ¼ d0 k
(37)
The interfacial height and time are non-dimensionalized as H ¼ h=d0 and t ¼ t=tc , respectively. The jump conditions and constitutive equation, Eqs. (16), (23)e(26), at the interface, Y ¼ H ðX; tÞ, then come as:
J¼
ReðV H X UÞ FH t ReF
(41)
2
J 2 Rem z ¼
2z2 z2 H X ðVX H X UX Þ þ H X UY þ UX 2
F
ðP Pv Þ
z4 H XX Gz þ ReH 3 CaF3
J 2 Rem z2 ðP Pv Þ ¼
z4 H XX Gz þ ReH 3 CaF3
(56)
Pe F2 2 UY ¼ MazFðqX þ H X qY Þ
(57)
ЕJ F qY ¼ 0
(58)
(42) 2 PeUY F2 2 þ Pez 4H X UX þ VX F2 2
2.3. The evolution equation
¼ MazFðH X qY þ qX Þ
(43)
JЕF qY þ z H X qX ¼ 0
2
(44)
J ¼1q
(45)
where,
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F ¼ 1 þ ðzH X Þ2
The system of differential equations, Eqs. (52)e(55), together with boundary conditions described above are solved herein so as to obtain a single differential equation describing the evolution of an evaporating interface. We begin by integrating Eq. (55) twice, subject to the wall boundary and interfacial jump conditions, Eq. (40) and Eq. (58), respectively. The dimensionless temperature is then found to be:
q ¼ ЕJ FY (46)
(59)
To obtain a relation for the mass flux, we can substitute Eq. (45) into Eq. (59) and evaluate at the interface, Y ¼ H ðX; tÞ, thereby obtaining the dimensionless mass flux relation as:
214
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
2.4. Linear stability analysis
1 J¼ ЕFH þ 1
(60)
The substitution of Eq. (60) into Eq. (59) yields a new relation for the dimensionless temperature as:
q¼
Е FY Е FH þ 1
(61)
Next, we solve the Navier-Stokes equation in the Y-direction by integrating Eq. (54) and obtaining the constant of integration from the normal stress jump condition, Eq. (56), and substituting the previously obtained mass flux relation, Eq. (60), yielding:
P ¼ Pv þ
Rem z2 ðЕFH þ 1Þ
2
z4 H XX Gz ReH 3 CaF3
(62)
Given the pressure, we can now solve Navier-Stokes in the Xdirection by differentiating Eq. (62) with respect to X and substituting into Eq. (53). The resulting substitution yields a result such that the right hand side is only a function of X and t. Equation (53) is then integrated twice, subject to no slip at the wall, Eq. (38), and the shear stress jump condition, Eq. (57), at the interface. Finally, the X-component of the liquid velocity is obtained as:
2EMazFFX 1 ЕF2 H X L Y 2 þ EMazF2 H X ðЕ FH þ 1ÞY U¼ Pe F2 2 ðЕ FH þ 1Þ2 (63) where
L¼
z4 ð3FX H XX FH XXX Þ 2Rem z2 ðЕFH X þ ЕFX H Þ CaF4 ðЕFH þ 1Þ3 3GzH X þ
ReH 4
(64)
To analyze the stability of an evaporating film, we perturb the base state film thickness H by a small perturbation hp :
H ðX; tÞ ¼ HðX; tÞ þ hp ðX; tÞ
(66)
The perturbation is assumed to be of the form:
hp ðX; tÞ ¼ dHðtÞ expðikXÞ
(67)
Substituting Eq. (66) into Eq. (A1) and removing the base state by way of binomial series and simplifying, we obtain the resulting leading-order differential equation of the following form:
hpt ¼ A k4 þ B k2 þ C hp
(68)
Equation (68) describes the growth rate of a perturbation with the coefficients provided in Appendix B. By using Eq. (68), it is then possible to describe the range of disturbance wavenumbers for which a small disturbance will grow (instability) or decay (stability). 3. Analysis and discussion The results of the linear stability analysis for the thin film are presented in this section. All the results are obtained for n-pentane as the working fluid. The base profile is that of a steady-state evaporating curved thin film, which can be computed from the steady-state evolution equation, Eq. (A1). Given the non-linearity, no closed form solution is possible to be used as a base profile HðX; tÞ, making the linear stability analysis highly complex. Therefore, we first consider an evaporating flat thin film and then a modified Young-Laplace profile as the base states. Then, we consider the more realistic, though more complex, curved evaporating thin film for the base state. 3.1. Flat thin film
We can now solve the continuity equation, by differentiating Eq. (63) with respect to X and substituting into Eq. (52) and integrate once with respect to Y subject to the no penetration boundary condition, Eq. (39). Using the transverse velocity component V in conjunction with Eq. (60) and Eqs. (63) and (64), we can solve the mass jump condition, Eq. (41), at the interface, thereby obtaining the evolution equation. We note that the evolution equation contains F, which is a function of X and t, as denoted in Eq. (46). It is known that the thin film slope and characteristic ratio are sufficiently small such that, zH x ≪1 to make the following approximation:
A flat thin film base state was considered as an approximation to study the stability of an evaporating curved thin film. The quasiparallel theory was also used in previous works, for example in Unsal [16] for the stability of condensing flows and in Pratt, Brown and Hallinan [17] for the evaporating meniscus. For an evaporating thin film to be at steady state, a flow of fluid into the thin film is needed to replenish that which is lost to evaporation. If we consider an evaporating flat film of a steady height, the base state will be of the following form:
FðX; tÞy1
H¼1
(69)
HX ¼ HXX ¼ HXXX ¼ HXXXX ¼ 0
(70)
(65)
Therefore, we simplify the evolution equation, subject to Eq. (65) thereby yielding a thin-film evolution equation. Thus, we obtain a single non-linear differential equation which describes how the time rate of change of the film thickness is related to mass loss due to phase change, long-range molecular forces, surface tension, vapor recoil, thermocapillarity and viscous forces. The resulting evolution equation is sizable and presented in Appendix A. Because of the large size of the resulting equation, all parametric derivation is carried out in Maple and the resulting differential equation is solved using MATLAB.
Substituting the base state conditions, Eqs. (69) and (70), into the growth rate equation, Eq. (68), and after some manipulation, we obtain the coefficients as:
z4 Re
A ¼
6Ca
(71)
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
B ¼
zEð2PeReRem z 9PeG 3MaReÞ 6PeðE þ 1Þ3
C ¼
zE2 ð3PeGðE þ 3Þ þ MaReðE þ 4ÞÞ 6PeðE þ 1Þ3
3PeGz 6PeðE þ 1Þ3
ERe ðЕ þ 1Þ2
H¼ (72)
(73)
Here we note that the stability of the meniscus is not spatially dependent as all of the base state derivatives with respect to X are zero, thereby leaving the solution as a function of the perturbation (disturbance) wavenumber k, and the applied superheat DT. Fig. 2 shows a plot of the growth rate as a function of the disturbance wavenumber for a range of wall superheats.
Fig. 2. Growth rate as a function of disturbance wavenumber.
The negative growth rates indicate stability. As the superheat increases, the range of wavenumbers capable of destabilizing the flat film decreases. In the unstable region, there is a cross over after which the higher wall superheats lead to larger growth rates. For a given wall superheat, the growth rate continuously increases toward a maximum as k/0.
W 2
215
ffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 W þ d0 ðXxc Þ 2
d0
(74)
Substituting the base state profile, Eq. (74), into the growth rate equation, Eq. (68), we obtain the fourth-order polynomial growth rate equation for the modified Young-Laplace base case. As the modified Young-Laplace base state is a function of position and is continuous and differentiable, the stability analysis becomes a function of position, unlike that of a flat thin-film case. The resulting growth rate equation is notably more complex than that of the flat film case and not reproduced here as it is comprised of 875 operands. For the modified Young-Laplace base state we plot the growth rate for a typical case, DT ¼ 0:05, W ¼ 10 mm and k ¼ 0:5; as it is shown in Fig. 3. It should be noted that in the proximity of the symmetry point at x ¼ W=2, HX /∞, and the assumption zHX ≪1 breaks down. For the typical case presented in Fig. 3, the solutions should be sufficiently reliable for the meniscus positions less than 4.5 mm where zHX is of an order of magnitude of 102 or smaller. It can be seen from Fig. 3 that the film is increasingly stable as one moves away from the evaporating thin-film region towards the intrinsic meniscus. The positive growth rates observed in the evaporating thin-film region indicates instability as illustrated in Fig. 3. The growth rate rapidly changes sign as a function of the position in the evaporating thin-film region. As the wavenumber decreases asymptotically toward zero, k/0, the scatter in the thin-film region reduces. For a range of wavenumbers less than some critical value, take for example kc ¼ 0:5, the growth-rate curve becomes smooth without scatter and the thin-film remains stable for DT < 105 K in a channel of W ¼ 10 mm. For the same critical wavenumber, as the channel width increases, the interfacial temperature jump at which the instability is first observed also decreases suggesting that the large channels are more susceptible to instability. While the growth rates increase with increasing wavenumber, the intrinsic meniscus region remains generally unaffected. Thus, the growthrate analysis indicates an inherent stability in the intrinsic meniscus region while the evaporating thin-film region is susceptible to instability. Considering the local heat transfer and evaporation rates peak in the evaporating thin-film region [18], the finding that the thin-film region could be the origin of instability is not surprising. This point will be further discussed later.
3.2. Modified Young-Laplace profile
3.3. Curved thin film
A base state described by the Young-Laplace equation is of interest as it has a closed form and is easy to implement for the linear stability analysis. The modified Young-Laplace equation represents a more realistic base state than that of a flat thin film. The modified Young-Laplace base state is described by the segment of a circle, centered on the symmetry plane of the channel and reduced in width by twice the non-evaporating thinfilm thickness. While the non-evaporating thin-film thickness is a function of the wall superheat, its scale as compared to the channel width places emphasis on the geometry over the applied superheat. Thus, we define the modified Young-Laplace base state as:
As there is no closed form solution for the steady-state evaporating curved thin-film base state, the governing equations need to be solved numerically to obtain the thin-film thickness and the associated derivatives at discrete locations along the extended meniscus. Then, these values can be substituted into the growthrate equation, Eq. (68), at each corresponding location to obtain the growth rate along the entire meniscus. Note that the base state is timeeindependent as the evaporated fluid will be replenished by the flow in the capillary channel. To obtain the thin film profile, the steady-state form of the evolution equation, Eq. (A1), needs to be solved. Using a slightly modified version of Eq. (A1), we first solve the governing non-
216
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
Fig. 3. Growth rate as function of position for the modified Young-Laplace base state, DT ¼ 0:05 K, W ¼ 10 mm and k ¼ 0:5.
linear equations using a fourth order Runge-Kutta method. Solving the evaporating thin-film equation requires a set of initial conditions to achieve the target curvature of 2=W at the far field with prescribed constant curvature tolerance of ± 0:02=W. As the sensitivity of the non-linear equation lies with the first derivative, hx , a shooting technique is employed to find the required initial slope to obtain the desired far field boundary condition. The results are found to be extremely sensitive to the initial conditions. The details of the solution methodology were presented in Ball,
Polansky and Kaya [19]. As the thin-film profile is a function of wall superheat, the profiles must be obtained over a range of wall superheats for a given channel. To ensure the validity of the stability model in the thin-film region, the slopes of the thin-film profiles are checked against the assumption: zHX ≪1. For the cases tested, we found the assumption to hold. The obtained thin-film thickness and its derivatives are shown in Fig. 4. Given the strong dependence on superheat, the thin-film curvature profiles transition from an asymptotic rise to 2=W for small superheats to an overshoot before achieving the far field curvature as illustrated in Fig. 4(b). The location of the overshoot peak is noted to shift spatially toward the non-evaporating thinfilm region with increasing superheats. The higher order derivatives also change sign in response to the change in curvature as shown in Fig. 4(c) and (d). These changes in sign on the higher derivative values along the thin film directly affect the growth rate results through Eq. (68). The resulting disturbance cutoff wavenumbers kc obtained using the curved thin-film solutions are shown in Fig. 5 for a range of wall superheats and meniscus locations in a channel 100 mm wide. The cutoff wavenumber is calculated by assigning the growth rate to zero, hpt =hp ¼ 0. As the thin-film solution is dependent upon the fluid properties, superheat and position so too is the cutoff wavenumber. The broken appearance of the surface as the cutoff wavenumber approaches to zero, kc /0; is a consequence of the sharp gradients coupled with the numerical resolution of thin-film solutions in this region. The stable region corresponds to the negative growth values and occupies the top portion of the curved plane.
Fig. 4. (a) Thin film thickness (b) local curvature (c) third derivative (d) fourth derivative for wall superheats ranging from 0.01 to 0.20 K, W ¼ 100 mm.
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
217
Looking at the most critical characteristics, we consider the locations where the maximum evaporative mass flux, interface temperature gradient, liquid pressure gradient and thin-film curvature occur in the base case profile. Fig. 7 shows the difference between the locations where the growth rate peaks ðxpeak Þ and a thin-film characteristic reaches a maximum ðxmax Þ, normalized by the length of the extended thin-film ðL Þ, as D ¼ ðxpeak xmax Þ=L . The local curvature has the smallest degree of offset and the differences continuously decrease with increasing wall superheat. The local curvature appears to be the most influential parameter. However, the differences shown in Fig. 7 are insufficiently varied to draw a sound conclusion. It is clear that all these thin-film characteristics reach a maximum in the evaporating thin-film region. As it was stated previously, the evaporating thin-film region seems to be the most likely region of the onset of instabilities.
Fig. 5. The disturbance cutoff wavenumber as a function of wall superheat and meniscus location, W ¼ 100 mm.
4. Conclusions
The cutoff wavenumber surface illustrated in Fig. 5 suggests that there are two regions of an evaporating thin film with a positive growth rate. The first occurs nearest the non-evaporating thin-film region, and the second occurs for a small portion of the thin film before tending back to negative growth rates towards the intrinsic meniscus region. The shape of the curve suggests that there is a peak cutoff wavenumber above which meniscus is always stable. Fig. 6 shows the growth rates for a range of wall superheats (0.01e0.20 K) for a channel width of 100 mm and a disturbance wavenumber of 0.15. The growth rates for the thin-film base case are observed to be smooth and continuous as compared to the sporadic nature found in the modified Young-Laplace base case. As the applied superheat increases, so too does the peak value of the growth rate. This suggests that higher superheats are associated with strong instability and will drive the meniscus to rupture faster than lower superheats for a given perturbation. Similar to the modified YoungLaplace base case, the growth rate starts off positive (destabilizing) in the thin film region with a peak positive growth rate and transitions to a negative (stabilizing) further along the meniscus. The position of the maximum positive growth rate shifts towards the nonevaporating thin-film region with increasing superheats. The location where the growth rate peaks should drive the instability. Thus, it would be interesting to investigate what other thin-film characteristics reach a maximum around this location.
The stability of an evaporating meniscus is studied using a linear stability analysis. For the base state, three different profiles are considered; a flat replenished thin film, modified Young-Laplace profile and a curved thin film. The flat thin-film analysis shows that the range of destabilizing wavenumbers decreases for increasing wall superheats. For the modified Young-Laplace base profile, the instability is principally driven by the position along the meniscus and the perturbation wavenumber. Though the solution is subject to the assumption zHX ≪1, the validity extends well into the intrinsic meniscus where the growth rate is found to be negative thereby indicating strong stability. The growth rate analysis also indicates that the evaporating thin-film region is susceptible to instability, where the local heat transfer and evaporation rates peak. The first two base profiles have the advantage of leading to the closed form solutions; however, they are not a very accurate representation of an evaporating meniscus. For a more accurate analysis, the steady-state thin-film profiles are numerically solved by using a fourth order Runge-Kutta method. It is found that there exists a cutoff wavenumber above which the meniscus always remains stable. The results predict strong instability for higher superheats. The growth rates are found to be positive (destabilizing) in the thin film region where they reach a maximum and become negative (stabilizing) towards the intrinsic meniscus. As the superheat increases, the position of the maximum positive growth rate shifts towards the non-evaporating thin-film region. No direct
Fig. 6. Growth rates calculated as a function of meniscus location for a range of wall superheats; W ¼ 100 mm and k ¼ 0:15.
Fig. 7. Comparison of the locations where the growth rate and thin-film characteristics peak for a channel of W ¼ 100 mm.
218
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
correlation is found between the location of the peak growth rate and individual thin-film characteristics. The results obtained using modified Young-Laplace and curved thin-film base profiles both indicate that the evaporating thin-film region is the most likely origin of the instability.
Appendix B. The perturbation growth rate equation coefficients
4
4
Acknowledgements
The support from Natural Sciences and Engineering Research Council of Canada under Discovery Grants Program is greatly appreciated.
5E2 Rez4 H5 3CaðEH þ 1Þ
B ¼
Appendix A. The thin film evolution equation
19E2 MaRezFFX H 3 H X 4E3 MaRezF2 FX H 4 H X 2 2 6Pe F 2 ðEFH þ 1Þ 3Pe F2 2 ðEFH þ 1Þ3
4
þ
2E2 MaRezF3 FX H 3 H X Rez H 3 H XXXX 2 6CaF4 3Pe F2 2 ðEFH þ 1Þ2 E2 MaRezF2 H 3 H
2
XX 2
3Pe F 2 ðEFH þ 1Þ
2H
EMaRezFH XX 2Pe F2 2 ðEFH þ 1Þ
þ
EMaRezFH H 2XX E2 MaRezFFXX H 4 2 2 3Pe F 2 ðEFH þ 1Þ Pe F2 2 ðEFH þ 1Þ
þ
2E3 MaRezF3 H 3 H 2X 2EMaRezFF2X H 3 2 2 3 3Pe F2 2 ðEFH þ 1Þ 3Pe F 2 ðEFH þ 1Þ
2E2 MaRezF2 F2X H 4 2E3 MaRezFF2X H 5 2 2 3 3Pe F2 2 ðEFH þ 1Þ2 3Pe F 2 ðEFH þ 1Þ 2EMaRezFX H 2 H X EMaRezF2 FX H 2 H X þ 2 2 Pe F 2 ðEFH þ 1Þ Pe F2 2 ðEFH þ 1Þ
þ
2
þ
þ
2CaF5 E2 ReRem z F2X H 5
FðEFH þ 1Þ
4
3
3FðEFH þ 1Þ
þ
GzH XX GzH 2X 2FH 2FH 2
2
2
5EReRem z FX H 3 H X 3
3FðEFH þ 1Þ
2
þ þ þ
E2 ReRem z FH 3 H 2X ðEFH þ 1Þ4 Rez4 FX H 3 H XXX 5
CaF
Rez4 FXX H 3 H XX 2CaF
5
2H 2 X 3
EReRem z H Re E FH þ 1 ðEFH þ 1Þ
5
2ðEH þ 1Þ
PeðEH þ 1Þ
5E4 GzH3 2ðEH þ 1Þ
ðEH þ 1Þ
5
þ
þ
Rez4 H 2 H X H XXX
E3 MaRezH2 HX2
5
PeðEH þ 1Þ
2CaF
(A1)
5E2 MaRezH 3 3PeðEH þ 1Þ5
5ERez4 H2 HX HXXX CaðEH þ 1Þ
CaðEH þ 1Þ5
10E2 Rez H3 HX HXXX CaðEH þ 1Þ EMaRezHHXX
5
þ
5
PeðEH þ 1Þ
2PeðEH þ 1Þ5 E4 ReH3 ðEH þ 1Þ5
þ
þ
E5 GzH2 HX2 5 ðEH þ 1Þ
þ
13E3 MaRezH3 HXX 6PeðEH þ 1Þ5
6PeðEH þ 1Þ
þ
3E2 ReH ðEH þ 1Þ5
10E2 GzHX2 5
5E4 GzH2 HXX 5
CaðEH þ 1Þ
ERe ðEH þ 1Þ5
Rez H2 HXXXX 2CaðEH þ 1Þ5
10E3 GzHX2 ðEH þ 1Þ5
2CaðEH þ 1Þ5
5E3 Rez H5 HXXXX CaðEH þ 1Þ5 5E2 GzHXX 5
ðEH þ 1Þ
5EGzHXX 5
þ
5E4 Rez H6 HXXXX
4
ðEH þ 1Þ5
4
5
2CaðEH þ 1Þ
þ
EReRem z2 H2 HXX
ðEH þ 1Þ5 H2
Rez HHX HXXX
5
5EGzHX2
þ
ðEH þ 1Þ H
E5 Rez H7 HXXXX
6PeðEH þ 1Þ5
4
2ðEH þ 1Þ5 H 2
2CaðEH þ 1Þ5
E5 MaRezH 5 HXX
ðEH þ 1Þ5
4
þ
E2 ReRem z2 H 3 HXX
GzHXX
5ERez H3 HXXXX
2ðEH þ 1Þ
CaðEH þ 1Þ5
XX 5
ðEH þ 1Þ5
ðEH þ 1Þ5 H 3 ðEH þ 1Þ
þ
3E3 ReH 2
GzHX2
CaðEH þ 1Þ5
10E3 Rez H4 HX HXXX
5E4 MaRezH4 H
5E2 MaRezH 2 HXX
E5 Rez4 H6 HX HXXX
5
4
5E4 Rez H 5 HX HXXX
þ
4
ðEH þ 1Þ5
2
5E4 GzHHX2 5
3
5E2 GzH
2ðEH þ 1Þ5 H
4
EReRem z2 H 3 H XX
6PeðEH þ 1Þ5
Gz
4
CaF6
E5 MaRezH 6
5
EMaRezH2 2PeðEH þ 1Þ
3ðEH þ 1Þ5
5E3 GzH 2
5
4
2Rez4 F2X H 3 H XX
3ðEFH þ 1Þ
5
4
þ
EMaRezFXX H 3 E2 MaRezF2X H 4 þ 3Pe F2 2 ðEFH þ 1Þ Pe F2 2 ðEFH þ 1Þ2
þ
ðEFH þ 1Þ4
EReRem z FXX H 4
5
5EGz
2
2EReRem z HHX2 2E2 ReRem z H 2 HX2 EMaRezHX2 C ¼ þ þ 5 5 ðEH þ 1Þ ðEH þ 1Þ PeðEH þ 1Þ5
2E2 ReRem z FX H 4 H X 2
E5 GzH4
E4 MaRezH5
(B1)
(B2)
2
3Rez FX H 2 H X H XX
PeðEH þ 1Þ
5
6CaðEH þ 1Þ5
EReRem z H3
þ
5
2
3E2 MaRezF2 H 2 H 2X EMaRezF2X H 3 þ 2 3FPe F 2 ðEFH þ 1Þ 2Pe F2 2 ðEFH þ 1Þ2
4
3ðEH þ 1Þ
2E3 MaRezH4
2ðEH þ 1Þ
6CaðEH þ 1Þ
2E2 ReRem z H4
þ
5
Rez4 H 3
5
2
E3 ReRem z H 5 3ðEH þ 1Þ
5ERez4
5
2
Ht¼
4
E5 Rez H 8 5ReE4 z H7 10E3 Rez H6 A ¼ 5 5 6CaðEH þ 1Þ 6CaðEH þ 1Þ 6CaðEH þ 1Þ5
2ðEH þ 1Þ H
4
5E2 Rez H4 HXXXX CaðEH þ 1Þ5
E5 GzH3 HXX 2ðEH þ 1Þ5 5E3 GzHHXX ðEH þ 1Þ5 (B3)
J. Polansky, T. Kaya / International Journal of Thermal Sciences 107 (2016) 209e219
References [1] Sheludko A. Thin liquid films. Adv Colloid Interface Sci 1967;1(4):391e464. [2] Vrij A. Possible mechanism for the spontaneous rupture of thin, free liquid films. Discuss Faraday Soc 1966;42:23e33. [3] Ruckenstein E, Jain R. Spontaneous rupture of thin liquid films. J Chem Soc Faraday Trans 2 Mol Chem Phys 1974;70:132e47. [4] Williams M, Davis S. Nonlinear theory of film rupture. J Colloid Interface Sci 1982;90(1):220e8. [5] Sharma A, Ruckenstein E. An analytical nonlinear theory of thin film rupture and its application to wetting films. J Colloid Interface Sci 1986;113(2):456e79. [6] Davis S. Thermocapillary instabilities. Annu Rev Fluid Mech 1987;19:403e35. [7] Bankoff S. Stability of liquid flow down a heated inclined plane. Int J Heat Mass Transf 1971;14:377e85. [8] Burelbach J, Bankoff S, Davis S. Nonlinear stability of evaporating/condensing liquid films. J Fluid Mech 1988;195:463e94. [9] Shklyaev O, Fried E. Stability of an evaporating thin liquid film. J Fluid Mech 2007;584:157e83. [10] Benselama A, Harmand S, Sefiane K. Thermocapillary effects on steadily evaporating contact line: a perturbative local analysis. Phys Fluids 2012;24(7).
219
pp. 072105e1-38. [11] Avramenko A, Shevchuk I, Harmand S, Tyrinov A. Thermocapillary instability in an evaporating two-dimensional thin layer film. Int J Heat Mass Transf 2015;91:77e88. [12] Oron A, Davis S, Bankoff S. Long-scale evolution of thin liquid films. Rev Mod Phys 1997;69(3):931e80. [13] Craster R, Matar O. Dynamics and stability of thin liquid films. Rev Mod Phys 2009;81(3):1131e98. [14] Wayner P. Intermolecular forces in phase-change heat transfer: 1998 Kern award review. Am Inst Chem Eng J 1999;45(10):2055e68. [15] Wanyer P. The effect of interfacial mass transport on flow in thin liquid films. Colloids Surf. 1991;52:71e84. [16] Unsal M. Non-linear stability of forced vapour flow condensation. Int J Heat Mass Transf 1988;31(8):1619e26. [17] Pratt D, Brown J, Hallinan K. Thermocapillary effects on the stability of a heated, curved meniscus. ASME J Heat Transf 1998;120(1):220e6. [18] Carey VP. Liquid-Vapor Phase-Change Phenomena. Taylor & Francis; 1992. [19] Ball G, Polansky J, Kaya T. Investigation of particular features of the numerical solution of an evaporating thin film in a channel. Front Heat Mass Transf 2013;4: 1e9.