- Email: [email protected]
Simulation and validation of the dynamics of liquid films evaporating on horizontal heater surfaces
Applied Thermal Engineering 48 (2012) 486e494
Contents lists available at SciVerse ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Simulation and validation of the dynamics of liquid films evaporating on horizontal heater surfaces Shengjie Gong, Weimin Ma*, Truc-Nam Dinh 1 Department of Physics, Royal Institute of Technology (KTH), Roslagstullsbacken 21, 106 91 Stockholm, Sweden
h i g h l i g h t s < Liquid film dynamics under thermal influence is simulated. < Long wave theory and minimum free energy theory are applied. < Impacts of heat flux and surface wettability on film dynamics are investigated. < A correlation is developed to predict the critical thickness of film rupture.
a r t i c l e i n f o
a b s t r a c t
Article history: Received 11 August 2011 Accepted 13 May 2012 Available online 23 May 2012
In this study a non-linear governing equation based on lubrication theory is employed to model the thinning process of an evaporating liquid film and ultimately predict the critical thickness of the film rupture under impacts of various forces resulting from mass loss, surface tension, gravity, vapor recoil and thermo-capillary. It is found that the thinning process in the experiment is well reproduced by the simulation. The film rupture is caught by the simulation as well, but it underestimates the measured critical thickness at the film rupture. The reason may be that the water wettability of the heater surfaces is not taken into account in the model. Thus, the minimum free energy criterion is used to obtain a correlation which combines the contact angle (reflection of wettability) with the critical thickness from the simulation. The critical thicknesses predicted by the correlation have a good agreement with the experimental data (the maximum deviation is less than 10%). Ó 2012 Elsevier Ltd. All rights reserved.
Keywords: Liquid film Evaporation Rupture Critical thickness
1. Introduction Liquid films spreading over and evaporating on solid surfaces are commonly encountered in daily life and industrial applications that involve processes such as cooling (e.g., condensation), heating (e.g., boiling), coating, cleaning and lubrication. Studies on liquid film dynamics are of paramount importance to achieve an efficient and safe operation of the related equipments. For instance, the instability of a liquid film will affect the coating quality, while the rupture of a near-wall liquid film (e.g., macrolayer in boiling) will deteriorate heat transfer and lead to safety issues (destruction of heater surface due to burnout) in power plants. In our previous work, a novel experimental method using a confocal optical sensor was developed to measure the thickness evolution of an evaporating liquid film, and first-of-a-kind data about the thinning and rupture process of a thin liquid film on * Corresponding author. Tel.: þ46 8 5537 8821; fax: þ46 8 5537 8830. E-mail address: [email protected] (W. Ma). 1 Present address: Idaho National Laboratory, Idaho Falls, ID 83415, USA. 1359-4311/$ e see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.applthermaleng.2012.05.021
a horizontal surface heated from below was obtained [1,2]. It was observed that a liquid film ruptures when it was thinning to a threshold value (called critical thickness). The critical thickness of the liquid films performed in the tests was found to increase with increasing heat flux, and within the range of 60e150 mm for water films and 150e250 mm for ethanol films. In addition, the film instability was greatly affected by the properties of heater surfaces. Wettability of the surfaces proves to have a profound influence on film stability, and the critical thickness decreases with the reduction in contact angle. The wettability (contact angle) alone, however, could not explain all the effects of material properties and surface conditions on film instability. Thermalephysical properties of the heater surfaces also manifest themselves through the onset of liquid film instability. For instance, given two surfaces with a comparable contact angle, the one (e.g., copper or aluminum) with a high thermal conductivity exhibited a thinner critical film than the one with a low thermal conductivity (e.g., stainless steel or titanium). For the essentially wetting case (contact angle near zero) of a hexane film formed on an aged titanium surface, the liquid film was found resilient to rupture, and instability was not detected
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
until the thickness was less than 15 mm at low heat flux. No instability was observed for high-evaporation-rate cases. This paper is regarding a companion work to the experimental study, with the focus on theoretical analysis of film dynamics. The objectives are two-fold: The high-quality data obtained in the experiment serve a rich benchmark to assess analytical models of film instability, while the interpretation of the experiment and further understanding of the phenomena require a comprehensive analysis which can discern the effects of competing forces that govern the film behavior at diminishing thickness and rupture. Over the past three decades, a large number of theoretical studies have been carried out to investigate the film instability. Notably, the minimum thickness of liquid film flowing down a vertical or inclined adiabatic solid surface has been predicted theoretically on basis of minimum total energy criteria [3e6], resulting in a reasonable agreement with the limited experimental data available under isothermal conditions. For the liquid films on horizontal solid surfaces, Sharma et al. [7,8] developed a theory to predict the critical thickness at which liquid film rupture occurs by hole formation when the free energy of the film-solid system becomes equal to the free energy of the hole-liquid-solid system. The model predicted film rupture thickness is several hundred micrometers which are in good agreement with the experimental data obtained for liquid films on selected non-wetting solid surfaces (e.g., Teflon, polyethylene, PMMA, wax). The limitations of those models are that they are neither inclusive of all key instability sources in the model nor able to capture the dynamic process of film rupture. The long wave theory [9,10] is capable to predict film dynamics and rupture process. Based on long wave theory Tan et al. [11] theoretically analyzed steady thermo-capillary flows of thin liquid layers without consideration of evaporating mass loss. The predicted results have a good agreement with the experimental data [12] obtained from thermo-capillary flows of silicon oil (nonvolatile liquid) film with thickness from 0.125 mm to 1.684 mm on a heated horizontal plate. Oron et al. [13] provided perhaps the most comprehensive review on the multifaceted subject of thin film dynamics modeling. They presented a unified mathematical system to predict the long-scale evolution of thin liquid films based on the long wave theory. The set of mathematical evolution equations has its root in the work of [10], taking into account the influential factors such as van der Waal forces, surface tension, gravity, thermo-capillary, mass loss and vapor recoil force. In their review, Oron et al. emphasized the importance of experimentally assessing and validating models derived from the long wave theory. Craster & Matar [14] presented a comprehensive review of the work carried out on thin films flows, focusing attention on the studies undertaken after the review by Oron et al. [13]. They pointed out that for the modeling of thin film dynamics the lubrication approximation was still a useful approach to elucidate a wide variety of flows in which films have small aspect ratios. Simulations of lattice Boltzmann method (LBM) and molecular dynamics (MD) also see their role in micro fluid dynamics. While the application of the MD simulation is still limited to ultra-thin liquid films [15], the LBM simulation incorporates microscopic physics into macroscopic model, and is able to simulating twophase flow near heater walls and liquid (droplet) behavior on partial wetting surfaces [16]. Yan et al. also gave an overview of the recent development of the lattice Boltzmann method for modeling engineering encountered single-phase and multiphase flow problems to demonstrate that the LBM is a useful tool for mesoscale modeling [16]. In the present study, the non-linear stability analysis approach developed by Burelbach et al. [10] based on long wave principle is employed to simulate the liquid film evolution on a heated planar
487
solid surface. Due to exclusion of contact angle in the formulation, the influence of wettability on film dynamics is not (at least not explicitly) modeled in the long wave theory. To account for this, the predicted critical thickness of the liquid is corrected by a correlation from the minimum total free energy criterion [7] which relates the critical thickness with contact angle. In other words, the present study is trying to take the advantages of both the long wave theory and the free energy theory, so that the modeling approach does not only consider all key instability sources, but also capture the dynamic process of film rupture. 2. Mathematical formulation 2.1. Governing equation of an evaporating liquid film A liquid film evaporating on the horizontal surface of a planar heater is as illustrated in Fig. 1 where the liquid film is bounded above by its vapor and below by the heater surface, but unbounded in the horizontal direction. To describe the evolution and dynamics of the liquid film, two-dimensional Cartesian coordinates (x-z) are used to describe the system. The variation of liquid film in the direction perpendicular to the x-z plane is assumed negligible, since the test section of the experiment [2] to be simulated here was designed in such a way that the film dynamics in the y-direction (perpendicular to the x-z plane) is suppressed by a short distance between two side walls, and thus the interfacial wave propagation is dominant in the x-direction. The heater surface is located at z ¼ 0 and the vaporeliquid interface at z ¼ h(x,t) where z designates vertical coordinate and h the film thickness that is the function of horizontal coordinate x and time t. The liquid is assumed to be a Newtonian fluid of constant thermo-physical properties (kinematic viscosity n, density r and thermal conductivity k). This is a free-boundary problem. In order to remove the complexity of treating the free boundaries (interfaces between liquid and gas/vapor) in the governing equations, one can resort to the long wave theory if variations along the liquid film are much more gradual than those normal to it [13]. The theory is applied in a variety of areas in classical physics, such as lubrication theory in viscous liquid films. In the lubrication theory or long-wave-theory approach an asymptotic reduction of the governing equations and boundary conditions was implemented [10], and consequently a single nonlinear partial differential equation formulated in terms of the local thickness of the liquid film can be obtained:
H E Hs þA X þS H3 HXXX G H 3 HX þ X X H þK H X 3 2 E2 H KM H þ HX þ HX ¼ 0 D H þK Pr HþK X X
(1)
where X and H are the dimensionless coordinate and the film thickness by scaling x and h to the initial averaged film thickness of
Fig. 1. The sketch of an evaporating liquid film on a heater surface.
488
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
h0; s is non-dimensional time related to the viscous time tn ¼ h20 =n; A is the scaled dimensionless Hamaker constant; S is the dimensionless surface tension constant; G is the dimensionless gravity constant; E is evaporation number which is the ratio of the viscous timescale to evaporative timescale; K is a constant measuring the degree of non-equilibrium at the evaporating interface; D is proportional to the ratio of vapor to liquid densities; M is Marangoni number; and Pr is Prandtl number. All these dimensionless parameters can be expressed as follows:
X ¼
x ; h0
s h S ¼ 0 20 ; 3rn
H ¼
h ; h0
s¼
3 rv D ¼ ; 2 r
t ; tv
A ¼
A0 ; 6ph0 rn2
gDTh0 M ¼ ; 2rnk
n Pr ¼ ; k
E ¼
configuration with a hole (Fig. 2b) is less than that of the flat film (Fig. 2a). Otherwise, the liquid film will rewet the surface and its integrity is kept. Thus, the threshold thickness for initiating liquid film instability can be obtained if the free energy of the liquid film in Fig. 2a is equal to that of the configuration in Fig. 2b:
s Al pr22 þ ssl Al pr12 þ Ah s þ pr12 ssa Al s Al ssl ¼ 0 (4)
tn kDT ¼ rnL tE 3
K ¼
kTs2 bh0 rv L2
In the above expressions, s0 is the surface tension at the saturation temperature Ts, A0 is Hamaker constant, b is the accommodation coefficient for evaporation with the value of w1 for water at 100 C [17], Rg is universal gas constant, Mw is the liquid molecular weight, and evaporative timescale tE ¼ rh20 L=kDT where L is the latent heat and DT ¼ TiTs is the temperature difference between the heater surface temperature Ti and the saturation temperature Ts. The value of M is related to
ds g ¼ dT
(2)
!
2pRg MW
(3)
where C0 ¼ 0.07583 N/m and C1 ¼ 1.477 104 N/(m C) for water [18]. The second term on the left-hand side of Eq. (1) is representing the van der Waals molecular force which is negligible if the liquid thickness is at the magnitude more than some micrometers [13]. The sixth term is concerning the non-dimensional recoil force, and the last term (with the product of K, M and Pr1) is the nondimensional thermo-capillary force. Periodic boundary conditions are assumed for Eq. (1). 2.2. Minimum free energy criterion for liquid film rupture According to minimum free energy principle, the horizontal liquid film will be unstable and rupture if the free energy of the
2
Applying the Young’s equation ssa ¼ scos q þ ssl, one can have
s Ah þ pr12 cosq pr22 ¼ 0
(5)
where Al is the liquid film area and Ah is the hole surface area. If the external energy (due to forces) used to create the hole with radius of r2 is assumed as,
DE ¼ 2pr2 hc s
(6)
the critical (threshold) thickness for liquid film rupture is resolved as [7]
The surface tension s is a linear function of the liquid temperature T (in C) as the following:
s ¼ C0 C1 $T
1
hc ¼
sffiffiffiffiffiffiffiffiffirffiffiffiffiffiffiffiffiffi f ðqÞ DE 4ðqÞ 2ps
ð0 < q < p=2Þ
(7)
where q is the contact angle, and
f ðqÞ ¼
hc r1
¼
" # 2 1 ð1 cosqÞ sinq 1 þ 2 2 esin2 q þ ð1 cosqÞ
4ðqÞ ¼ coshðf ðqÞ=sinqÞ þ cosqsinhðf ðqÞ=sinqÞ
(8)
(9)
This indicates p the critical thickness of liquid film rupture is ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi proportional to f ðqÞ=4ðqÞ, which represents the impact of wettability on liquid film instability. Since the above-mentioned lubrication theory accounts for the effects of various external forces on film dynamics and instability, Eq. (7) can be rewritten as
Fig. 2. Schematic of hole formation in a liquid film (0 < q < p/2).
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
sffiffiffiffiffiffiffiffiffi f ðqÞ $h hc ¼ C 4ðqÞ c;cal
(10)
where hc,cal is the critical liquid film thickness for rupture predicted by the lubrication theory, and C is a constant to be determined from the experimental data. 3. Experimental setup and its modeling
489
The first step is to see how the different forces affect the liquid film instability under uniform heat flux over the heater surface. The approach is then used to predict the critical thickness of the liquid film under the experimental conditions [2]. Finally, the minimum energy criterion is adopted to account for wettability effect on the critical thickness, and a correlation is developed to predict the critical thickness of an evaporating liquid film under non-uniform temperature profile of heater surface and various wettability conditions.
3.1. Test section As shown in Fig. 3, the heater surface is a rectangular metal sheet mounted on a heating system. The metal sheet is 30 mm long, 20 mm wide and 0.42 mm thick. The materials of the sheet selected in the present study are silicon wafer, copper, aluminum, stainless steel and titanium. The metallic sheets are glued to the copper block using silver epoxy to minimize thermal resistance, with the contact region of 10 mm 15 mm. Cartridge heaters are embedded in the copper block and powered by the programmable DC power supply. Copper temperature under the metallic sheet is monitored by thermocouples. A glass vessel (60 mm 20 mm) is designed to surround the test surface, to supply water and avoid edge effect due to capillary. A stereo microscope is employed to visualize film rupture. A confocal optical sensor with the measuring range of 300 mm and the nominal resolution of 0.012 mm is used to measure the dynamic thickness of the liquid film. More details of the experimental setup can be found in [2]. 3.2. Simulation model Based on the geometry as shown in Fig. 3, it is assumed that the liquid film is symmetric around the centerline perpendicular to the heater surface. Thus, the question is simplified as a twodimensional problem as shown in Fig. 4. Water film is considered in the present study. In the lubrication theory, the length of the liquid film on the metallic surface is defined to be one wave length l ¼ 30 mm, and the ratio of the initial thickness to the length of the liquid film is h0/l ¼ 0.01 which is small enough to make the long wave theory applicable. The heat flux over the heating region is obtained through evaporative mass loss rate, and further verified by balance between power input and heat loss. 4. Results and discussions This section presents the dynamics of the evaporating liquid film using the lubrication theory approach represented by Eq. (1).
4.1. Liquid films on heater surfaces with uniform heat flux In this study it is assumed that the entire heater surface (30 mm 20 mm) has a uniform (identical) heat flux, i.e., to ignore the localized actual heating region (10 mm 15 mm) in the experiment. The purpose is to have ideas on the governing forces affecting liquid film instability at given conditions. According to Eq. (1), gravity force mainly helps to stabilize the liquid film, while high order differential surface tension force contributes to smooth the liquid film surface. The other forces (evaporation mass loss item, evaporation recoil force and thermocapillary force) will destabilize the liquid film. In order to check the effect of these forces on liquid film evolution, a liquid film with initial averaged thickness of h0 ¼ 300 mm, 30 mm and 3 mm are analyzed, respectively. The initial profile of the liquid film due to perturbation is assumed to be hinitial ¼ h0[1e0.1 cos(2px/l)], and the heat flux profile and temperature profile are assumed to be uniform over the heater surface. Thus the value of DT (temperature difference between the heater surface temperature and the saturation temperature) is the same everywhere, and equal to 5 C according to the nominal heat flux of 36.8 kW/m2 and experimental conditions [2]. The actual heat flux employed in the simulation is 32.4 kW/m2 by taking heat loss into account. For the liquid films of h0 ¼ 300 mm and 30 mm, the perturbation is damped quickly (s < 10) and a flat surface forms as shown in Fig. 5a and b. This is because gravity is dominant during liquid film evolution. However, with further decrease in initial liquid film thickness as shown in Fig. 5c, the stabilizing forces (especially the gravity force) are diminishing but the destabilizing forces are getting pronounced, leading to zero thickness (rupture) in the center of the film. Since most of the terms in Eq. (1) are the functions of initial liquid film thickness, it is not surprising that after the film thickness decreases to a certain value the destabilizing forces exceeds the stabilizing ones, so that the perturbation is no longer suppressed.
Fig. 3. Test section of liquid film evaporating on a locally heated solid surface.
490
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
Fig. 4. Schematic of the two-dimensional evaporating liquid film for simulation.
Fig. 6 shows the thinning processes of the liquid films along the centerline of heater surface, where the thinning curves with absence of certain forces are also plotted. From the figures one can conclude: i) the thinning curve is similar to the experimental data [2]; ii) the
Fig. 6. Absence of forces vs. the evolution of a water film at DT ¼ 5 C under uniform heating. (a): h0 ¼ 300 mm; (b): h0 ¼ 30 mm; (c): h0 ¼ 3 mm (G e gravity; S e surface tension; TCethermo-capillary; E evaporation number).
mass loss and vapor recoil (related to evaporation number E) are dominant forces responsible for the thinning of liquid films, i.e., the main destabilizing factor comes from evaporation; iii) for an initially thick liquid film, the gravity force suppresses the perturbation rapidly (cf. Fig. 6a), rendering a stable liquid film; iv) for an initially thinner liquid film (h0 ¼ 30 mm as shown in Fig. 6b), the surface tension (capillary force) overtakes the gravity as the dominant force for film stabilization; v) the thermo-capillary acting as destabilizing force plays a negligible role in all cases. 4.2. Liquid films on heater surfaces with non-uniform heat flux Fig. 5. Liquid film evolution on a uniform-temperature heater surface, DT ¼ 5 C. (a) h0 ¼ 300 mm; (b) h0 ¼ 30 mm; (c) h0 ¼ 3 mm.
In this case a realistic configuration as shown in Fig. 4 is considered, where uniform heat flux is applied to the region of
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
491
Fig. 7. Temperature profile along the different heater surfaces (q ¼ 32.4 kW/m2).
10 mm 15 mm in direct contact with the heated copper block, instead of the entire heat surface (30 mm 20 mm). The temperature profile of the heater interface is calculated by solving heat conduction problem for the metal sheet which is 30 mm long, 20 mm wide and 0.42 mm thick (cf. Figs. 3 and 4). The temperature in the wide direction is assumed to be uniform, which results in a two-dimensional heat conduction problem for the heater surface. The downward boundary at the heating region (10 mm 15 mm) is set at a constant heat flux, while the remaining downward surface is assumed to have zero heat flux (insulated). The upward surface and side surface in contact with water are assigned a convective boundary condition with a constant heat transfer coefficient calculated from the heat flux and measured heater temperature. The FLUENT code is employed as the computational vehicle. Fig. 7 shows the calculated temperature profiles of different heater surfaces which have different thermal conductivities as listed in Table 1 at the heat flux of q ¼ 32.4 kW/m2 over the heating zone. It is clear that there is a steep temperature gradient near the edge of the heating zone. In other words, the heater surface temperature is a function of the x coordinate, i.e. DT(x) ¼ Ti(x)Ts. Fig. 7 also tells that under the same heat flux, the heater surface with a higher thermal conductivity will have a smaller DT and a flatter temperature profile along the surface. Given the temperature profile on titanium heater surface as shown in Fig. 7, the film thickness profile calculated by Eq. (1) is as shown in Fig. 8 at the moment when the film thickness at the centerline reaches zero (i.e., film rupture begins). The film thickness at the periphery when the rupture occurs is called the critical thickness, to consistent with the definition in the experimental work [2]. The calculated critical thickness appears relevant to its initial thickness (see the cases with h0 ¼ 50, 75 and 100 mm), but remains unchanged after the initial thickness of liquid film is much greater than 95.9 mm as the cases with h0 ¼ 300 mm and h0 ¼ 150 mm in Fig. 8a. For the case with h0 ¼ 100 mm which is approaching 95.9 mm, the calculated critical thickness is 92.4 mm (the difference is less than 4% compared with the value 95.9 mm). Fig. 8b shows the time period till the film rupture for a liquid film
Fig. 8. Liquid film evaporating on a titanium heater surface with non-uniform temperature profile (q ¼ 32.4 kW/m2); h0 ¼ 300 mm. (a) Liquid film thickness profile at rupture for cases with varied initial thicknesses; (b) Rupture time vs. initial film thickness.
Table 1 Thermal conductivity of heater surface materials. Material
Thermal conductivity (100 C)
Copper Aluminum Silicon Titanium SS304
394.5 W/m-K 238.5 W/m-K 104 W/m-K 20.7 W/m-K 16.2 W/m-K
Fig. 9. Liquid film profiles on different heater surfaces with non-uniform temperature profiles (q ¼ 32.4 kW/m2); h0 ¼ 300 mm. (a) Liquid thickness profiles at rupture; (b) Calculated critical thickness at liquid film rupture.
492
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
Fig. 11. The measured critical thickness vs. the predicted by Eq. (11) for different Ti surfaces.
Fig. 10. Impact of heat flux on the critical thickness of liquid film on non-uniformtemperature heater surfaces (Ti_400 refers to the Ti surface aged at 400 C).
with different initial thickness. To avoid the influence of the initial condition, all the following calculations use the same initial thickness of h0 ¼ 300 mm, identical to the initial value in the experiment [2] obtaining the critical thickness at film rupture. Fig. 9a presents the liquid film profiles on different heater surfaces (with the temperature profiles as calculated in Fig. 7) when the film rupture occurs. Accordingly, the calculated rupture thicknesses hc,cal are shown in Fig. 9b. Apparently the critical thickness decreases with increasing thermal conductivity. This is consistent with the experimental finding, and also explains why the film stability benefits from a good thermal conduction of the heater surface [2]. It should be noted it is a challenge to obtain the precise temperature profile under the evaporating film. The numerical results of the conduction problem presented here is just a first-cut estimation of temperature profiles affected the substrate’s thermal conductivity. However, the thermal resistance of the heater substrate does affect the evaporation, as demonstrated by Sefiane & Bennacer [19] for sessile droplets’ evaporation under thermal influence.
Based on more calculations by varying heat flux, the calculated critical thickness of the liquid film rupture on either titanium or copper surface is found to increase with increasing heat flux, as illustrated in Fig. 10. This tendency is consistent with the experimental data (the solid triangles and squares), although the measured critical thicknesses are underestimated by the calculations. The calculations also confirms the finding of [2] that the critical thickness of liquid film rupture on the copper surface is much smaller than that on the titanium surface, and the main reason for this is that the higher thermal conductivity helps the better uniformity in temperature of the heater surface. If we assume that the discrepancy between the calculated critical thickness and the measured critical thickness for film rupture is due to wettability which is not included in Eq. (1), one can use Eq. (10) to predict the critical thickness of liquid film rupture. Based on the experimental data [2], the constant C is correlated as in Table 2, taking the values of 2.4, 2.24, 2.33 and 2.46 for different titanium surfaces (with different contact angles). The averaged value of the C constant is 2.34, with which Eq. (10) can be written as follow for the titanium surfaces.
Table 2 Contact angle of water on titanium surfaces. Material
Aging temp. ( C)
Contact angle ( )
C in Eq. (10)
Average C
Ti_20 Ti_200 Ti_400 Ti_800
e 200 400 800
51.2 55.4 43.5 39.2
2.40 2.24 2.33 2.46
2.34
Fig. 12. Critical thickness vs. hear flux for copper surface.
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
sffiffiffiffiffiffiffiffiffi f ðqÞ $h hc ¼ 2:34 4ðqÞ c;cal
(11)
Fig. 11 is the comparison of the prediction of Eq. (11) with the critical thickness of liquid film rupture obtained in the experiment [2]. One can see the deviation is within 10% for all data points. Using Eq. (11) the critical thickness of water film on the copper surface (with average contact angle of 45.4 ) is predicted as in Fig. 12, where the experimental values and the calculated one by Eq. (1) are also plotted. The agreement implies the applicability of Eq. (11) to the surfaces made of other materials, particularly at a low heat flux. 5. Conclusions In this work the lubrication theory is employed to analyze liquid film dynamics/evolution evaporating on heater surfaces, with the objectives to interpret our experimental findings on liquid film evolution and its critical thickness at rupture by using a confocal optical measurement system [2], and to examine the capability of the theory and the influence of acting forces in liquid film dynamics. While the influences of liquid properties, heat flux, and thermal conductivity of heater surface are captured by the simulation, influence of the wettability is considered via a minimum free energy criterion. Generally speaking, the thinning processes of the liquid films are well captured by the simulation of the lubrication theory. For the case with ideally uniform heat flux over the heater surface, at the thickness level of tens micrometers for the liquid film, the forces due to gravity and surface tension (working as stabilizing) are comparable with the forces due to evaporative mass loss and recoil force (working as destabilizing), and thus instability of the liquid film results. For the case of non-uniform heating, a temperature profile of the heater surface is obtained by solving a simplified two-dimensional model. Given the temperature profile, the critical thickness for the film rupture can be predicted by the lubrication theory simulation, although its value is lower than the experimental data. By introducing a coefficient (which is a function of contact angle) based on the minimum free energy criterion, the predicted critical thickness can be modified to reflect the influence of surface wettability. The so-obtained critical thicknesses have a good agreement with the experimental ones for both titanium and copper surfaces, with a maximum deviation less than 10%. The results also explain why the critical thickness on a copper surface is thinner than that on a titanium surface, given the conditions with similar other parameters. It is because the good thermal conductivity of copper surface helps to decrease the nonuniformity of temperature distribution on the heat surface, and therefore increases the resilience of the liquid film to rupture. Acknowledgements This study is made possible by research grant VR-2005-5729 from Vetenskapsrådets (Swedish Research Council) and supported by the NORTHNET RM1 and APRI7 research programs. The authors thank the staff at the Nuclear Power Safety Laboratory for their technical support in experimental setup. Nomenclature A A0 Al Ah C0
Dimensionless Hamaker constant Hamaker constant Area of liquid film Surface area of the vapor hole Constant 75.83$103 N/m
C1 C D E DE g G h h0 hc hc,cal hc,exp H k K L M Mw Pr r r1 r2 Rg S t tn tE T Ti Ts DT x X
493
Constant 0.1477$103 N/m C Constant, Eq. (10) Ratio of vapor density and liquid density Evaporative number External energy for creating a hole Gravitational acceleration 9.8 m/s2 Non-dimensional gravity Liquid film thickness m Initial liquid film thickness m Critical thickness m Calculated critical thickness by Eq. (1) m Critical thickness of film rupture in experiment m Non-dimensional thickness Thermal conductivity W/(m$K) Degree of non-equilibrium at the interface Latent heat 2257 kJ/kg Marangoni number Molecular weight of water 0.018 kg/mol Prandtl number Radius of vapor hole m Radius of dry hole on the solid surface m Maximum radius of vapor hole at the liquidevapor interface m Universal gas constant 8.314 J/mol K Non-dimensional surface tension Time s Viscous time s Evaporative time s Temperature C Heater surface temperature C Saturated temperature C Temperature difference, TiTs C Horizontal coordinate Non-dimensional horizontal coordinate
Greek letters s Non dimensional time l Wave length m r Liquid density 958 kg/m3 rv Vapor density 0.6 kg/m3 n Kinematic viscosity 0.294 106 m2/s s Liquid-air surface tension N/m ssl Solid-liquid surface tension N/m ss Solid-air surface tension N/m s0 Surface tension at saturated temperature 0.0588 N/m g Surface tension gradient k Thermal diffusivity m2/s b Accommodation coefficient for evaporation 1 q Contact angle a Convective heat transfer coefficient W/(m2 K) References [1] S.J. Gong, W.M. Ma, T.N. Dinh, Diagnostic techniques for the dynamics of a thin liquid film under forced flow and evaporating conditions, Microfluidics and Nanofluidics 9 (2010) 1077e1089. [2] S.J. Gong, W.M. Ma, T.N. Dinh, An experimental study of rupture dynamics of evaporating liquid films on different heater surfaces, International Journal of Heat and Mass Transfer 54 (2011) 1538e1547. [3] S.G. Bankoff, Minimum thickness of a draining liquid film, International Journal of Heat Mass Transfer 14 (1970) 2143e2146. [4] J. Mikielewicz, J.R. Moszynski, Minimum thickness of a liquid film flowing vertically down a solid surface, International Journal of Heat Mass Transfer 19 (1970) 771e776. [5] A.B. Ponter, K.M. Åswald, Minimum thickness of a liquid film flowing down a vertical surface-validity of Mikielewicz and Moszynski’s equation, International Journal of Heat Mass Transfer 20 (1977) 575e576.
494
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
[6] A. Doniec, Flow of a laminar liquid film down a vertical surface, Chemical Engineering Science 43 (1988) 847e854. [7] A. Sharma, E. Ruckenstein, Dewetting of solids by the formation of holes in macroscopic liquid films, Journal of Colloid and Interface Science 133 (1989) 358e368. [8] A. Sharma, E. Ruckenstein, Energetic criteria for the breakup of liquid films on non-wetting solid surface, Journal of Colloid and Interface Science 137 (1990) 433e445. [9] M.B. Williams, S.H. Davis, Nonlinear theory of film rupture, Journal of Colloid Interface Science 90 (1982) 220e228. [10] J.P. Burelbach, S.G. Bankoff, S.H. Davis, Nonlinear stability of evaporating/ condensing liquid films, Journal of Fluid Mechanics 195 (1988) 463e494. [11] M.J. Tan, S.G. Bankoff, S.H. Davis, Steady thermocapillary flows of thin liquid layers, I. Theory, Physics of Fluids A 2 (1990) 313e321. [12] J.P. Burelbach, S.G. Bankoff, S.H. Davis, Steady thermocapillary flows of thin liquid layers, II. Experiment, Physics of Fluids A 2 (1990) 322e333.
[13] A. Oron, S.H. Davis, S.G. Bankoff, Long-scale evolution of thin liquid films, Reviews of Modern Physics 69 (1997) 931e980. [14] R.V. Craster, O.K. Matar, Dynamics and stability of thin liquid films, Reviews of Modern Physics 81 (2009) 1131e1198. [15] C.Y. Ji, Y.Y. Yan, A molecular dynamics simulation of liquidevapouresolid system near triple-phase contact line of flow boiling in a microchannel, Applied Thermal Engineering 28 (2008) 195e202. [16] Y.Y. Yan, Y.Q. Zu, B. Dong, LBM: a useful tool for mesoscale modeling of single-phase and multiphase flow, Applied Thermal Engineering 31 (2011) 649e655. [17] M.S. Plesset, A. Prosperetti, Flow of vapor in a liquid enclosure, Journal of Fluid Mechanics 78 (1976) 433e444. [18] V.P. Carey, Liquid-vapor Phase-change Phenomena, Tayor and Francis Group, LLC, 2008, p. 55. [19] K. Sefiane, R. Bennacer, An expression for droplet evaporation incorporating thermal effects, Journal of Fluid Mechanics 667 (2011) 260e271.
Contents lists available at SciVerse ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Simulation and validation of the dynamics of liquid films evaporating on horizontal heater surfaces Shengjie Gong, Weimin Ma*, Truc-Nam Dinh 1 Department of Physics, Royal Institute of Technology (KTH), Roslagstullsbacken 21, 106 91 Stockholm, Sweden
h i g h l i g h t s < Liquid film dynamics under thermal influence is simulated. < Long wave theory and minimum free energy theory are applied. < Impacts of heat flux and surface wettability on film dynamics are investigated. < A correlation is developed to predict the critical thickness of film rupture.
a r t i c l e i n f o
a b s t r a c t
Article history: Received 11 August 2011 Accepted 13 May 2012 Available online 23 May 2012
In this study a non-linear governing equation based on lubrication theory is employed to model the thinning process of an evaporating liquid film and ultimately predict the critical thickness of the film rupture under impacts of various forces resulting from mass loss, surface tension, gravity, vapor recoil and thermo-capillary. It is found that the thinning process in the experiment is well reproduced by the simulation. The film rupture is caught by the simulation as well, but it underestimates the measured critical thickness at the film rupture. The reason may be that the water wettability of the heater surfaces is not taken into account in the model. Thus, the minimum free energy criterion is used to obtain a correlation which combines the contact angle (reflection of wettability) with the critical thickness from the simulation. The critical thicknesses predicted by the correlation have a good agreement with the experimental data (the maximum deviation is less than 10%). Ó 2012 Elsevier Ltd. All rights reserved.
Keywords: Liquid film Evaporation Rupture Critical thickness
1. Introduction Liquid films spreading over and evaporating on solid surfaces are commonly encountered in daily life and industrial applications that involve processes such as cooling (e.g., condensation), heating (e.g., boiling), coating, cleaning and lubrication. Studies on liquid film dynamics are of paramount importance to achieve an efficient and safe operation of the related equipments. For instance, the instability of a liquid film will affect the coating quality, while the rupture of a near-wall liquid film (e.g., macrolayer in boiling) will deteriorate heat transfer and lead to safety issues (destruction of heater surface due to burnout) in power plants. In our previous work, a novel experimental method using a confocal optical sensor was developed to measure the thickness evolution of an evaporating liquid film, and first-of-a-kind data about the thinning and rupture process of a thin liquid film on * Corresponding author. Tel.: þ46 8 5537 8821; fax: þ46 8 5537 8830. E-mail address: [email protected] (W. Ma). 1 Present address: Idaho National Laboratory, Idaho Falls, ID 83415, USA. 1359-4311/$ e see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.applthermaleng.2012.05.021
a horizontal surface heated from below was obtained [1,2]. It was observed that a liquid film ruptures when it was thinning to a threshold value (called critical thickness). The critical thickness of the liquid films performed in the tests was found to increase with increasing heat flux, and within the range of 60e150 mm for water films and 150e250 mm for ethanol films. In addition, the film instability was greatly affected by the properties of heater surfaces. Wettability of the surfaces proves to have a profound influence on film stability, and the critical thickness decreases with the reduction in contact angle. The wettability (contact angle) alone, however, could not explain all the effects of material properties and surface conditions on film instability. Thermalephysical properties of the heater surfaces also manifest themselves through the onset of liquid film instability. For instance, given two surfaces with a comparable contact angle, the one (e.g., copper or aluminum) with a high thermal conductivity exhibited a thinner critical film than the one with a low thermal conductivity (e.g., stainless steel or titanium). For the essentially wetting case (contact angle near zero) of a hexane film formed on an aged titanium surface, the liquid film was found resilient to rupture, and instability was not detected
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
until the thickness was less than 15 mm at low heat flux. No instability was observed for high-evaporation-rate cases. This paper is regarding a companion work to the experimental study, with the focus on theoretical analysis of film dynamics. The objectives are two-fold: The high-quality data obtained in the experiment serve a rich benchmark to assess analytical models of film instability, while the interpretation of the experiment and further understanding of the phenomena require a comprehensive analysis which can discern the effects of competing forces that govern the film behavior at diminishing thickness and rupture. Over the past three decades, a large number of theoretical studies have been carried out to investigate the film instability. Notably, the minimum thickness of liquid film flowing down a vertical or inclined adiabatic solid surface has been predicted theoretically on basis of minimum total energy criteria [3e6], resulting in a reasonable agreement with the limited experimental data available under isothermal conditions. For the liquid films on horizontal solid surfaces, Sharma et al. [7,8] developed a theory to predict the critical thickness at which liquid film rupture occurs by hole formation when the free energy of the film-solid system becomes equal to the free energy of the hole-liquid-solid system. The model predicted film rupture thickness is several hundred micrometers which are in good agreement with the experimental data obtained for liquid films on selected non-wetting solid surfaces (e.g., Teflon, polyethylene, PMMA, wax). The limitations of those models are that they are neither inclusive of all key instability sources in the model nor able to capture the dynamic process of film rupture. The long wave theory [9,10] is capable to predict film dynamics and rupture process. Based on long wave theory Tan et al. [11] theoretically analyzed steady thermo-capillary flows of thin liquid layers without consideration of evaporating mass loss. The predicted results have a good agreement with the experimental data [12] obtained from thermo-capillary flows of silicon oil (nonvolatile liquid) film with thickness from 0.125 mm to 1.684 mm on a heated horizontal plate. Oron et al. [13] provided perhaps the most comprehensive review on the multifaceted subject of thin film dynamics modeling. They presented a unified mathematical system to predict the long-scale evolution of thin liquid films based on the long wave theory. The set of mathematical evolution equations has its root in the work of [10], taking into account the influential factors such as van der Waal forces, surface tension, gravity, thermo-capillary, mass loss and vapor recoil force. In their review, Oron et al. emphasized the importance of experimentally assessing and validating models derived from the long wave theory. Craster & Matar [14] presented a comprehensive review of the work carried out on thin films flows, focusing attention on the studies undertaken after the review by Oron et al. [13]. They pointed out that for the modeling of thin film dynamics the lubrication approximation was still a useful approach to elucidate a wide variety of flows in which films have small aspect ratios. Simulations of lattice Boltzmann method (LBM) and molecular dynamics (MD) also see their role in micro fluid dynamics. While the application of the MD simulation is still limited to ultra-thin liquid films [15], the LBM simulation incorporates microscopic physics into macroscopic model, and is able to simulating twophase flow near heater walls and liquid (droplet) behavior on partial wetting surfaces [16]. Yan et al. also gave an overview of the recent development of the lattice Boltzmann method for modeling engineering encountered single-phase and multiphase flow problems to demonstrate that the LBM is a useful tool for mesoscale modeling [16]. In the present study, the non-linear stability analysis approach developed by Burelbach et al. [10] based on long wave principle is employed to simulate the liquid film evolution on a heated planar
487
solid surface. Due to exclusion of contact angle in the formulation, the influence of wettability on film dynamics is not (at least not explicitly) modeled in the long wave theory. To account for this, the predicted critical thickness of the liquid is corrected by a correlation from the minimum total free energy criterion [7] which relates the critical thickness with contact angle. In other words, the present study is trying to take the advantages of both the long wave theory and the free energy theory, so that the modeling approach does not only consider all key instability sources, but also capture the dynamic process of film rupture. 2. Mathematical formulation 2.1. Governing equation of an evaporating liquid film A liquid film evaporating on the horizontal surface of a planar heater is as illustrated in Fig. 1 where the liquid film is bounded above by its vapor and below by the heater surface, but unbounded in the horizontal direction. To describe the evolution and dynamics of the liquid film, two-dimensional Cartesian coordinates (x-z) are used to describe the system. The variation of liquid film in the direction perpendicular to the x-z plane is assumed negligible, since the test section of the experiment [2] to be simulated here was designed in such a way that the film dynamics in the y-direction (perpendicular to the x-z plane) is suppressed by a short distance between two side walls, and thus the interfacial wave propagation is dominant in the x-direction. The heater surface is located at z ¼ 0 and the vaporeliquid interface at z ¼ h(x,t) where z designates vertical coordinate and h the film thickness that is the function of horizontal coordinate x and time t. The liquid is assumed to be a Newtonian fluid of constant thermo-physical properties (kinematic viscosity n, density r and thermal conductivity k). This is a free-boundary problem. In order to remove the complexity of treating the free boundaries (interfaces between liquid and gas/vapor) in the governing equations, one can resort to the long wave theory if variations along the liquid film are much more gradual than those normal to it [13]. The theory is applied in a variety of areas in classical physics, such as lubrication theory in viscous liquid films. In the lubrication theory or long-wave-theory approach an asymptotic reduction of the governing equations and boundary conditions was implemented [10], and consequently a single nonlinear partial differential equation formulated in terms of the local thickness of the liquid film can be obtained:
H E Hs þA X þS H3 HXXX G H 3 HX þ X X H þK H X 3 2 E2 H KM H þ HX þ HX ¼ 0 D H þK Pr HþK X X
(1)
where X and H are the dimensionless coordinate and the film thickness by scaling x and h to the initial averaged film thickness of
Fig. 1. The sketch of an evaporating liquid film on a heater surface.
488
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
h0; s is non-dimensional time related to the viscous time tn ¼ h20 =n; A is the scaled dimensionless Hamaker constant; S is the dimensionless surface tension constant; G is the dimensionless gravity constant; E is evaporation number which is the ratio of the viscous timescale to evaporative timescale; K is a constant measuring the degree of non-equilibrium at the evaporating interface; D is proportional to the ratio of vapor to liquid densities; M is Marangoni number; and Pr is Prandtl number. All these dimensionless parameters can be expressed as follows:
X ¼
x ; h0
s h S ¼ 0 20 ; 3rn
H ¼
h ; h0
s¼
3 rv D ¼ ; 2 r
t ; tv
A ¼
A0 ; 6ph0 rn2
gDTh0 M ¼ ; 2rnk
n Pr ¼ ; k
E ¼
configuration with a hole (Fig. 2b) is less than that of the flat film (Fig. 2a). Otherwise, the liquid film will rewet the surface and its integrity is kept. Thus, the threshold thickness for initiating liquid film instability can be obtained if the free energy of the liquid film in Fig. 2a is equal to that of the configuration in Fig. 2b:
s Al pr22 þ ssl Al pr12 þ Ah s þ pr12 ssa Al s Al ssl ¼ 0 (4)
tn kDT ¼ rnL tE 3
K ¼
kTs2 bh0 rv L2
In the above expressions, s0 is the surface tension at the saturation temperature Ts, A0 is Hamaker constant, b is the accommodation coefficient for evaporation with the value of w1 for water at 100 C [17], Rg is universal gas constant, Mw is the liquid molecular weight, and evaporative timescale tE ¼ rh20 L=kDT where L is the latent heat and DT ¼ TiTs is the temperature difference between the heater surface temperature Ti and the saturation temperature Ts. The value of M is related to
ds g ¼ dT
(2)
!
2pRg MW
(3)
where C0 ¼ 0.07583 N/m and C1 ¼ 1.477 104 N/(m C) for water [18]. The second term on the left-hand side of Eq. (1) is representing the van der Waals molecular force which is negligible if the liquid thickness is at the magnitude more than some micrometers [13]. The sixth term is concerning the non-dimensional recoil force, and the last term (with the product of K, M and Pr1) is the nondimensional thermo-capillary force. Periodic boundary conditions are assumed for Eq. (1). 2.2. Minimum free energy criterion for liquid film rupture According to minimum free energy principle, the horizontal liquid film will be unstable and rupture if the free energy of the
2
Applying the Young’s equation ssa ¼ scos q þ ssl, one can have
s Ah þ pr12 cosq pr22 ¼ 0
(5)
where Al is the liquid film area and Ah is the hole surface area. If the external energy (due to forces) used to create the hole with radius of r2 is assumed as,
DE ¼ 2pr2 hc s
(6)
the critical (threshold) thickness for liquid film rupture is resolved as [7]
The surface tension s is a linear function of the liquid temperature T (in C) as the following:
s ¼ C0 C1 $T
1
hc ¼
sffiffiffiffiffiffiffiffiffirffiffiffiffiffiffiffiffiffi f ðqÞ DE 4ðqÞ 2ps
ð0 < q < p=2Þ
(7)
where q is the contact angle, and
f ðqÞ ¼
hc r1
¼
" # 2 1 ð1 cosqÞ sinq 1 þ 2 2 esin2 q þ ð1 cosqÞ
4ðqÞ ¼ coshðf ðqÞ=sinqÞ þ cosqsinhðf ðqÞ=sinqÞ
(8)
(9)
This indicates p the critical thickness of liquid film rupture is ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi proportional to f ðqÞ=4ðqÞ, which represents the impact of wettability on liquid film instability. Since the above-mentioned lubrication theory accounts for the effects of various external forces on film dynamics and instability, Eq. (7) can be rewritten as
Fig. 2. Schematic of hole formation in a liquid film (0 < q < p/2).
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
sffiffiffiffiffiffiffiffiffi f ðqÞ $h hc ¼ C 4ðqÞ c;cal
(10)
where hc,cal is the critical liquid film thickness for rupture predicted by the lubrication theory, and C is a constant to be determined from the experimental data. 3. Experimental setup and its modeling
489
The first step is to see how the different forces affect the liquid film instability under uniform heat flux over the heater surface. The approach is then used to predict the critical thickness of the liquid film under the experimental conditions [2]. Finally, the minimum energy criterion is adopted to account for wettability effect on the critical thickness, and a correlation is developed to predict the critical thickness of an evaporating liquid film under non-uniform temperature profile of heater surface and various wettability conditions.
3.1. Test section As shown in Fig. 3, the heater surface is a rectangular metal sheet mounted on a heating system. The metal sheet is 30 mm long, 20 mm wide and 0.42 mm thick. The materials of the sheet selected in the present study are silicon wafer, copper, aluminum, stainless steel and titanium. The metallic sheets are glued to the copper block using silver epoxy to minimize thermal resistance, with the contact region of 10 mm 15 mm. Cartridge heaters are embedded in the copper block and powered by the programmable DC power supply. Copper temperature under the metallic sheet is monitored by thermocouples. A glass vessel (60 mm 20 mm) is designed to surround the test surface, to supply water and avoid edge effect due to capillary. A stereo microscope is employed to visualize film rupture. A confocal optical sensor with the measuring range of 300 mm and the nominal resolution of 0.012 mm is used to measure the dynamic thickness of the liquid film. More details of the experimental setup can be found in [2]. 3.2. Simulation model Based on the geometry as shown in Fig. 3, it is assumed that the liquid film is symmetric around the centerline perpendicular to the heater surface. Thus, the question is simplified as a twodimensional problem as shown in Fig. 4. Water film is considered in the present study. In the lubrication theory, the length of the liquid film on the metallic surface is defined to be one wave length l ¼ 30 mm, and the ratio of the initial thickness to the length of the liquid film is h0/l ¼ 0.01 which is small enough to make the long wave theory applicable. The heat flux over the heating region is obtained through evaporative mass loss rate, and further verified by balance between power input and heat loss. 4. Results and discussions This section presents the dynamics of the evaporating liquid film using the lubrication theory approach represented by Eq. (1).
4.1. Liquid films on heater surfaces with uniform heat flux In this study it is assumed that the entire heater surface (30 mm 20 mm) has a uniform (identical) heat flux, i.e., to ignore the localized actual heating region (10 mm 15 mm) in the experiment. The purpose is to have ideas on the governing forces affecting liquid film instability at given conditions. According to Eq. (1), gravity force mainly helps to stabilize the liquid film, while high order differential surface tension force contributes to smooth the liquid film surface. The other forces (evaporation mass loss item, evaporation recoil force and thermocapillary force) will destabilize the liquid film. In order to check the effect of these forces on liquid film evolution, a liquid film with initial averaged thickness of h0 ¼ 300 mm, 30 mm and 3 mm are analyzed, respectively. The initial profile of the liquid film due to perturbation is assumed to be hinitial ¼ h0[1e0.1 cos(2px/l)], and the heat flux profile and temperature profile are assumed to be uniform over the heater surface. Thus the value of DT (temperature difference between the heater surface temperature and the saturation temperature) is the same everywhere, and equal to 5 C according to the nominal heat flux of 36.8 kW/m2 and experimental conditions [2]. The actual heat flux employed in the simulation is 32.4 kW/m2 by taking heat loss into account. For the liquid films of h0 ¼ 300 mm and 30 mm, the perturbation is damped quickly (s < 10) and a flat surface forms as shown in Fig. 5a and b. This is because gravity is dominant during liquid film evolution. However, with further decrease in initial liquid film thickness as shown in Fig. 5c, the stabilizing forces (especially the gravity force) are diminishing but the destabilizing forces are getting pronounced, leading to zero thickness (rupture) in the center of the film. Since most of the terms in Eq. (1) are the functions of initial liquid film thickness, it is not surprising that after the film thickness decreases to a certain value the destabilizing forces exceeds the stabilizing ones, so that the perturbation is no longer suppressed.
Fig. 3. Test section of liquid film evaporating on a locally heated solid surface.
490
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
Fig. 4. Schematic of the two-dimensional evaporating liquid film for simulation.
Fig. 6 shows the thinning processes of the liquid films along the centerline of heater surface, where the thinning curves with absence of certain forces are also plotted. From the figures one can conclude: i) the thinning curve is similar to the experimental data [2]; ii) the
Fig. 6. Absence of forces vs. the evolution of a water film at DT ¼ 5 C under uniform heating. (a): h0 ¼ 300 mm; (b): h0 ¼ 30 mm; (c): h0 ¼ 3 mm (G e gravity; S e surface tension; TCethermo-capillary; E evaporation number).
mass loss and vapor recoil (related to evaporation number E) are dominant forces responsible for the thinning of liquid films, i.e., the main destabilizing factor comes from evaporation; iii) for an initially thick liquid film, the gravity force suppresses the perturbation rapidly (cf. Fig. 6a), rendering a stable liquid film; iv) for an initially thinner liquid film (h0 ¼ 30 mm as shown in Fig. 6b), the surface tension (capillary force) overtakes the gravity as the dominant force for film stabilization; v) the thermo-capillary acting as destabilizing force plays a negligible role in all cases. 4.2. Liquid films on heater surfaces with non-uniform heat flux Fig. 5. Liquid film evolution on a uniform-temperature heater surface, DT ¼ 5 C. (a) h0 ¼ 300 mm; (b) h0 ¼ 30 mm; (c) h0 ¼ 3 mm.
In this case a realistic configuration as shown in Fig. 4 is considered, where uniform heat flux is applied to the region of
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
491
Fig. 7. Temperature profile along the different heater surfaces (q ¼ 32.4 kW/m2).
10 mm 15 mm in direct contact with the heated copper block, instead of the entire heat surface (30 mm 20 mm). The temperature profile of the heater interface is calculated by solving heat conduction problem for the metal sheet which is 30 mm long, 20 mm wide and 0.42 mm thick (cf. Figs. 3 and 4). The temperature in the wide direction is assumed to be uniform, which results in a two-dimensional heat conduction problem for the heater surface. The downward boundary at the heating region (10 mm 15 mm) is set at a constant heat flux, while the remaining downward surface is assumed to have zero heat flux (insulated). The upward surface and side surface in contact with water are assigned a convective boundary condition with a constant heat transfer coefficient calculated from the heat flux and measured heater temperature. The FLUENT code is employed as the computational vehicle. Fig. 7 shows the calculated temperature profiles of different heater surfaces which have different thermal conductivities as listed in Table 1 at the heat flux of q ¼ 32.4 kW/m2 over the heating zone. It is clear that there is a steep temperature gradient near the edge of the heating zone. In other words, the heater surface temperature is a function of the x coordinate, i.e. DT(x) ¼ Ti(x)Ts. Fig. 7 also tells that under the same heat flux, the heater surface with a higher thermal conductivity will have a smaller DT and a flatter temperature profile along the surface. Given the temperature profile on titanium heater surface as shown in Fig. 7, the film thickness profile calculated by Eq. (1) is as shown in Fig. 8 at the moment when the film thickness at the centerline reaches zero (i.e., film rupture begins). The film thickness at the periphery when the rupture occurs is called the critical thickness, to consistent with the definition in the experimental work [2]. The calculated critical thickness appears relevant to its initial thickness (see the cases with h0 ¼ 50, 75 and 100 mm), but remains unchanged after the initial thickness of liquid film is much greater than 95.9 mm as the cases with h0 ¼ 300 mm and h0 ¼ 150 mm in Fig. 8a. For the case with h0 ¼ 100 mm which is approaching 95.9 mm, the calculated critical thickness is 92.4 mm (the difference is less than 4% compared with the value 95.9 mm). Fig. 8b shows the time period till the film rupture for a liquid film
Fig. 8. Liquid film evaporating on a titanium heater surface with non-uniform temperature profile (q ¼ 32.4 kW/m2); h0 ¼ 300 mm. (a) Liquid film thickness profile at rupture for cases with varied initial thicknesses; (b) Rupture time vs. initial film thickness.
Table 1 Thermal conductivity of heater surface materials. Material
Thermal conductivity (100 C)
Copper Aluminum Silicon Titanium SS304
394.5 W/m-K 238.5 W/m-K 104 W/m-K 20.7 W/m-K 16.2 W/m-K
Fig. 9. Liquid film profiles on different heater surfaces with non-uniform temperature profiles (q ¼ 32.4 kW/m2); h0 ¼ 300 mm. (a) Liquid thickness profiles at rupture; (b) Calculated critical thickness at liquid film rupture.
492
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
Fig. 11. The measured critical thickness vs. the predicted by Eq. (11) for different Ti surfaces.
Fig. 10. Impact of heat flux on the critical thickness of liquid film on non-uniformtemperature heater surfaces (Ti_400 refers to the Ti surface aged at 400 C).
with different initial thickness. To avoid the influence of the initial condition, all the following calculations use the same initial thickness of h0 ¼ 300 mm, identical to the initial value in the experiment [2] obtaining the critical thickness at film rupture. Fig. 9a presents the liquid film profiles on different heater surfaces (with the temperature profiles as calculated in Fig. 7) when the film rupture occurs. Accordingly, the calculated rupture thicknesses hc,cal are shown in Fig. 9b. Apparently the critical thickness decreases with increasing thermal conductivity. This is consistent with the experimental finding, and also explains why the film stability benefits from a good thermal conduction of the heater surface [2]. It should be noted it is a challenge to obtain the precise temperature profile under the evaporating film. The numerical results of the conduction problem presented here is just a first-cut estimation of temperature profiles affected the substrate’s thermal conductivity. However, the thermal resistance of the heater substrate does affect the evaporation, as demonstrated by Sefiane & Bennacer [19] for sessile droplets’ evaporation under thermal influence.
Based on more calculations by varying heat flux, the calculated critical thickness of the liquid film rupture on either titanium or copper surface is found to increase with increasing heat flux, as illustrated in Fig. 10. This tendency is consistent with the experimental data (the solid triangles and squares), although the measured critical thicknesses are underestimated by the calculations. The calculations also confirms the finding of [2] that the critical thickness of liquid film rupture on the copper surface is much smaller than that on the titanium surface, and the main reason for this is that the higher thermal conductivity helps the better uniformity in temperature of the heater surface. If we assume that the discrepancy between the calculated critical thickness and the measured critical thickness for film rupture is due to wettability which is not included in Eq. (1), one can use Eq. (10) to predict the critical thickness of liquid film rupture. Based on the experimental data [2], the constant C is correlated as in Table 2, taking the values of 2.4, 2.24, 2.33 and 2.46 for different titanium surfaces (with different contact angles). The averaged value of the C constant is 2.34, with which Eq. (10) can be written as follow for the titanium surfaces.
Table 2 Contact angle of water on titanium surfaces. Material
Aging temp. ( C)
Contact angle ( )
C in Eq. (10)
Average C
Ti_20 Ti_200 Ti_400 Ti_800
e 200 400 800
51.2 55.4 43.5 39.2
2.40 2.24 2.33 2.46
2.34
Fig. 12. Critical thickness vs. hear flux for copper surface.
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
sffiffiffiffiffiffiffiffiffi f ðqÞ $h hc ¼ 2:34 4ðqÞ c;cal
(11)
Fig. 11 is the comparison of the prediction of Eq. (11) with the critical thickness of liquid film rupture obtained in the experiment [2]. One can see the deviation is within 10% for all data points. Using Eq. (11) the critical thickness of water film on the copper surface (with average contact angle of 45.4 ) is predicted as in Fig. 12, where the experimental values and the calculated one by Eq. (1) are also plotted. The agreement implies the applicability of Eq. (11) to the surfaces made of other materials, particularly at a low heat flux. 5. Conclusions In this work the lubrication theory is employed to analyze liquid film dynamics/evolution evaporating on heater surfaces, with the objectives to interpret our experimental findings on liquid film evolution and its critical thickness at rupture by using a confocal optical measurement system [2], and to examine the capability of the theory and the influence of acting forces in liquid film dynamics. While the influences of liquid properties, heat flux, and thermal conductivity of heater surface are captured by the simulation, influence of the wettability is considered via a minimum free energy criterion. Generally speaking, the thinning processes of the liquid films are well captured by the simulation of the lubrication theory. For the case with ideally uniform heat flux over the heater surface, at the thickness level of tens micrometers for the liquid film, the forces due to gravity and surface tension (working as stabilizing) are comparable with the forces due to evaporative mass loss and recoil force (working as destabilizing), and thus instability of the liquid film results. For the case of non-uniform heating, a temperature profile of the heater surface is obtained by solving a simplified two-dimensional model. Given the temperature profile, the critical thickness for the film rupture can be predicted by the lubrication theory simulation, although its value is lower than the experimental data. By introducing a coefficient (which is a function of contact angle) based on the minimum free energy criterion, the predicted critical thickness can be modified to reflect the influence of surface wettability. The so-obtained critical thicknesses have a good agreement with the experimental ones for both titanium and copper surfaces, with a maximum deviation less than 10%. The results also explain why the critical thickness on a copper surface is thinner than that on a titanium surface, given the conditions with similar other parameters. It is because the good thermal conductivity of copper surface helps to decrease the nonuniformity of temperature distribution on the heat surface, and therefore increases the resilience of the liquid film to rupture. Acknowledgements This study is made possible by research grant VR-2005-5729 from Vetenskapsrådets (Swedish Research Council) and supported by the NORTHNET RM1 and APRI7 research programs. The authors thank the staff at the Nuclear Power Safety Laboratory for their technical support in experimental setup. Nomenclature A A0 Al Ah C0
Dimensionless Hamaker constant Hamaker constant Area of liquid film Surface area of the vapor hole Constant 75.83$103 N/m
C1 C D E DE g G h h0 hc hc,cal hc,exp H k K L M Mw Pr r r1 r2 Rg S t tn tE T Ti Ts DT x X
493
Constant 0.1477$103 N/m C Constant, Eq. (10) Ratio of vapor density and liquid density Evaporative number External energy for creating a hole Gravitational acceleration 9.8 m/s2 Non-dimensional gravity Liquid film thickness m Initial liquid film thickness m Critical thickness m Calculated critical thickness by Eq. (1) m Critical thickness of film rupture in experiment m Non-dimensional thickness Thermal conductivity W/(m$K) Degree of non-equilibrium at the interface Latent heat 2257 kJ/kg Marangoni number Molecular weight of water 0.018 kg/mol Prandtl number Radius of vapor hole m Radius of dry hole on the solid surface m Maximum radius of vapor hole at the liquidevapor interface m Universal gas constant 8.314 J/mol K Non-dimensional surface tension Time s Viscous time s Evaporative time s Temperature C Heater surface temperature C Saturated temperature C Temperature difference, TiTs C Horizontal coordinate Non-dimensional horizontal coordinate
Greek letters s Non dimensional time l Wave length m r Liquid density 958 kg/m3 rv Vapor density 0.6 kg/m3 n Kinematic viscosity 0.294 106 m2/s s Liquid-air surface tension N/m ssl Solid-liquid surface tension N/m ss Solid-air surface tension N/m s0 Surface tension at saturated temperature 0.0588 N/m g Surface tension gradient k Thermal diffusivity m2/s b Accommodation coefficient for evaporation 1 q Contact angle a Convective heat transfer coefficient W/(m2 K) References [1] S.J. Gong, W.M. Ma, T.N. Dinh, Diagnostic techniques for the dynamics of a thin liquid film under forced flow and evaporating conditions, Microfluidics and Nanofluidics 9 (2010) 1077e1089. [2] S.J. Gong, W.M. Ma, T.N. Dinh, An experimental study of rupture dynamics of evaporating liquid films on different heater surfaces, International Journal of Heat and Mass Transfer 54 (2011) 1538e1547. [3] S.G. Bankoff, Minimum thickness of a draining liquid film, International Journal of Heat Mass Transfer 14 (1970) 2143e2146. [4] J. Mikielewicz, J.R. Moszynski, Minimum thickness of a liquid film flowing vertically down a solid surface, International Journal of Heat Mass Transfer 19 (1970) 771e776. [5] A.B. Ponter, K.M. Åswald, Minimum thickness of a liquid film flowing down a vertical surface-validity of Mikielewicz and Moszynski’s equation, International Journal of Heat Mass Transfer 20 (1977) 575e576.
494
S. Gong et al. / Applied Thermal Engineering 48 (2012) 486e494
[6] A. Doniec, Flow of a laminar liquid film down a vertical surface, Chemical Engineering Science 43 (1988) 847e854. [7] A. Sharma, E. Ruckenstein, Dewetting of solids by the formation of holes in macroscopic liquid films, Journal of Colloid and Interface Science 133 (1989) 358e368. [8] A. Sharma, E. Ruckenstein, Energetic criteria for the breakup of liquid films on non-wetting solid surface, Journal of Colloid and Interface Science 137 (1990) 433e445. [9] M.B. Williams, S.H. Davis, Nonlinear theory of film rupture, Journal of Colloid Interface Science 90 (1982) 220e228. [10] J.P. Burelbach, S.G. Bankoff, S.H. Davis, Nonlinear stability of evaporating/ condensing liquid films, Journal of Fluid Mechanics 195 (1988) 463e494. [11] M.J. Tan, S.G. Bankoff, S.H. Davis, Steady thermocapillary flows of thin liquid layers, I. Theory, Physics of Fluids A 2 (1990) 313e321. [12] J.P. Burelbach, S.G. Bankoff, S.H. Davis, Steady thermocapillary flows of thin liquid layers, II. Experiment, Physics of Fluids A 2 (1990) 322e333.
[13] A. Oron, S.H. Davis, S.G. Bankoff, Long-scale evolution of thin liquid films, Reviews of Modern Physics 69 (1997) 931e980. [14] R.V. Craster, O.K. Matar, Dynamics and stability of thin liquid films, Reviews of Modern Physics 81 (2009) 1131e1198. [15] C.Y. Ji, Y.Y. Yan, A molecular dynamics simulation of liquidevapouresolid system near triple-phase contact line of flow boiling in a microchannel, Applied Thermal Engineering 28 (2008) 195e202. [16] Y.Y. Yan, Y.Q. Zu, B. Dong, LBM: a useful tool for mesoscale modeling of single-phase and multiphase flow, Applied Thermal Engineering 31 (2011) 649e655. [17] M.S. Plesset, A. Prosperetti, Flow of vapor in a liquid enclosure, Journal of Fluid Mechanics 78 (1976) 433e444. [18] V.P. Carey, Liquid-vapor Phase-change Phenomena, Tayor and Francis Group, LLC, 2008, p. 55. [19] K. Sefiane, R. Bennacer, An expression for droplet evaporation incorporating thermal effects, Journal of Fluid Mechanics 667 (2011) 260e271.