pubs.acs.org/Langmuir © 2009 American Chemical Society
Structure of Circulation Flows in Polymer Solution Droplets Receding on Flat Surfaces Yu Yoshitake, Shohei Yasumatsu, Masayuki Kaneda, Koichi Nakaso, and Jun Fukai* Department of Chemical Engineering, Graduate School of Engineering, Kyushu University, Motooka 744, Fukuoka 819-0395 Japan Received August 31, 2009. Revised Manuscript Received October 24, 2009 In a previous report where internal flows were experimentally visualized in polymer solution droplets receding on a lyophobic surface [Kaneda et al., Langmuir 2008, 24, 9102-9109], the direction of the circulation flow was found to depend on solvent and solute concentration. To identify the reason for this finding, the internal flow in the droplet is investigated numerically. A mathematical model predicts that double circulation flows initiate after a single flow develops at high Marangoni numbers, while only a single circulation flow develops at low Marangoni numbers. The dependencies of the calculated velocities on the solvent and the initial solute concentration agree qualitatively with experiment. It is concluded that the difference of the flow directions that were investigated experimentally is due to such a change in the flow structures. The effects of the contact angle and dimensions on transport phenomena in a droplet are also discussed.
Introduction Fundamental knowledge of fluid dynamics and heat and mass transfer in microscopic droplets on substrates has been of recent interest for controlling the self-assembly of nanoparticles and proteins, formation of thin films, etc. Deegan et al.1 developed a simple mathematical model that considers mass conservation in the radial direction under the assumptions of uniform solute concentration and velocity fields in the vertical direction. They succeeded in using the model to predict formation of a ring stain. The improved models have been frequently used to predict the morphology of the films formed.2,3 Subsequently, the equations of continuity and momentum in the cylindrical coordinate system were numerically solved to discuss the fluid dynamics in droplets.4,5 In addition, the effect of thermal Marangoni forces on the fluid dynamics was investigated numerically.6,7 Hu et al.8 reported the Marangoni effect on local particles density in droplets. Though it is an experimental study, Kajiya et al.9 successfully measured the transient distribution of the mean solute concentration in the vertical direction and investigated the development of the ring stain. The aforementioned studies focused on the fluid flows after the contact line of the droplet was pinned. In contrast to these studies, there are only a few studies targeting the flows before self-pinning. Kaneda et al.10 experimentally visualized the circulation flows in polymer solution droplets receding on a lyophobic surface. They reported that (a) the circulation flow is not observed with no solute present; (b) the velocity of the circulation flow increases *To whom correspondence should be addressed. E-mail: jfukai@ chem-eng.kyushu-u.ac.jp. Tel.: þ81-92-802-2744. (1) Deegan, R. D.; Bakajin, O.; Dupont, T. F.; Huber, G.; Nagel, S. F.; Witten, T. A. Phys. Rev. E. 2000, 62, 756–765. (2) Morozumi, Y.; Ishizuka, H.; Fukai, J. J. Chem. Eng. Jpn. 2004, 37, 778–784. (3) Ozawa, K.; Nishitani E., E.; Doi, M. Jpn. J. Appl. Phys. 2005, 44, 4229–4234. (4) Fischer, B. J. Langmuir 2002, 18, 60–67. (5) Hu, H.; Larson, R. G. Langmuir 2005, 21, 3963–3971. (6) Hu, H.; Larson, R. G. Langmuir 2005, 21, 3972–3980. (7) Girard, F.; Antoni, M.; Faure, S.; Steinchen, A. Langmuir 2006, 22, 11085– 11091. (8) Hu, H.; Larson, R. G. J. Phys. Chem. B 2006, 110, 7090–7094. (9) Kajiya, T.; Kaneko, D.; Doi, M. Langmuir 2008, 24, 12369–12374. (10) Kaneda, M.; Hyakuta, K.; Takao, Y.; Ishizuka, H; Fukai, J. Langmuir 2008, 24, 9102–9109.
Langmuir 2010, 26(6), 3923–3928
with increasing initial solute concentration; (c) the direction of the circulation flow is upward on the axis of symmetry in most cases; (d) the direction of the flow on the axis of symmetry is reversed at high solute concentration. Observations a-c suggest that the circulation flow is dominated by solutal forces rather than thermal forces. Kaneda et al. numerically concluded that items a-c are derived from the solutal Marangoni and Rayleigh forces, not from thermal forces. However, no reasonable explanation for observation (d) has yet been found. A main purpose of the present study is to identify the reason for observation (d) using a numerical technique. In the previous study10 the calculation domain was assumed to be a hemisphere in the spherical coordinate system, to restrict the contact angle of the droplet to 90°. In the current study, the governing equations are given in the cylindrical coordinate system to observe the droplets at various contact angles. The modeling for the local evaporation rate on the vapor-liquid interface is also improved. The calculations are performed for polystyrene-anisole and polystyreneacetophenone solution droplets receding on silicon wafer. A goal of this study is to quantitatively investigate the selfpinning mechanism because it is one of the key factors to dominate the film morphologies. The self-pinning is known to be controlled by contact angle. On the other hand, the droplets used in the previous study are much larger than those used by industrial processes because it is difficult to visualize the velocity fields in the inkjet droplets. Accordingly the effects of contact angle and dimensions of the droplet on transport phenomena are also discussed. Numerical Model. Consider a sessile droplet evaporating on a semifinite substrate on the plane r-z in the cylindrical coordinate system. The solute concentration dependencies of the fluid density F, fluid viscosity μ, surface tension σ, and mass diffusion coefficient DAB are considered. The subscript A and B indicate solute and solvent, respectively. The other physical properties are assumed to be constants. μ and DAB are represented by μ = μl0fμ(cA) and DAB = DAB0fD (cA), respectively, where cA is the solute concentration. The subscript 0 indicates the initial condition. fμ(cA) and fD(cA) represent the solute-concentration dependencies of the viscosity and diffusion coefficient. The computation domain for the fluid flow and solute concentration
Published on Web 11/23/2009
DOI: 10.1021/la903245m
3923
Article
Yoshitake et al. Table 1. Physical Properties of Liquids and Silicon Substrate
acetophenone anisole substrate (Si)
density F [kg/m3]
heat capacity Cp [J kg/K]
thermal conductivity k [W/(m 3 K)]
latent heat L [J/kg]
diffusion coefficient of solute in vapor Dv [m2/s]
1.023 103 0.989 103 2.326 103
1.704 103 1.930 103 0.810 103
0.1471 0.156 162
4.45 105 4.16 105
6.84 10-6 7.05 10-6
Table 2. Nondimensional Numbers Corresponding to d0 = 1.05 mm solute
cA0
Pr
Sc
Ras
Mas
anisole anisole anisole acetophenone acetophenone acetophenone
0.005 0.01 0.03 0.005 0.01 0.03
15.2 19.0 46.0 22.4 26.7 53.0
1.25 103 1.95 103 1.15 104 2.85 103 4.02 103 1.59 104
6.06 12.1 36.6 0.354 0.711 2.13
2.16 8.65 77.9 0.206 0.774 6.64
is the droplet, while that for the temperature is the droplet and substrate. During the receding process, the droplet is assumed to similarly shrink, keeping the contact angle at ψc. The instantaneous wetting radius rc(t) thus decreases from the initial radius rc0 with the lapse of time t. rc is selected as the characteristic length of the dimensionless variables to conveniently consider receding of the contact line. Solutal Marangoni and Rayleigh convections are considered. Boussinesq approximation is applied in the conservation equation of momentum (eq 2). The dimensionless governing equations for the dimensionless velocity vector U, temperature Θ, and solute concentration CA are r3U ¼ 0
ð1Þ
DU rc ¼ -rP þ Pr½r 3 fμ frU þ ðrUÞT g þ Dτ rc0
3 Pr Ras CA ð2Þ
DΘi Ri 2 ¼ r Θi Dτ Ri0
ði ¼ l, wÞ
DCA fD ¼ r 2 CA Dτ Le
ð3Þ ð4Þ
where D/Dτ (= ∂/∂τ þ UR∂/∂R þ UZ∂/∂Z) is the substantial derivative operation, and subscripts l and w indicate the droplet and substrate, respectively. Uniform temperature and solute concentration are assumed as the initial conditions, that is, U = 0, T = T0, and CA = CA0. Considering the surface tension gradient on the free surface, the boundary condition for the fluid dynamics on the free surface is given by DU S rc DCA Mas ¼ ð5Þ U N ¼ 0, DN rc0 DS where N and S are the dimensionless coordinates normal and tangential to the free surface, respectively. The local evaporation rate of the solvent is given by the equation proposed by Hu et al.8 Considering loss of the solvent due to evaporation, the other boundary conditions on the free surface are given by DΘl rc BiΘl ¼ ð1 -R2 Þ -λ þ ð6Þ DN rc0
fD0
DCA Δc ¼ ð1 -R2 Þ -λ 1 þ CA c0 DN
3924 DOI: 10.1021/la903245m
ð7Þ
λ is a function of the contact angle ψc:8 λ ¼ ðπ -2Ψc Þ=ð2π -2Ψc Þ
ð8Þ
The boundary condition on the exposed surface of the substrate is given by the same equation as eq 8, but the first term on the right side is neglected. The boundary conditions at the liquid-solid interface (Z = 0, 0 e R e 1) are given by U R ¼ U Z ¼ 0,
kl
DΘ DΘ ¼ kw , DZ DZ
Θl ¼ Θw
ð9Þ
The dimensionless parameters in the above equations are defined as U = u/u0, P = p/(Fu02), τ = u0t/rc, R = r/rc, r* = rrc, Θ = (T - T0)/ΔT and CA = (cA - cA0)/ΔcA, where U is the velocity vector, P is pressure, and r is the deviation operation. The characteristic quantities are given by u0 = ν0/rc, ΔT = Lmr=0rc/ kl0, ΔcA = mr=0cA0rc/(DAB0Fl0), where ν is dynamic viscosity, and k is thermal conductivity. It should be noted that rc vanishes from ΔT and Δc if eq 10 is substituted into their definitions. The dimensionless numbers are defined by Pr ¼ ¼
ν Rl hrc0 , Le ¼ , Bi ¼ , Ras Rl0 DAB0 kl gð -1=FÞðDF=DcA Þr3c0 ΔcA ðDσ=DcA Þrc0 ΔcA , Mas ¼ ð10Þ ν0 Rl0 Rl0 μl0
Considering that the partial vapor pressures of the solvents in the bulk are nearly zero, the local evaporation rate mr=0 at r = 0 is given by8 2 ( )3 Dv cv π 2 4 mr ¼0 ¼ ð0:27Ψc þ 1:3Þ 1 - 0:2339 Ψc þ 0:3619 5 4 rc ð11Þ where Dv is the diffusion coefficient of the solvent in vapor, and cν is saturated vapor concentration of the solvent on the free surface. The dependencies of the solute concentration on the viscosities and surface tensions of the solutions are measured by a corn-andplate viscometer and Wilhelmy plate method, respectively. According to the measured results, they are approximated by ( expð34:41cA -6:419Þ ½Pa 3 s at cA < 0:10 kg=kg ð12Þ μ ¼ expð22:20cA -5:162Þ ½Pa 3 s at cA g0:10 kg=kg Dσ=DcA ¼ 5:82 10 -2 cA
at cA < 0:20 kg=kg
ð13Þ
for anisole solutions and ( expð34:33cA -6:923Þ ½Pa 3 s at cA < 0:05 kg=kg μ ¼ ð14Þ expð26:54cA -6:041Þ ½Pa 3 s at cA g0:05 kg=kg Dσ=DcA ¼ 5:60 10 -2 cA þ 4:00 10 -5 at cA < 0:20 kg=kg
ð15Þ
Langmuir 2010, 26(6), 3923–3928
Yoshitake et al.
Article
Figure 1. Time variation of the calculated temperature distribution for anisole solution at c0 = 0.005 kg/kg. The heat transfer is dominated by the conduction rather than the convection at least at cA0 < 0.03 kg/kg.
for acetophenone solutions. DAB is estimated using Wilke and Chan’s equation.11 The values of other physical properties are shown in Table 1. The calculations continue until 200 s. This time duration is much shorter than the drying process. However, this is an allowable limit because of the capacity of the computer used. The numerical procedure is described in the Appendix.
Results Temperature Distribution. According to the experimental condition,10 the initial droplet diameter d0 before deposit and the contact angle are given as 1.05 mm and 70°, respectively. These values correspond to rc0 = 0.76 mm. The dimensionless numbers observed are shown in Table 2. A typical result for the time variation of the temperature distribution for anisole solution is shown in Figure 1. The droplet is cooled from the free surface mainly due to the evaporation of solvent at the initial stage. On the other hand, the heat flow from the substrate prevents the solid-liquid interface from cooling. As a result the isothermal lines become nearly parallel to the substrate surface. The temperature distributions almost fully develop at 5.0 s, changing very little thereafter. The temperature differences along the axis of symmetry are 0.312 K for anisole solutions and 0.363 K for acetophenone solutions. These values are independent of cA0 even when the velocity increases with increasing cA0 (see Figures 2 and 3). Heat transfer is accordingly dominated by conduction heat transfer rather than convection. (11) Geankoplis, C. J. Transport Process and Separation Process Principles; Prentice Hall Professional Technical Reference: Upper Saddle River, NJ, 2003.
Langmuir 2010, 26(6), 3923–3928
Effect of Initial Solute Concentration. Figure 2 shows the time variations of isoconcentration contours and velocity vector for anisole solutions at cA0 = 0.005 kg/kg. A circulation flow has already developed at 1.0 s. The velocity increases with lapse of time. The circulation flow is initiated by the Rayleigh effect and then progressed by the Marangoni effect.10 It should be noted that the contact line in the figure is retained at R = 1 because the radial direction is dimensionalized with the instantaneous wetting radius rc as aforementioned. The change of the droplet volume is small even at 200 s. The ratio of the wetting diameter to the initial one (=rc/r0) at 200 s is 0.908 in this case. Figure 3 shows the results for cA0 = 0.03 kg/kg. A circulation flow also develops at 1.0 s. The velocities are larger than those for cA0 = 0.005 kg/kg because ∂σ/∂cA increases with increasing cA as shown in eqs 13 and 15. The solute concentration dependence of the fluid velocity was discussed in detail in the precious report.10 In this case, the velocity of the circulation flow is large enough to move the internal fluids with low solute concentrations toward the free surface at the inclination angle of about 45° as shown in Figure 3b. As a result, the solute concentration on the free surface, or the surface tension, is locally minimized at this point, resulting in double circulations as t > 10.0 s. The double circulation continues at least in the present time duration. It should be noted that the flow direction on the axis of symmetry changes from upward to downward when the double circulation is developed. Comparison of Figures 2 and 3 shows that the concentration distributions depend on the velocity fields. This fact suggests that mass transfer is dominated by not only molecular diffusion but also by convection. DOI: 10.1021/la903245m
3925
Article
Yoshitake et al.
Figure 2. Isoconcentration contours (left side) and velocity vector (right side). Anisole solution cA0 = 0.005 kg/kg; initial droplet diameter = 1.05 mm; rc/r0 = 0.908 at t = 200 s. A single circulation flow develops at low initial concentrations. The mass transfer rate is dominated by diffusion and convection.
Figure 3. Isoconcentration contours (left side) and velocity vector (right side). Anisole solution cA0 = 0.03 kg/kg; initial droplet diameter = 1.05 mm, rc/r0 = 0.908 at t = 200 s. A single circulation flow develops at low initial concentrations. The mass transfer rate is dominated by diffusion and convection.
The fluid velocity at Z = 0.3 on the axis of symmetry is plotted with time in Figure 4. This position corresponds to where the fluid velocities were experimentally measured.10 The positive velocity indicates upward flow on the axis of symmetry, while the negative velocity indicates downward flow. First, the maxima of experimental velocities were of the order 10-6-10-5 m/s for anisole solutions and 10-6 m/s for acetophenone solutions. The orders of the calculated velocities agree with the experiments. Second, the present model predicts that the velocities changed from positive to negative under some conditions. The negative velocity means the development of double circulation flows. Consequently, the observed change from positive to negative velocity does not show the change in the direction of a single circulation flow, but it shows the development of double circulation flows. Actually the experimental domain of the visualization was limited because of
refraction on the free surface. Thus the velocity field where another circulation appears could not be investigated if double circulation flow developed. When the double circulations develop, the absolute value of the velocity rapidly increases to a local maximum value, and then decreases. In Figure 3, the area of the circulation flow on the side of the contact line tends to extend toward the axis of symmetry. This change is more clearly found in the results for Ψc = 30° (see Figure 6). The model thus predicts that the double circulations finally return to a single circulation. This result agrees qualitatively with the experimental result that the negative velocity changes to positive velocity at a later stage of the drying process. The experimental critical values of cA0 between the single and double circulations were in the range 0.03-0.08 kg/kg for anisole
3926 DOI: 10.1021/la903245m
Langmuir 2010, 26(6), 3923–3928
Yoshitake et al.
Article
Figure 4. Time variation of the fluid velocity at Z = 0.3 on the axis of symmetry. The positive sign of the velocity indicates an upward flow while the negative sign indicates downward flow. The negative velocity shows the development of double circulation flows. When the double circulations develop, the flow is experimentally observed just as the direction of the circulation flow might be changed.
solutions, and >0.20 kg/kg for the acetophenone solutions. The calculations underestimate the critical value of cA0. This suggests that they underestimate the diffusion effect on the mass transfer, which must result in inaccuracies of the diffusion coefficients. In Figure 3, DAB0 = 3 10-10 m2/s is estimated by Wilke and Chan’s equation.11 However, no double circulation flows developed if DAB0 = 3 10-9 m2/s were assumed. Effect of Contact Angle. It is worthwhile investigating transport phenomena at small contact angles because it is more difficult to experimentally visualize the internal flow as contact angle decreases. Figure 5 shows the results for anisole solution at Ψc = 30° and cA0 = 0.02 kg/kg. The initial droplet diameter d0 is the same as that for Ψc = 70°. It is emphasized that a violent circulation flow develops near the contact line. This is because, according to eq 9, the local evaporation rate more rapidly increases near the contact line, resulting in a larger gradient of the solute concentration, or the surface tension. By investigating the results in detail, it is found that there are weak circulation flows on both sides of the strong flow at the initial stage. The strong circulation flow extends its area toward the axis of symmetry, finally to eliminate a weak circulation as schematically shown in Figure 6. Though the calculations were attempted for anisole solution qat cA0 = 0.03 kg/kg, the solutions of the nonlinear algebraic equations were not convergent, probably because of instability of the fluid flows. Effect of Droplet Dimension. Figure 7 shows the calculated results for anisole solution at cA0 = 0.03 kg/kg. The calculation conditions are the same as for Figure 2 except for the droplet diameter d0 = 70 μm. Mas = 5.35 and Ras = 1.17 10-2 in this case. The model predicts the development of the double Langmuir 2010, 26(6), 3923–3928
Figure 5. Isoconcentration contours (left side) and velocity vector (right side) for 0.02 kg/kg anisole solution. Contact angle=30°. rc/ r0 = 0.905 at t = 200 s. A violent circulation flow develops near the contact line. This is because, according to eq 11, the local evaporation rate increases more rapidly near the contact line, resulting in a larger gradient of the solute concentration.
Figure 6. Change in circulation flows with time. (a) The circulation flow in the middle is strong while those on both sides are weak. (b) The strong circulation extends its area toward the axis of symmetry, to vanish a weak circulation on the left side.
Figure 7. Isoconcentration contours (left side) and velocity vector (right side) for 0.03 kg/kg anisole solutions at τ = 30.3. rc/r0 = 0.904. The maximum velocities have the same order of magnitude as the large droplet in Figure 3d because the evaporation rate given by eq 11 increases with decreasing droplet size. DOI: 10.1021/la903245m
3927
Article
Yoshitake et al.
Figure 8. The computation domain descretized by finite elements. The numbers of elements are (a) 3400 and (b) 3544.
circulation flows. If the evaporation rate were independent of the droplet dimension, the fluid velocity would have decreased with decreasing droplet size. However, the local evaporation rate given by eq 11 increases with decreasing the droplet dimension. As a result, the maximum velocities for Figure 3d and Figure 7, which are the results at similar dimensionless time, have the same order of magnitude. Figure 7 also demonstrates that convection mass transfer still plays an important role in mass transfer. No double circulation flow is predicted to develop for a small droplet of acetophenone solution at cA0 = 0.03 kg/kg: Mas = 0.455 and Ras = 6.81 10-4. As numerically demonstrated in the previous study,10 the Rayleigh effect is smaller than the Marangoni effect except in the period of time before the circulation flow fully develops. Considering the results for the large and small droplets, the critical value of Ma is of the order of magnitude 100.
Conclusion The significant findings in this study are summarized as follows: (1) The reason for the change in the directions of the circulation flow experimentally observed at high initial solute concentrations is probably the development of double circulation flows. (2) The present model predicts the experimental results that the fluid velocity increases with increase of the initial solute concentration, and the sign of the characteristic velocity changes from positive to negative at a critical initial solute concentration. (3) The present model seems to underestimate the diffusion effect on the mass transfer. This is probably due to inaccuracy of the estimated diffusion coefficient. (4) The critical Marangoni number between single and double circulation flows has the order of magnitude 100 regardless of the inaccuracy of the diffusion coefficient. (5) At a low contact angle, three circulation flows are developed at the initial stage. The middle one furthermore develops while the others disappear. (6) The order of the velocity field does not depend strongly on the droplet size. This is because the evaporation rate increases in inverse proportion to the droplet diameter.
3928 DOI: 10.1021/la903245m
Appendix The corner between the solid surface and the free surface becomes sharp with decreasing the contact angle. A quadrangular finite element is largely distorted if it is adjusted on such a corner, which becomes a cause of the calculation errors. The mathematical model is solved numerically by invoking the artificial compressibility method,12 which allows triangular finite elements for resolving the fluid dynamics. An implicit method is utilized for the numerical integration in time. At each time step, nonlinear algebraic equations for the velocity and pressure fields are solved using the Newton-Raphson method, then the linear algebraic equations for the temperature and concentration fields are solved numerically. These procedures are iterated until the velocities converge with an allowable error (=10-5). The temperature and concentration fields were found to sufficiently converge when the velocities converge under the present conditions. Then the instantaneous values of rc are calculated by subtracting the amount of the evaporating solvent from the droplet volume. Typical grids used are shown in Figure 8. The finite elements near the interface are much smaller than those in the internal regions to satisfy the boundary conditions as accurately as possible. Especially, fine finite elements (=10-3) are used near the free surface.13 When the grid size is not sufficiently small, the conservation of mass becomes worse. The total mass of the solute was checked in each time step, and found to be changed by only 2% at 200 s. The number of the elements could be reduced to about two-third, while maintaining the same level of accuracy, by using fine elements near the free surfaces. The time step in the calculations is about 10-3 s. Even if it was twice, the calculated results did not change because the implicate method is utilized. (12) Morozumi, Y.; Nakamoto, M.; Fukai, J.; Miyatake, O. Kagaku Kogaku Ronbunshu 2002, 28, 55–60. (13) Kawahara, M; Hirano, H. Int. J. Numer. Methods Fluids 1983, 3, 137–163.
Langmuir 2010, 26(6), 3923–3928
