Marangoni Effects Control Flow in a Heated Evaporating Sessile Droplet

1. Introduction/Literature Review

The evaporation of sessile droplets is a basic interfacial process [1] that appears in many technologies, including microelectronics cooling [2], inkjet printing [3], and advanced surface patterning [4,5]. The performance of these systems depends strongly on how heat and mass are transported inside the droplet [6], where laminar convection can either enhance or limit evaporation [7]. Understanding how buoyancy and thermocapillary (Marangoni) forces compete is therefore important, and numerical modeling is now a standard tool to resolve these coupled mechanisms [8]. This work builds on previous multiphysics models of sessile droplet evaporation [9,10] and connects them to a broader context that includes nanoporous evaporation and microscale heat sinks [33,36].

Recent studies have examined how engineered fluids and surfaces can be used to control heat transfer. For example, Hussain et al. studied nanofluids in solar collectors and showed that nanoparticle suspensions can increase thermal efficiency in solar energy systems [32]. Yesudasan demonstrated that, in hydrophilic nanopores, the critical diameter for continuous evaporation lies between 3 nm and 4 nm, which sets a geometric scale for stable phase change in nanoscale confinements [33]. More recently, Yesudasan analyzed the thermal dynamics of heat pipes with sub-critical nanopores and showed that sub-critical diameters can strongly influence capillary flow and local evaporation balance [34]. Mohammed and Yesudasan used molecular dynamics simulations to examine how water structure and transport properties change inside nanopores [35], and related work on MEMS heat sinks has shown how small-scale channel and fin designs can enhance or limit overall thermo-hydraulic performance [36]. Together, these studies highlight that evaporation and convection are very sensitive to geometry and confinement at micro- and nanoscales.

The foundational physics of internal flow can be understood through a hierarchy of complexity, starting with buoyancy-driven natural convection. Akkuş et al. showed that a steadily fed hemispherical droplet can sustain a single recirculating cell when surface-tension gradients are absent, revealing a simple but limited buoyancy-controlled regime [11,12]. More detailed models add additional physical effects. Barmi and Meinhart [13] and Xu and Luo [14] reported that evaporation can evolve from an early Marangoni-dominated stage to a later capillary-dominated stage [15,16,17], although some models neglect buoyancy in the surrounding vapor [18]. Demirci et al. [19], together with earlier work by Erbil [20] and Ajaev [21], emphasized that thermocapillarity can dominate internal flow when droplets evaporate into their own vapor, but they also highlighted the uncertainty in kinetic parameters such as the mass accommodation coefficient.

At smaller length scales, additional transitions appear. Semenov et al. used simulations to show that pinned droplets can switch between diffusion-limited and kinetic-limited regimes depending on size, and Brutin summarized related phenomena in a comprehensive monograph on droplet wetting and evaporation [22,23]. Integral evaporation metrics have also been studied: Picknett and Bexon [24] and Nguyen and Nguyen [25] derived expressions for droplet lifetime under different evaporation modes, showing a nonlinear dependence on contact angle but assuming purely diffusive, isothermal conditions. More complete coupling between fluid flow and heat transfer was reported by Kang et al. [26], who included evaporative cooling, internal thermocapillarity, and external natural convection. Several recent studies have reported complementary insights into thermal transport and evaporation phenomena in micro- and nanoscale systems using both numerical and molecular approaches [32,33,34,35,36]. Christy et al. examined binary-mixture droplets and found rich flow transitions as composition and evaporation progressed [27]. These studies clarify how convection can both promote evaporation (by mixing cooler liquid with warmer regions) and impede it (by enhancing evaporative cooling at the surface).

Across these works, a consistent picture emerges. Pure buoyancy-driven flow, as in the reference case by Akkuş et al. [11], is rarely realized in practical systems that usually have some degree of Marangoni forcing. More recent experiments and simulations, such as those by Rocha et al. [28], show that solutal Marangoni effects can even dominate over thermal Marangoni effects in some mixtures. Many studies [13,19,26] identify thermocapillarity as a first-order mechanism controlling internal circulation and overall evaporation rates [29]. The dominant physics depends on contact angle, substrate superheat, and ambient composition [30,31,37], and the relative importance of each mechanism can change over time as the droplet shape evolves [13].

Even with this progress, important gaps remain. Flow-regime maps that clearly separate buoyancy-, thermocapillary-, and capillary-dominated regimes are still incomplete [38]. The quantitative contributions of each mechanism to heat and mass transfer are also not fully isolated across parameter space [39]. There is still uncertainty in key parameters such as the surface-tension temperature coefficient, which can significantly change predicted Marangoni strengths [40]. The present COMSOL model is designed to address some of these gaps by starting from a validated buoyancy-only baseline [11] and then adding Marangoni stresses and more detailed evaporation boundary conditions. This framework allows a controlled comparison between the two convective regimes and supports systematic studies of wall temperature and contact angle.

2. Methods and Materials

COMSOL Multiphysics 6.2 is used to model steady-state evaporation from a sessile water droplet on a heated wall in a two-dimensional axisymmetric domain. The axisymmetric formulation allows three-dimensional effects to be represented in a simpler 2D geometry by revolving the solution around the central axis. The computational domain is a hemispherical region subdivided into 18 subdomains, 55 boundaries, and 38 vertices, with each wedge sector spanning $5^{\circ }$. The central region represents the liquid droplet. The geometry is built step by step: a primary circular sector of ${90}^{\circ }$ is created, Boolean operations are applied to subtract the inlet, mirror operations extend the domain, and a final union operation assembles the full axisymmetric configuration.

A small circular inlet with radius $0.1$ mm at the center of the substrate provides a continuous mass supply at a controlled temperature. This inlet is subtracted from the main circle during the geometry construction. The key geometric features are summarized in Figure 1.

The mesh is generated using a free quadrilateral scheme with boundary-layer refinement near the inlet and along the free surface. The final mesh contains 6888 elements, with 6024 triangles in transition regions and 864 quadrilaterals in more regular zones. The boundary layers are refined to capture sharp gradients in velocity and temperature along the curved interface and near the inlet, as illustrated in Figure 2, Figure 3 and Figure 4. A mesh-refinement study with an “extra fine” setting increases the element count to 21,458 triangles and 1440 quadrilaterals. The difference in key output quantities between the finest two meshes remains below 1 %, showing that the chosen resolution is sufficient.

The fluid flow inside the droplet is modeled with the Laminar Flow interface. The governing equations are the stationary, incompressible Navier–Stokes equations with gravity $g=9.81m/s^2$ applied in the vertical z-direction. This captures natural convection due to density variations as well as forced convection due to the inlet flow. Heat transfer is modeled with the Heat Transfer in Fluids interface, accounting for both conduction and convection within the liquid. A reference temperature of $293.15$ K is used.

Marangoni stresses at the free surface are introduced through a temperature-dependent surface-tension relation

$$\sigma ={\sigma }^0-{\mathrm{Y}}_{T_{\mathrm{amb}}}\left(T-T_{\mathrm{amb}}\right),$$

(1)

where ${\sigma }^0=0.072N/m$ is the surface tension at ambient temperature, ${\mathrm{Y}}_{T_{\mathrm{amb}}}=1.6\mathrm{e}-4N/\left(m·K\right)$ is the surface-tension temperature coefficient, and the equilibrium contact angle is set to ${90}^{\circ }$. This relation is applied only at the interface, and COMSOL automatically converts the gradient of $\sigma$ into a shear stress along the free surface.

The boundary conditions are chosen to match the physical system, as illustrated in Figure 5. At the inlet, a normal inflow velocity is applied,

$$u_{\mathrm{in}}={\displaystyle \frac{{\dot{m}}_{\mathrm{total}}}{\rho \left(T_w\right)\pi r_{\mathrm{in}}^2}},$$

(2)

with the total mass flow adjusted to balance evaporation. Here, $\rho \left(T_w\right)$ is the liquid density evaluated at the wall temperature. At the free surface, a normal outflow velocity condition is used,

$$u={\displaystyle \frac{{\dot{m}}_{\mathrm{surf}}}{\rho \left(T_w\right)}},$$

(3)

together with the Marangoni shear stress generated by the temperature-dependent surface tension. The substrate is modeled as a no-slip wall with fixed temperature $T=T_w$, and the symmetry axis at $r=0$ enforces zero normal velocity and zero normal heat flux.

Thermally, the heated wall is maintained at $T_w$, while the free surface experiences both convective heat loss and evaporative cooling:

$$q=h_{nc}\left(T-T_{\mathrm{amb}}\right)+{\dot{m}}_{\mathrm{surf}}{h}_{fg}.$$

(4)

The axis of symmetry is adiabatic. All correlations for $h_{nc}$, $h_m$, and ${\dot{m}}_{\mathrm{surf}}$ are implemented via user-defined variables, and a small regularization parameter ${\epsilon }_{\mathrm{reg}}$ is introduced in denominators to avoid numerical singularities.

The liquid is water with temperature-dependent density, viscosity, thermal conductivity, and specific heat defined via interpolation from ${20}^{\circ }\mathrm{C}$ to ${70}^{\circ }\mathrm{C}$. Ambient air properties (density, thermal conductivity, viscosity, and specific heat) are treated as constants. Key derived variables are listed in Table 1, and constants are listed in Table 2.

The solver uses a stationary study with a temperature sweep from $T_w=303.15K$ to $313.15$ K with 11 equally spaced values. A direct solver with Newton–Raphson iterations, adaptive damping, and solution-based convergence is employed, solving about $1.6\times {10}^5$ degrees of freedom. After convergence, the solution is post-processed to compute local evaporative flux, total evaporation rate, and global circulation patterns. These numerical procedures are similar in spirit to those used in earlier droplet-evaporation studies [41,42].

3. Results and Discussion

3.1. Flow-Regime Transition: Buoyancy to Marangoni Dominance

Figure 6 shows the velocity magnitude and streamlines for the buoyancy-only case at $T_w=313.15K$ (${40}^{\circ }\mathrm{C}$). The flow is driven by density variations: warm fluid near the base rises along the axis, cools near the free surface, and then moves downward along the interface, forming a single toroidal vortex. The peak velocity is on the order of ${10}^{-4}\mathrm{m}/\mathrm{s}$, which is typical of buoyancy-dominated convection in small droplets [41,42].

When Marangoni stresses are added, the flow field changes dramatically, as shown in Figure 7. The interfacial temperature gradient generates a surface-tension gradient ${\mathrm{Y}}_{T_{\mathrm{amb}}}(T-T_{\mathrm{amb}})$ that pulls fluid along the surface from the warm apex (where surface tension is lower) to the cooler edge (where surface tension is higher). This surface motion drives a set of counter-rotating vortices inside the droplet. The characteristic velocity rises to $O\left({10}^{-3}\right)\mathrm{m}/\mathrm{s}$, consistent with a Marangoni-dominated regime [43,44]. In this regime, buoyancy still contributes, but the main circulation pattern is set by the interfacial shear.

3.2. Temperature Distribution and Thermal Boundary Layers

The different flow structures lead to different temperature fields. In the buoyancy-only case (Figure 8 and Figure 9), the isotherms are almost horizontal and stacked, which means that conduction is the main heat-transfer mechanism, with convection playing a secondary role. The single vortex slightly distorts the isotherms but does not strongly mix the entire droplet, so temperature gradients remain relatively thick and smooth [41,42].

In the combined buoyancy–Marangoni case, the temperature distribution becomes more complex. As shown in Figure 10 and Figure 11, the isotherms bend strongly towards the free surface. Warm fluid is rapidly pulled along the interface and redistributed, which thins the thermal boundary layer and promotes stronger mixing inside the droplet. This pattern is consistent with previous observations of thermocapillary-driven convection in evaporating droplets [10,28]. The resulting field is far from stratified, and the interior temperature becomes more uniform due to the multiple vortices.

3.3. Evaporative Flux Distribution and Mass Transfer

The local evaporative mass flux ${\dot{m}}_{\mathrm{surf}}$ depends on both surface temperature and vapor concentration above the droplet. Because Marangoni convection changes the surface temperature distribution, it also modifies the evaporation pattern. Over the studied range $T_w=303.15K$ to $313.15$ K, the average evaporative mass flux and total axisymmetric mass-loss rate are summarized in Table 3 and Table 4.

Figure 12 compares the average evaporative flux as a function of $T_w$ for the two cases. Both show a gradual increase with wall temperature, but the combined buoyancy–Marangoni case consistently has a higher flux at the same $T_w$. This indicates that stronger internal circulation refreshes the interface and sustains a larger driving force for evaporation, in line with observations in nanofluid-based collectors [32] and nanopore evaporation studies [33,34].

Figure 13 shows the total evaporative mass rate. The trend is similar: both curves rise with temperature, but the Marangoni-enabled configuration yields substantially higher total evaporation rates. The relative increase is on the order of 40–50 % over the range of temperatures studied. This behavior is analogous to the enhanced heat and mass transfer that has been reported in MEMS heat sinks with optimized microchannel and fin designs [36] and in nanoporous systems that support continuous evaporation [33].

4. Conclusions

This study compares buoyancy-driven and Marangoni-driven convection in a steadily fed sessile water droplet on a heated substrate using a two-dimensional axisymmetric COMSOL model. In the buoyancy-only case, the droplet develops a single toroidal vortex with peak velocities on the order of ${10}^{-4}\mathrm{m}/\mathrm{s}$ and a nearly stratified temperature field. When Marangoni stresses are added, the internal flow reorganizes into multiple counter-rotating vortices with velocities on the order of ${10}^{-3}\mathrm{m}/\mathrm{s}$. The temperature field becomes strongly distorted near the free surface, thermal boundary layers thin, and both the average evaporative mass flux and the total evaporation rate increase by roughly 40–50 % over the studied wall-temperature range. These results show that Marangoni stresses can dominate buoyancy even for modest superheats.

The present findings are consistent with broader trends in microscale and nanoscale heat-transfer research. Nanofluid-based solar collectors use enhanced thermal properties of suspensions to improve heat transfer in solar systems [32], while nanopore studies identify critical diameters for continuous evaporation and stable phase change [33,34]. Molecular dynamics simulations of confined water [35] and thermo-hydraulic analyses of MEMS heat sinks [36] also show that small changes in geometry and surface conditions can cause large changes in flow and heat-transfer behavior. In the same spirit, the current droplet model demonstrates that adding realistic Marangoni stresses is essential when predicting evaporation and cooling performance in practical applications such as microfluidic cooling, additive manufacturing, and thermal management of compact electronic devices.