Fluid flow structures in an evaporating sessile droplet depending on the droplet size and properties of liquid and substrate

We investigate numerically quasi-steady internal flows in an axially symmetrical evaporating sessile droplet depending on the ratio of substrate to fluid thermal conductivities, fluid volatility, contact angle and droplet size. Temperature distributions and vortex structures are obtained for droplets of 1-hexanol, 1-butanol and ethanol. To this purpose, the hydrodynamics of an evaporating sessile drop, effects of the thermal conduction in the droplet and substrate and diffusion of vapor in air have been jointly taken into account. The equations have been solved by finite element method using ANSYS Fluent. The phase diagrams demonstrating the number and orientation of the vortices as functions of the contact angle and the ratio of substrate to fluid thermal conductivities, are obtained and analyzed for various values of parameters. In particular, influence of gravity on the droplet shape and the effect of droplet size have been considered. We have found that the phase diagrams of highly volatile droplets do not contain a subregion corresponding to a reversed single vortex, and their single-vortex subregion becomes more complex. The phase diagrams for droplets of larger size do not contian subregions corresponding to a regular single vortex and to three vortices. We demonstrate how the single-vortex subregion disappears with a gradual increase of the droplet size.

on a heat transfer in the immediate vicinity of the symmetry axis piercing the apex. The transition between the opposite circulation directions, taking place with a variation of the relative substrateliquid thermal conductivity, has been identified in [13] under the conditions that differ from those obtained in Ref. [12]. Hu and Larson demonstrated that fluid circulation in the vortex can reverse its sign at a critical contact angle for a drop placed onto substrates with finite thicknesses [14].
The influence of substrate temperature has been studied in a number of experiments [15][16][17].
The thermal conduction processes throughout the droplet can generally result in a nonmonotonic spatial dependence of the surface temperature and in more complicated convection patterns inside a drop. The temperature distribution can be qualitatively understood as a result of matching the heat transfer through the solid-liquid and liquid-vapor interfaces. With varying relative substrate-liquid thermal conductivity and/or contact angle, transitions between regimes with different numbers of vortices and/or circulation directions take place [18,19].
In order to form the vortex structure inside the droplet in a controlled manner, one should know the dependence of the flow patterns on such externally varying parameters as thermal conductivities of liquid and substrate, substrate thickness, liquid volatility, contact angle and droplet size. Here we study such a dependence by performing the numerical simulation of Marangoni convection in a drop. The obtained results extend those of Ref. [19], which mainly focus on the substrate thickness but did not consider the dependence of fluid flows on the liquid volatility and droplet size.

II. BASIC EQUATIONS
We consider axially symmetrical evaporating sessile droplet with a contact line pinned to a hydrophilic solid substrate ( Figure 1). The droplet shape can be described with good accuracy within a spherical cap approximation when the Bond number Bo = ρgh 0 R/(2σ sin θ) and the capillary number Ca = ηv/σ are much smaller than unity (see the notations in Table I).
The axial symmetry of thermocapillary fluid flow can break down for a strong evaporation of a substantially heated droplet and a resulting strong fluid velocity [20][21][22]. However, axisymmetric flows can usually be observed at room temperature and under atmospheric pressure [12].
The dynamics of vapor concentration u(r, z) in the surrounding atmosphere is described by the diffusion equation The boundary conditions are: u = u s on the drop surface, u = 0 far away from the drop, ∂u/∂r = 0 and ∂u/∂z = 0 on the axes r = 0 and z = 0, correspondingly.
In the quasi-stationary approximation, the diffusion equation (1) can be replaced by the Laplace's equation ∆u = 0. The approximation can be employed if the vapor concentration adjusts on time scales much less than the total drying time (i.e., R 2 /D is much smaller than t f ≈ 0.2ρRh 0 /(Du s )), and if the Stefan flow can be neglected. The latter condition holds at room temperature and under atmospheric pressure, when the saturated vapor concentration is much smaller than the density of air [23]. Within the quasi-stationary approximation and the spherical cap approximation, the analytical solution of the diffusion equation can be obtained which gives the inhomogeneous evaporation flux from the surface of the evaporating droplet [24] J s (r) = Du s R where x(r) = r 2 cos θ/R 2 + 1 − r 2 sin 2 θ/R 2 and P −1/2+iτ (x) is the Legendre polynomial. Expr.
(2) can be conveniently approximated with high accuracy as [24] where J 0 (θ) can be determined by the following expressions [25]: Λ(θ) = 0.2239(θ − π/4) 2 + 0.3619, J 0 (π/2) = Du s /R. parameters Saturated vapor density u s g/cm 3 1.46 · 10 −4 2.76 · 10 −5 6.55 · 10 −6 The above model of diffusion-limited evaporation assuming a constant vapor pressure along the interface is valid for a relatively small temperature gradients and/or only a minor slope in the vapor pressure chart. For a particular cases considered in the present study, the temperature difference along the interface is found to be of the order of 1 K, and there is ∼ 5% relative deviation of saturated vapor density as a result of the 1 K change in temperature. Therefore, the above model can be employed to a first approximation.
It follows from (4) that the evaporating flux density is inhomogeneous along the surface and increases substantially on approach to the pinned contact line. A resulting inhomogeneous mass flow and corresponding heat transfer modify the temperature distribution along the droplet surface and, therefore, result in appearence of Marangoni forces associated with the temperature-dependent surface tension, which, in turn, generate convection inside the droplet.
The basic hydrodynamic equations inside the drop are the Navier-Stokes equations and the continuity equation for the incompressible fluid div v = 0.
The buoyancy force has been neglected in Eq. (8), because the buoyancy-induced convection is much weaker than Marangoni flow when ρgh 2 β/(7σ ) 1, where β is thermal expansion coefficient [27]. The latter condition holds for all the cases considered in the present study, including droplets of larger size in Sec. IV, where we always have h < 0.1 cm.
The calculation of temperature distribution is carried out using the thermal conduction equation For relatively small and slowly evaporating droplets, when P e = vR/κ 1, the convective heat transfer is much smaller than conductive heat transfer, and, hence, the velocity field does not influence the thermal conduction: This is the case, for example, for the droplet of 1-hexanol (see the parameters in Table I). The thermal conduction inside the substrate is also described by (11). Considering a fixed geometric shape, we also disregard the capillary flow, which is the coffeestain outward flow [24]. This is justified in presence of a strong recirculating flow due to the Marangoni effect. Indeed, the velocity of the capillary flow is of the order of R/t f when the contact angle is not too small (see, e.g., the expressions in [28][29][30]  The simulation includes several steps.

• Geometry and mesh generation
We use two-dimensional geometry of the axially symmetrical surface according to Table I.
The surface curvature radius and the droplet height are equal to R/ sin θ and R(1/ sin θ − cot θ), correspondingly. In ANSYS Fluent, the symmetry axis should be denoted as x-axis (as distinct from z-axis in Fig. 1). The next step is to divide the simulation region into small computational cells.

• Defining model and setting properties
We choose the pressure-based solver type, 2D axisymmetric geometry, steady mode for time, viscous (laminar) flow and enable the calculation of energy in the model. We specify all physical properties of the fluid such as viscosity and density. The values are listed in Table I including the temperature derivative of the surface tension σ T = −∂σ/∂T which allows us to specify the properties of Marangoni stress.

• Setting boundary conditions and solution method
The boundary conditions are set according to Section II. Some of the boundary conditions, such as at the symmetry axis, are automatically taken into account in ANSYS Fluent. The rate of heat loss Q 0 = LJ s (r) is specified with user-defined function (UDF) written in the C programming language. The solution algorithm is SIMPLE (Semi-Implicit Method for Pressure Linked Equations). It uses a relationship between velocity and pressure corrections to enforce mass conservation and to obtain the pressure field. The algorithm is written in such a way that the continuity equation is automatically satisfied. See [31] for detailed description of the SIMPLE algorithm.
• Post-processing We display the simulation results: velocity field, absolute values of velocity and temperature distribution, and analyse the vortex structure.
We note that our mesh contained at least 70000 cells for each contact angle. We used at least We have determined the phase diagram regions via the numerically obtained fluid flow fields.
The surface temperature distrubition usually follows the vortex structure in the droplet, although there are some exceptions to this rule due to the inertia of the fluid flow.
In our phase diagrams, the blue curve and black curve correspond to the change of sign of the tangential component of temperature gradient at the liquid-vapor interface, at the droplet apex and near the contact line, correspondingly. As shown in Ref. [19], with an increase of the substrate thickness, the subregion II 2 becomes dominating over II 3 . Further increase of the substrate thickness results in shifting the regions I and II to larger contact angles. The blue curve is situated below the black curve at small values of h R = h S /R, while for h R 0.05 the blue curve is above the black curve [19]. The phase diagrams in Fig. 2 correspond to the latter case. The main features of the phase diagram can be qualitatively understood as a result of matching the heat transfer across the solid-liquid interface and the vaporization heat through the liquidvapor interface (see details in [18,19]). As a result of the matching and competition of these two effects, vortex structures are formed which depend on the values of k R and θ.
Consider now the volatility effect on the phase diagram. Fig. 2 shows the phase diagrams obtained for the three droplets of different liquids, where the parameters of the droplets are listed in Table I. Since the saturated vapor density of 1-hexanol, 1-butanol and ethanol differ by orders of magnitude as seen in Table I, the evaporation rates (volatilities) of these liquids described by Eqs.(4)-(7) also significantly differ.
As seen in Fig. 2, with an increase of the fluid volatility, the subregion III disappears; the subregion II 2 increases at the expense of the subregion III; the subregion II 3 becomes larger at the expense of the subregion I; the subregion I diminishes for small values of k R and becomes more complex for highly volatile liquids.
Domination of the evaporative cooling as compared to the heat flow at the substrate-liquid θ R(mm) R s (mm) h s (mm) Bond number k R = 0.2 k R = 0.5 k R = 1 k R = 2 k R = 5 k R = 7 interface leads to the reversed single-vortex regime for small values of k R and at small contact angles, for weakly volatile liquids [19]. For highly volatile liquids, however, the cooling of the droplet bulk is a faster process resulting in a larger temperature difference between the droplet bulk and the substrate, which, in turn, leads to increasing the heat flow at the substrate-liquid interface close to the contact line even for moderately low values of k R . This is in accordance with the disappearance of the subregion III for highly volatile liquids as observed in Fig. 2.
Increasing the fluid volatility results in larger chances of domination of the vaporization heat for relatively large contact angles. This is in accordance with the shift of the subregion I at relatively small values of k R to larger contact angles. Understanding of the complex structure of subregion I shown in Fig. 2(c) for relatively large values of k R requires an additional study.
Consider now the effect of droplet size on the phase diagram. The fluid flow structure has been calculated for droplet size exceeding that in Table I Nonspherical droplet free surface was obtained with the Runge-Kutta numerical method using the approach described in [26].
We find that the phase diagrams for droplets of larger size do not contain the subregions I and II 3 . Table II shows numbers of vortices in the droplet of 1-butanol of different size for θ ≈ 30 • .
It demonstrates that droplets of capillary size contain a single-vortex for θ ≈ 30 • , while larger droplets contain two vortices in this case. Although the existence of subregion I depends on k R , the subregion I disappears with gradually increasing the droplet radius R for each k R , as shown in Table II.
The absence of subregion I for large droplets can be understood as follows. In order for the subregion I to occur, the heat transfer through the solid-liquid interface should dominate compared to the vaporization heat. For droplets of larger size, the former heat flow becomes relatively weak for the regions far away from the substrate. Therefore, the vaporization heat dominates there, and the vortex splits into two vortices.

V. CONCLUSION
We have performed numerical simulations of the fluid flow structure in evaporating sessile droplets depending on thermal conductivities of liquid and substrate, liquid volatility, contact angle and droplet size, and we have analyzed the corresponding phase diagrams in the k R -θ plane.
We find that with an increase of the fluid volatility, the subregion I becomes more complex and the subregion III disappears; the subregion II 2 increases at the expense of the subregion III; the subregion II 3 becomes larger at the expense of the subregion I. Also, we have shown that the phase diagrams for droplets of larger size do not contain the subregions I and II 3 .
The results may provide a better understanding of the Marangoni effect of drying droplets and may be useful for preparing the fluid flow vortex structure in a controlled manner by selecting substrates and liquids with appropriate properties. Also, the results may potentially be useful to better understand evaporative deposition patterns and evaporative self-assembly of nanostructures.