MODELING THE EVAPORATION OF DRYING SESSILE DROPLETS WITH BUOYANCY DRIVEN INTERNAL CONVECTION

Droplet evaporation is a fundamental phenomenon encountered in diverse applications such as inkjet printing, DNA mapping, film coating, and electronics cooling. Modeling the evaporation process of a sessile droplet is complicated because of the coupling of several physical phenomena occurring in different phases and various magnitudes such as the buoyant convection of the liquid in millimeter size droplets and that of the surrounding air/water vapor mixture, in the order of meters. In this study, the theoretical framework presented previously for the steadily fed droplets [Int J Therm Sci, 158 (2020) 106529] is extended to resolve the evaporation of drying droplets with a pinned contact line. Based on the quasi-steady-state assumption, buoyant convection inside the droplet and diffusive-convective transport of vapor in the gas domain are modeled. As a test case, drying process of a water droplet with a 68◦ initial contact angle on a heated substrate is simulated and the predictions of the model are interpreted.

evaporation process of a sessile droplet is complicated because of the coupling of several physical phenomena occurring in different phases and various magnitudes such as the buoyant convection of the liquid in millimeter size droplets and that of the surrounding air/water vapor mixture, in the order of meters. In this study, the theoretical framework presented previously for the steadily fed droplets [Int J Therm Sci, 158 (2020) 106529] is extended to resolve the evaporation of drying droplets with a pinned contact line. Based on the quasi-steady-state assumption, buoyant convection inside the droplet and diffusive-convective transport of vapor in the gas domain are modeled. As a test case, dry-ing process of a water droplet with a 68 • initial contact angle on a heated substrate is simulated and the predictions of the model are interpreted.

INTRODUCTION
Droplet evaporation has been extensively studied by the researchers due not only to the diversity of its applications but also its complexity. Modeling of this multiphase problem is not straightforward; therefore, past studies highlighted only restricted aspects of the problem. There is an extensive body of literature regarding droplet evaporation. Here we highlight several of them.
Despite its small size, energy transport is not restricted with the conduction in evaporating droplets [1]. Consequently, convection mechanisms should be taken into consideration. Outside the droplet, vapor diffusion is not the sole transport mechanism, which necessitates the inclusion of the gas convection mechanisms (e.g. gas buoyancy) [2,3]. Substrate conductivity can affect the internal flow by reversing the direction of circulation [4]. Non-monotonic evaporation distribution can lead to non-monotonic interfacial temperature distribution for droplets with buoyancy-driven internal convection [3,5]. Stefan flow of air substantially contributes to the evaporation. Omission of it results in underestimated evaporation rates [6,7].
The presence of thermocapillary flow in water droplet is controversial in the literature [8,9]. Reason of the conflicting experimental results was attributed to the high sensitive of water droplet surfaces to contamination [10]. Thus, Marangoni flow do not always initiate. However, buoyant forces are always present and shape the velocity field inside the droplet. Therefore, in this work, the effect of buoyancy on the evaporation of drying sessile droplets is investigated by developing a model that solves full compressible Navier-Stokes equations by utilizing temperature dependent thermophysi-cal properties in both phases.

MODELING
Mass and energy transfer is inherently coupled at the interface. Therefore, an iterative modeling approach that is adopted by past studies [3,5,11], is utilized. Quasi-steady-state successive simulations are carried out for droplets with reducing contact angles to simulate the drying stages of a sessile droplet. The connections between the droplets with reducing contact angle are established by assigning the deformation velocity at the interface.
The model is utilized to simulate an experiment conducted in [12], where a pinned sessile droplet placed on a flat heated hydrophilic substrate evaporates into the air. Due to its small size, droplet has a shape of spherical cap. An axisymmetric modeling approach is adopted and the corresponding problem domain is shown in Figure 1. Gas volume is selected sufficiently large to eliminate any boundary effect on the problem.
Conservation equations for mass, linear momentum, and energy are solved in both phases. Additionally, vapor transport equation is solved in the gas phase. The governing equations are summarized below: Symmetry boundary conditions are applied at the axis. At the far field, in the case of inflow, constant ambient temperature and pressure values are assigned, while, in the case of outflow, At the liquid-vapor interface, mass and energy balances are secured. Evaporative mass flux determines the normal component of the interfacial liquid velocity as follows: Interfacial gas velocity is determined based on the Stefan flow of air: Tangential velocity is the same at the interface for both phases. Evaporative mass flux is calculated based on the diffusive and conductive transfer of vapor at the interface:ṁ Although the majority of energy is transferred by the evaporation at the interface, conduction to gas and radiation to the ambient still contributes to the interfacial energy transfer: Governing equations are solved using Finite Element Method based solver of COMSOL Multiphysics software. Embedded grid generator of COMSOL is utilized to create the solution mesh. Mesh independency of the solution is secured by conducting mesh iterations.

RESULTS AND DISCUSSION
Geometrical parameters and far field conditions utilized in the simulation are summarized in Table 1. Droplet lifetime is predicted as 636 s by the simulation. This duration was 638 s in the experiment of Sobac and Brutin [12]. Consequently, the deviation is only 0.3%, suggesting the accuracy of the modeling.
Resultant outer flow of the gas and the temperature distribution are shown in Figure 2, where the droplet is not apparent because of its small size with respect to the solution domain. It should be noted that a vortex structure forms near the substrate, which has a certain influence on the evaporation rate of the droplet. At the initial stages of the evaporation (e.g. at relatively higher contact angles), a buoyancy induced flow is observed. Heated liquid raises along the centerline, during which it cools down. Then the liquid turns and moves along the interface, during which it gets warm. This circulation forms a convection cell as shown in Figure 3. Flow of the vapor stemming from the interface is shown by white streamlines in the same figure. The vapor gets away from the interface under the influence of outer gas flow.
When the droplet becomes flatter (e.g. at relatively smaller contact angles), radial flow becomes apparent as shown in Figure 4. This mechanism is responsible for the well-known coffee ring effect [13]. It arises in order to replenish the liquid at the contact line, where the evaporation rate is divergent due to thinning film.  Distributions of the gas velocity magnitude are shown in Figure 5 and Figure 6 for the contact angles of 68 • and 20 • , respectively. Bound-ary layer forms on the substrate due to no slip boundary conditions. Near the contact line, jetlike flow of the gas indicates the intensified evaporation rates. Above the interface, gas flow field is shaped based on the interaction of outer flow and Stefan flow. In the case at θ = 68 • , Stefan flow is in the same orientation with the outer gas flow. Hence, a smooth upward gas flow forms near the interface. In the case at θ = 20 • , Stefan flow opposes the outer gas flow. As a result, a stagnation point forms above the apex. It should be noted that gas flow field is highly influential on the evaporation rates. Therefore, proper coupling of the phases is a must in the modeling of droplet evaporation. This the reason behind the underestimated evaporation rates predicted by early studies [14].
To conclude, predictions of the model excellently match with the experimental results in terms of global evaporation parameters, e.g. droplet lifetime. While buoyancy induced convection governs the internal liquid flow at the higher contact angles, radial flow is effective at the smaller contact angles. Future work will focus on the transition between these flow patterns. Moreover, thermocapillarity will be included into the model.