Regular Article 常规文章

Marangoni circulation in evaporating droplets in the presence of soluble surfactants
在可溶性表面活性剂存在下蒸发液滴中的马兰戈尼循环

2.1. Drop evolution

A droplet that is deposited on a substrate is considered. The contact angle $\theta$ is assumed ‘small’ ($\theta ⩽{40}^{\circ }$ as shown by Hu and Larson [50], [51]) and the density $\rho$ and dynamic viscosity $\mu$ are taken constant. Also, the typical drop height $H\le {10}^{-3}$ m, which implies that the Bond number $\mathit{Bo}=\frac{\rho {\mathit{gH}}^2}{\sigma }\ll 1$ (with gravitational acceleration g and surface tension $\sigma$). Thus, the effects of gravity can be disregarded [52], [53]. Note that this assumption only holds, because there are no possible effects of buoyancy in a single component droplet. For multicomponent droplets, gravity cannot be neglected as shown by Edwards et al. [54] and Li et al. [55]. The Reynolds number is much smaller than unity and a cylindrical coordinate system $(r,\alpha ,z)$ is used with the assumptions of no swirl (angular velocity $U_{\alpha }=0$) and axisymmetry ($\frac{\partial }{\partial \alpha }=0$).

Given this case, the Navier–Stokes equations can be rewritten into a 1D evolution equation for the height $h(r,t)$ as a function of radial coordinate r and time t [21]:(1)$\frac{\partial h}{\partial t}=\frac{1}{r\mu }\frac{\partial }{\partial r}\left(\frac{{\mathit{rh}}^3}{3}\frac{\partial p}{\partial r}-\frac{{\mathit{rh}}^2}{2}\frac{\partial \sigma }{\partial r}\right)+w_e\text{.}$Here, p denotes the pressure inside the droplet and $w_e$ the evaporative volume flux, which is negative. From this equation it can be recognized that fluid in a higher pressure region will tend to flow towards a lower pressure region, while fluid in a lower surface tension region will tend to flow towards a higher surface tension region (the Marangoni effect). Furthermore, any evaporation of fluid will accordingly result in a decrease in local drop height. For a derivation of Eq. (1) and the velocity field ($u,w$), see Section S1 of the Supplementary material.

The pressure in the droplet is dominated by surface tension effects. Therefore, the pressure p can be given by the Laplace pressure:(2)$p=-\frac{1}{r}\frac{\partial }{\partial r}\left(\sigma \frac{r{\partial }_{r}h}{\sqrt{1+{({\partial }_{r}h)}^2}}\right)\text{.}$Substituting the Laplace pressure in Eq. (1) shows that the droplet will tend toward a spherical cap shape. If there is no evaporation ($w_e=0$) a droplet will tend towards an equilibrium, constant curvature shape ($\frac{\partial p}{\partial r}=\frac{\partial \sigma }{\partial r}=0$). Note that $\sigma$ is also part of the derivative in Eq. (2), because it is not necessarily constant in the presence of surfactants. Thus, rather than writing $\sigma$ outside the derivative as is usually done (e.g. see [12]), it should be included in it as shown by Thiele et al. [56]. Note that Eq. (2) is beyond the lubrication limit, because of the square root in the denominator. In traditional lubrication theory, this term is approximated as unity.

The boundary and initial conditions that $h(r,t)$ is subjected to are given by:(3)${\left(\frac{\partial h}{\partial r}\right)}_{r=0}=0,$(4)${\left(\frac{{\partial }^{3}h}{\partial r^3}\right)}_{r=0}=0,$(5)$h(R,t)=0,$(6)$h(r,0)=h_0(r)\text{.}$Here, R is the drop radius and $h_0$ is the initial drop profile, given by the previously mentioned spherical cap shape. The considered cases involve a contact line that is pinned rather than one that moves [57], [58]. This means the drop radius will remain constant during evaporation, while the contact angle decreases, as opposed to a decreasing radius with a constant contact angle. Contact line pinning typically occurs for relatively small contact angles and rough substrates and is also promoted by the presence of surfactants [8], [38], [41], [43].

2.2. Interfacial surfactant transport equation

At the liquid–air interface of the droplet a surfactant concentration $\mathrm{\Gamma }(r,t)$ can be defined, which can be used to describe the transport of adsorbed surfactant molecules along the interface. As shown by Wong et al. [59], the surfactant transport equation is given by:(7)$\begin{array}{cc}\hfill \frac{\partial \mathrm{\Gamma }}{\partial t}=-& \frac{1}{r}\frac{\partial ({\mathit{rU}}_t\mathrm{\Gamma })}{\partial s}+\frac{\mathrm{\Gamma }{\partial }_{t}h}{1+{({\partial }_{r}h)}^2}\left(\frac{{({\partial }_{r}h)}^2}{1+{({\partial }_{r}h)}^2}+\frac{1}{r}\frac{\partial h}{\partial r}\right)+\hfill \\ \hfill & \frac{D_{\mathrm{\Gamma }}}{r}\frac{\partial (r{\partial }_s\mathrm{\Gamma })}{\partial s}+\frac{\partial h}{\partial t}\frac{{\partial }_{r}h}{1+{({\partial }_{r}h)}^2}\frac{\partial \mathrm{\Gamma }}{\partial r}+J_{\mathrm{\Gamma }\varphi }\text{.}\hfill \end{array}$Here, $U_t$ is the fluid velocity tangential to the liquid–air interface, $D_{\mathrm{\Gamma }}$ the surface diffusion coefficient and $J_{\mathrm{\Gamma }\varphi }$ is the rate with which surfactant is exchanged with the bulk through adsorption and desorption. Note that $J_{\mathrm{\Gamma }\varphi }$ can be either positive and negative. The derivative $\frac{\partial }{\partial s}$ is the surface derivative which can be written as $\frac{\partial }{\partial s}=\frac{1}{\sqrt{1+{({\partial }_{r}h)}^2}}\frac{\partial }{\partial r}$.

The first term on the right-hand side of Eq. (7) accounts for convective transport tangential to the interface and the second term for transport normal to the interface, which in practice boils down to a source- or sink-like effect as the interfacial curvature decreases or increases respectively. The third term denotes diffusion, the fourth term corrects for the displacement of the surface coordinates that move along the normal of the surface and the last term is the adsorption/desorption of surfactant from the bulk onto the interface and vice versa. For a derivation of Eq. (7), see Section S.2 in the Supplementary material.

The surfactant concentration $\mathrm{\Gamma }(r,t)$ is subjected to the following boundary conditions and initial condition:(8)${\left(\frac{\partial \mathrm{\Gamma }}{\partial r}\right)}_{r=0}=0,$(9)${\left(\frac{\partial \mathrm{\Gamma }}{\partial r}\right)}_{r=R}=0,$(10)$\mathrm{\Gamma }(r,0)={\mathrm{\Gamma }}_0\text{.}$These denote the symmetry condition, no-flux condition and initial (homogeneous) surfactant concentration ${\mathrm{\Gamma }}_0$ respectively.

Since surfactants adsorbed at the interface decrease the surface tension, an equation of state for $\sigma$ is required to close the problem. In this work, the Frumkin equation is used, which is analogous to the Langmuir isotherm and considers that surfactant molecules occupy a finite amount of space at the interface [60]. The Frumkin equation is given by:(11)$\sigma ={\sigma }_0+R_{u}T{\mathrm{\Gamma }}_{\infty }\ln \left(1-\frac{\mathrm{\Gamma }}{{\mathrm{\Gamma }}_{\infty }}\right)\text{.}$Here, ${\sigma }_0$ is the surface tension of the pure liquid, $R_u$ is the universal gas constant, T the temperature (which is assumed constant) and ${\mathrm{\Gamma }}_{\infty }$ the maximum possible surfactant concentration. Note that for $\mathrm{\Gamma }\ll {\mathrm{\Gamma }}_{\infty }$ the equation can be approximated by a linear, dilute equation of state.

Given this equation of state for surface tension and the radial velocity (as defined by Equation (S.7) in the Supplementary material) a typical Marangoni velocity can be defined as $u_{\mathit{Mar}}=\frac{{\mathit{HR}}_{u}T{\mathrm{\Gamma }}_0}{\mu R}$. This typical velocity is used in SubSection 2.6 to define several relevant dimensionless numbers.

2.3. Bulk surfactant transport equation

Besides surfactant being adsorbed onto the interface, there is also surfactant dissolved in the bulk of the droplet. The bulk monomer concentration $\varphi$ is considered to be a function of r and t and independent of z, which is allowed as rapid vertical diffusion is assumed (as outlined by [35], [61]). In order to let the concentration be independent of $h(r,t)$, the evolution of the bulk monomer concentration is described as a function of $\psi (r,t)=\varphi (r,t)h(r,t)$, as introduced by Thiele et al. [62]. The transport equation can then be written as:(12)$\frac{\partial \psi }{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(\frac{{\mathit{rh}}^2\psi }{3\mu }\frac{\partial p}{\partial r}-\frac{\mathit{rh}\psi }{2\mu }\frac{\partial \sigma }{\partial r}+D_{\varphi }h\frac{\partial \varphi }{\partial r}\right)-J_{\mathrm{\Gamma }\varphi }\sqrt{1+{({\partial }_{r}h)}^2}-J_{M\varphi }N\text{.}$In this equation, three terms can be distinguished. The first term is a transport term (with $D_{\varphi }$ the diffusion coefficient of surfactant monomers in the bulk), which takes into account convective and diffusive transport, the second term is the adsorption/desorption of surfactant from the bulk onto the interface and vice versa, including a factor that compensates for the interface geometry, and the third term is the micelle formation rate $J_{M\varphi }$ multiplied with the preferred number of surfactant monomers N which form a single micelle.

When the bulk surfactant concentration exceeds a certain threshold, it becomes energetically more favorable for the molecules to cluster together and form micelles. This threshold is called the critical micelle concentration (CMC). In practice, the number of surfactant monomers that form a single micelle tends to strongly prefer a single value N, which is also assumed in this work [63].

Similarly to the monomer bulk concentration $\varphi$, the micelle bulk concentration M is given in terms of $\zeta (r,t)=M(r,t)h(r,t)$:(13)$\frac{\partial \zeta }{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(\frac{{\mathit{rh}}^2\zeta }{3\mu }\frac{\partial p}{\partial r}-\frac{\mathit{rh}\zeta }{2\mu }\frac{\partial \sigma }{\partial r}+D_{M}h\frac{\partial M}{\partial r}\right)+J_{M\varphi }\text{.}$Here, $D_M$ is the diffusion coefficient of micelles. Note that it is assumed here that micelles do not adsorb directly onto the interface. This is a common assumption in literature [30], [33], [35], that follows from the fact that surfactant monomers need to dissociate from micelles before they can be adsorbed onto the interface [64]. Since the micelle formation/decomposition process is modelled as a single step, it is consistent to only allow individual monomers to adsorb onto the interface. Future models, however, could include multi-step models.

The bulk concentrations $\psi (r,t)$ and $\zeta (r,t)$ are subject to the following boundary conditions and initial conditions:(14)${\left(\frac{\partial \psi }{\partial r}\right)}_{r=0}={\left(\frac{\partial \zeta }{\partial r}\right)}_{r=0}=0,$(15)${\left(\frac{\partial \psi }{\partial r}\right)}_{r=R}={\left(\frac{\partial \zeta }{\partial r}\right)}_{r=R}=0,$(16)$\psi (r,0)={\varphi }_{0}h(r,t),$(17)$\zeta (r,0)=M_{0}h(r,t)\text{.}$Similarly to the boundary and initial conditions of $\mathrm{\Gamma }$, these denote the symmetry condition, the no-flux condition at the contact line and the initial, constant bulk concentrations ${\varphi }_0$ and $M_0$ respectively. All initial surfactant concentrations (${\mathrm{\Gamma }}_0,{\varphi }_0$ and $M_0$) are always chosen to be in equilibrium, so initially $J_{\mathrm{\Gamma }\varphi }=J_{M\varphi }=0$. A derivation of both Eqs. (12), (13) is given in Section S.3 of the Supplementary material.

2.4. Surfactant adsorption

The transport between the interface and the bulk is a continuous adsorption and desorption of molecules that overall tends towards a dynamic equilibrium. Furthermore, there is a limited amount of space available at the interface for surfactants, so as the interfacial concentration increases the adsorption rate tends to decrease, while the desorption rate increases. This behavior can be described by the following reaction equation [33], [35]:(18)$S+\varphi \underset{k_a^{\mathrm{\Gamma }}}{\overset{k_d^{\mathrm{\Gamma }}}{\rightleftharpoons }}\mathrm{\Gamma }\text{.}$Here, S ($=1-\frac{\mathrm{\Gamma }}{{\mathrm{\Gamma }}_{\infty }}$) indicates the fraction of available space at the interface, $k_a^{\mathrm{\Gamma }}$ is the interfacial adsorption coefficient, and $k_d^{\mathrm{\Gamma }}$ the interfacial desorption coefficient. Thus, the interfacial adsorption flux is given by:(19)$J_{\mathrm{\Gamma }\varphi }=k_a^{\mathrm{\Gamma }}\varphi \left(1-\frac{\mathrm{\Gamma }}{{\mathrm{\Gamma }}_{\infty }}\right)-k_d^{\mathrm{\Gamma }}\mathrm{\Gamma }\text{.}$In this equation it can indeed be recognized that the first, adsorption term increases as the bulk concentration increases and approaches zero as $\mathrm{\Gamma }\to {\mathrm{\Gamma }}_{\infty }$, while the second, desorption term becomes more negative as $\mathrm{\Gamma }$ increases. This behavior corresponds indeed to Eq. (18).

Note that surfactant adsorption onto the substrate is not taken under consideration here, because this will not have a direct effect on the internal flow. Marangoni flow can only occur as a result of a surface tension gradient at a free interface. Any indirect influences of substrate sorption on the flow – surfactant being removed or added to the bulk – are considered less significant.

Similar to Eq. (18), a reaction equation can be written for the formation of micelles:(20)$N\varphi \underset{k_a^M}{\overset{k_d^M}{\rightleftharpoons }}M\text{.}$Here, $k_a^M$ is the micelle formation coefficient and $k_d^M$ the micelle decomposition coefficient. This reaction equation can subsequently be written into a micelle formation rate $J_{M\varphi }$ as well [30], [33], [35]:(21)$J_{M\varphi }=k_a^M{\varphi }^N-k_d^{M}M\text{.}$Given this equation, it is possible to approximate the reaction constants in terms of the CMC and concentrations. It can be seen that the equilibrium concentration is given by $\frac{k_a^M}{k_d^M}=\frac{M}{{\varphi }^N}$. Now, if $\mathrm{\Phi }(=\mathit{NM}+\varphi )$ is the total concentration of surfactant monomers in the bulk, substitution of the initial total concentration ${\mathrm{\Phi }}_0$ results in:(22)$\frac{k_a^M}{k_d^M}=\frac{{\mathrm{\Phi }}_0-\mathit{CMC}}{N{(\mathit{CMC})}^N}\text{.}$Here, it should be noted that in equilibrium $\varphi =\mathit{CMC}$, given that there are micelles present. Of course, it is still required to choose one of the reaction constants to define the time scales of the reactions.

2.5. Evaporation

The evaporative volume flux $w_e$ is given by:(23)$w_e=\frac{D_v{M}_l{p}_{\mathit{sat},l}}{\rho R_{u}T}\frac{\partial {\widehat{p}}_l}{\partial n}\sqrt{1+{\left(\frac{\partial h}{\partial r}\right)}^2}\text{.}$Here, $D_v$ is the vapor diffusivity in air, $M_l$ the molar mass of the liquid, $p_{\mathit{sat},l}$ the liquid saturation pressure, which is assumed constant, and ${\widehat{p}}_l=\frac{p_l}{p_{\mathit{sat},l}}$, with $p_l$ the local vapor pressure. The derivative $\frac{\partial }{\partial n}$ is the spatial derivative along the normal vector (pointing outwards) [21].

In order to calculate $\frac{\partial {\widehat{p}}_l}{\partial n}$, it is required to describe the vapor field. As shown by Deegan [11] and Hu and Larson [65], vapor diffusion can be considered instantaneous. Together with the assumption of no convection, this implies that the vapor field can be described by the Laplace equation:(24)${\nabla }^2{\widehat{p}}_l=0\text{.}$The corresponding boundary conditions are given by:(25)${\widehat{p}}_l{|}_{z=h}=1\mathrm{for}r<R\text{;}$(26)$\frac{\partial {\widehat{p}}_l}{\partial z}{|}_{z=0}=0\mathrm{for}r>R\text{;}$(27)${\widehat{p}}_l=\mathit{RH}\mathrm{for}(r,z)\to \infty \text{.}$Here RH is the relative humidity. These equations represent saturated vapor at the droplet surface, no vapor penetration at the substrate and ambient relative humidity far away from the droplet. An analytical solution to Eq. (24) is derived in Section S.4 of the Supplementary material.

2.6. Dimensionless parameters

One may note that the traditional capillarity number $\mathit{Ca}=\mu w_e/{\sigma }_0$ is not listed in Table 1. This dimensionless number is omitted, because for any physically relevant case $\mathit{Ca}\ll 1$. Varying Ca therefore has negligible effect on the drop evolution and internal flow field, because surface tension forces always dominate over viscous drag forces, resulting in little deviations from the spherical cap shape. Simulations show that changing Ca by a factor ${10}^3$ only results in a 0.0067% change in total pressure gradient from the drop center toward the contact line.

3.1. Below CMC without diffusion

Simulations are carried out for various ranges of the dimensionless parameters. The numerical procedure is outlined in Section S.5 of the Supplementary material. The initial contact angle ${\theta }_0$ is always set to ${25}^{\circ }$ and the contact line is pinned. During each simulation, at least 200 time steps are calculated ($~$5% of the drying time), after which the velocity field is analyzed.

Only the initial stage of the evaporation process is considered, because this allows one to capture a ‘snapshot’ of the flow dynamics. It is still required to allow the internal flow dynamics to evolve sufficiently, ensuring that the initial state has no significant effect on the flow field. If the full evaporation process would be taken into account, this would give a skewed perspective on the flow dynamics, because the considered dimensionless numbers change over time. For example, the height decreases over time and the interfacial surfactant concentration increases over time. Also, the bulk concentration may exceed the CMC at some point in time. Analysis of the entire evaporation process is outside the scope of this paper and can be considered in future work. Nevertheless, the results of this work can still be used to predict flow patterns during the entire evaporation process, as long as the relevant dimensionless numbers can be estimated.

Visual analysis of the velocity field has led to the definition of three separate regimes: the ‘coffee-ring regime’, where the flow field looks similar to pure droplet evaporation, the ‘Marangoni regime’, where a clear Marangoni eddy can be distinguished, and the ‘transition regime’, which shows behavior of both the coffee-ring and Marangoni regime. Representative velocity plots of the three regimes are shown in Fig. 1. In all three velocity regimes an outward, capillary flow can be distinguished, that is caused by selective evaporation at the contact line. However, for the Marangoni regime (and in some degree for the transition regime) there is a backflow close to the interface that convects adsorbed surfactant towards the center, where it desorbs into the bulk again. All three regimes have a zero-fluid velocity at both the liquid–air and liquid–solid interface at r = 0, which is in accordance with the ‘Hairy ball theorem’ (Poincaré-Brouwer theorem). This theorem predicts that there necessarily exists at least one zero-velocity point on the surfaces of compressible liquids and interfaces allowing mass, energy and momentum transport [66].

1. Download : Download high-res image (524KB)
2. Download : Download full-size image

Fig. 1. Flow regime classifications: (a) Marangoni regime, (b) Transition regime and (c) Coffee-ring regime. The red line indicates the position of the interface, given by $h(r)$.

Calculating the mean, dimensionless vorticity $\omega$ over the middle area where the Marangoni vortex tends to appear (around $R/2<r<5/6R$) shows that $\omega ⩽-1.2·{10}^{-4}$ corresponds to the Marangoni regime, $-1.2·{10}^{-4}<\omega ⩽-0.8·{10}^{-4}$ to the transition regime and $\omega >-0.8·{10}^{-4}$ to the coffee-ring regime. Here, the time scale $t_d=R^2/D_v$ is used to make the vorticity dimensionless.

Alternatively, as a cross-validation the resulting velocity fields can be classified by calculating the fraction of velocity vectors opposing the typical coffee-ring flow. Here, numerical analysis shows that the Marangoni regime corresponds to more than 24% of the radial velocity vectors and more than 9% of the axial velocity vectors pointing to the center and upward respectively. Similar, the transition regime corresponds to velocity fields that are not in the Marangoni regime with more than 22% of the radial velocity vectors and more than 6% of the axial velocity vectors pointing to the center and upward respectively. If a lower proportion of the velocity vectors points centerwise or upwards, the velocity field is considered to be in the coffee-ring regime. This gives similar results to the classification with vorticity and small variations of the values do not yield significantly different classification results. Further classification of regimes in this paper is made through the mean vorticity with frequent, random visual checks for extra verification.

Fig. 2 shows the regime charts for simulations without diffusion and without micelles. As can be seen in Fig. 2a, for $\mathrm{Tr}<1$, both $1/\mathrm{Ev}$ and $\mathrm{Tr}$ have very similar effects on the occurence of Marangoni flow. This makes sense, because as $1/\mathrm{Ev}$ decreases, the evaporative velocity starts to dominate over the Marangoni velocity. The evaporative effects are thus too strong for the Marangoni effect to counter, resulting in a coffee-ring regime flow. Similarly, as $\mathrm{Tr}$ decreases, the effects of adsorption and desorption become less significant (given that De remains constant). This will result in behavior as if the surfactant is insoluble, which is described by the coffee-ring regime [38]. Because the relative strength of the Marangoni effect increases, the interfacial concentration is kept homogeneous since the fluid velocity close to the interface is reduced and the concentration increase at the drop apex due to curvature effects becomes more dominant. At the same time less adsorption is taking place, so the increased bulk concentration close to the contact line will not result in enough interfacial adsorption to cause a positive concentration gradient. That increasing $\mathrm{Tr}$ tends to promote Marangoni circulation is also found by Jung et al., who used a lattice-gas model [39].

1. Download : Download high-res image (432KB)
2. Download : Download full-size image

Fig. 2. Marangoni circulation regimes without diffusion or micelles: (a) Evaporation number vs. Transport number ($\mathrm{De}=3.08\times {10}^3$), (b) Desorption number vs. Transport number ($1/\mathrm{Ev}=2.82\times {10}^6$) and (c) Desorption number vs. Evaporation number ($\mathit{Tr}=5.0\times {10}^{-4}$).

For $\mathrm{Tr}>1$, the occurrence of Marangoni circulation becomes mostly limited by $1/\mathrm{Ev}$. Only if a certain threshold value of $1/\mathrm{Ev}$ is exceeded, Marangoni vortices can be formed. At the same time, $\mathrm{Tr}$ is largely irrelevant, because the desorption kinetics are no longer dominated by the Marangoni velocity and thus are not a limiting factor anymore. That $1/\mathrm{Ev}$ still is a limiting factor however, is not surprising. Decreasing surfactant concentration ${\mathrm{\Gamma }}_0$ or increasing the evaporation rate (e.g. through an increase in $D_v$) will eventually always result in a coffee-ring flow, because either the droplet can be considered pure or the evaporation becomes too dominant.

It is important to note that a correct interpretation of $\mathrm{Tr}$ involves not only the desorption rate, but also the adsorption rate. If $\mathrm{Tr}$ is varied by changing $k_d$ the value of $k_a$ changes accordingly, given that $\mathrm{De}$ should remain constant. Increasing $\mathrm{Tr}$ thus does not only mean that interfacial desorption of surfactants becomes more dominant with respect to the Marangoni velocity, but interfacial sorption in general becomes more dominant.