- Email: [email protected]
Characteristics of an evaporating thin film of a highly wetting liquid in a groove
International Journal of Heat and Mass Transfer 79 (2014) 829–837
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Characteristics of an evaporating thin film of a highly wetting liquid in a groove H. Honda a,⇑, O. Makishi b a b
Kyushu University, Kasuyamachi, Kasuyagun, Fukuoka 811-2307, Japan Department of Mechanical Systems Engineering, Okinawa National College of Technology, Nago, Okinawa 905-2192, Japan
a r t i c l e
i n f o
Article history: Received 11 April 2014 Received in revised form 13 June 2014 Accepted 3 August 2014
Keywords: Evaporation Thin film Highly wetting liquid Marangoni effect Groove Numerical Analysis
a b s t r a c t The behavior of an evaporating thin film of a highly wetting liquid filled in a rectangular groove was studied theoretically. Following the previous studies, the liquid film was divided into an adsorbed film region, a micro region and an intrinsic meniscus region. A theoretically based method was employed to estimate the Hamaker constant. The Marangoni effect due to the thickness change of a very thin liquid film was included in the basic equation. The fourth order differential equation for the liquid film thickness was solved by the finite difference method. Three sets of the combination of boundary conditions at both ends of the micro region were examined and the set of boundary conditions which gave the most realistic solution was decided. The solutions were obtained for a wide range of the micro region length. The principle of maximum entropy production was introduced to estimate the most probable solution. The results were compared with previous results obtained by the shooting method and the effect of the solution method on the numerical results was discussed. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction Heat transfer during evaporation from a very thin liquid film plays a very important role in thermal equipments such as grooved heat pipes, microchannel evaporators and microfin tube evaporators, and in nucleate boiling. The behavior of thin liquid films is considerably different from that of thick films due to the effects of long-range intermolecular forces. A comprehensive review of relevant literature is given by Wayner [1]. Potash and Wayner [2] presented a theoretical model of evaporation from the extended meniscus of a highly wetting liquid on a vertical wall, which consisted of an adsorbed film region with no evaporation, a micro region with a very thin liquid film and an intrinsic meniscus region. A number of investigators extended this model to the analysis of heat transfer in the grooved heat pipes [3–7], microchannel evaporators [8] and microfin-tube evaporators [9]. This model was also used for the analysis of heat transfer during nucleate boiling [10–12] and evaporation of liquid droplets [13]. The basic equations for the evaporating thin liquid film are reduced to the fourth order ordinary differential equation for the liquid film thickness d with respect to the coordinate x measured along the heated surface. At the boundary between the adsorbed ⇑ Corresponding author. Tel./fax: +81 92 938 0570. E-mail address: [email protected] (H. Honda). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.08.044 0017-9310/Ó 2014 Elsevier Ltd. All rights reserved.
film region and the micro region, x ¼ 0, where the evaporation rate is zero, the boundary conditions are given by [7,8]
d ¼ d0
ð1Þ
dd=dx ¼ 0
ð2Þ
2
2
ð3Þ
3
3
ð4Þ
d d=dx ¼ 0 d d=dx ¼ 0
where d0 is the thickness of the adsorbed film. According to Wayner [1], d0 is given by
d0 ¼
AT v 6pql hlv DT
1=3 for d0 < 10 nm
ð5Þ
where A is the Hamaker constant of the system consisting of vapor, liquid film and solid substrate. Solution of the above mentioned fourth order differential equation subject to the boundary conditions given by Eqs. (1)–(4) results in the trivial solution of d ¼ d0 all over the surface [7,8]. In order to avoid this difficulty, Xu and Carey [3] and Ha and Peterson [6] adopted a simplified model in which the capillary effect on the liquid pressure was neglected. Khrustalev and Faghri [5] adopted a simplified model in which the radius of curvature of the evaporating liquid film r was assumed to be equal to that of the intrinsic meniscus rb all over the micro
830
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
Nomenclature A B1, B2 hlv m P Pd Pm DP Q q r T DT w x, y
Hamaker constant of vapor–liquid film-substrate system, J parameters in Eq. (29) latent heat of evaporation, J/kg mass flow rate of liquid, kg/m pressure, Pa disjoining pressure, Pa recoil pressure, Pa total pressure difference (¼ Pv Pl ), Pa heat flow rate, W/m heat flux, W/m2 radius of curvature of liquid film, m temperature, K wall superheat (¼ T w T v ), K groove width, m coordinates along and normal to the groove wall, respectively, Fig. 1, m
e
angle, Fig. 1, ° thermal conductivity, W/m K dynamic viscosity, Pa s density, kg/m3 surface free energy, J/m2
k
l q r
Subscripts b connecting point between micro region and intrinsic meniscus region c bottom of intrinsic meniscus i vapor–liquid interface l liquid max maximum value new new value old old value v vapor w wall 0 adsorbed film
Greek symbols d liquid film thickness, m
region. Stephan and Busse [4], Wang et al. [8] and Stephan and Hammer [10] solved the set of basic equations by the shooting method, in which one of the boundary values at x ¼ 0 was perturbed so that r approached rb at a large x. Do et al. [7] solved the basic equations with the boundary conditions of
r ¼ rb
ð6Þ 0
dd=dx ¼ tan e
ð7Þ
at x ¼ 0, where tan e was the slope of the liquid film at the boundary between the micro region and the intrinsic meniscus region and x0 was the coordinate measured in the x direction from the boundary. The value of r b was assumed to be a function of the liquid film thickness at x0 ¼ 0. The end of the micro region was determined by the condition that the evaporation rate there was equal to 0. Among the solution methods described above, that of Stephan and Hammmer [10] is considered to be most accurate, because it does not include additional assumptions for the boundary values except for the choice of the perturbed parameter. Honda and Wang [9] solved the fourth order differential equation for the liquid film thickness by the finite difference method. The boundary conditions at x ¼ 0 were Eq. (2) and 0
d ¼ d0r
ð8Þ
where d0r (> d0 ) was the thickness of the liquid film considering the curvature effect. The geometrical condition that the liquid film in the micro region was smoothly connected with the intrinsic meniscus region at x ¼ xb was given by Eq. (6) and
dd=dx ¼ tan e
ð9Þ
This solution method is attractive because the boundary conditions at both ends of the micro region are satisfied correctly. However, the obtained solution was not realistic on the point that, at x ¼ 0, the liquid film thickness was not continuous and the curvature of the liquid film took the maximum value. Thus more work is required to decide the combination of boundary conditions at x ¼ 0 and x ¼ xb which gives the most realistic solution. As described in [1], the value of the Hamaker constant depends on the combination of vapor, liquid and solid substrate which constitutes the evaporation system. In the previous studies, however,
sufficient information on the physical properties required to calculate A was not available. In the present study, the value of A is estimated using a theoretically based method. Another problem that has not been addressed in the previous studies is the Marangoni effect which was caused by the surface energy change of a very thin liquid film along the surface as a result of the change in the liquid film thickness. In the present study, the interfacial shear stress due to the Marangoni effect is taken into account in the basic equation. The fourth order differential equation for the liquid film thickness is solved by the finite difference method developed by Honda and Wang [9]. Numerical calculations are conducted for three sets of the combination of boundary conditions at x ¼ 0 and x ¼ xb and the numerical results are compared with each other. Numerical solutions are obtained for the evaporation of FC-72 filled in a small rectangular groove made of copper. Comparison is made between the cases in which the Marangoni effect is taken into account and neglected. Comparison is also made with the previous results of Raj et al. [13] for the evaporation of a FC-72 droplet on a heated surface in which the shooting method of Stephan and Hammer [10] was adopted. 2. Theoretical model Consider evaporation of a highly wetting apolar liquid filled in a rectangular groove. Fig. 1 shows the physical model and the coordinates. The coordinate x is measured downward along the groove from the boundary between the adsorbed film region and the micro region, and the coordinate y is measured vertically outward from the groove surface. Depending on the amount of liquid filled in the groove, the coordinate x ¼ 0 is located at the edge of the groove or on the groove wall. The coordinate xc denotes the position of liquid meniscus at the center of the groove. For the micro region, the momentum equation in the x-direction is written as 2
ll
d u dy
2
¼
dP l dx
ð10Þ
The boundary conditions are
u ¼ 0 at y ¼ 0
ð11Þ
831
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
where r1 denotes the value of r at d ¼ 1. The value of r decreases monotonically, from that of the substrate to that of the bulk liquid, with increasing liquid film thickness [17]. Thus Eq. (19) can be written as
r ¼ rl
A
ð20Þ
12pd2
where rl is the surface free energy of the bulk liquid (i.e., surface tension). In the previous studies, the dependence of r on d was overlooked and r was assumed to be equal to rl . It is relevant to estimate here the magnitude of the second term on the right hand side of Eq. (20) relative to the first term. The second term takes the largest value at the smallest value of the liquid film thickness; i.e., d ¼ d0 . For one of the conditions of numerical calculation discussed in Section 3, the values of rl and A=12pd20 are calculated to be 9.83 mN/m and 1.07 mN/m, respectively. Thus the second term on the right hand side of Eq. (20) is not negligible at a small value of d. The interfacial shear stress due to the Marangoni effect is given by
dr A dd ¼ dx 6pd3 dx
si ¼
Substitution of Pl and into Eq. (13) yields
Fig. 1. Physical model and coordinates.
ll
at y ¼ d
ð12Þ
where si denotes the interfacial shear stress due to the Marangoni effect. Integration of Eq. (10) subject to the boundary conditions given by Eqs. (11) and (12) yields the solution for u. Then integration of u over the liquid film thickness yields the expression for the liquid mass flow rate m as follows:
m ¼ ql
q d3 dPl d2 udy ¼ l þ si ll 3 dx 2
d
0
! ð13Þ
The P l is expressed in terms of P v , the disjoining pressure Pd , the 2 capillary pressure r=r and the recoil pressure P m ¼ q2 =qv hlv as follows:
Pl ¼ Pv Pd
r r
þ Pm
ð14Þ
where r is the surface free energy of the system consisted of vapor, liquid film and heated wall. Following Wayner [1], P d is approximated by
Pd ¼
A
pffiffiffiffiffiffiffiffiffi Al As
ð16Þ
where Al and As are the Hamaker constants of liquid and solid, respectively. The liquid curvature 1=r is related to d by 2
ql d3 d A r þ Pm 3ll dx 4pd3 r
m¼
The P d is related to
dr Pd ¼ dd
ð17Þ
r by the following equation [16]:
q ¼ hlv
12pd2
ð19Þ
ð23Þ
dm dx
ð24Þ
Substitution of m given by Eq. (22) into Eq. (24) yields
q¼
q l h lv d 3 d A r d þ Pm dx 4pd3 r 3ll dx
ð25Þ
Assuming a linear temperature drop across the liquid film, q is given by
kl ðT w T i Þ d
ð26Þ
where Ti is the surface temperature of the liquid film. Following Wayner et al. [18] and assuming that the accommodation coefficient is equal to unity, q is also written as
q¼
2 pRv T v
1=2
qv h2lv
Tv
Ti Tv þ
Tv A r þ P m ql hlv 6pd3 r
ð27Þ
Eliminating T i from Eq. (26) and Eq. (27) yields
ð18Þ
A
Thus, for the present case in which the Marangoni effect is included, the disjoining pressure term is virtually larger than that for the case without the Marangoni effect by 50%. The heat flux q is related to m by
q¼
Substitution of P d given by Eq. (15) into Eq. (18) and then integration the resulting equation from d ¼ d to d ¼ 1 yield
ð22Þ
ql d3 dDP ql d3 d A rl ¼ þ Pm 3ll dx 3ll dx r 6pd3
2
1 d d=dx ¼ r f1 þ ðdd=dxÞ2 g3=2
If the Marangoni effect is neglected, m is given by
q¼
Neglecting cohesion in the vapor, A is given by the following equation [1]:
r ¼ r1
m¼
ð15Þ
6pd3
A ¼ Al
si given by Eqs. (14) and (21) respectively,
du ¼ si dy
Z
ð21Þ
d þ qkl Thv2
v lv
kl T A r DT þ v þ Pm 3 pRv T v 1=2 ql hlv 6pd r
ð28Þ
2
2
Substitution of Pm ¼ q2 =qv hlv into Eq. (28) and then solution of the resulting equation for q yields
q¼
1 ðB1 2
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi B21 4B2 Þ
ð29Þ
832
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
1=2 q q h3 q q h3 where B1 ¼ vk Tl v lv d þ qkl Thv2 pR2v T v and B2 ¼ v T vl lv DTþ l v lv qv h2lv 6pAd3 rr . Equating q given by Eqs. (25) and (28), the fourth order differential equation for d is derived as follows:
ql hl d 3 d A r d þ Pm dx 4pd3 r 3ll dx ¼
d þ qkl Thv2
kl T A r DT þ v þ Pm 3 pRv T v 1=2 ql hlv 6pd r
ð30Þ
2
v lv
It is relevant to note here that if the Marangoni effect is neglected, the terms A=4pd3 and r=r in Eq. (30) are replaced by A=6pd3 and rl =r, respectively. Thus, for the case without the Marangoni effect, the differential equation for d is given by
ql hl d 3 d A rl d þ Pm dx 6pd3 3ll dx r ¼
d þ qkl Thv2
v lv
kl T A rl DT þ v þ P m pRv T v 1=2 ql hlv 6pd3 r
d þ qkl Thv2
kl T A r DT þ v pRv T v 1=2 ql hlv 6pd3 r
ð34Þ
2
Substitution of q ¼ 0, d ¼ d0r and r ¼ r 0 at x ¼ 0 yields
ð31Þ
3. Numerical calculation Since one of the objectives of the present study is to compare the numerical results obtained by the finite difference method with the previous results obtained by the shooting method, numerical calculations are conducted for the conditions adopted by Raj et al. [13]; i.e., evaporation of FC-72 at 32 °C on the copper surface with the wall superheats of 1 K, 2 K, 4 K and 8 K. For apolar liquids, Al can be estimated by the following approximation [14]:
ð32Þ
Thus the value of Al for FC-72 at 32 °C is calculated to be 2:06 1020 J. The value of As for copper is obtained from the table of measured data given by Eichenlaub et al. [15] as 2:73 1019 J. Substitution of these values in Eq. (16) yields A ¼ 5:40 1020 J. This value is somewhat smaller than A ¼ 3:96 1020 J assumed by Raj et al. [13]. It should be mentioned that Raj et al. [13] defined the Hamaker constant by A=6p instead of A. Eq. (30) is a highly non-linear differential equation. This equation was solved by the implicit finite difference method developed by Honda and Wang [9]. The unsteady term ql @d=@t was added to the left hand side of Eq. (30) and the solution was obtained as the steady state of an unsteady liquid flow subject to the assumed initial distribution of d and the boundary conditions at both ends of the micro region. The liquid in the thin film is driven by the combined pressure gradient and interfacial shear stress and the effect of the assumed initial distribution disappears as the calculation proceeds. By choosing an appropriate value for the time step, the steady state solution is obtained relatively easily. As described in Section 2, the boundary conditions at x ¼ xb are given by Eqs. (6) and (9). For the rectangular groove, the geometrical condition that the liquid film in the micro region is smoothly connected with the intrinsic meniscus region at x ¼ xb is given by the following equation:
r b ¼ ðw=2 db Þ= cos e
q¼
v lv
2
Al ¼ 2:1 1018 rl J
Case B: d ¼ d0r and dd=dx ¼ 0 at x ¼ 0; r ¼ rb and dd=dx ¼ e at x ¼ xb . 2 2 Case C: d ¼ d0 , dd=dx ¼ 0 and d d=dx ¼ 0 at x ¼ 0; r ¼ rb and dd=dx ¼ e at x ¼ xb . Among the three sets of the boundary conditions, Case B was assumed by Honda and Wang [9]. For Case A and Case C, it follows from Eq. (17) that the capillary pressure at x ¼ 0 is equal to 0. For Case B, the capillary pressure at x ¼ 0 is not zero [9]. In order to obtain the capillary pressure and the liquid film thickness at x ¼ 0, an iterative procedure was adopted as follows: For a low heat flux where the recoil pressure is negligible as compared with the other pressure components, Eq. (28) can be approximated by
ð33Þ
Thus the boundary values db and rb must satisfy Eq. (33). In order to obtain the solution, two more boundary conditions are required at x ¼ 0. This indicates that only two of the four physically reasonable boundary conditions at x ¼ 0, given by Eqs. (1)–(4), are satisfied by the solution. In the present study, three sets of the boundary conditions were assumed. These were 2 2 Case A: d ¼ d0 and d d=dx ¼ 0 at x ¼ 0; r ¼ rb and dd=dx ¼ e at x ¼ xb .
1=3 6p ql hlv DT r d0r ¼ Tv A r0
ð35Þ
Since the value of r 0 was not known a priori, an iterative calculation was required to obtain the solution. At the beginning of the numerical calculation, the value of 1=r0 was assumed to be 0. Then, after the solution was obtained, the calculated value of 1=r0 was substituted into Eq. (35) to obtain the new assumption of d0r . The numerical calculation was repeated until a converged solution was obtained. For Case C, it was estimated that a converged solution could not be obtained, because the number of the boundary conditions was larger than the required number. In order to prove this estimate, the numerical calculation was also conducted for Case C. Since the values of the geometrical parameters xb and e were not given a priori, numerical calculations were conducted assuming a range of values for these parameters. Then the numerical results were compared with each other to determine the most appropriate values of these parameters. The groove width w was changed in the range of 0.1–1.0 mm. The grid size was assumed to be 5 1010 m or less for 0 6 x 6 5 108 m depending on the condition and it was increased gradually with increasing x. The starting value of xb was assumed to be 109 —108 m depending on the condition. At the beginning of the numerical calculation, a very small value was assumed for the time step Ds and it was increased with time. The local film thickness changed with time and approached the steady state value. The numerical calculation was continued until the convergence criteria j1 dold =dnew j < 104 and jql hlv ðdnew dold Þ=Dsqj < 104 were satisfied at all grid points. Since the finite difference solution obtained by the present method satisfies the boundary conditions at both ends of the calculation domain correctly, the errors in the calculated values of d, q, r etc. are estimated to be less than 0.01%. After a converged solution was obtained, the calculation was repeated increasing the value of xb . The numerical results revealed that the continuous distribution of d was not obtained when xb was increased beyond some value. The maximum value of xb for which the converged solution was obtained, xb;max , was determined by using the bisection method as follows: (1) Denote the largest value of xb for which the converged solution was obtained as xb;old , and the next assumption of xb for which the converged solution was not obtained as xb;old . (2) Conduct the numerical calculation using xb ¼ ðxb;old þ xb;old Þ=2 as the new assumption of the micro region length. (3) Depending on the condition whether the converged solution was obtained or not obtained, denote xb as either xb;old or xb;old . (4) Return to step (2) and repeat the numerical calculation until the convergence criterion j1 xb;old =xb;old j < 104 was satisfied. (5) Denote the last value of xb;old as xb;max .
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
833
4. Numerical results and discussion Converged solutions were obtained for all sets of the boundary conditions; i.e., Case A, Case B and Case C. For Case C, where the assumed value of the capillary pressure at x ¼ 0 was equal to 0, the calculated value of the capillary pressure took the maximum value at the first nodal point x ¼ x1 next to x ¼ 0 and the distribution of the capillary pressure in the region x P x1 was similar to that for Case B shown in Fig. 7. Thus the distribution of the capil2 2 lary pressure (and so the distribution of d d=dx ) was discontinuous between x ¼ 0 and x ¼ x1 . This indicated that the solution satisfying the set of boundary conditions of Case C did not exist. As described in Section 3, the numerical solutions were obtained for a range of the micro region length. Figs. 2–4, respectively, show the distributions of the liquid film thickness, the heat flux and the meniscus profiles for different values of the micro region lengths obtained for the set of boundary conditions of Case A. The conditions of numerical calculation are DT ¼ 1 K, w ¼ 0:2 mm and e ¼ 30 . The largest value of the micro region length in these figures (i.e., xb ¼ 17:7 lm) corresponds to xb;max . The liquid film thickness is smaller for larger xb . The slope of the liquid film thickness curve is close to 0 near x ¼ 0 but it increases rapidly with increasing x. It is seen from Fig. 3 that the heat flux takes the maximum value in the region x 6 0:065 lm. The maximum heat flux is almost unchanged for different xb , whereas the coordinate at the maximum heat flux point increases as xb
Fig. 4. Comparison of meniscus profiles for different micro region lengths.
increases. Comparison of the heat flux distributions for different xb reveals that the highest total heat transfer rate is obtained by the case of xb ¼ xb;max . The principle of maximum entropy production has been found to give successful predictions for the steadystate properties of a variety of dynamic systems [19]. Thus, if no geometrical constraint is imposed on the meniscus profile, the most probable solution is estimated to be that for xb ¼ xb;max . It is seen from the meniscus profiles shown in Fig. 4 that the apparent contact angle of the meniscus decreases as xb increases. As a result, the coordinate xc at the bottom of the meniscus is larger for larger xb . In Fig. 4, the smallest value and the largest value of xc are 0.058 mm and 0.071 mm, respectively. Thus, if the liquid level in the groove is smaller than 0.058 mm, the steady state solution of the evaporating thin film does not exist. For the liquid level in between 0.058 mm and 0.071 mm, the steady state solution will be that for which the value of xc is equal to the liquid level. For the liquid level larger than 0.071 mm, the steady state solution will be that for xb ¼ xb;max . In the discussion that follows, only the results for xb ¼ xb;max are presented. Figs. 5–7, respectively, compare the distributions of the heat flux, the total pressure difference and the capillary pressure between the solutions obtained for the set of boundary conditions of Case A and Case B. The conditions of the numerical calculation are DT ¼ 4 K, w ¼ 0:2 mm and e ¼ 30 . The heat flux curve has a positive slope at x ¼ 0 for both cases. The heat flux increases rapidly with increasing x, takes the maximum value at x 0:037 lm, and
Fig. 2. Comparison of liquid film thicknesses for different micro region lengths.
Fig. 3. Comparison of heat fluxes for different micro region lengths.
Fig. 5. Comparison of heat fluxes between the solutions obtained for the set of boundary conditions of Case A and Case B.
834
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
numerical results reveals that dd=dx ¼ 0:0128 at x ¼ 0. For the numerical results presented in Figs. 5–7, the ratio of d0r (given by Eq. (35)) to d0 (given by Eq. (5)) is calculated to be 1.020. Thus, for Case B, both the liquid film thickness and the liquid curvature are not continuous at x ¼ 0. From the comparison of the numerical results for Case A, Case B and Case C described above, it may be concluded that the most realistic solution is obtained by the set of boundary conditions of Case A. In Figs. 5–7, the solution of Raj et al. [13] is also shown for comparison. The solution was obtained for the micro region of an evaporating liquid droplet on a flat surface The boundary conditions at x ¼ 0 were d ¼ d0 , dd=dx ¼ 0, DP ¼ A=6pd30 and Q ¼ Q 0 , where Q was the perturbed parameter whose physical meaning was the heat flow rate through the liquid film at x ¼ 0. The first boundary condition is the same as that assumed in Case A, whereas the second boundary condition is the same as that assumed in Case B. Since DP is defined by DP ¼ 6pAd3 þ rrl P m , the third boundary Fig. 6. Comparison of total pressure differences between the solutions obtained for the set of boundary conditions of Case A and Case B.
then decreases monotonically with further increasing x. The difference in the heat flux between Case A and Case B is most significant in the region 0 < x < 0:004 lm. At x ¼ 0:004 lm, the value of q for Case B is larger than that for Case A by 14% and the difference increases as x decreases. The maximum heat flux and the coordinate at the maximum heat flux point are slightly smaller for Case A. It is seen from Fig. 6 that the slope of the total pressure difference curve is negative all over the surface for both Case A and Case B. This indicates that the mass flow rate through the liquid film, given by Eq. (22), is negative all over the surface. For Case A, the amount of liquid outflow at x ¼ 0 is calculated to be 0.67% of the total evaporation rate. However, it is not probable that such an outflow actually occurs. This may indicate the limitation of the macroscopic formulation in dealing with the evaporation phenomena from a very thin liquid film. Most significant difference between Case A and Case B is seen in the distributions of the capillary pressure shown in Fig. 7. For Case A, the capillary pressure is zero at x ¼ 0. It increases as x increases, then takes a maximum value, and then decreases monotonically with further increasing x. For Case B, on the other hand, the capillary pressure takes a maximum value at x ¼ 0. It decreases sharply as x increases. Then it takes a minimum value, and then a maximum value, and then decreases monotonically with further increasing x. The above results show that, for Case A, both the liquid film thickness and the liquid curvature are continuous at x ¼ 0. However, the slope of the liquid film is not continuous at x ¼ 0. In the adsorbed film region x < 0 where d ¼ d0 , it is apparent that dd=dx ¼ 0. On the other hand, examination of the
condition is correct only for the case of 6pAd3 rr0 . The value of 0
Q 0 was determined so that r ¼ r b was satisfied at x ¼ xb . According to Raj et al. [13], the exact values of r b and xb did not have a noticeable influence on the numerical solution. This is in contrast to the present results in which the solution was considerably influenced by the assumed value of xb . As seen from Figs. 5–7, a large difference is observed between the present results and the results of Raj et al. [13]. Most significant difference is seen in the distributions of calculated values near x ¼ 0. For the present results, the calculated values change sharply in the region 0 6 x 6 0:015 lm. For Raj et al. [13], on the other hand, the changes are very small in the region 0 6 x 6 0:02 lm. This is due to the difference in the solution method between the two cases. Raj et al. [13] assumed an infinitesimal perturbation Q 0 for Q . Thus the initial variations of the calculated values from the boundary values at x ¼ 0 were very small. It is relevant to note here that when Q 0 ¼ 0 is assumed, no change occurs for the calculated values all over the surface. The calculated values change first very slowly, and then more and more rapidly, with increasing x. This indicates that the effect of the assumed value of Q 0 decreases as x increases. In Fig. 7 for the capillary pressure, the calculated value by Raj et al. [13] is close to zero in the region 0 < x < 0:01 lm. Raj et al. [13] have shown a figure in which the distributions of the pressure components and the total pressure difference were plotted on the semi-logarithmic scale. According to the figure, the capillary pressure shows a trend similar to the present solution for Case B in the region close to x ¼ 0. This is probably due to the fact that the same boundary condition of dd=dx ¼ 0 at x ¼ 0 were adopted in both cases. However, the capillary pressure obtained by Raj et al. [13] is three orders of magnitude smaller than that for Case B. This was due to the fact that an infinitesimal perturbation Q 0 was assumed for Q . For this case the boundary values satisfy the inequality 6pAd3 rr0 and the boundary 0
condition DP ¼ A=6pd30 holds correctly. However, for the present numerical results for Case B, r=r0 is not negligible as compared to
Fig. 7. Comparison of capillary pressures between the solutions obtained for the set of boundary conditions of Case A and Case B.
A=6pd30r (see Figs. 6 and 7). If the shooting method is used to obtain the solution for this case, the numerical calculation must be conducted with both Q 0 and 1=r 0 as the perturbed parameters. Thus, it will not be easy to obtain the solution satisfying the boundary conditions at x ¼ xb . It may be concluded from the above discussion that the finite difference method is superior to the shooting method on the point that the boundary conditions at both ends of the micro region are satisfied correctly. In the discussion that follows, only the results for xb ¼ xb max obtained for the set of boundary conditions of Case A are presented. Figs. 8–10, respectively, compare the distributions of the liquid film thickness, the heat flux and the total pressure difference between the cases in which the Marangoni effect is included and
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
Fig. 8. Comparison of liquid film thicknesses between the solutions in which the Marangoni effect is considered and neglected.
Fig. 9. Comparison of heat fluxes between the solutions in which the Marangoni effect is considered and neglected.
835
Marangoni effect in the x direction and a smaller liquid film thickness than the case without the Marangoni effect is required to maintain the liquid flow. Comparison of the results between the two cases reveals that the largest difference in the liquid film thickness occurs at x ¼ 0:029 lm. At this point, the liquid film thickness for the case without the Marangoni effect is larger than that for the case with the Marangoni effect by 34.5%. It is seen from Fig. 9 that the largest difference in the heat flux between the two cases occurs at x ¼ 0:013 lm, where the heat flux for the case without the Marangoni effect is higher than that for the case with the Marangoni effect by 17.7%. It is interesting to note here that the film thickness at x ¼ 0:013 lm is a little bit larger for the case without the Marangoni effect (see Fig. 8). Thus, the above results indicate that the temperature drop across the liquid film, T w T i , is larger for the case without the Marangoni effect than the case with the Marangoni effect by more than 17.7%. The two heat flux curves cross at x ¼ 0:035 lm. In the region of x > 0:035 lm, the heat flux for the case with the Marangoni effect is lager than that for the case without the Marangoni effect. As a result, the difference in the total heat transfer rate between the two cases is small. For Fig. 9, the calculated values of the total heat transfer rates obtained by integrating the heat flux distributions over the x direction are 1.585 W/m and 1.564 W/m for the cases with and without the Marangoni effect. It is seen from Fig. 10 that the slope of the total pressure difference curve near x ¼ 0 is steeper for the case without the Marangoni effect. This indicates that a higher pressure drop than the case with the Marangoni effect is required to maintain the liquid flow. Figs. 11 and 12, respectively, compare the distributions of the liquid film thickness and the heat flux for different slopes of the liquid film at the boundary between the micro region and the intrinsic meniscus region. The other conditions of the numerical calculation are DT ¼ 4 K and w ¼ 0:2 mm. Due to the geometrical condition at x ¼ xb given by Eq. (33), the numerical solution was not obtained for e 6 15 . For 20 6 e, the liquid film thickness is larger for larger e. However, in the range of 30 6 e 6 55 , the liquid film thickness distribution is almost unchanged. The values of xb;max corresponding to the assumed values of e are also shown in Figs. 11 and 12. For 20 6 e 6 40 , xb;max increases rapidly as e increases. Then xb;max increases gradually with further increasing e. It is interesting to note here that the values of xb;max are considerably different for the solutions in the range of e of 30 6 e 6 55 which show almost the same liquid film thicknesses distribution. This indicates that the liquid film profiles for these solutions are close to the same circular arc in the region of x P 2:56 lm. It is
Fig. 10. Comparison of total pressure differences between the solutions in which the Marangoni effect is considered and neglected.
neglected in the basic equation; i.e., solutions of Eq. (30) and Eq. (31), respectively. The conditions of the numerical calculation are DT ¼ 4 K, w ¼ 0:2 mm and e ¼ 30 . The liquid film thickness is smaller for the case with the Marangoni effect. For this case the liquid is pulled by the interfacial shear stress due to the
Fig. 11. Comparison of liquid film thicknesses for different slopes of liquid film at the boundary between micro region and intrinsic meniscus region.
836
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
Fig. 12. Comparison of heat fluxes for different slopes of liquid film at the boundary between micro region and intrinsic meniscus region.
Fig. 14. Comparison of meniscus profiles for different groove widths.
seen from Fig. 12 that the maximum heat flux point is almost unchanged for 20 6 e 6 55 and it moves toward a smaller x with further increasing e. The total heat transfer rate is larger for smaller e. Since the solutions are obtained for a wide range of e, a physically acceptable criterion is required to determine the most appropriate value of e. In the region close to x ¼ xb , the twodimensionality of the heat conduction and the liquid flow becomes more and more significant as e increases. Thus the assumption of one-dimensional heat conduction and liquid flow adopted in the theoretical model is not appropriate for large values of e. On the other hand, for small values of e, xb;max decreases rapidly as e decreases. As a result, the liquid film thickness at x ¼ xb;max , db , also decreases rapidly and a very high evaporation heat transfer occurs near x ¼ xb;max . It is seen from Fig. 11 that for e ¼ 20 and e ¼ 25 , the liquid film thickness at x ¼ xb;max is less than 0.06 lm. Since the evaporating mass in the micro region is supplied from the intrinsic meniscus region, a high liquid flow velocity occurs in the intrinsic meniscus region near x ¼ xb;max . Thus the assumption of constant liquid pressure in the intrinsic meniscus region does not hold any more. These results indicate that too large values and too small values of e are not appropriate for the boundary condition. It may be concluded from the above discussion that the appropriate value of e to be assumed in the numerical calculation is 30°. Figs. 13 and 14, respectively, compare the distributions of the heat flux and the meniscus profile for different groove widths.
The other conditions of numerical calculation are DT ¼ 4 K and e ¼ 30 . It is seen from Fig. 13 that the heat flux distribution in the region 0 6 x 6 0:5lm is almost unchanged irrespective of the groove width. As expected from the heat flux distribution, the meniscus profiles for different w, shown in Fig. 14, are almost unchanged in the region close to x ¼ 0. However, due to the geometrical condition at x ¼ xb , the meniscus profiles for different groove widths tend to circular arcs with the radii of curvature rb given by Eq. (33) as x increases. Fig. 15 compares the heat flux distributions for different wall superheats. The other conditions of numerical calculation are w ¼ 0:2 mm and e ¼ 30 . It is seen from Fig. 15 that the maximum heat fluxes for DT ¼ 1 K, 2 K, 4 K and 8 K are 1.02 MW/m2, 2.10 MW/m2, 4.28 MW/m2 and 8.64 MW/m2, respectively. Thus the maximum heat flux is roughly proportional to the wall superheat. However, this does not mean that the liquid film thickness at the maximum heat flux point is unchanged for different DT. As seen from Eq. (29), q is a complex function of DT, d and r. The maximum heat flux point moves toward a smaller x as the wall superheat increases. In Fig. 15, the results of Raj et al. [13] are also shown for comparison. The general trend of the data are similar to the present case except that the maximum heat flux point is shifted to a larger x and the maximum heat fluxes are about 15% higher than the present values. In order to check the possibility that the difference in the maximum heat fluxes between the present results and Raj et al. [13] was due the difference in the assumed value of A, numerical calculations were also conducted using the A
Fig. 13. Comparison of heat fluxes for different groove widths.
Fig. 15. Comparison of heat fluxes for different wall superheats.
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
value assumed by Raj et al. [13]. The numerical results showed that the two solutions were quite close to each other, indicating that the solution was rather insensitive to the value of A. Thus the difference in the maximum heat flux between the present results and Raj et al. [13] was probably due to the difference in the method of numerical calculation, which resulted in different distributions of the liquid film thickness.
837
given by Eq. (33) as x increased. Comparison of the present results obtained by the finite difference method and the previous results obtained by the shooting method revealed that the finite difference method was superior to the shooting method on the point that the boundary conditions at both ends of the micro region were satisfied correctly. Conflict of interest
5. Conclusions None declared. The basic equations describing the behavior of an evaporating thin film of a highly wetting liquid in a groove is reduced to the fourth order ordinary differential equation for the liquid film thickness. At the boundary between the adsorbed film region and the micro region, there are four physically reasonable boundary conditions given by Eqs. (1)–(4). On the other hand, at the boundary between the micro region and the intrinsic meniscus region, there are two geometrical conditions given by Eqs. (6) and (9). Thus there are a number of possibilities for the selection of the combination of boundary conditions at the boundary between the adsorbed film region and the micro region. In the present study, three sets of the combination of boundary conditions, Case A, Case B and Case C, were examined. The Hamaker constant was estimated by using a theoretically based method. The numerical solutions were obtained by the finite difference method developed by Honda and Wang [9]. Comparison of the numerical results revealed that the most realistic solution was obtained by the set of boundary conditions of Case A. The numerical results showed a small amount of liquid outflow at the boundary between the adsorbed film region and the micro region. This seemed to indicate the limitation of the macroscopic formulation in describing the evaporation phenomena from a very thin liquid film. The solutions were obtained for a wide range of the micro region length. As a result, the meniscus profile was considerably different depending on the micro region length. On the basis of the principle of maximum entropy production, the most probable solution, for the case where no geometrical constraint was imposed on the meniscus profile, was estimated to be that for the largest micro region length. The numerical results showed that the Marangoni effect due to the thickness change of a very thin liquid film was not negligible. The solutions were obtained for a wide range of the assumed value of the slope of liquid film at the boundary between the micro region and the intrinsic meniscus region. On the basis of physical consideration on the numerical results, the appropriate value of the slope to be assumed in the numerical calculation was concluded to be 30°. For the groove widths of 0.1–1.0 mm, the numerical solutions in the small x region were almost unchanged. But the meniscus profiles tended to circular arcs with different radii of curvatures r b
References [1] P.C. Wayner Jr., Intermolecular forces in phase-change heat transfer: 1998 Kern award review, AIChE J. 45 (1999) (1998) 2055–2068. [2] M. Potash Jr., P.C. Wayner Jr., Evaporation from a two-dimensional extended meniscus, Int. J. Heat Mass Transfer 15 (1972) 1851–1863. [3] X. Xu, V.P. Carey, Film evaporation from a micro-grooved surface – an approximate heat transfer model and its comparison with experimental data, J. Thermophys. Heat Transfer 4 (1990) 512–520. [4] P.C. Stephan, C.A. Busse, Analysis of the heat transfer coefficient of grooved heat pipe evaporator walls, Int. J. Heat Mass Transfer 35 (1992) 383–391. [5] D. Khrustalev, A. Faghri, Heat transfer during evaporation on capillary-grooved structures of heat pipes, Trans. ASME J. Heat Transfer 117 (1995) 740–747. [6] J.M. Ha, G.P. Peterson, The interline heat transfer of evaporating thin films along a micro grooved surface, Trans. ASME J. Heat Transfer 118 (1996) 747– 755. [7] K.H. Do, S.J. Kim, S.V. Garimella, A mathematical model for analyzing the thermal characteristics of a flat micro heat pipe with a grooved wick, Int. J. Heat Mass Transfer 51 (2008) 4637–4650. [8] H. Wang, S.V. Garimella, J.Y. Murthy, Characteristics of an evaporating thin film in a microchannel, Int. J. Heat Mass Transfer 50 (2007) 3933–3942. [9] H. Honda, Y.S. Wang, Theoretical study of evaporation heat transfer in horizontal microfin tubes: stratified flow model, Int. J. Heat Mass Transfer 47 (2004) 3971–3983. [10] P. Stephan, J. Hammer, A new model for nucleate boiling heat transfer, Wärme und Stoffübertragung 30 (1994) 119–125. [11] J.H. Lay, V.K. Dhir, Shape of a vapor stem during nucleate boiling of saturated liquids, ASME J. Heat Transfer 117 (1995) 394–401. [12] C. Kunkelmann, P. Stephan, CFD simulation of boiling flows using the volumeof-fluid method within openFOAM, Numer. Heat Transfer, Part A 56 (2009) 631–646. [13] R. Raj, C. Kunkelmann, P. Stephan, J. Plawsky, J. Kim, Contact line behavior for a highly wetting fluid under superheated conditions, Int. J. Heat Mass Transfer 55 (2012) 2664–2675. [14] J.N. Israelachvili, Intermolecular and Surface Forces, second ed., Academic Press, London, 1992. [15] S. Eichenlaub, C. Chan, S.P. Beaudoin, Hamaker constants in integrated circuit metallization, J. Colloid Interface Sci. 248 (2002) 389–397. [16] B.V. Derjaguin, N.V. Churaev, Structural component of disjoining pressure, J. Colloid Interface Sci. 49 (1974) 249–255. [17] R.J. Waltman, A. Khurshudov, G.W. Tyndall, Autophobic dewetting of perfluoropolyether films on amorphous-nitrogenated carbon surfaces, Tribol. Lett. 12 (2002) 163–169. [18] P.C. Wayner Jr., Y.K. Kao, L.V. LaCroix, The interline heat transfer coefficient of an evaporating wetting film, Int. J. Heat Mass Transfer 19 (1976). 487-49. [19] R.K. Niven, Steady state of a dissipative flow-controlled system and the maximum entropy production principle, Phys. Rev. E 80 (021113) (2009) 1–15.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Characteristics of an evaporating thin film of a highly wetting liquid in a groove H. Honda a,⇑, O. Makishi b a b
Kyushu University, Kasuyamachi, Kasuyagun, Fukuoka 811-2307, Japan Department of Mechanical Systems Engineering, Okinawa National College of Technology, Nago, Okinawa 905-2192, Japan
a r t i c l e
i n f o
Article history: Received 11 April 2014 Received in revised form 13 June 2014 Accepted 3 August 2014
Keywords: Evaporation Thin film Highly wetting liquid Marangoni effect Groove Numerical Analysis
a b s t r a c t The behavior of an evaporating thin film of a highly wetting liquid filled in a rectangular groove was studied theoretically. Following the previous studies, the liquid film was divided into an adsorbed film region, a micro region and an intrinsic meniscus region. A theoretically based method was employed to estimate the Hamaker constant. The Marangoni effect due to the thickness change of a very thin liquid film was included in the basic equation. The fourth order differential equation for the liquid film thickness was solved by the finite difference method. Three sets of the combination of boundary conditions at both ends of the micro region were examined and the set of boundary conditions which gave the most realistic solution was decided. The solutions were obtained for a wide range of the micro region length. The principle of maximum entropy production was introduced to estimate the most probable solution. The results were compared with previous results obtained by the shooting method and the effect of the solution method on the numerical results was discussed. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction Heat transfer during evaporation from a very thin liquid film plays a very important role in thermal equipments such as grooved heat pipes, microchannel evaporators and microfin tube evaporators, and in nucleate boiling. The behavior of thin liquid films is considerably different from that of thick films due to the effects of long-range intermolecular forces. A comprehensive review of relevant literature is given by Wayner [1]. Potash and Wayner [2] presented a theoretical model of evaporation from the extended meniscus of a highly wetting liquid on a vertical wall, which consisted of an adsorbed film region with no evaporation, a micro region with a very thin liquid film and an intrinsic meniscus region. A number of investigators extended this model to the analysis of heat transfer in the grooved heat pipes [3–7], microchannel evaporators [8] and microfin-tube evaporators [9]. This model was also used for the analysis of heat transfer during nucleate boiling [10–12] and evaporation of liquid droplets [13]. The basic equations for the evaporating thin liquid film are reduced to the fourth order ordinary differential equation for the liquid film thickness d with respect to the coordinate x measured along the heated surface. At the boundary between the adsorbed ⇑ Corresponding author. Tel./fax: +81 92 938 0570. E-mail address: [email protected] (H. Honda). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.08.044 0017-9310/Ó 2014 Elsevier Ltd. All rights reserved.
film region and the micro region, x ¼ 0, where the evaporation rate is zero, the boundary conditions are given by [7,8]
d ¼ d0
ð1Þ
dd=dx ¼ 0
ð2Þ
2
2
ð3Þ
3
3
ð4Þ
d d=dx ¼ 0 d d=dx ¼ 0
where d0 is the thickness of the adsorbed film. According to Wayner [1], d0 is given by
d0 ¼
AT v 6pql hlv DT
1=3 for d0 < 10 nm
ð5Þ
where A is the Hamaker constant of the system consisting of vapor, liquid film and solid substrate. Solution of the above mentioned fourth order differential equation subject to the boundary conditions given by Eqs. (1)–(4) results in the trivial solution of d ¼ d0 all over the surface [7,8]. In order to avoid this difficulty, Xu and Carey [3] and Ha and Peterson [6] adopted a simplified model in which the capillary effect on the liquid pressure was neglected. Khrustalev and Faghri [5] adopted a simplified model in which the radius of curvature of the evaporating liquid film r was assumed to be equal to that of the intrinsic meniscus rb all over the micro
830
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
Nomenclature A B1, B2 hlv m P Pd Pm DP Q q r T DT w x, y
Hamaker constant of vapor–liquid film-substrate system, J parameters in Eq. (29) latent heat of evaporation, J/kg mass flow rate of liquid, kg/m pressure, Pa disjoining pressure, Pa recoil pressure, Pa total pressure difference (¼ Pv Pl ), Pa heat flow rate, W/m heat flux, W/m2 radius of curvature of liquid film, m temperature, K wall superheat (¼ T w T v ), K groove width, m coordinates along and normal to the groove wall, respectively, Fig. 1, m
e
angle, Fig. 1, ° thermal conductivity, W/m K dynamic viscosity, Pa s density, kg/m3 surface free energy, J/m2
k
l q r
Subscripts b connecting point between micro region and intrinsic meniscus region c bottom of intrinsic meniscus i vapor–liquid interface l liquid max maximum value new new value old old value v vapor w wall 0 adsorbed film
Greek symbols d liquid film thickness, m
region. Stephan and Busse [4], Wang et al. [8] and Stephan and Hammer [10] solved the set of basic equations by the shooting method, in which one of the boundary values at x ¼ 0 was perturbed so that r approached rb at a large x. Do et al. [7] solved the basic equations with the boundary conditions of
r ¼ rb
ð6Þ 0
dd=dx ¼ tan e
ð7Þ
at x ¼ 0, where tan e was the slope of the liquid film at the boundary between the micro region and the intrinsic meniscus region and x0 was the coordinate measured in the x direction from the boundary. The value of r b was assumed to be a function of the liquid film thickness at x0 ¼ 0. The end of the micro region was determined by the condition that the evaporation rate there was equal to 0. Among the solution methods described above, that of Stephan and Hammmer [10] is considered to be most accurate, because it does not include additional assumptions for the boundary values except for the choice of the perturbed parameter. Honda and Wang [9] solved the fourth order differential equation for the liquid film thickness by the finite difference method. The boundary conditions at x ¼ 0 were Eq. (2) and 0
d ¼ d0r
ð8Þ
where d0r (> d0 ) was the thickness of the liquid film considering the curvature effect. The geometrical condition that the liquid film in the micro region was smoothly connected with the intrinsic meniscus region at x ¼ xb was given by Eq. (6) and
dd=dx ¼ tan e
ð9Þ
This solution method is attractive because the boundary conditions at both ends of the micro region are satisfied correctly. However, the obtained solution was not realistic on the point that, at x ¼ 0, the liquid film thickness was not continuous and the curvature of the liquid film took the maximum value. Thus more work is required to decide the combination of boundary conditions at x ¼ 0 and x ¼ xb which gives the most realistic solution. As described in [1], the value of the Hamaker constant depends on the combination of vapor, liquid and solid substrate which constitutes the evaporation system. In the previous studies, however,
sufficient information on the physical properties required to calculate A was not available. In the present study, the value of A is estimated using a theoretically based method. Another problem that has not been addressed in the previous studies is the Marangoni effect which was caused by the surface energy change of a very thin liquid film along the surface as a result of the change in the liquid film thickness. In the present study, the interfacial shear stress due to the Marangoni effect is taken into account in the basic equation. The fourth order differential equation for the liquid film thickness is solved by the finite difference method developed by Honda and Wang [9]. Numerical calculations are conducted for three sets of the combination of boundary conditions at x ¼ 0 and x ¼ xb and the numerical results are compared with each other. Numerical solutions are obtained for the evaporation of FC-72 filled in a small rectangular groove made of copper. Comparison is made between the cases in which the Marangoni effect is taken into account and neglected. Comparison is also made with the previous results of Raj et al. [13] for the evaporation of a FC-72 droplet on a heated surface in which the shooting method of Stephan and Hammer [10] was adopted. 2. Theoretical model Consider evaporation of a highly wetting apolar liquid filled in a rectangular groove. Fig. 1 shows the physical model and the coordinates. The coordinate x is measured downward along the groove from the boundary between the adsorbed film region and the micro region, and the coordinate y is measured vertically outward from the groove surface. Depending on the amount of liquid filled in the groove, the coordinate x ¼ 0 is located at the edge of the groove or on the groove wall. The coordinate xc denotes the position of liquid meniscus at the center of the groove. For the micro region, the momentum equation in the x-direction is written as 2
ll
d u dy
2
¼
dP l dx
ð10Þ
The boundary conditions are
u ¼ 0 at y ¼ 0
ð11Þ
831
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
where r1 denotes the value of r at d ¼ 1. The value of r decreases monotonically, from that of the substrate to that of the bulk liquid, with increasing liquid film thickness [17]. Thus Eq. (19) can be written as
r ¼ rl
A
ð20Þ
12pd2
where rl is the surface free energy of the bulk liquid (i.e., surface tension). In the previous studies, the dependence of r on d was overlooked and r was assumed to be equal to rl . It is relevant to estimate here the magnitude of the second term on the right hand side of Eq. (20) relative to the first term. The second term takes the largest value at the smallest value of the liquid film thickness; i.e., d ¼ d0 . For one of the conditions of numerical calculation discussed in Section 3, the values of rl and A=12pd20 are calculated to be 9.83 mN/m and 1.07 mN/m, respectively. Thus the second term on the right hand side of Eq. (20) is not negligible at a small value of d. The interfacial shear stress due to the Marangoni effect is given by
dr A dd ¼ dx 6pd3 dx
si ¼
Substitution of Pl and into Eq. (13) yields
Fig. 1. Physical model and coordinates.
ll
at y ¼ d
ð12Þ
where si denotes the interfacial shear stress due to the Marangoni effect. Integration of Eq. (10) subject to the boundary conditions given by Eqs. (11) and (12) yields the solution for u. Then integration of u over the liquid film thickness yields the expression for the liquid mass flow rate m as follows:
m ¼ ql
q d3 dPl d2 udy ¼ l þ si ll 3 dx 2
d
0
! ð13Þ
The P l is expressed in terms of P v , the disjoining pressure Pd , the 2 capillary pressure r=r and the recoil pressure P m ¼ q2 =qv hlv as follows:
Pl ¼ Pv Pd
r r
þ Pm
ð14Þ
where r is the surface free energy of the system consisted of vapor, liquid film and heated wall. Following Wayner [1], P d is approximated by
Pd ¼
A
pffiffiffiffiffiffiffiffiffi Al As
ð16Þ
where Al and As are the Hamaker constants of liquid and solid, respectively. The liquid curvature 1=r is related to d by 2
ql d3 d A r þ Pm 3ll dx 4pd3 r
m¼
The P d is related to
dr Pd ¼ dd
ð17Þ
r by the following equation [16]:
q ¼ hlv
12pd2
ð19Þ
ð23Þ
dm dx
ð24Þ
Substitution of m given by Eq. (22) into Eq. (24) yields
q¼
q l h lv d 3 d A r d þ Pm dx 4pd3 r 3ll dx
ð25Þ
Assuming a linear temperature drop across the liquid film, q is given by
kl ðT w T i Þ d
ð26Þ
where Ti is the surface temperature of the liquid film. Following Wayner et al. [18] and assuming that the accommodation coefficient is equal to unity, q is also written as
q¼
2 pRv T v
1=2
qv h2lv
Tv
Ti Tv þ
Tv A r þ P m ql hlv 6pd3 r
ð27Þ
Eliminating T i from Eq. (26) and Eq. (27) yields
ð18Þ
A
Thus, for the present case in which the Marangoni effect is included, the disjoining pressure term is virtually larger than that for the case without the Marangoni effect by 50%. The heat flux q is related to m by
q¼
Substitution of P d given by Eq. (15) into Eq. (18) and then integration the resulting equation from d ¼ d to d ¼ 1 yield
ð22Þ
ql d3 dDP ql d3 d A rl ¼ þ Pm 3ll dx 3ll dx r 6pd3
2
1 d d=dx ¼ r f1 þ ðdd=dxÞ2 g3=2
If the Marangoni effect is neglected, m is given by
q¼
Neglecting cohesion in the vapor, A is given by the following equation [1]:
r ¼ r1
m¼
ð15Þ
6pd3
A ¼ Al
si given by Eqs. (14) and (21) respectively,
du ¼ si dy
Z
ð21Þ
d þ qkl Thv2
v lv
kl T A r DT þ v þ Pm 3 pRv T v 1=2 ql hlv 6pd r
ð28Þ
2
2
Substitution of Pm ¼ q2 =qv hlv into Eq. (28) and then solution of the resulting equation for q yields
q¼
1 ðB1 2
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi B21 4B2 Þ
ð29Þ
832
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
1=2 q q h3 q q h3 where B1 ¼ vk Tl v lv d þ qkl Thv2 pR2v T v and B2 ¼ v T vl lv DTþ l v lv qv h2lv 6pAd3 rr . Equating q given by Eqs. (25) and (28), the fourth order differential equation for d is derived as follows:
ql hl d 3 d A r d þ Pm dx 4pd3 r 3ll dx ¼
d þ qkl Thv2
kl T A r DT þ v þ Pm 3 pRv T v 1=2 ql hlv 6pd r
ð30Þ
2
v lv
It is relevant to note here that if the Marangoni effect is neglected, the terms A=4pd3 and r=r in Eq. (30) are replaced by A=6pd3 and rl =r, respectively. Thus, for the case without the Marangoni effect, the differential equation for d is given by
ql hl d 3 d A rl d þ Pm dx 6pd3 3ll dx r ¼
d þ qkl Thv2
v lv
kl T A rl DT þ v þ P m pRv T v 1=2 ql hlv 6pd3 r
d þ qkl Thv2
kl T A r DT þ v pRv T v 1=2 ql hlv 6pd3 r
ð34Þ
2
Substitution of q ¼ 0, d ¼ d0r and r ¼ r 0 at x ¼ 0 yields
ð31Þ
3. Numerical calculation Since one of the objectives of the present study is to compare the numerical results obtained by the finite difference method with the previous results obtained by the shooting method, numerical calculations are conducted for the conditions adopted by Raj et al. [13]; i.e., evaporation of FC-72 at 32 °C on the copper surface with the wall superheats of 1 K, 2 K, 4 K and 8 K. For apolar liquids, Al can be estimated by the following approximation [14]:
ð32Þ
Thus the value of Al for FC-72 at 32 °C is calculated to be 2:06 1020 J. The value of As for copper is obtained from the table of measured data given by Eichenlaub et al. [15] as 2:73 1019 J. Substitution of these values in Eq. (16) yields A ¼ 5:40 1020 J. This value is somewhat smaller than A ¼ 3:96 1020 J assumed by Raj et al. [13]. It should be mentioned that Raj et al. [13] defined the Hamaker constant by A=6p instead of A. Eq. (30) is a highly non-linear differential equation. This equation was solved by the implicit finite difference method developed by Honda and Wang [9]. The unsteady term ql @d=@t was added to the left hand side of Eq. (30) and the solution was obtained as the steady state of an unsteady liquid flow subject to the assumed initial distribution of d and the boundary conditions at both ends of the micro region. The liquid in the thin film is driven by the combined pressure gradient and interfacial shear stress and the effect of the assumed initial distribution disappears as the calculation proceeds. By choosing an appropriate value for the time step, the steady state solution is obtained relatively easily. As described in Section 2, the boundary conditions at x ¼ xb are given by Eqs. (6) and (9). For the rectangular groove, the geometrical condition that the liquid film in the micro region is smoothly connected with the intrinsic meniscus region at x ¼ xb is given by the following equation:
r b ¼ ðw=2 db Þ= cos e
q¼
v lv
2
Al ¼ 2:1 1018 rl J
Case B: d ¼ d0r and dd=dx ¼ 0 at x ¼ 0; r ¼ rb and dd=dx ¼ e at x ¼ xb . 2 2 Case C: d ¼ d0 , dd=dx ¼ 0 and d d=dx ¼ 0 at x ¼ 0; r ¼ rb and dd=dx ¼ e at x ¼ xb . Among the three sets of the boundary conditions, Case B was assumed by Honda and Wang [9]. For Case A and Case C, it follows from Eq. (17) that the capillary pressure at x ¼ 0 is equal to 0. For Case B, the capillary pressure at x ¼ 0 is not zero [9]. In order to obtain the capillary pressure and the liquid film thickness at x ¼ 0, an iterative procedure was adopted as follows: For a low heat flux where the recoil pressure is negligible as compared with the other pressure components, Eq. (28) can be approximated by
ð33Þ
Thus the boundary values db and rb must satisfy Eq. (33). In order to obtain the solution, two more boundary conditions are required at x ¼ 0. This indicates that only two of the four physically reasonable boundary conditions at x ¼ 0, given by Eqs. (1)–(4), are satisfied by the solution. In the present study, three sets of the boundary conditions were assumed. These were 2 2 Case A: d ¼ d0 and d d=dx ¼ 0 at x ¼ 0; r ¼ rb and dd=dx ¼ e at x ¼ xb .
1=3 6p ql hlv DT r d0r ¼ Tv A r0
ð35Þ
Since the value of r 0 was not known a priori, an iterative calculation was required to obtain the solution. At the beginning of the numerical calculation, the value of 1=r0 was assumed to be 0. Then, after the solution was obtained, the calculated value of 1=r0 was substituted into Eq. (35) to obtain the new assumption of d0r . The numerical calculation was repeated until a converged solution was obtained. For Case C, it was estimated that a converged solution could not be obtained, because the number of the boundary conditions was larger than the required number. In order to prove this estimate, the numerical calculation was also conducted for Case C. Since the values of the geometrical parameters xb and e were not given a priori, numerical calculations were conducted assuming a range of values for these parameters. Then the numerical results were compared with each other to determine the most appropriate values of these parameters. The groove width w was changed in the range of 0.1–1.0 mm. The grid size was assumed to be 5 1010 m or less for 0 6 x 6 5 108 m depending on the condition and it was increased gradually with increasing x. The starting value of xb was assumed to be 109 —108 m depending on the condition. At the beginning of the numerical calculation, a very small value was assumed for the time step Ds and it was increased with time. The local film thickness changed with time and approached the steady state value. The numerical calculation was continued until the convergence criteria j1 dold =dnew j < 104 and jql hlv ðdnew dold Þ=Dsqj < 104 were satisfied at all grid points. Since the finite difference solution obtained by the present method satisfies the boundary conditions at both ends of the calculation domain correctly, the errors in the calculated values of d, q, r etc. are estimated to be less than 0.01%. After a converged solution was obtained, the calculation was repeated increasing the value of xb . The numerical results revealed that the continuous distribution of d was not obtained when xb was increased beyond some value. The maximum value of xb for which the converged solution was obtained, xb;max , was determined by using the bisection method as follows: (1) Denote the largest value of xb for which the converged solution was obtained as xb;old , and the next assumption of xb for which the converged solution was not obtained as xb;old . (2) Conduct the numerical calculation using xb ¼ ðxb;old þ xb;old Þ=2 as the new assumption of the micro region length. (3) Depending on the condition whether the converged solution was obtained or not obtained, denote xb as either xb;old or xb;old . (4) Return to step (2) and repeat the numerical calculation until the convergence criterion j1 xb;old =xb;old j < 104 was satisfied. (5) Denote the last value of xb;old as xb;max .
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
833
4. Numerical results and discussion Converged solutions were obtained for all sets of the boundary conditions; i.e., Case A, Case B and Case C. For Case C, where the assumed value of the capillary pressure at x ¼ 0 was equal to 0, the calculated value of the capillary pressure took the maximum value at the first nodal point x ¼ x1 next to x ¼ 0 and the distribution of the capillary pressure in the region x P x1 was similar to that for Case B shown in Fig. 7. Thus the distribution of the capil2 2 lary pressure (and so the distribution of d d=dx ) was discontinuous between x ¼ 0 and x ¼ x1 . This indicated that the solution satisfying the set of boundary conditions of Case C did not exist. As described in Section 3, the numerical solutions were obtained for a range of the micro region length. Figs. 2–4, respectively, show the distributions of the liquid film thickness, the heat flux and the meniscus profiles for different values of the micro region lengths obtained for the set of boundary conditions of Case A. The conditions of numerical calculation are DT ¼ 1 K, w ¼ 0:2 mm and e ¼ 30 . The largest value of the micro region length in these figures (i.e., xb ¼ 17:7 lm) corresponds to xb;max . The liquid film thickness is smaller for larger xb . The slope of the liquid film thickness curve is close to 0 near x ¼ 0 but it increases rapidly with increasing x. It is seen from Fig. 3 that the heat flux takes the maximum value in the region x 6 0:065 lm. The maximum heat flux is almost unchanged for different xb , whereas the coordinate at the maximum heat flux point increases as xb
Fig. 4. Comparison of meniscus profiles for different micro region lengths.
increases. Comparison of the heat flux distributions for different xb reveals that the highest total heat transfer rate is obtained by the case of xb ¼ xb;max . The principle of maximum entropy production has been found to give successful predictions for the steadystate properties of a variety of dynamic systems [19]. Thus, if no geometrical constraint is imposed on the meniscus profile, the most probable solution is estimated to be that for xb ¼ xb;max . It is seen from the meniscus profiles shown in Fig. 4 that the apparent contact angle of the meniscus decreases as xb increases. As a result, the coordinate xc at the bottom of the meniscus is larger for larger xb . In Fig. 4, the smallest value and the largest value of xc are 0.058 mm and 0.071 mm, respectively. Thus, if the liquid level in the groove is smaller than 0.058 mm, the steady state solution of the evaporating thin film does not exist. For the liquid level in between 0.058 mm and 0.071 mm, the steady state solution will be that for which the value of xc is equal to the liquid level. For the liquid level larger than 0.071 mm, the steady state solution will be that for xb ¼ xb;max . In the discussion that follows, only the results for xb ¼ xb;max are presented. Figs. 5–7, respectively, compare the distributions of the heat flux, the total pressure difference and the capillary pressure between the solutions obtained for the set of boundary conditions of Case A and Case B. The conditions of the numerical calculation are DT ¼ 4 K, w ¼ 0:2 mm and e ¼ 30 . The heat flux curve has a positive slope at x ¼ 0 for both cases. The heat flux increases rapidly with increasing x, takes the maximum value at x 0:037 lm, and
Fig. 2. Comparison of liquid film thicknesses for different micro region lengths.
Fig. 3. Comparison of heat fluxes for different micro region lengths.
Fig. 5. Comparison of heat fluxes between the solutions obtained for the set of boundary conditions of Case A and Case B.
834
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
numerical results reveals that dd=dx ¼ 0:0128 at x ¼ 0. For the numerical results presented in Figs. 5–7, the ratio of d0r (given by Eq. (35)) to d0 (given by Eq. (5)) is calculated to be 1.020. Thus, for Case B, both the liquid film thickness and the liquid curvature are not continuous at x ¼ 0. From the comparison of the numerical results for Case A, Case B and Case C described above, it may be concluded that the most realistic solution is obtained by the set of boundary conditions of Case A. In Figs. 5–7, the solution of Raj et al. [13] is also shown for comparison. The solution was obtained for the micro region of an evaporating liquid droplet on a flat surface The boundary conditions at x ¼ 0 were d ¼ d0 , dd=dx ¼ 0, DP ¼ A=6pd30 and Q ¼ Q 0 , where Q was the perturbed parameter whose physical meaning was the heat flow rate through the liquid film at x ¼ 0. The first boundary condition is the same as that assumed in Case A, whereas the second boundary condition is the same as that assumed in Case B. Since DP is defined by DP ¼ 6pAd3 þ rrl P m , the third boundary Fig. 6. Comparison of total pressure differences between the solutions obtained for the set of boundary conditions of Case A and Case B.
then decreases monotonically with further increasing x. The difference in the heat flux between Case A and Case B is most significant in the region 0 < x < 0:004 lm. At x ¼ 0:004 lm, the value of q for Case B is larger than that for Case A by 14% and the difference increases as x decreases. The maximum heat flux and the coordinate at the maximum heat flux point are slightly smaller for Case A. It is seen from Fig. 6 that the slope of the total pressure difference curve is negative all over the surface for both Case A and Case B. This indicates that the mass flow rate through the liquid film, given by Eq. (22), is negative all over the surface. For Case A, the amount of liquid outflow at x ¼ 0 is calculated to be 0.67% of the total evaporation rate. However, it is not probable that such an outflow actually occurs. This may indicate the limitation of the macroscopic formulation in dealing with the evaporation phenomena from a very thin liquid film. Most significant difference between Case A and Case B is seen in the distributions of the capillary pressure shown in Fig. 7. For Case A, the capillary pressure is zero at x ¼ 0. It increases as x increases, then takes a maximum value, and then decreases monotonically with further increasing x. For Case B, on the other hand, the capillary pressure takes a maximum value at x ¼ 0. It decreases sharply as x increases. Then it takes a minimum value, and then a maximum value, and then decreases monotonically with further increasing x. The above results show that, for Case A, both the liquid film thickness and the liquid curvature are continuous at x ¼ 0. However, the slope of the liquid film is not continuous at x ¼ 0. In the adsorbed film region x < 0 where d ¼ d0 , it is apparent that dd=dx ¼ 0. On the other hand, examination of the
condition is correct only for the case of 6pAd3 rr0 . The value of 0
Q 0 was determined so that r ¼ r b was satisfied at x ¼ xb . According to Raj et al. [13], the exact values of r b and xb did not have a noticeable influence on the numerical solution. This is in contrast to the present results in which the solution was considerably influenced by the assumed value of xb . As seen from Figs. 5–7, a large difference is observed between the present results and the results of Raj et al. [13]. Most significant difference is seen in the distributions of calculated values near x ¼ 0. For the present results, the calculated values change sharply in the region 0 6 x 6 0:015 lm. For Raj et al. [13], on the other hand, the changes are very small in the region 0 6 x 6 0:02 lm. This is due to the difference in the solution method between the two cases. Raj et al. [13] assumed an infinitesimal perturbation Q 0 for Q . Thus the initial variations of the calculated values from the boundary values at x ¼ 0 were very small. It is relevant to note here that when Q 0 ¼ 0 is assumed, no change occurs for the calculated values all over the surface. The calculated values change first very slowly, and then more and more rapidly, with increasing x. This indicates that the effect of the assumed value of Q 0 decreases as x increases. In Fig. 7 for the capillary pressure, the calculated value by Raj et al. [13] is close to zero in the region 0 < x < 0:01 lm. Raj et al. [13] have shown a figure in which the distributions of the pressure components and the total pressure difference were plotted on the semi-logarithmic scale. According to the figure, the capillary pressure shows a trend similar to the present solution for Case B in the region close to x ¼ 0. This is probably due to the fact that the same boundary condition of dd=dx ¼ 0 at x ¼ 0 were adopted in both cases. However, the capillary pressure obtained by Raj et al. [13] is three orders of magnitude smaller than that for Case B. This was due to the fact that an infinitesimal perturbation Q 0 was assumed for Q . For this case the boundary values satisfy the inequality 6pAd3 rr0 and the boundary 0
condition DP ¼ A=6pd30 holds correctly. However, for the present numerical results for Case B, r=r0 is not negligible as compared to
Fig. 7. Comparison of capillary pressures between the solutions obtained for the set of boundary conditions of Case A and Case B.
A=6pd30r (see Figs. 6 and 7). If the shooting method is used to obtain the solution for this case, the numerical calculation must be conducted with both Q 0 and 1=r 0 as the perturbed parameters. Thus, it will not be easy to obtain the solution satisfying the boundary conditions at x ¼ xb . It may be concluded from the above discussion that the finite difference method is superior to the shooting method on the point that the boundary conditions at both ends of the micro region are satisfied correctly. In the discussion that follows, only the results for xb ¼ xb max obtained for the set of boundary conditions of Case A are presented. Figs. 8–10, respectively, compare the distributions of the liquid film thickness, the heat flux and the total pressure difference between the cases in which the Marangoni effect is included and
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
Fig. 8. Comparison of liquid film thicknesses between the solutions in which the Marangoni effect is considered and neglected.
Fig. 9. Comparison of heat fluxes between the solutions in which the Marangoni effect is considered and neglected.
835
Marangoni effect in the x direction and a smaller liquid film thickness than the case without the Marangoni effect is required to maintain the liquid flow. Comparison of the results between the two cases reveals that the largest difference in the liquid film thickness occurs at x ¼ 0:029 lm. At this point, the liquid film thickness for the case without the Marangoni effect is larger than that for the case with the Marangoni effect by 34.5%. It is seen from Fig. 9 that the largest difference in the heat flux between the two cases occurs at x ¼ 0:013 lm, where the heat flux for the case without the Marangoni effect is higher than that for the case with the Marangoni effect by 17.7%. It is interesting to note here that the film thickness at x ¼ 0:013 lm is a little bit larger for the case without the Marangoni effect (see Fig. 8). Thus, the above results indicate that the temperature drop across the liquid film, T w T i , is larger for the case without the Marangoni effect than the case with the Marangoni effect by more than 17.7%. The two heat flux curves cross at x ¼ 0:035 lm. In the region of x > 0:035 lm, the heat flux for the case with the Marangoni effect is lager than that for the case without the Marangoni effect. As a result, the difference in the total heat transfer rate between the two cases is small. For Fig. 9, the calculated values of the total heat transfer rates obtained by integrating the heat flux distributions over the x direction are 1.585 W/m and 1.564 W/m for the cases with and without the Marangoni effect. It is seen from Fig. 10 that the slope of the total pressure difference curve near x ¼ 0 is steeper for the case without the Marangoni effect. This indicates that a higher pressure drop than the case with the Marangoni effect is required to maintain the liquid flow. Figs. 11 and 12, respectively, compare the distributions of the liquid film thickness and the heat flux for different slopes of the liquid film at the boundary between the micro region and the intrinsic meniscus region. The other conditions of the numerical calculation are DT ¼ 4 K and w ¼ 0:2 mm. Due to the geometrical condition at x ¼ xb given by Eq. (33), the numerical solution was not obtained for e 6 15 . For 20 6 e, the liquid film thickness is larger for larger e. However, in the range of 30 6 e 6 55 , the liquid film thickness distribution is almost unchanged. The values of xb;max corresponding to the assumed values of e are also shown in Figs. 11 and 12. For 20 6 e 6 40 , xb;max increases rapidly as e increases. Then xb;max increases gradually with further increasing e. It is interesting to note here that the values of xb;max are considerably different for the solutions in the range of e of 30 6 e 6 55 which show almost the same liquid film thicknesses distribution. This indicates that the liquid film profiles for these solutions are close to the same circular arc in the region of x P 2:56 lm. It is
Fig. 10. Comparison of total pressure differences between the solutions in which the Marangoni effect is considered and neglected.
neglected in the basic equation; i.e., solutions of Eq. (30) and Eq. (31), respectively. The conditions of the numerical calculation are DT ¼ 4 K, w ¼ 0:2 mm and e ¼ 30 . The liquid film thickness is smaller for the case with the Marangoni effect. For this case the liquid is pulled by the interfacial shear stress due to the
Fig. 11. Comparison of liquid film thicknesses for different slopes of liquid film at the boundary between micro region and intrinsic meniscus region.
836
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
Fig. 12. Comparison of heat fluxes for different slopes of liquid film at the boundary between micro region and intrinsic meniscus region.
Fig. 14. Comparison of meniscus profiles for different groove widths.
seen from Fig. 12 that the maximum heat flux point is almost unchanged for 20 6 e 6 55 and it moves toward a smaller x with further increasing e. The total heat transfer rate is larger for smaller e. Since the solutions are obtained for a wide range of e, a physically acceptable criterion is required to determine the most appropriate value of e. In the region close to x ¼ xb , the twodimensionality of the heat conduction and the liquid flow becomes more and more significant as e increases. Thus the assumption of one-dimensional heat conduction and liquid flow adopted in the theoretical model is not appropriate for large values of e. On the other hand, for small values of e, xb;max decreases rapidly as e decreases. As a result, the liquid film thickness at x ¼ xb;max , db , also decreases rapidly and a very high evaporation heat transfer occurs near x ¼ xb;max . It is seen from Fig. 11 that for e ¼ 20 and e ¼ 25 , the liquid film thickness at x ¼ xb;max is less than 0.06 lm. Since the evaporating mass in the micro region is supplied from the intrinsic meniscus region, a high liquid flow velocity occurs in the intrinsic meniscus region near x ¼ xb;max . Thus the assumption of constant liquid pressure in the intrinsic meniscus region does not hold any more. These results indicate that too large values and too small values of e are not appropriate for the boundary condition. It may be concluded from the above discussion that the appropriate value of e to be assumed in the numerical calculation is 30°. Figs. 13 and 14, respectively, compare the distributions of the heat flux and the meniscus profile for different groove widths.
The other conditions of numerical calculation are DT ¼ 4 K and e ¼ 30 . It is seen from Fig. 13 that the heat flux distribution in the region 0 6 x 6 0:5lm is almost unchanged irrespective of the groove width. As expected from the heat flux distribution, the meniscus profiles for different w, shown in Fig. 14, are almost unchanged in the region close to x ¼ 0. However, due to the geometrical condition at x ¼ xb , the meniscus profiles for different groove widths tend to circular arcs with the radii of curvature rb given by Eq. (33) as x increases. Fig. 15 compares the heat flux distributions for different wall superheats. The other conditions of numerical calculation are w ¼ 0:2 mm and e ¼ 30 . It is seen from Fig. 15 that the maximum heat fluxes for DT ¼ 1 K, 2 K, 4 K and 8 K are 1.02 MW/m2, 2.10 MW/m2, 4.28 MW/m2 and 8.64 MW/m2, respectively. Thus the maximum heat flux is roughly proportional to the wall superheat. However, this does not mean that the liquid film thickness at the maximum heat flux point is unchanged for different DT. As seen from Eq. (29), q is a complex function of DT, d and r. The maximum heat flux point moves toward a smaller x as the wall superheat increases. In Fig. 15, the results of Raj et al. [13] are also shown for comparison. The general trend of the data are similar to the present case except that the maximum heat flux point is shifted to a larger x and the maximum heat fluxes are about 15% higher than the present values. In order to check the possibility that the difference in the maximum heat fluxes between the present results and Raj et al. [13] was due the difference in the assumed value of A, numerical calculations were also conducted using the A
Fig. 13. Comparison of heat fluxes for different groove widths.
Fig. 15. Comparison of heat fluxes for different wall superheats.
H. Honda, O. Makishi / International Journal of Heat and Mass Transfer 79 (2014) 829–837
value assumed by Raj et al. [13]. The numerical results showed that the two solutions were quite close to each other, indicating that the solution was rather insensitive to the value of A. Thus the difference in the maximum heat flux between the present results and Raj et al. [13] was probably due to the difference in the method of numerical calculation, which resulted in different distributions of the liquid film thickness.
837
given by Eq. (33) as x increased. Comparison of the present results obtained by the finite difference method and the previous results obtained by the shooting method revealed that the finite difference method was superior to the shooting method on the point that the boundary conditions at both ends of the micro region were satisfied correctly. Conflict of interest
5. Conclusions None declared. The basic equations describing the behavior of an evaporating thin film of a highly wetting liquid in a groove is reduced to the fourth order ordinary differential equation for the liquid film thickness. At the boundary between the adsorbed film region and the micro region, there are four physically reasonable boundary conditions given by Eqs. (1)–(4). On the other hand, at the boundary between the micro region and the intrinsic meniscus region, there are two geometrical conditions given by Eqs. (6) and (9). Thus there are a number of possibilities for the selection of the combination of boundary conditions at the boundary between the adsorbed film region and the micro region. In the present study, three sets of the combination of boundary conditions, Case A, Case B and Case C, were examined. The Hamaker constant was estimated by using a theoretically based method. The numerical solutions were obtained by the finite difference method developed by Honda and Wang [9]. Comparison of the numerical results revealed that the most realistic solution was obtained by the set of boundary conditions of Case A. The numerical results showed a small amount of liquid outflow at the boundary between the adsorbed film region and the micro region. This seemed to indicate the limitation of the macroscopic formulation in describing the evaporation phenomena from a very thin liquid film. The solutions were obtained for a wide range of the micro region length. As a result, the meniscus profile was considerably different depending on the micro region length. On the basis of the principle of maximum entropy production, the most probable solution, for the case where no geometrical constraint was imposed on the meniscus profile, was estimated to be that for the largest micro region length. The numerical results showed that the Marangoni effect due to the thickness change of a very thin liquid film was not negligible. The solutions were obtained for a wide range of the assumed value of the slope of liquid film at the boundary between the micro region and the intrinsic meniscus region. On the basis of physical consideration on the numerical results, the appropriate value of the slope to be assumed in the numerical calculation was concluded to be 30°. For the groove widths of 0.1–1.0 mm, the numerical solutions in the small x region were almost unchanged. But the meniscus profiles tended to circular arcs with different radii of curvatures r b
References [1] P.C. Wayner Jr., Intermolecular forces in phase-change heat transfer: 1998 Kern award review, AIChE J. 45 (1999) (1998) 2055–2068. [2] M. Potash Jr., P.C. Wayner Jr., Evaporation from a two-dimensional extended meniscus, Int. J. Heat Mass Transfer 15 (1972) 1851–1863. [3] X. Xu, V.P. Carey, Film evaporation from a micro-grooved surface – an approximate heat transfer model and its comparison with experimental data, J. Thermophys. Heat Transfer 4 (1990) 512–520. [4] P.C. Stephan, C.A. Busse, Analysis of the heat transfer coefficient of grooved heat pipe evaporator walls, Int. J. Heat Mass Transfer 35 (1992) 383–391. [5] D. Khrustalev, A. Faghri, Heat transfer during evaporation on capillary-grooved structures of heat pipes, Trans. ASME J. Heat Transfer 117 (1995) 740–747. [6] J.M. Ha, G.P. Peterson, The interline heat transfer of evaporating thin films along a micro grooved surface, Trans. ASME J. Heat Transfer 118 (1996) 747– 755. [7] K.H. Do, S.J. Kim, S.V. Garimella, A mathematical model for analyzing the thermal characteristics of a flat micro heat pipe with a grooved wick, Int. J. Heat Mass Transfer 51 (2008) 4637–4650. [8] H. Wang, S.V. Garimella, J.Y. Murthy, Characteristics of an evaporating thin film in a microchannel, Int. J. Heat Mass Transfer 50 (2007) 3933–3942. [9] H. Honda, Y.S. Wang, Theoretical study of evaporation heat transfer in horizontal microfin tubes: stratified flow model, Int. J. Heat Mass Transfer 47 (2004) 3971–3983. [10] P. Stephan, J. Hammer, A new model for nucleate boiling heat transfer, Wärme und Stoffübertragung 30 (1994) 119–125. [11] J.H. Lay, V.K. Dhir, Shape of a vapor stem during nucleate boiling of saturated liquids, ASME J. Heat Transfer 117 (1995) 394–401. [12] C. Kunkelmann, P. Stephan, CFD simulation of boiling flows using the volumeof-fluid method within openFOAM, Numer. Heat Transfer, Part A 56 (2009) 631–646. [13] R. Raj, C. Kunkelmann, P. Stephan, J. Plawsky, J. Kim, Contact line behavior for a highly wetting fluid under superheated conditions, Int. J. Heat Mass Transfer 55 (2012) 2664–2675. [14] J.N. Israelachvili, Intermolecular and Surface Forces, second ed., Academic Press, London, 1992. [15] S. Eichenlaub, C. Chan, S.P. Beaudoin, Hamaker constants in integrated circuit metallization, J. Colloid Interface Sci. 248 (2002) 389–397. [16] B.V. Derjaguin, N.V. Churaev, Structural component of disjoining pressure, J. Colloid Interface Sci. 49 (1974) 249–255. [17] R.J. Waltman, A. Khurshudov, G.W. Tyndall, Autophobic dewetting of perfluoropolyether films on amorphous-nitrogenated carbon surfaces, Tribol. Lett. 12 (2002) 163–169. [18] P.C. Wayner Jr., Y.K. Kao, L.V. LaCroix, The interline heat transfer coefficient of an evaporating wetting film, Int. J. Heat Mass Transfer 19 (1976). 487-49. [19] R.K. Niven, Steady state of a dissipative flow-controlled system and the maximum entropy production principle, Phys. Rev. E 80 (021113) (2009) 1–15.