Stagnation Point of Surface Flow during Drop Evaporation - Langmuir

May 2, 2018 - At the start of evaporation, surface temperature was monotonous with no extrema (region A). As the contact angle θ decreased near regio...
2 downloads 3 Views 2MB Size
Subscriber access provided by Kaohsiung Medical University

New Concepts at the Interface: Novel Viewpoints and Interpretations, Theory and Computations

Stagnation Point of Surface Flow During Drop Evaporation Lihui Wang, and Michael T Harris Langmuir, Just Accepted Manuscript • DOI: 10.1021/acs.langmuir.8b00627 • Publication Date (Web): 02 May 2018 Downloaded from http://pubs.acs.org on May 2, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Stagnation point of surface flow during drop evaporation Lihui Wang and Michael T. Harris∗ Davidson School of Chemical Engineering, Purdue University 480 Stadium Mall Drive, West Lafayette, IN 47907-2100 E-mail: [email protected] Abstract Capillary flow and Marangoni flow influence flow patterns of an evaporating liquid drop. While it is obvious that Marangoni stress on the drop surface affects the surface flow direction, we found that capillary flow also has an impact. The numerical results of this study showed a stagnation point near the contact line, which was further explained by the lubrication theory. The stagnation point is produced by the competing effects of Marangoni flow and capillary flow, and emerges when the contact angle is small, because the divergence of the capillary flow near the contact line increases as the contact angle decreases. The radial position of the stagnation point from the numerical results (rnumerical ≈ 0.995) agreed with the experimentally observed stagnation point (rexperimental > 0.992).

Introduction Microfluidic flow in sessile drop evaporation is worth studying because it not only affects common natural processes like coffee stains, but also affects various applications like inkjet printing, 1,2 deposition of DNA sequencing micoarrays, 3,4 manufacture of novel electronic 1

ACS Paragon Plus Environment

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

materials, 5,6 nanochromatography, 7 and other evaporative self-assembly techniques. 8–10 In these processes, flow patterns in evaporating drops affect the colloidal deposition. Therefore, understanding the development of these flow patterns is crucial to elucidating and controlling deposition patterns of colloidal particles in evaporating drops. The evaporation flux along the spherical sessile drop surface has been widely studied and is considered controlled by diffusion. The evaporation flux has been shown to diverge near a pinned contact line. 11,12 The divergent evaporation flux and the pinned contact line together result in the capillary flows (Figure 1a) toward the contact line replenishing the evaporating liquid at the contact line. 13–15 Marangoni flow is caused by the surface tension change due to temperature change. The evaporative cooling flux along the drop surface is divergent near the contact line because it is proportional to the evaporation mass flux. Apart from evaporative cooling, the temperature at the drop-air interface is also affected by the heat diffusion from the substrate to the drop surface, because the drop surface temperature is lowered by evaporative cooling. The evaporative cooling and the heat diffusion together lead to a gradient of the temperature along the drop surface which causes a change of the surface tension leading to a Marangoni flow. 14 Because the surface tension is negatively correlated with the temperature and the Marangoni flow is dragged to the high surface tension region, temperature-induced Marangoni flow moves toward the low temperature region. A circulation flow profile (Figure 1b) is formed from the two flows, the capillary flow and the Marangoni flow. This has been observed in experiments with water drops 16–18 as well as other organic drops. 19,20 This circulation flow, compared to pure capillary flow, changes the deposition time 21 and the deposition patterns 20 of colloidal drops. The direction of this circulation changes under various conditions. Hu and Larson did both numerical and analytical calculations of velocity profiles, and showed that the circulation flow changes from counterclockwise to clockwise as the drop evaporates. 22 Later Ristenpart et al. gave an analytical result showing that the direction of the circulation de-

2

ACS Paragon Plus Environment

Page 2 of 23

Page 3 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

(a)

(b)

Figure 1: Flow profiles of (a) capillary flow only and (b) capillary flow and Marangoni flow pends on the relative thermal conductivity kR , which is a ratio of the substrate conductivity to the liquid conductivity. 8 Furthermore, Xu et al. improved the criteria by including the thickness of the substrate (d) on which the drop rests. They proved that the circulation direction is determined by a new parameter RN =

1 d , kR R

where R is the radius of the contact

line. 23 However, these results assumed a monotonous change in the temperature on the drop surface from the apex of the drop to the contact line. In contrast, Zhang et al. numerically solved the temperature field with a commercial package and indicated a nonmonotonic distribution of the surface temperature. 24 Similarly, Barash used the finite difference method to solve the transport equations to compute both the flow profile and the surface temperature of the drop. Their results manifested the existence of multi-vortices aroused by the nonmonotonic surface temperature change. 25,26 This implies that stagnation points exist on the drop surface where the surface flow is zero between two vortices. Stagnation points were also observed in the experiments of water drop evaporation by different groups, 16,17,27 but the cause of the stagnation points was either not mentioned or not validated. One proposed explanation was that a stagnation point exists where the direction of the Marangoni stress changes on the drop surface. 16,25 However, despite the fact that the stagnation point existed within 20µm away from the contact line in the experiments, 16,27 the analytical results and the numerical results with the Marangoni-stress explanation all

3

ACS Paragon Plus Environment

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

predicted that the stagnation points were much farther away from the contact line. Another explanation was proposed by Xu et al. that a stagnation point would appear when the Marangoni effect was strong enough which was indicated by a large Marangoni number (a ratio of the Marangoni stress to the viscous stress). 28 The location of the stagnation point predicted by this paper was closer to the actual experimentally observed location, but still not within 20µm away from the contact line. Therefore, the previous theories and the experimental results did not agree in terms of the location of the stagnation point. The numerical results in the current paper clearly show the effect of the capillary flow on stagnation points at the drop surface near the contact line. Thus, the surface flow direction cannot be determined solely based on the Marangoni stress. In comparison to previous research, 21,22,26 our simulations gave more detailed information of the fluid motion and temperature variations near the contact line in the later stages of drop evaporation (when the contact angle is less than 15o ). Elucidation of the fluid flow profile near the contact line at later stages is significant because, according to the experiments of Devlin et al., 21 particle deposition mainly happens at the late stage of evaporation. Furthermore, particle accumulation near the contact line involves more interesting phenomena, such as particle separation when particles of different sizes are involved. 29,30

Methods We set up a transient, axisymmetric model consisting of a set of dimensionless equations which are explained as follows. The variables used to calculate the dimensionless groups and the characteristic values in the model are listed in Table 1. The dimensionless equations were transformed from dimensional equations by choosing characteristic values as Table 2. The transformation introduced the dimensionless groups shown in Table 3. Inside the drop, four dimensionless variables were computed, which were velocity v (including radial velocity u and axial velocity v), pressure p and temperature T . The Navier-

4

ACS Paragon Plus Environment

Page 4 of 23

Page 5 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Table 1: Simulation Parameters Variable Definition R drop base radius ρ liquid density µ viscosity kl liquid thermal conductivity α thermal diffusivity T0 initial temperature σ0 initial surface tension β d˜ σ /dT˜ cv vapor saturation concentration Hum relative humidity D diffusivity of vapor in air Hv latent heat of vaporization ks substrate thermal conductivity ls half length of the substrate d substrate thickness

Value 1 × 10−3 9.97 × 102 8.90 × 10−4 0.58 1.39 × 10−7 25 71.91 -0.1657 2.32 × 10−2 0.5 26.1 × 10−5 2264.76 0.80 1.3 × 10−3 0.15 × 10−3

Units m kg/m3 Pa·s W/m/K m2 /s ◦ C dyn/cm dyn/cm/K kg/m3 m2 /s kJ/kg W/m/K m m

Table 2: Characteristic Variables Variable length velocity surface tension

Definition lc = R c vc = Dc ρlc σc = σ 0

Value 1 × 10−3 6.5 × 10−7 71.91

temperature stress vapor concentration

Tc = Hv Dcc /k τc = µvc /lc cc = cv

2.5 5.8 × 10−7 2.32 × 10−2

Units Dimensionless variable m l = ˜l/lc m/s u = u˜/vc , v = v˜/vc dyn/cm σ = (˜ /σc σ − σ0 )  K T = T˜ − T0 /Tc Pa kg/m3

τ = τ˜/τc c = c˜/cc

Table 3: Dimensionless Groups Variable Definition Re Re = ρvc lc /µ Ca Ca = µvc /σc Pe P e = vc lc /α c REH REH = HkvlDc Tc kR kR = ks /kl ˜ P eM a P eM a = vM a h/α

5

Value 7.4 × 10−4 8.1 × 10−9 4.7 × 10−3 2.5 1.4 ˜ change with h

ACS Paragon Plus Environment

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 23

Stokes equation and the continuity equation were used to solve for the motion of the incompressible fluid in the drop:  Re ·

∂v + v · Ov ∂t

 =O·T

(1)

O·v =0

(2)

where T = −pI + [Ov + (Ov)T ]. Since the model was axisymmetric, the azimuth coordinate could be neglected. Here the force balance included the inertia and the viscous stress, as well as the capillary force and the Marangoni stress incorporated through the boundary condition on the liquid-gas interface

Can · T = −2Hσn + Os σ,

(3)

where n is the normal vector along the interface, 2H is twice the mean curvature, and σ is the dimensionless surface tension. The Marangoni stress term Os σ implies the change of surface tension, and this change comes from the temperature change as

dσ dT

= β·

Tc , σc

where β is constant for a small range of temperature change, which was satisfied in this problem. 8 Gravity was neglected because it was small compared to the capillary force 15 or the Marangoni stress. 22 Another boundary condition at the interface, the kinematic boundary condition, was given by

n · (v − vs ) = J,

(4)

where J is the evaporation flux solved from quasi-steady state vapor diffusion O2 c = 0, and is given by

J = J0(θ) (1 − r2 )−λ(θ) ,

6

ACS Paragon Plus Environment

(5)

Page 7 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

where r is the dimensionless distance from the axis of symmetry (0 6 r 6 1) and θ is the contact angle. The form of Equation (5) was given by Hu and Larson, 12 and was transformed into dimensionless form here as J0(θ) = (1 − Hum)(0.27θ2 + 1.30)[0.6381 − 0.2239(θ − π/4)2 ], where Hum is the relative humidity, and λ(θ) = 0.5 − θ/π. The thermal energy in the drop was modeled by the heat convection-diffusion equation:  Pe

∂T + v · OT ∂t



= O2 T.

(6)

The boundary conditions were given at the liquid-vapor interface for evaporative cooling as (7) and at the liquid-substrate interface for heat flux balance as (8), where relative vaporization heat REH =

Hv Dcc k l Tc

is the ratio of evaporative cooling to heat conduction, and relative

heat conductivity kR = ks /kl is the ratio of substrate heat conductivity to the liquid heat conductivity.

REH · J = −n · OT   1 n · OTs − OT = 0 kR

(7) (8)

In the boundary condition (8), Ts is the substrate temperature solved from the substrate heat conduction. The lower boundary of the substrate was assumed to be constant at the initial temperature T0 . There assumed to be no heat flux (n · OT = 0) at the substrate-vapor interface and at the end of the substrate (r = ls where ls is the half length of the substrate). These boundary conditions are also shown in Figure 2. The initial conditions were v = 0, θ = 50◦ and T = 0 (viz. T˜ = T0 , where T˜ is the dimensional temperature). The model was solved numerically by the Galerkin finite element method (G/FEM) with the same discretization method used in the work of Widjaja and Harris 31 and Devlin et al.. 21 However, we further improved the mesh distribution to pack more elements near the contact line so that not only does the calculation continued to converge as the contact angle decreased till θ ≈ 3o , but also the velocity profiles were tracked as close as 0.1µm near the contact line. 7

ACS Paragon Plus Environment

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

mass balance: n ⋅(v − vs ) = J (r,θ )

drop-air : interface

momentum balance: Can ⋅T = −2H σ n + ∇ sσ

the axis of symmetry u=0 ∂v / ∂r = 0 n ⋅∇T = 0 n

Page 8 of 23

heat balance: REH ⋅ J(r,θ ) = −n ⋅∇T

v(u,v), p, T

n

drop-substrate interface:

n ⋅(∇Ts − 1 / kR ⋅∇T ) = 0

z

v=0

r

n ⋅∇Ts = 0 n n

n ⋅∇Ts = 0

Ts

n ⋅∇Ts = 0 Ts = 0 (T!s = T0 )

Figure 2: Boundary conditions (More details on mesh distribution, numerical convergence and accuracy can be referred to in supporting information)

Results and Discussion Surface Temperature Effect of Heat Convection Both heat convection and diffusion influence surface temperature. As shown in Figure 3, the numerical solution of the heat convection-diffusion equation and the solution of heat diffusion equation were compared. The figure showed the effect of heat convection on the drop surface temperature profile as the contact angle θ changes. When θ > 15◦ , the surface temperature is significantly different with or without heat convection. The difference made by the heat convection comes from Marangoni flow. If the Marangoni flow is not incorporated to the calculation, the Peclet number (ratio of heat convection to diffusion) P e =

vc lc α

∼ 10−2 is small, thus heat convection is negligible. 22 In contrast, if the

Marangoni flow is incorporated, a new dimensionless group, P eM a , needs to be defined to

8

ACS Paragon Plus Environment

Page 9 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

24.8

24.7

24.6

T(◦ C) 24.5

θ = 40◦ θ = 30◦ θ = 20◦ θ = 16◦ θ = 10◦ θ = 5◦

24.4

With convection Only diffusion

24.3 0

0.2

0.4

r

0.6

0.8

1

Figure 3: Surface temperature under various θ changing as evaporation goes. Heat convection is negligible when θ < 15◦ . The explanation to the sharp increase near r = 1 is attached in the supporting information.

compare heat convection and diffusion, instead of P e. The velocity scale in P eM a is a Marangoni flow velocity, not vc (=

Dcc ) ρlc

which is essentially

the characteristic velocity of the capillary flow. This was because the Marangoni flow significantly increases the fluid velocity, compared to having only capillary flow. When Marangoni flow exists, the numerical solution of Navier-Stokes equation showed that the dimensionless axial velocity v = v˜/vc (with v˜ representing the dimensional axial velocity) was as large as 103 ; this tells that the Marangoni velocity was as large as 103 times the characteristic velocity vc . The radial Marangoni velocity uM a was derived from a balance of the Marangoni stress and the viscous stress (Figure 4a) as in (9); the axial Marangoni velocity vM a was derived in (10) based on incompressibility. uM a Tc ∼µ R H vM a uM a ∼ H R

β

(9) (10)

The length scale in P eM a was the drop thickness H instead of lc = R, the radius of the drop base. This is because heat transport in the axial direction is more significant than

9

ACS Paragon Plus Environment

Langmuir

τM a ∼ β TRc

Ma τviscous ∼ µ uH

0.4

0.3

(a) The balance of the Marangoni stress and the viscous stress

H 1 − cos(θ) = R sin(θ)

H/R

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 23

0.2

0.1

0 5

10

15

20

25

30

35

40

θ (°)

(b) Temperature profile in the drop when θ = 38◦ : temperature difference is larger in the axial direction than in the radial direction.

(c) The dependence of the aspect ratio on θ

Figure 4: (a) a diagram of stress balance; (b)temperature profile; (c)the aspect ratio of the drop in the radial direction. As shown by Figure 4b, temperature difference is larger between the drop surface and the drop base than that of the radial direction. Because of this same reason, the velocity scale should be the axial Marangoni velocity vM a instead of the radial Marangoni velocity uM a . So the new dimensionless number characterizing the influence of heat convection is P eM a =

vM a H . α

It scales as P eM a ∼

βTc R µα

· (H )3 ∼ 104 ( H )3 after substituting Equation R R

(9) and (10), whose value depends on the aspect ratio H/R. In this case, when H/R > 0.1, P eM a > 10. So H/R > 0.1 is the condition when heat convection is important. The dependence of the aspect ratio on θ was derived from the geometry of a spherical cap as (11), and was shown by Figure 4c.

H/R = (1 − cosθ) /sinθ

10

ACS Paragon Plus Environment

(11)

Page 11 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

When θ > 15◦ , H/R > 0.1. This is a situation when heat convection should make a difference. Validated by Figure 3, surface temperatures with and without heat convection are obviously different for θ > 15◦ . At the late evaporation stage when θ < 10◦ , H/R < 0.1. This is a situation when heat convection can be neglected. It is again supported by Figure 3 that surface temperatures were almost the same with and without heat convection. Therefore, heat convection plays an important role in determining the surface temperature distribution when the contact angle is large (θ > 15◦ for water drop with R = 1mm), and can be neglected only at the later stage of evaporation.

Extremum Points Surface temperature change is not monotonous when the drop is flat (at a small contact angle) as shown in Figure 3. Extrema (dT /ds = 0) exist when the temperature is observed from the apex of the drop to the contact line along the liquid-vapor interface. These are also where dσ/ds = 0. The direction of Marangoni stress on the liquid-gas interface is opposite on either side of an extremum. As shown in Figure 5a, Marangoni stress τM a =

σc Oσ lc s

is in

the same direction as s (called outward in the following context) if dσ/ds > 0, and inward if dσ/ds < 0. Since dσ/ds is negatively corelated to dT /ds, Marangoni stress is outward when dT /ds < 0. The location of the extrema under different θ are shown in Figure 5b. At the start of evaporation, surface temperature was monotonous with no extrema (region A). As the contact angle θ decreased near Region B during evaporation, an extremum emerged. Then the extremum bifurcated into two extrema at different locations r for a certain θ in region B. dσ/ds < 0 between the two extrema. Later as θ decreased to region C, only one extremum still existed for the surface temperature. Given a pair of (θ, r), dσ/ds < 0 if it is on the right side of the extremum line, and dσ/ds > 0 if it is on the left side.

11

ACS Paragon Plus Environment

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 23

(a) The direction of Marangoni stress

1

0.8

r

dσ/ds > 0 C

0.6

B

A

0.4

dσ/ds < 0 0.2

0

5

10

θ(◦ )

15

(b) Locations of extremum points (dσ/ds = 0) As drop evaporates, θ decreases from the right to the left on the abscissa axis. Given a pair of (θ, r), dσ/ds < 0 if it is on the right side of the extremum line, and dσ/ds > 0 if on the left side.

Figure 5: Extremum points (dσ/ds = 0) vs. Marangoni stress direction

12

ACS Paragon Plus Environment

Page 13 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Flow Profiles and Stagnation Points Flow profiles change during drop evaporation (Figure 6). At first with a large contact angle, a perfect single counter-clockwise circulation flow was formed as Figure 6a. This circulation flow maintained till extremum points emerged on the surface. Then it changed to a multicirculation flow (shown with θ = 14◦ in Figure 6b). Later, a clockwise circulation flow was developed as shown in Figure 6c. Under multi circulations, stagnation points of surface flow (t · v = 0, with t as the tangential vector shown in Figure 5a) exist. The location of these points are shown in Figure 7a. Similar to Figure 5b, the surface flow is inward on the right side of the line and outward on the left side. When θ = 20◦ , (r, θ) is on the right side of the line for any r, which means that the surface flow is inward at any location of the drop surface, implying the counter-clockwise circulation. When θ = 14◦ , there were three stagnation points. The surface velocity changed direction across each stagnation point. As in Figure 6b, from the top to the contact line along the surface, the flow changed from inward to outward and to inward again and to outward again. The stagnation points at r < 0.95 were caused by the change of Marangoni stress direction. The comparison of the plot of the extremum points and of the stagnation points in Figure 8a showed that the locations of the two points coincided for the r < 0.95 portion. Since the extremum points are the points where Marangoni stress changes directions, these stagnation points can be called Marangoni-induced stagnation points. For example, when θ = 15◦ , there were two Marangoni-induced stagnation points according to Figure 8a (the points where red line and blue line coincide). The flow profile near one of the Marangoniinduced stagnation points is shown in Figure 7b. However, some of the stagnation points can not be explained by the directional change of the Marangoni stress as the two lines do not coincide for r > 0.99 in Figure 8a. These are the stagnation points close to the contact line. For example, when θ = 10◦ , there were two stagnation points, but only one extremum point. The stagnation point around r ≈ 0.95 13

ACS Paragon Plus Environment

Langmuir

Page 14 of 23

20.016

14.003

contact angle = 20.016

contact angle = 14.003

(a)

(b)

0.2

θ = 20◦

z 0.1

0

10.004

0.2

0.4

r

(c)

θ = 14◦

0.1

0.2

0

0.2

0.6

0.8

0

1

0

0.2

0.4

0.6

r

0.8

1

0.1

θ = 10◦ z 0

0

0.2

0.4

r

0.6

0.8

1

Figure 6: Flow profiles of a drop semi-section under various θ (not to scale) The flow pattern is still a single counter-clockwise circulation when θ = 20◦ . Later as contact angle = 15.000 Marangoni stagnation points emerge, it changes to a multi-circulation flow. At the late stage, a clockwise circulation flow is developed.

1

θ = 15◦

z

0.1

0.8

0.05

Outward flow (t · v < 0)

0

0.6

0.5

0.6 contact angle = 10.003

0.7

r

r (b) flow profile near the Marangoniinduced stagnation point

0.4

Inward flow (t · v < 0)

0.004

0.2

θ = 10◦ z

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0.3

z

0.3

0.002

0

8

10

12

14

16

τM a

18

θ(◦ ) 0

(a) the locations of the stagnation points (t · v = 0) r as θ decreases with evaporation The surface flow is inward on the right side of the line and outward on the left side.

0.98

0.99

(c) flow profile near the capillaryinduced stagnation point (not to scale)

Figure 7: Locations of stagnation points and flow profiles near them

14

1

r

ACS Paragon Plus Environment

Page 15 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

coincided with the extremum point, thus it was Marangoni-induced. The other one around r ≈ 0.995 was an extra stagnation point. The flow field around this stagnation point is shown by Figure 7c. Such a stagnation point was also observed experimentally with an evaporating water drop with a contact angle of 10◦ by Xu and Luo, 16 and existed less than 17µm away from the contact line, giving r > 0.992 (in a dimensionless form). Therefore, the value of r in the simulation for this stagnation point agreed well with experiments. capillary-induced stagnation points 1

1

extremum points (dσ/ds = 0) stagnation points (t · v = 0)

0.8

0.6

r 0.4

dσ/ds > 0 Outward flow (t · v > 0)

0.2

0

5

extremum points Stagnation points

0.8

r

0.6

offset of Marangoni-induced stagnation points

0.4

dσ/ds < 0 Inward flow (t · v < 0)

10

θ(◦ )

0.2

5

15

(a) Water

10

θ(◦ )

15

20

(b) Octane

Figure 8: Comparison between locations of stagnation points and extremum points The locations of stagnation points are not exactly where τM a = 0 due to capillary flow.

Capillary-Induced Stagnation Point The stagnation points at r > 0.99 are caused by the divergence of capillary flow near the contact line. To study the influence of this divergence on stagnation points, a different boundary condition was used to eliminate the divergence. Widjaja and Harris 31 studied capillary flow with no Marangoni stress incorporated, and found that radial velocities were much less with a constant evaporation flux imposed across the drop-gas interface. Thus, we imposed a constant evaporation flux similar to Widjaja and Harris’s work but incorporated Marangoni stress, and found that the surface flow was inward with no stagnation point 15

ACS Paragon Plus Environment

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 23

near the contact line at θ = 10◦ (Figure S6). Therefore, only if the capillary flow velocity increases rapidly (caused by the divergent evaporation flux) near the contact line does the extra stagnation point exist. So it is named a capillary-induced stagnation point. The essence of the capillary-induced stagnation point is that the surface velocity is in the opposite direction to the Marangoni stress between the stagnation point and the contact line (Figure 7c). This was explained by the lubrication theory near the contact line (Figure 9). The lubrication theory can be used because the contact angle was small (θ < 15◦ ) when stagnation points emerged near the contact line. The height-averaged radial velocity for capillary flow was derived from mass balance as Equation (12), 15 where A(θ) is a coefficient that changes as the drop evaporates. The mass balance is the same with or without Marangoni stress; therefore, the height-averaged radial velocity u¯r(r,θ) is essentially the same as (12) when Marangoni flow exists.

u¯r(r,θ) = A(θ)

−λ(θ) i 1h 1 − r2 − 1 − r2 r

(12)

The radial velocity profile at each value of r is parabolic in z, as ur(r,θ) = Az 2 + Bz + C, from the lubrication theory. The following boundary conditions (dimensionless form) were r = τM a with τM a as the Marangoni stress in s direction. applied: vr |z=0 = 0 and ∂v ∂z z=h(r,t) The surface velocity in s direction vs derived from these conditions is 3 1 vs = vr |z=h(r,t) = u¯r(r,θ) + τM a h(r, t). 2 4

(13)

lim u ¯r(r,θ) = ∞ with the physical divergent evaporation flux condition; τM a h(r, t) From (12), r→1 lim vs > 0 for whichever value of τM a . That is to say, the surface flow is bounded. Therefore, r→1

near the contact line is moving toward the contact line even when τM a is negative (away from the contact line). Such a case did not happen for a constant evaporation flux condition (as shown by the numerical results of Figure S6 in the supporting information), because u¯r(r,θ) was also bounded in (13) under this condition. 16

ACS Paragon Plus Environment

Page 17 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Figure 9: Lubrication theory: When τM a is negative, vs can be positive or negative based on the value of u¯r(r,θ)

Offset of Marangoni-Induced Stagnation Points The phenomenon that the surface flow is in the opposite direction to Marangoni stress does not only happen where capillary flow velocity diverges (between capillary-induced stagnation points and contact lines), but also where Marangoni stress is small (close to extremum points). It is demonstrated by the offset of the line of Marangoni-induced stagnation points from the line of extremum points as shown in Figure 8b. Between the two lines are the conditions when dσ/ds < 0 yet the surface flow is outward (t·v > 0). For example, when the contact angle of the octane drop was 7◦ , the locations of the Marangoni-induced stagnation point and the corresponding extremum point differ by more than 0.1R. τM a and the surface flow were in opposite directions between the two points (Figure 10a). This indicates that the surface flow is not dominated by Marangoni stress. This offset can also be explained by the lubrication theory. According to (13), vs is positive when τM a > − 38

u ¯r(r,θ) . h(r,t)

In other words, for a negative τM a ,when

|τM a |