Modeling unsteady B\'enard-Marangoni instabilities in drying volatile droplets on a heated substrate

We study unsteady internal flows in a sessile droplet of capillary size evaporating in constant contact line mode on a heated substrate. Three-dimensional simulations of internal flows in evaporating droplets of ethanol and silicone oil have been carried out. For describing the Marangoni flows we find it necessary to account for the diffusion of vapor in air, the thermal conduction in all three phases and thermal radiation. The equations have been solved numerically by finite element method using ANSYS Fluent. As a result of the simulations, the nonstationary behavior of B\'enard-Marangoni (BM) instabilities is obtained. At the first stage, a flower structure of BM cells near the triple line emerge. For smaller contact angles, the cells grow in size and occupy the central region of the droplet surface. Being closely connected with recent experimental and theoretical studies, the results obtained help to analyze and resolve the associated issues.


I. INTRODUCTION
The issues of droplet evaporation have attracted significant attention both from the theory point of view and in connection with the development of new applications, including pattern formation on surfaces, preparation of ultra-clean surfaces, deposition of microarrays of DNA and RNA molecules, protein crystallography, combustion of fuel droplets in engine, diagnosis of diseases, development of modern printing methods for inkjet printers, manufacture of new electronic and optical devices and a number of other areas [1][2][3][4][5][6].
The nonuniform mass flow and heat transfer during a drop evaporation have a crucial influence on the evaporation kinetics. In particular, they change the temperature distribution along the droplet surface and can cause thermocapillary convection inside the droplet due to the temperaturedependent surface tension. Axisymmetric Marangoni flows can often be observed at room temperature and under atmospheric pressure [7][8][9][10].
A thermocapillary fluid flow with the broken axial symmetry can spontaneously appear inside a substantially heated droplet with an intense evaporation and sufficiently high fluid velocity. In a heated droplet, a Marangoni instability due to non-uniform temperature distribution can lead to arXiv:2008.10287v1 [physics.flu-dyn] 24 Aug 2020 the fluid motion of the particular form of longitudinal rolls or hydrothermal waves (HTWs), as was first demonstrated in [11] and confirmed in numerous consequent studies [12][13][14][15][16][17][18][19]. HTWs are timeperiodic traveling waves in the liquid, induced by instabilities of the steady Marangoni convection, which were originally predicted for the liquid layer on a rigid plate by Smith and Davis [20,21] and later observed in [22,23]. According to Smith and Davis [20,21], the disturbance propagates perpendicular to the surface temperature gradient at small Prandtl numbers and almost exactly in the upstream direction at large Prandtl numbers [20,21].
The HTWs in droplets were originally identified by Sefiane et al. in 2008 [11], when the evaporation of sessile droplets of ethanol, methanol, water, and Fluorinert FC-72 was investigated.
Later more detailed observations of HTWs have been reported for the ethanol droplets on four different substrates with different thermal conductivities [12]. HTWs moving in azimuthal direction have been observed for volatile droplets of ethanol and methanol. For FC-72, convective cells emerged near the droplet apex and moved towards the contact line. The number of waves has been found to increase with the increase of substrate thermal conductivity and its temperature, and to decrease with the decrease of the droplet height. On the other hand, no hydrothermal instability has been found in evaporating water droplets.
Three stages of the evaporation process with thermoconvective instabilities have been distinguished in [13] in the droplets of methanol, ethanol, and FC-72 on a polytetrafluoroethylene substrate: heating of the drop (which reaches almost the substrate temperature), evaporation with thermoconvective instabilities (the heat flux is maximal at this stage), and evaporation without thermoconvective patterns (the droplet is more like a film at this stage). The main phase is characterized by the presence of convection cells, and the origin of the cells' movement is associated with the temperature difference between the heated substrate and the ambient temperature. Temporal evolution of the HTWs was also described by means of simultaneous thermal and fluid motion observations in a drying ethanol droplet on a heated substrate [14]. HTWs traveling in azimuthal direction in a droplet of ethanol were demonstrated to take place both under normal gravity and microgravity conditions, so that the gravity effects on the behavior of HTWs are small [15].
Sefiane et al. [16] demonstrated that HTWs in FC-72 droplets are bulk waves that extend across the entire droplet volume. Also, the influence of the HTWs on the heat transfer in the solid substrate and on evaporation rate spatial distribution was studied in [16]. A specific kind of HTWs was determined by Karapetsas et al. in [17], where three-dimensional spiral-like timeperiodic traveling waves were observed that are organized radially and are uniform along the azimuthal angle. Shi et al. [24] determined critical Marangoni numbers for the emergence of the Marangoni convection instabilities in a sessile droplet with low volatility on a heated substrate by a series of three-dimensional computer simulations.
New features of the problem in question have been revealed by Semenov et al. [25], who have found that the cell patterns they determined in the droplet with a large Prandtl number represent the BM cell structure rather than HTWs. The authors have observed Marangoni instability in an ethanol droplet of capillary size on a heated and highly conductive substrate. They also performed numerical simulations by one-sided model described in [26], where they took into account, in particular, an unsteady diffusion-limited evaporation for non-isothermal pinned sessile droplet and The authors found that the number of BM cells increased with the decrease of the droplet height.
The interplay of the BM convection cells and HTWs has been observed in the recent work [27] in a volatile sessile droplet of silicone oil on a copper substrate at constant contact angle mode over a wide range of substrate temperatures. When the substrate temperature was higher than room temperature, the BM quasistationary convection cell structure was observed and found to evolve into an oscillating state, as the evaporation progresses. When the substrate was colder than the surrounding air, the oscillating BM cells were observed only when the temperature of the substrate is higher than 10.6 • C. The cell size became larger and their number was smaller when the temperature of the substrate goes up above the temperature of the surrounding gas. When the substrate temperature is above 14.3 • C, the coexistence of advancing hydrothermal waves and BM vibrational convection was observed. The paper reports critical values of the contact angle at which the transition between different types of convection occured.
In particular, the authors [27] considered the case when the substrate temperature is 7.  [25,27]. At the same time, the spatial structures of BM cells at sufficiently small contact angles, which were observed in [27] and numerically obtained in [25], differ substantially. The similarity of the results obtained in [25,27] is expected as the droplet parameters considered in the two papers are close to each other. Also the relevant physical properties of 0.65 cSt silicone oil and properties of the ethanol are similar, and highly conductive heated substrates at similar temperatures are employed in both works. It is therefore surprising that the results are quite different at a later stage of the droplet evaporation.
The purpose of this paper is to analyze the behavior of unsteady BM cells in volatile sessile droplets on heated substrates including the stage of comparatively small contact angles. It is also of interest to identify the main physical mechanisms, which are responsible for the behavior studied and have to be taken into account in the corresponding simulations.
We have theoretically investigated the Marangoni instability patterns in evaporating sessile droplets of ethanol and silicone oil for the particular parameters which are close to each of the experiments [25,27]. Based on the three-dimensional numerical simulations, the fluid flows and geometric characteristics of the BM patterns in the droplet have been described. The mathematical model was developed that takes account of heat transfer in all three phases, hydrodynamics, vapor diffusion in air, Stefan flow in the gas and thermal radiation. Using ANSYS Fluent software, threedimensional computer simulations were carried out and the unsteady temperature distributions and velocity fields in the drop were obtained. We have theoretically recovered both the flower structure of BM cells near triple line for sufficiently large contact angles and the BM cell pattern where larger BM cells occupy the fluid-gas interface for smaller contact angles.
The paper is organized as follows. Sec. II describes mathematical model and method of numerical solution. Relevance of various physical effects for the problems under consideration are also analyzed in Sec. II. In Sec. III we present our main results. Sec. IV contains discussion.

II. MATHEMATICAL MODEL
We consider three-dimensional nonstationary problem for an axially symmetrical sessile droplet on a heated solid substrate. The surface is taken as a spherical cap. This approximation is valid under the condition Bo 1, i.e., when the influence of gravitation on the droplet shape is small.
Here Bo = ρgh 0 R/(2σ sin θ) is the Bond number, and h 0 is the droplet height. We use the notations in Table I.

A. Diffusion of vapor in air
The dynamics of the vapor concentration u in the surrounding atmosphere is described by the diffusion equation The boundary conditions are as follows: 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, respectively.
Semenov et. al derive in [26] and employ in [25] analytical equations for the unsteady diffusionlimited evaporation from sessile droplets, which agree with the numerical results of [28]. However, these unsteady effects are small on time scales larger than the time required for the Brownian particle to pass the characteristic length R. i.e., for t R 2 /D ∼ 0.7 s., for the problems under consideration. Therefore, the evaporation kinetics can be considered as a steady state process and one can use the quasi-stationary approximation ∆u = 0 (2) with the same boundary conditions together with the spherical cap approximation, and then find the analytical solution for the evaporation flux [29,30]. It is demonstrated in [29] the the analitical solution can be approximated with high accuracy by the following relation where λ(θ) = 1/2 − θ/π. (3) can be obtained as follows [31]: Λ(θ) = 0.2239(θ − π/4) 2 + 0.3619, where u s is the saturated vapor density at the droplet free surface. Saturated vapor pressure is determined by the Antoine equation (with p s in (Pa) and T in ( • K)): where A, B, and C are empirical constants. We use whereR is the universal gas constant and M is molar mass of the vapor.
It follows from (3) that the mass loss resulting from vaporization occurs nonuniformly along the free surface of the liquid layer, significantly increasing near the pinned contact line. A nonuniform mass flow during evaporation and the corresponding heat transfer change the temperature distribution along the droplet surface and cause Marangoni forces due to temperature-dependent surface tension. These forces produce thermocapillary convection inside a droplet.

B. Stefan flow in the gas
Because the saturated vapor concentration is not much smaller than the density of air (see Table I), it is natural to consider the contribution of the Stefan flow in the gas for the problems under consideration [33].
As a result of the convective mass transport in the gas, the evaporation rate increases as [26]: where u g = M g p g /(RT ), p g is the air pressure, It follows that the Stefan flow in the gas increases the evaporation rate by approximately 10% and 18% for droplets of ethanol and 0.65 cSt silicone oil, respectively.
The results below in Sec. III do not account for the Stefan flow in the gas. However, our preliminary numerical calculations for the 0.65 cSt silicone oil droplet with parameters specified in Table I and contact angle 20.49 • show that the influence of the Stefan flow on the resulting surface temperature distribution is quite small.

C. Hydrodynamics
The governing equations for the fluid dynamics inside the drop are the Navier-Stokes equations and the continuity equation for the incompressible fluid ∂v ∂t Here ν = η/ρ is kinematic viscosity, v is fluid velocity, p is pressure. The boundary conditions for the equations (10-11) take the form: v = 0 at the substrate-fluid interface and ∂σ/∂τ = η(∂v τ /∂n − v τ ∂ϕ/∂τ ) at the fluid-gas interface [28].

D. Heat transfer in the drop
The temperature distribution is obtained with the heat conduction equation The heat conduction inside the substrate is obtained with The boundary conditions for the thermal conduction are as follows: T = T s at the lower bound of substrate, ∂T s /∂n = 0 at the substrate-gas interface, k s ∂T s /∂z = k∂T L /∂z at the substratefluid interface. Here T L and T S correspond to the liquid and substrate temperature, k and k S are thermal conductivities of liquid and substrate, respectively.
Following the experiments [25,27], the calculations in the present work involve a highly conducting substrate. In this case, the boundary conditions can be written simply as T = T s at the substrate-fluid interface.
The boundary condition at the fluid-gas interface is discussed in Sec. II E.

E. Heat flux at the fluid-gas interface
The boundary condition for the heat transfer at the droplet surface is Here, is the total heat flux from the droplet surface, which includes latent heat of evaporation, effects of the thermal conduction in liquid and gas and thermal radiation. The estimates below show that all these effects are important for droplets evaporating from heated substrates. Also, n is a normal vector to the droplet upper surface, k air ∂T /∂n air is the heat flow associated with non-uniform temperature distribution in an ambient gas, k air is the thermal conductivity of air, J (r) is the local evaporation rate determined by (4), is the net radiation flux between the drop and the environment; here σ B is the Stefan-Boltzman constant, ε is the emissivity; T and T a is the droplet surface temperature and temperature of the ambient gas, respectively.
The relation (16) can be employed under the following condition. First, it is important that the droplet size is much larger than the characteristic wavelengths of the thermal radiation, which are determined by the Planck distribution. Under this condition the emissivity almost does not depend on the frequency [34]. And vise versa, for objects with dimensions much smaller than the thermal wavelength, super-Planckian far-field radiative heat transfer is possible [35,36]. Secondly, it is important that the temperature difference inside the drop is much smaller than the temperature difference between the droplet and the ambient gas.
Let's obtain the relation for (∂T /∂n) air .
The quasi-stationary temperature distribution in the gas phase is determined by the equation with the boundary conditions: T = T s at the droplet surface (under the assumption that the temperature difference inside the drop is much smaller than the temperature difference between the droplet and the ambient gas), T = T a away from the droplet. This problem is mathematically equivalent to Eq. (2) with the corresponding boundary conditions, therefore, one can employ the same formulas for the analytical solution. It follows from (3)-(6) that the quantity J(r)/(Du s ) depends only on the geometry of the problem. It follows from the mathematical equivalence that Therefore, the main boundary condition at the droplet surface takes the form (15), where (∂T /∂n) air is calculated similarly to Eqns.(3)-(6):

F. Heat transfer in the gas phase
Here, we estimate the relative strength of the heat fluxes LJ(r) and k air (∂T /∂n) air , which are summands in (15). If follows from (18) that the heat fluxes are equal when Hence, the contribution of the effect of heat conduction in air matches the evaporation cooling heat flux for The present work considers the case T s −T a ∼ 10 K for droplets of ethanol and silicone oil, so taking into account the effect of heat transfer in the gas phase is worthwhile for obtaining quantitative results.

G. Thermal radiation
Here, we estimate the relative strength of the evaporative cooling heat flux LJ(r) and the thermal radiation heat flux (16) between the droplet and the environment.
It follows from Eqns.(3)-(6) that the heat fluxes are roughly equal for Hence, the contribution of the effect of thermal radiation between the droplet and the environment matches the evaporative cooling heat flux for The present work considers the case T s − T a ∼ 10 K for droplets of ethanol and silicone oil, so taking into account the effect of thermal radiation is worthwhile for obtaining quantitative results.

H. Numerical simulation
We solve the thermal conduction equation and the Navier-Stokes equations by the finite element method (FEM) using computational fluid dynamics simulation package ANSYS Fluent. The simulation includes several steps.

• Geometry and mesh generation
We create three-dimensional geometry of the axially symmetrical droplet. The surface curvature radius and the droplet height are equal to R/ sin θ and R/(1/ sin θ − cot θ), respectively.
The next step is to divide the simulation region into small computational cells. The number of cells is between 200 and 500 thousand depending on contact angle.

• Defining model and setting properties
We choose the pressure-based solver type, transient 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 density, viscosity, the temperature derivative of the surface tension σ T , which allows us to specify the properties of Marangoni stress, etc.

• Setting boundary conditions and solution method
The boundary conditions are set according to Sec. II. The rate of heat loss (15) is specified with a user-defined function (UDF) written in the C programming language. The solution algorithm is SIMPLE (Semi-Implicit Method for Pressure Linked Equations) [37]. The main feature of this method is usage 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.
• Post-processing We display the simulation results: velocity field, absolute values of velocity, and temperature distribution and then analyze the vortex structure and wave patterns. HTW instability mechanism [20,21] tends to move the BM cells from the droplet apex towards the contact line. Fig. 2 shows that the oscillation of the pattern is very slow compared to the characteristic frequences of HTWs in [11,12]. Therefore, we confirm the existence of irregular oscillatory BM convection rather than HTWs in this case. The Marangoni convection cells are sufficiently deep and are extended in the droplet bulk, as seen in Fig. 3.   Interestingly, we have found that the well-known analytical solution for the diffusion-limited evaporation rate [29], that significantly simplifies the calculation of the heat transfer in the droplet, can also be employed for calculation of thermal conduction in the gas phase.
Our results agree well with the experimental results of [25,27]. We confirm the conclusion of these works that the instability in this case is a nonstationary BM convection rather than HTWs.
We also clarify the behavior of the BM instabilities with decrease of the contact angle. We find that a decrease of the contact angle results in increasing BM cells, eventually transforming the flower pattern into the spatial structure, where the BM cells occupy the whole droplet. Such a behavior has been observed experimentally in [27] but has not been obtained numerically until now, to the best of our knowledge. We find that this behavior applies to the drying ethanol droplet of [25] as well. This resolves the discrepancy between the BM cell behavior in [27] and in [25].
We believe that our results represent a useful step towards a better understanding of the thermocapillary instabilities in evaporating droplets.