2.1. Problem statement

A sketch of the problem is provided in figure 1(a). A sessile drop of initial height ${\hat{H}}_{0}$ and contact radius $\hat{R}_{0}$ resting on a heated substrate evaporates into a non-saturated surrounding gas at atmospheric pressure; the caret denotes dimensional variables. The liquid is a pure substance (distilled water in the experiments) while the gas is modelled as a two-component mixture of variable composition: a non-condensable gas (nitrogen in the experiments) and the liquid’s vapour. All three pure substances are regarded as Newtonian fluids. The density, dynamic viscosity, thermal conductivity, and specific heat capacity for the liquid are denoted by $\hat{{\it\rho}}_{l}$ , $\hat{{\it\mu}}_{l}$ , $\hat{k}_{l}$ , and ${\hat{c}}_{pl}$ , respectively. The equivalent physical properties for the non-condensable gas (designated by subscript ‘1’) and vapour (designated by subscript ‘2’) are represented by $\hat{{\it\rho}}_{g1}$ , $\hat{{\it\mu}}_{g1}$ , $\hat{k}_{g1}$ , ${\hat{c}}_{pg1}$ , and $\hat{{\it\rho}}_{g2}$ , $\hat{{\it\mu}}_{g2}$ , $\hat{k}_{g2}$ , ${\hat{c}}_{pg2}$ , respectively. The relative amount of the gaseous components is represented by the vapour mass fraction, ${\it\omega}$ , defined as the mass of vapour per unit mass of gas mixture. The coefficient of binary molecular diffusion is $\hat{D}$ , and the specific latent heat is ${\rm\Delta}{\hat{h}}_{v}$ . The liquid–gas interface is endowed with a surface tension, $\hat{{\it\sigma}}$ , which decreases monotonically with temperature, $\hat{T}$ , according to $\hat{{\it\sigma}}=\hat{{\it\sigma}}_{0}-\hat{{\it\gamma}}(\hat{T}-\hat{T}_{a})$ ; here, $\hat{{\it\sigma}}_{0}$ is the surface tension at the reference temperature $\hat{T}_{a}$ , and $\hat{{\it\gamma}}=-\partial \hat{{\it\sigma}}/\partial \hat{T}$ denotes its temperature-dependence coefficient. The substrate is considered to remain isothermal at temperature $\hat{T}_{w}$ while the gas far from the drop is at temperature $\hat{T}_{a}$ .

Figure 1. Schematic of (a) the mathematical model and (b) the experimental set-up.

The following scaling is introduced to render the flow variables dimensionless:

Here, $\boldsymbol{x}=(x,\,y,\,z)$ and $\boldsymbol{u}=(u,\,v,\,w)$ are the coordinate and velocity vectors, $\hat{U} _{0}=\hat{{\it\gamma}}{\rm\Delta}\hat{T}/\hat{{\it\mu}}_{l}$ represents the characteristic thermocapillary velocity with ${\rm\Delta}\hat{T}=\hat{T}_{w}-\hat{T}_{a}$ , $p$ is the pressure, and $t$ denotes the time. The physical properties of the non-condensable gas and vapour are taken into consideration via the following ratios:

The gas-mixture properties, denoted by $\hat{{\it\rho}}_{g}$ , $\hat{{\it\mu}}_{g}$ , $\hat{k}_{g}$ and ${\hat{c}}_{pg}$ , depend on the relative amount of each component as follows:

The initial drop geometry ( $H_{0}$ , $R_{0}$ ) and the dimensions of the computational domain ( $L_{x}$ , $L_{y}$ , $L_{z}$ ) are

where $\hat{l}_{x}$ , $\hat{l}_{y}$ and $\hat{l}_{z}$ are the width, depth and height of the computational region.

2.2. Governing equations

The fully coupled two-phase dynamics of the drop and surrounding gas is modelled via DNS with a DI method (see Anderson, McFadden & Wheeler Reference Anderson, McFadden and Wheeler1998 for a review of this methodology). According to this method, the mathematically sharp interface is replaced with an interface of finite thickness, which can be prescribed (Jacqmin Reference Jacqmin1999). We use a DI method based on that developed by Ding, Spelt & Shu (Reference Ding, Spelt and Shu2007) for pure hydrodynamic problems (following Jacqmin Reference Jacqmin1999, and Badalassi, Ceniceros & Banerjee Reference Badalassi, Ceniceros and Banerjee2003), but with a number of extensions to address non-isothermal conditions and phase change.

DI methods are grounded in two fundamental principles. The first was given by van der Waals (Reference van der Waals1979). If the interface between two fluids is regarded as a narrow layer of finite thickness, the volumetric free energy density $\hat{f}$ across it is dependent on the local composition (represented here by the order parameter $c$ ) and their gradients,

where $\hat{{\it\varepsilon}}$ is a measure of the interface thickness, ${\it\Psi}(c)$ is the bulk energy density with minimum levels corresponding to the fluids’ two stable phases, and ${\it\alpha}$ is a dimensionless adjustment parameter. In addition, van der Waals (Reference van der Waals1979) hypothesized that equilibrium interface profiles are those that minimize the total free energy $\hat{\mathscr{F}}=\int _{{\it\Omega}}\hat{f}\,\text{d}\hat{V}$ , where ${\it\Omega}$ represents the interfacial domain. From the calculus of variations (Jacqmin Reference Jacqmin1999), these profiles must satisfy

where $\hat{{\it\phi}}=\partial \hat{\mathscr{F}}/\partial c$ is the chemical potential. The second key element of DI methods was given by Cahn & Hilliard (Reference Cahn and Hilliard1958, Reference Cahn and Hilliard1959) and Cahn (Reference Cahn1961), who extended the van der Waals hypothesis to time-dependent situations, postulating that irregular interface profiles recover equilibrium by virtue of diffusion fluxes which are proportional to chemical potential gradients. On these premises, and selecting the liquid volume fraction as an order parameter ( $c=1$ in the liquid, $c=0$ in the gas and $0<c<1$ in the interface region), it can be shown that the spatiotemporal evolution of $c$ is governed by the phase-field, advection–diffusion Cahn–Hilliard equation (Ding et al. Reference Ding, Spelt and Shu2007)

where $\mathit{Pe}=\hat{U} _{0}{\hat{H}}_{0}/(\hat{M}_{0}\hat{{\it\phi}}_{0})$ is the Péclet number, $M=c(1-c)$ is a diffusive coefficient referred to as ‘mobility’, ${\it\phi}={\it\varepsilon}^{-1}[{\it\Psi}^{\prime }(c)-{\it\varepsilon}^{2}{\rm\nabla}^{2}c]$ is the chemical potential, and ${\it\Psi}(c)={\textstyle \frac{1}{4}}c^{2}(1-c)^{2}$ . With appropriate boundary conditions, (2.7) conforms to a mass-conservative interface-capturing method, in which the interface profile is kept regular by means of the diffusive term. Some recent examples in which this equation has been used include droplet spreading (Ding & Spelt Reference Ding and Spelt2007a ), droplet motion under shear flow (Ding & Spelt Reference Ding and Spelt2008; Ding, Gilani & Spelt Reference Ding, Gilani and Spelt2010), and two-layer flows (Valluri et al. Reference Valluri, Naraigh, Ding and Spelt2010).

For an evaporating drop, the total mass of liquid and vapour is globally conserved, but the mass of liquid is not; this decreases in time due to the liquid-to-gas phase change. Thus, accounting for interface mass transfer requires adding an interface sink/source term in (2.7), which leads to the phase change version of the Cahn–Hilliard equation,

where $S$ is the volumetric interface mass-transfer rate (only computed in the interface region). The negative sign on the right-hand side of (2.8) results from our sign convention, namely $S>0~(S<0)$ for evaporation (condensation). Following the traditional assumption that the gas at the interface is (quasi-) saturated with vapour (Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten1997; Picknett & Bexon Reference Picknett and Bexon1977; Hu & Larson Reference Hu and Larson2002; Cazabat & Guena Reference Cazabat and Guena2010), which is valid for slow evaporation cases like the one considered here, $S$ is calculated as follows:

where $p_{s}$ denotes the saturation vapour pressure, $\mathscr{M}=\hat{M}_{g2}/\hat{M}_{g1}$ is the vapour-gas molar weight ratio, $\mathscr{P}=\hat{{\it\rho}}_{l}\hat{U} _{0}^{2}/\hat{p}_{r}$ , is a dimensionless group including the absolute pressure $\hat{p}_{r}$ , and $t_{s}$ is a measure of interface mass-transfer rate small enough to maintain saturation condition at the interface. The first term within the brackets represents the saturation vapour mass fraction calculated using Raoult’s law and Dalton’s law of partial pressure, i.e. it is assumed that the gaseous components form an ideal mixture. The detailed justification of this expression is given in Sáenz et al. (Reference Sáenz, Valluri, Sefiane, Karapetsas and Matar2014). The saturation pressure dependence on $T$ is approximated with a standard Antoine expression, $\log _{10}(\hat{p}_{s})=\hat{A}-\hat{B}/({\hat{C}}+\hat{T})$ , where $\hat{A}$ , $\hat{B}$ and ${\hat{C}}$ are the empirical coefficients.

The overall mass balance is completed with the continuity equation, while the conservation of momentum is governed by the Navier–Stokes equations

where ${\it\rho}=c+{\it\rho}_{g}(1-c)$ and ${\it\mu}=c+{\it\mu}_{g}(1-c)$ are the one-fluid form of the density and viscosity respectively, $\mathit{Re}=\hat{{\it\rho}}_{l}\hat{U} _{0}{\hat{H}}_{0}/\hat{{\it\mu}}_{l}$ is the Reynolds number, $\mathit{We}=\hat{{\it\rho}}_{l}\hat{U} _{0}^{2}{\hat{H}}_{0}/(6\sqrt{2}\hat{{\it\sigma}}_{0})$ denotes the Weber number, ${\it\phi}\boldsymbol{{\rm\nabla}}c$ is the DI interpretation of the mean surface tension stress enforcing (Jacqmin Reference Jacqmin1999; Ding et al. Reference Ding, Spelt and Shu2007), and $(\boldsymbol{{\rm\nabla}}_{~s}T){\it\delta}$ models the thermocapillarity effect, with $\boldsymbol{{\rm\nabla}}_{~s}=\boldsymbol{{\rm\nabla}}-\boldsymbol{n}(\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})$ being the gradient tangent to the interface, $\boldsymbol{n}=\boldsymbol{{\rm\nabla}}c/|\boldsymbol{{\rm\nabla}}c|$ the inwards unit vector normal to the interface, and ${\it\delta}=|\boldsymbol{{\rm\nabla}}c|$ . It should be noted that several authors (Jacqmin Reference Jacqmin2000; Ding & Spelt Reference Ding and Spelt2007a ; Yue, Zhou & Feng Reference Yue, Zhou and Feng2010) have already shown that the results from this DI approach converge to the sharp-interface limit in the limit of zero interfacial thickness.

The conservation of energy in both phases is computed via

Here, $k=c+k_{g}(1-c)$ and $c_{p}=(c+{\it\rho}_{g}c_{pg}(1-c))/(c+{\it\rho}_{g}(1-c))$ are the one-fluid form thermal conductivity and specific heat capacity, and $\mathit{Pr}=\hat{{\it\mu}}_{l}{\hat{c}}_{pl}/\hat{k}_{l}$ and $\mathit{Ja}={\hat{c}}_{pl}{\rm\Delta}\hat{T}/{\rm\Delta}{\hat{h}}_{v}$ are the Prandtl and Jakob numbers, respectively. Finally, the vapour mass fraction, ${\it\omega}$ , is governed by a general advection–diffusion transport equation only solved in the gas phase:

where $\mathit{Sc}=\hat{{\it\mu}}_{l}/(\hat{{\it\rho}}_{l}\hat{D})$ is the Schmidt number.

2.3. Boundary and initial conditions

The boundary conditions are selected to mimic the experiments presented in § 3. The lower boundary ( $z=0$ ) is modelled as an isothermal $(T=1)$ , no-slip $(\boldsymbol{u}=0)$ solid substrate with non-penetration of mass $(\boldsymbol{n}_{z}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\omega}=0)$ . Here, $\boldsymbol{n}_{x}$ , $\boldsymbol{n}_{y}$ , and $\boldsymbol{n}_{z}$ denote the unit vectors in the principal directions. The Cahn–Hilliard equation requires the Neumann boundary condition $\boldsymbol{n}_{z}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}=0$ to guarantee the volume conservation of each fluid (Ding & Spelt Reference Ding and Spelt2007b ). Another condition on $c$ is necessary, this being a Dirichlet or Neumann condition in pinned or MCL configurations, respectively. A number of models have been proposed to solve the dynamics of moving contact lines, e.g. Qian, Wang & Sheng (Reference Qian, Wang and Sheng2006), Ding & Spelt (Reference Ding and Spelt2007b ), Carlson, Do-Quang & Amberg (Reference Carlson, Do-Quang and Amberg2009) and Yue et al. (Reference Yue, Zhou and Feng2010). The approach suggested by Ding & Spelt (Reference Ding and Spelt2007b ) is adopted here. The contact angle is therefore prescribed by virtue of the boundary condition

where ${\it\theta}_{s}$ is the microscale contact angle, and $\boldsymbol{n}_{w}$ and $\boldsymbol{t}_{w}$ are the unit vectors normal and tangential to the solid wall, respectively (Ding & Spelt Reference Ding and Spelt2007b ). It is important to emphasize that prescribing a slip velocity is unnecessary, given that the slow contact-line motion is solved by virtue of finite diffusion (Sui et al. Reference Sui, Ding and Spelt2014). The temperature of the upper boundary $(z=L_{z})$ is set to the ambient temperature as recorded by the secondary thermocouple $(T=0)$ , while Neumann boundary conditions ( $\boldsymbol{n}_{x}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}T=0$ and $\boldsymbol{n}_{y}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}T=0$ ) are used at the vertical walls of the computational domain ( $x=\pm L_{x}/2$ and $y=\pm L_{y}/2$ ).

In the experimental set-up, the size of chamber was intentionally designed to be very large in comparison to that of the drop, so this would evaporate into a much larger volume of pure $\text{N}_{2}$ (dry gas). An approximately infinite amount of $\text{N}_{2}$ surrounding the drop makes it reasonable to assume that the phase change does not result in a noteworthy increment in the bulk humidity. Thus, the zero Dirichlet boundary condition is imposed in the vertical and top boundaries ( ${\it\omega}=0$ at $x=\pm L_{x}/2$ , $y=\pm L_{y}/2$ , and $z=L_{z}$ ). Provided that the computational domain is large enough, the vapour distribution around the drop becomes independent of the far-field boundary conditions. For the same reason, it is considered that the flow far away from the drop is at rest ( $\boldsymbol{u}=0$ at $x=\pm L_{x}/2$ , $y=\pm L_{y}/2$ , and $z=L_{z}$ ). Finally, the fact that the species equation is only solved in the gas phase requires the implicit imposition of the no-penetration condition ( $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\omega}=0$ , and $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{u}=0$ ) along the species-field boundary in the liquid–gas interface (contour $c=0.9$ ).

Initial conditions are also prescribed. At $t=0$ , the drop is a spherical cap of height $H_{0}=1$ and base radius $R_{0}$ ; the two-phase flow is at rest $(\boldsymbol{u}=0)$ , and the gas-mixture is completely dry $({\it\omega}=0)$ . The initial pressures in the gas, and within the drop are set to zero, and the Laplace pressure, respectively. In dimensionless form, the latter becomes ${\rm\Delta}p=1/(3\sqrt{2}WeR_{s})$ with $R_{s}=(H_{0}^{2}+R_{0}^{2})/(2H_{0})$ . From the experiments, it is known that the chamber was allowed to reach thermal equilibrium before the drop was deposited onto the substrate, and that the liquid dosed through the pipette was pumped from a reservoir located outside the chamber and was therefore at ambient temperature. These two considerations led us to select the initial drop temperature as $T=0$ , and the initial temperature of the gas as a linear distribution from the substrate to the ambient temperature, i.e.  $T=(L_{z}-z)/L_{z}$ . Both conditions are simultaneously met with $T=(1-c)(L_{z}-z)/L_{z}$ .

2.4. Numerical method

The governing equations (2.8), (2.10)–(2.13), coupled with the volumetric interface mass-transfer rate calculated with (2.9), constitute a system of seven partial differential equations (PDEs) with seven unknowns ( $c$ , $p$ , $u$ , $v$ , $w$ , $T$ and ${\it\omega}$ ). This system of equations is simultaneously solved (DNS) with a finite-volume (FV) discretization of the physical domain by means of a staggered marker-and-cell (MAC) grid: scalar variables ( $c$ , $p$ , $T$ , ${\it\omega}$ ) are stored in the cell centres while the velocity components ( $u$ , $v$ , $w$ ) are defined at the control-volume faces. Spatial derivatives are approximated by a centred scheme. Mesh refinement tests are carried out to ensure that the results are grid-independent.

The Cahn–Hilliard equation is solved numerically using the implicit–explicit strategy proposed by Badalassi et al. (Reference Badalassi, Ceniceros and Banerjee2003) with a second-order semi-backward difference formula (SBDF) scheme (see Ascher, Ruuth & Wetton Reference Ascher, Ruuth and Wetton1995). The convective term is explicitly computed approximating the fluxes at cell faces with an upwind (referred to the flow direction) fifth-order weighted essentially non-oscillatory (WENO) scheme (Liu, Osher & Chan Reference Liu, Osher and Chan1994) to increase the code’s stability to the shock-like nature of the $c$ profile across the interface (Ding et al. Reference Ding, Spelt and Shu2007). The resolution of the velocity and pressure fields is achieved by a standard projection method, wherein the intermediate velocity (without considering $p$ ) is found with the second-order Crank–Nicolson Adams–Bashforth (CNAB) scheme to approximate the diffusive (implicitly) and advective (explicitly) terms, respectively. The pressure is then added to update the velocity to its final value ensuring that the continuity condition is satisfied. The energy and species fields are advanced in time with similar semi-implicit CNAB techniques.

The solution strategy requires us to know the value of all the variables at the current time-step $(n)$ as well as $c$ at the previous time-step $(n-1)$ . The computation then proceeds by updating the volume fraction via (2.8) with $\boldsymbol{u}$ and $S$ evaluated at $(n)$ . Equations (2.13) and (2.12) are then solved sequentially. This order results naturally from the fact that the former involves $c^{n+1}$ (to evaluate ${\it\rho}_{g}$ ) while the latter additionally requires ${\it\omega}^{n+1}$ (to evaluate ${\it\rho}$ , $k$ and $c_{p}$ ). Finally, with $c^{n+1}$ , ${\it\omega}^{n+1}$ and $T^{n+1}$ , the velocity is updated solving (2.10) and (2.11) with the surface tension and physical properties evaluated at $(n+1/2)$ , as recommended by Ding et al. (Reference Ding, Spelt and Shu2007). Following the study conducted by Ding et al. (Reference Ding, Spelt and Shu2007), the Péclet number is set to $\mathit{Pe}=1/{\it\varepsilon}^{2}$ , choosing ${\it\varepsilon}$ to be proportional to the grid resolution, i.e.  ${\it\varepsilon}=0.5\min ({\rm\Delta}x,\,{\rm\Delta}y,\,{\rm\Delta}z)$ . This configuration leads to an interface thickness, i.e. a distance between contours of $c=0.1$ and 0.9, of around three times the minimum grid spacing. The relevant time scales in the system are the viscous, conductive, convective, molecular diffusivity and evaporative time scales. The marching time is selected to be smaller than all of them, i.e.  ${\rm\Delta}t<\min ({\rm\Delta}\hat{x}^{2}/\hat{{\it\nu}}_{l},\,{\rm\Delta}\hat{x}^{2}/\hat{{\it\alpha}}_{l},\,{\rm\Delta}\hat{x}/\hat{u} _{0},\,{\rm\Delta}\hat{x}^{2}/\hat{D},\,\hat{{\it\rho}}_{g}/{\hat{S}}_{max})/\hat{t}_{0}$ . This temporal scale is also used for  $t_{s}$ . The numerical method is coded in Fortran 90 for implementation on a shared-memory architecture, using Open Multi-Processing (OpenMP). The code is run on the supercomputer ARCHER (Cray XC30 Supercomputer) using a 2.7 GHz, 12-core E5-2697 v2 (Ivy Bridge) series processor for a typical simulation.

3.1. Experimental apparatus

An experimental apparatus was designed to measure the profile of evaporating drops using optical techniques. A schematic of the set-up is presented in figure 1(b). Very small drops $(\hat{V}_{0}\sim 10~{\rm\mu}\text{l})$ of distilled water were gently deposited on a heated substrate, allowed to evaporate, and simultaneously recorded from the side by a CCD camera and from the top by an infrared (IR) camera. The substrate was a copper cylinder of dimensions 3.8 cm diameter and 3.8 cm length. This was maintained at constant temperature, $\hat{T}_{w}$ , by means of a cartridge heater inserted within the cylinder from below (Omega CSS-10120/120V, 20W, 2.54 cm length and 0.64 cm diameter), a PID controller $\pm 0.2\,^{\circ }\text{C}$ (Omega CN77544-A2) and a K-type thermocouple located in a small hole 6 mm below the substrate’s upper surface. Copper and the cylinder dimensions were chosen to provide a substrate with very high thermal conductivity as well as very large thermal mass so that the surface remains isothermal throughout the evaporation process. To ensure a smooth finish, the copper substrate was progressively ground, and polished (with $0.05~{\rm\mu}\text{m}$ colloidal alumina suspension in the final stage). The emissivities for water and polished copper are ${\it\epsilon}=0.96$ and 0.03, respectively, which results in high-contrast IR images wherein the interface temperature and contact line are sharply defined.

The size of the drops was controlled via the pipette tip size, tube diameter, and peristaltic pump speed. In order to avoid imprecision in assessing the composition (physical properties), and relative humidity (driving force) of the air surrounding the drop, the substrate and dosing mechanism were placed within a closed, acrylic cubic chamber (dimensions $20\times 20\times 20~\text{cm}$ and 0.54 cm wall thickness) maintained at atmospheric pressure. A system of valves was attached to this chamber to allow the replacement of the ambient air by pure nitrogen ( $\text{N}_{2}$ ). The use of a closed chamber also guaranteed that the evaporation process was not distorted by externally induced convection currents and that the IR readings were not contaminated by radiation from the surroundings, since acrylic is opaque in the spectrum range of our IR camera. Inner reflections were also minimized by painting the interior of the chamber with black paint, except for two small gaps on opposite sides to allow the CCD recording and back-lighting of the drop, respectively.

The drop profiles were captured with a CCD camera ( $900~\text{pixel}\times 600~\text{pixel}$ , $9~{\rm\mu}\text{m}~\text{pixel}^{-1}$ ), capable of recording up to 30 frames per second (f.p.s.), connected to a video-digitizer board (frame grabber). LED back-lighting was employed to improve the video’s contrast without raising the temperature. The recorded side images were later post-processed using the Droplet Shape Analyser from Krüss (DSA v1.9, Krüss GmbH, Hamburg, Germany) in order to obtain the instantaneous height ${\hat{H}}$ , contact angle ${\it\theta}$ , base radius $\hat{R}$ and volume $\hat{V}$ throughout the evaporation process. Simultaneously, a midwave IR camera mounted directly above the substrate and facing vertically downwards onto the drop (FLIR Silver SC5600, spectrum range $3{-}5~{\rm\mu}\text{m}$ , $640~\text{pixel}\times 512~\text{pixel}$ , 100 f.p.s.) was used to record the interface temperature field and contact-line dynamics.

A cylindrical black body, as described in Kim et al. (Reference Kim, Kommer, Dessiatoun and Kim2012), was used to calibrate the IR readings. Before every experimental run, $\text{N}_{2}$ was blown through the chamber for a sufficiently long period of time to make sure that this was the only gas present within the test cell. The inlet and outlet valves were then closed and the system was allowed to reach equilibrium. Steady state was recognized by means of a secondary thermocouple located in the gas, 3.5 cm above the heated substrate, which also provided the real ‘bulk’ or ‘ambient’ temperature, $\hat{T}_{a}$ . Only after the system reached equilibrium, detected when $\hat{T}_{a}$ became steady, was a single drop released onto the heated surface and the recording process started. This procedure was repeated for every drop in this investigation to ensure consistency.

4.1. Comparisons with experiments

The validation of the numerical solutions is more challenging than in the non-heated configuration. This is not only because of the notable numerical difficulty associated with solving seven fully coupled non-linear PDEs, but also due to the large thermal gradients arising at the interface. Given that $p_{s}$ grows exponentially with $T$ , deviations as small as $\pm 2\,^{\circ }\text{C}$ at high temperatures lead to substantial changes in $m(t)$ .

Initially, the flow undergoes a transient process wherein the drop heats up and the vapour disperses in the gas. Simultaneously, the Marangoni effect drives fluid from the contact-line region (hot) towards the apex (cold) and therefore induces convection currents in both phases. This transient phase is examined in more detail in § 4.2. After this stage, the flow reaches a quasi-steady state wherein the two-phase flow characteristics evolve more gradually with $t$ in accordance with the evaporation time scale. The instantaneous evolution of ${\it\theta}$ , $H$ , $R$ and $V$ obtained with the simulations is compared with the experimental data in figure 4. Excellent agreement between experiments and simulations is observed in the initial CR stage. The average numerical evaporation rate is $m(t)=1.963\times 10^{-6}$ , 2.7 % lower than that experimentally observed, i.e.  $2.017\times 10^{-6}$ . The deviation increases slightly for $t>10.5\times 10^{5}$ , when the experimental drop evaporates according to the combined $\text{CA}+\text{CR}$ and SS modes. This is to be expected, as in the simulations the drop remains perfectly pinned and therefore the evaporation modes are not comparable at this stage. A more detailed discussion examining the instantaneous $m(t)$ is presented in § 6.

Figure 4. Comparison between experiments and numerical simulations for a water drop evaporating into $\text{N}_{2}$ placed on a heated substrate at $\hat{T}_{w}=40\,^{\circ }\text{C}$ with ${\rm\Delta}\hat{T}=14.6\,^{\circ }\text{C}$ . (a) Contact angle ${\it\theta}$ , (b) drop height $H$ , (c) contact radius $R$ , and (d) volume $V$ . The vertical dashed lines separate the evaporation modes observed experimentally: see § 3.2.

The second key point in the validation is the comparison between the experimental and numerical interface temperature, shown in figure 5(a). Overall, there is again good qualitative and quantitative agreement. The only discrepancies are found at the contact line, but these result from the fact that the IR emissions are distorted due to the transition from water (high emissivity) to polished copper (very low). As the contact line is approached, the accuracy of the $T$ readings diminishes due to the transition to a different emissivity; otherwise $T\simeq 1$ should be observed at the contact line with the IR camera. Figure 5(b) shows the distribution of $T$ in the bulk of the liquid and gas. The flow and liquid temperature is qualitatively similar to that previously presented by other authors (Mollaret et al. Reference Mollaret, Sefiane, Christy and Veyret2004; Hu & Larson Reference Hu and Larson2005a ; Girard et al. Reference Girard, Antoni, Faure and Steinchen2006; Yang et al. Reference Yang, Hong and Cheng2014). The gas temperature, not available in the literature, reveals a colder area above the drop due to the evaporative cooling. This thermal dip is progressively corrected as we move away from the drop until, at a certain distance $(r>6)$ , $T$ becomes $r$ -independent and decreases uniformly for increasing $z$ . The vapour distribution is presented in figure 5(c).

Figure 5. Water drop evaporating into $\text{N}_{2}$ with $\hat{T}_{w}=40\,^{\circ }\text{C}$ . (a) Comparison between the experimental and numerical interface temperature with $H=0.6$ , (b) liquid–gas temperature, and (c) vapour mass fraction in the initial stages of the evaporation process.

Figure 6. Initial evolution of the two-phase velocity for a water drop placed on a heated substrate at $\hat{T}_{w}=40\,^{\circ }\text{C}$ with ${\rm\Delta}\hat{T}=14.6\,^{\circ }\text{C}$ . (a) Initial conditions, $t=0$ , (b) $t=750$ (0.42 s), (c) $t=1.5\times 10^{3}$ (0.83 s), (d) $t=2.25\times 10^{3}$ (1.25 s), and (e) quasi-steady state, $t=3\times 10^{4}$ (16.7 s).

Figure 7. Transient heating of a water drop at $\hat{T}_{a}$ placed on a heated substrate at $\hat{T}_{w}$ . (a–e) $\hat{T}_{w}=40\,^{\circ }\text{C}$ with ${\rm\Delta}\hat{T}=14.6\,^{\circ }\text{C}$ : (a) initial conditions, $t=0$ , (b) $t=750$ (0.42 s), (c) $t=1.5\times 10^{3}$ (0.83 s), (d) $t=2.25\times 10^{3}$ (1.25 s), (e) quasi-steady state, $t=3\times 10^{4}$ (16.7 s). (f–j) $\hat{T}_{w}=70\,^{\circ }\text{C}$ with ${\rm\Delta}\hat{T}=39.8\,^{\circ }\text{C}$ : (f) initial conditions, $t=0$ , (g) $t=750$ (0.15 s), (h) $t=1.5\times 10^{3}$ (0.30 s), (i) $t=2.25\times 10^{3}$ (0.46 s), (j) quasi-steady state, $t=3\times 10^{4}$ (6.08 s).

4.2. Initial drop heating

To conclude this section on axisymmetric drops, we provide insights into such fundamental questions as how a water drop heats up when placed on a hot substrate. This initial stage has received very little attention in the literature. Note that Hu & Larson (Reference Hu and Larson2002, Reference Hu and Larson2005a ,Reference Hu and Larson b , Reference Hu and Larson2006), Girard et al. (Reference Girard, Antoni, Faure and Steinchen2006, Reference Girard, Antoni, Faure and Steinchen2008a ,Reference Girard, Antoni and Sefiane b ) and Girard & Antoni (Reference Girard and Antoni2008) described the transient drop evolution as a sum of steady-state solutions and therefore could not capture the dynamics of the initial stage.

The evolution in time of the two-phase flow velocity is presented in figure 6. During the initial transient phase, the liquid and gas velocities are comparable and significantly larger than those of the quasi-steady equilibrium; this is due to the severe thermal gradients resulting from considering an initially cold drop. As the flow converges to the equilibrium state, the gas velocity decreases more rapidly than that of the liquid and ultimately becomes significantly lower (figure 6 e). The coupled evolution of the thermal field is depicted in figure 7 for $\hat{T}_{w}=40$ and $70\,^{\circ }\text{C}$ to illustrate how the increasing role of the Marangoni effect influences the dynamics of the process. The strength of thermocapillarity is measured by the Marangoni number, $\mathit{Ma}=\hat{{\it\gamma}}\hat{R}_{0}{\rm\Delta}\hat{T}^{\ast }/(\hat{{\it\mu}}_{l}\hat{{\it\alpha}}_{l})$ , where $\hat{{\it\alpha}}_{l}$ is the liquid thermal diffusivity and ${\rm\Delta}\hat{T}^{\ast }$ is the temperature difference between the contact line and the apex. Calculation of this group shows that $\mathit{Ma}$ is 2.5 times larger for $\hat{T}_{w}=70\,^{\circ }\text{C}$ , i.e. the strength of thermocapillarity increases for increasing $\hat{T}_{w}$ . Figure 7(a, f) shows the initial conditions (cold drop) while figure 7(e,j) illustrates the quasi-steady $T$ field found after the initial abrupt heating. Between these states, the drop warms up very quickly due to the heat transported from the substrate, across the solid–liquid interface, and surrounding gas, across the liquid–gas interface (see figure 7 b,g). Some of the latter is absorbed by evaporation, which is strongest at the beginning given that ${\it\omega}_{t=0}=0$ , while the rest begins to increase the interface temperature. The temperature, $T$ , drops rapidly in the surrounding gas (its thermal diffusivity is two orders of magnitude larger than that of the liquid) and a thermal gradient arises along the interface triggering the Marangoni convection. Fluid from the contact-line region (at higher $T$ ) is driven along the interface to the apex, which is cooler due to evaporation, and then recirculated vertically downward towards the centre of the drop. At high $\hat{T}_{w}$ , the amount of heat associated with this Marangoni convection is significant and leads to the drop being heated through its core, as shown in figure 7(g–i). Figure 7(h,i) reveals that, rather than its centre, the coldest part of the drop is a toroidal-like region whose cross-section is approximately centred at $(r,\,z)=(0.8,\,0.3)$ . At high $\hat{T}_{w}$ , therefore, the drop reaches thermal equilibrium by means of interior heating due to the prominent role of advection. At lower $\hat{T}_{w}$ , on the other hand (see figure 7 b–d), the role played by the Marangoni effect is less significant and consequently the drop increases its temperature through its periphery via conduction. In this case the coldest region coincides with the drop’s core. Independently of $\hat{T}_{w}$ , the initial thermal adjustment takes place in ${<}3\,\%$ of the drop’s lifetime, after which the drop enters a quasi-steady phase wherein the flow characteristics evolve very slowly.

5.1. Elliptical contact area

Perfectly perpendicular dosing is key to generating spherical drops. However, in most real configurations the drop deposition is not exactly perpendicular to the substrate. A broad range of factors can be the cause of this, e.g. vertical dosing onto inclined substrates or drops released at a certain angle due to set-up configurations or even convection currents in the gas. In such cases, and assuming that the drop deposition is not violent, the resulting contact area is elliptical rather than circular. We therefore define the initial drop shape in Cartesian coordinates as follows:

where $a$ , $b$ , and $c$ are the semi-axis lengths along the $x$ , $y$ , and $z$ directions, respectively. The volume of a semi-ellipsoid $V_{e}$ is given by

while the base circumference $C_{e}$ can be estimated according to Ramanujan’s second approximation, i.e. 

where $h=(a-b)^{2}/(a+b)^{2}$ . It is known that the overall evaporation rate of a sessile drop is proportional to its volume (or contact angle) and to the length of its contact line, e.g. Birdi et al. (Reference Birdi, Vu and Winter1989), Rowan et al. (Reference Rowan, Newton and McHale1995) and Hu & Larson (Reference Hu and Larson2002). Thus, $a$ , $b$ and $c$ are selected so that the resulting ellipsoidal geometry provides the same volume and base perimeter as those for the spherical drop used as reference, i.e.  $V_{e}={\rm\pi}H(3R^{2}+H^{2})/6$ and $C_{e}=2{\rm\pi}R$ , where $H=0.98$ and $R=1.2$ . Equations (5.2) and (5.3) form a set of two relations with three unknowns. The missing expression necessary to complete an independent system is provided by the degree of deformation in the horizontal direction, ${\it\Lambda}=a/b$ , which is chosen to be ${\it\Lambda}=3/2$ . It follows that the resulting ellipsoidal geometry is defined by $a=1.4260$ , $b=0.9505$ , and $c=0.9815$ .

Our simulations show that an ellipsoid is not a possible shape for a sessile drop with an elliptical pinned contact line. The initial drop geometry (arbitrarily selected) undergoes a rapid change as surface tension enforces interface surface-area minimization for the given volume. Figure 8 illustrates the initial shape together with the final quasi-steady drop profile. The pinned contact line impedes the formation of a sphere but surface tension drives the interface shape to reach the closest possible configuration. In an attempt to even out the different axis lengths, there is interface pulling along the longest axis, $x$ (figure 8 a), while the interface expands in the other two directions, $y$ and $z$ (figure 8 b). The interface surface area of the resulting drop is smaller than the area of the ellipsoid but larger than the area of the spherical cap by 1.2 %, with the same volume and base perimeter. It should be noted that this initial shape adjustment is not caused by evaporation. The time interval in which this adjustment takes place is sufficiently short that the volume difference between the initial and the spherical drop shape is negligible $({\rm\Delta}V\approx 0.1\,\%)$ .

Figure 8. Initial semi-ellipsoidal interface (dashed line) and final drop shape (solid line) along the principal (a) $xz$ and (b) $yz$ planes for the drop with deformed contact area. This drop is defined by the semi-axis lengths of its elliptical contact area, $a=1.4260$ and $b=0.9505$ , and its volume, $V={\rm\pi}H(3R^{2}+H^{2})/6$ with $H=0.98$ and $R=1.2$ .

Figure 9. (a) Elliptical contact area of the deformed drop (solid line) along with that of the spherical drop (dashed line) employed for comparison. (b) Contact-angle ${\it\theta}$ and horizontal-curvature $1/R_{xy}$ distribution along the elliptical contact line. This is given in terms of the azimuthal angle ${\it\psi}$ defined as shown in (a). Here, ${\it\theta}_{max}=94^{\circ }$ and ${\it\theta}_{min}=59^{\circ }$ .

Figure 9(a) depicts how the resulting elliptical contact line of the deformed drop compares with the circular contact line of the reference drop, while figure 9(b) shows the contact-angle distribution along the triple line, which is not constant as in the spherical case. The position along the contact line is given in terms of the azimuthal angle, ${\it\psi}$ , defined as shown in figure 9(a). The contact angle varies between ${\it\theta}_{min}=59^{\circ }$ for the longest horizontal axis ( ${\it\psi}=0^{\circ }$ and $180^{\circ }$ ) and ${\it\theta}_{max}=94^{\circ }$ for the other ( ${\it\psi}=90^{\circ }$ and $270^{\circ }$ ). The average contact angle $\bar{{\it\theta}}=\int _{CL}{\it\theta}\,\text{d}l/\!\int _{CL}\,\text{d}l$ is $82^{\circ }$ . For a given $V$ , the range of contact angle variation throughout the base line, ${\rm\Delta}{\it\theta}={\it\theta}_{max}-{\it\theta}_{min}$ , increases for increasing deformation ${\it\Lambda}$ . Note that ${\rm\Delta}{\it\theta}=35^{\circ }$ with ${\it\Lambda}=3/2$ while ${\rm\Delta}{\it\theta}=0^{\circ }$ for ${\it\Lambda}=1$ (spherical cap). In addition, for a given ${\it\Lambda}$ , ${\rm\Delta}{\it\theta}$ decreases in time due to evaporation. For instance, when the drop has half of the original volume, i.e.  $V=V_{e}/2$ , we find that ${\rm\Delta}{\it\theta}=62^{\circ }-31^{\circ }=31^{\circ }$ and $\bar{{\it\theta}}=50^{\circ }$ . Note that in the limit $V\rightarrow 0$ then ${\rm\Delta}{\it\theta}\rightarrow 0$ . In other words, the conclusions are: (i) the more deformed the drop is the larger the range of ${\it\theta}$ , and (ii) evaporation makes the contact-angle distribution more homogeneous with time.

It is important to note that the contact-angle distribution is inversely related to the horizontal contact-line curvature $1/R_{xy}=|(\boldsymbol{{\rm\nabla}}_{H}\boldsymbol{\cdot }\boldsymbol{n})|_{z=0}$ , where $\boldsymbol{{\rm\nabla}}_{H}=\boldsymbol{{\rm\nabla}}-\boldsymbol{n}_{z}(\boldsymbol{n}_{z}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})$ is the horizontal gradient operator, and $\boldsymbol{n}$ and $\boldsymbol{n}_{z}$ are the unit vectors normal to the interface and to the lower boundary, respectively. Figure 9(b) illustrates the evolution of $1/R_{xy}$ between its maximum found at the ends of longest axis ( ${\it\psi}=0^{\circ }$ and $180^{\circ }$ ), where $1/R_{xy}=a/b^{2}=1.58$ , and minimum at the shortest ( ${\it\psi}=90^{\circ }$ and $270^{\circ }$ ), where $1/R_{xy}=b/a^{2}=0.47$ . The higher $1/R_{xy}$ at a certain point, the more this point behaves like a ‘corner’ and, as a result, experiences a greater degree of pulling, resulting in lower ${\it\theta}$ . To understand the physical mechanism behind this phenomenon, we need to examine the Young–Laplace equation, ${\rm\Delta}p={\it\sigma}(1/R_{1}+1/R_{2})$ . Since gravity has been neglected, the pressure in the liquid and gas are different but essentially constant and therefore ${\rm\Delta}p$ should be equal across any point along the interface. This imposes the constraint that if the contact-line curvature is large at one point ( $1/R_{1}$ ), the second principal curvature ( $1/R_{2}$ ), which dictates ${\it\theta}$ in this case, has to be smaller (larger ${\it\theta}$ ) in order to maintain a constant ${\rm\Delta}p$ . Note that, thanks to the drop symmetry, the principal radii of curvature of the interface at the extreme points A and B coincide with radii of curvature of the contact line and the interface contained in the principal vertical plane across the point (directly connected to ${\it\theta}$ ). This singularity may not occur for the intermediate points, a consideration which is revisited in the more complex case presented in § 5.2.

As was mentioned above, the evaporation rate $m(t)$ in spherical drops is known to be dependent on the length of the base circumference $C_{e}$ and volume $V$ , which is similar to saying that $m(t)$ depends on $C_{e}$ and ${\it\theta}$ . It is evident that the spherical condition together with $C_{e}$ and $V$ lead to a unique interface surface area and therefore to a characteristic $m(t)$ . Not surprisingly, this consideration does not hold for deformed drops given that the interface area is not uniquely defined with $C_{e}$ and $V$ only. It is found that $m(t)$ is 2.3 % larger than that for the reference spherical drop even though $C_{e}$ and $V$ are the same in both cases, which is in accordance with the order of surface area increment reported above. This finding is also consistent with the experimental observations by Jansen et al. (Reference Jansen, Zandvliet and Kooij2014), who recently reported higher evaporation rates for more elongated droplets.

Figure 10 depicts $T$ along the interface. The temperature ranges from $T=1$ at the triple line to $T=0.77$ at the apex. Given the elliptical contact area, one might have expected $T$ to follow rather closely a distribution of concentric ellipses of decreasing size for increasing $z$ . This is not what is observed in the simulations, however. The irregular drop shape leads to the thermal field being deformed in the $z$ direction, as shown in figure 10(b). Following the interface shape, the isotherms are stretched downwards along the longest axis. Hence, their vertical location varies depending on the azimuthal position. The curvature and path length towards the apex are also different at each point. This combination of factors gives rise to azimuthal thermal gradients along the interface.

Figure 10. Interface temperature for the drop with elliptical contact area. The solid black lines are contours of constant $z$ coordinate. (a) Front and (b) side views.

Figure 11 shows the three-dimensional (3D) flow within the drop by means of streamlines coloured by the local $T$ . The flow repeats itself periodically at each quadrant even though the data are from a three-dimensional simulation of the complete drop. Flow is driven from the contact line to the apex along the interface and then recirculated through the drop core back towards the contact line. However, the streamlines along the interface are bent towards the ends of the longest axis, as shown in figure 11(b), due to the azimuthal dependence of $\boldsymbol{{\rm\nabla}}_{~s}T$ ; the dashed arrows show the direction of the emerging azimuthal velocity components. This effect also induces azimuthal currents in the interior of the drop, which have been isolated for the first quadrant in figure 11 for clarity of presentation. At each quarter, the current’s origin coincides with the centre of the vortex at the vertical plane across the shortest axis of the ellipse, i.e. the $yz$ plane (figure 11 d), and flows into the centre of the second fundamental vortex generated along the longest axis, i.e. the $xz$ plane (figure 11 c). Each stream, therefore, is a bent toroidal-like feature, transporting liquid from areas located in perpendicular planes and at different heights.

Figure 11. Three-dimensional flow within the drop with elliptical contact area defined in figure 9. The interface and streamlines are coloured by temperature. (a) 3D, (b) top, (c) front and (d) side views. The solid and dashed arrows indicate the direction of the primary interface flow and secondary azimuthal currents, respectively.

In order to characterize the torroidal currents more quantitatively, the local volumetric helicity $\mathscr{H}=\boldsymbol{u}\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}\times \boldsymbol{u})$ is presented in figure 12. Note that $\mathscr{H}$ is large when the flow velocity $\boldsymbol{u}$ and the vorticity $(\boldsymbol{{\rm\nabla}}\times \boldsymbol{u})$ are not only large but also parallel, which makes $\mathscr{H}$ a very useful variable to provide a measure of the current’s azimuthal strength. At the vertical planes $\mathscr{H}=0$ because the $\boldsymbol{u}$ components perpendicular to the planes are 0. Note that for each plane the perpendicular velocity on one side is in the opposite direction to the perpendicular velocity from the other side. For instance, at the $yz$ plane, $(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}_{x})_{x\rightarrow 0}^{x<0}=-(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}_{x})_{x\rightarrow 0}^{x>0}$ , and therefore $(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}_{x})_{x=0}=0$ . Similarly for the $xz$ plane, $(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}_{y})_{y=0}=0$ given that $(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}_{y})_{y\rightarrow 0}^{y<0}=-(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}_{y})_{y\rightarrow 0}^{y>0}$ . As the fluid departs from the $yz$ plane, there is an increment in the flow’s $\mathscr{H}$ until it reaches a maximum and then decays to become 0 again at the $xz$ plane. Thus, a particle following this path will probably accelerate during the first stage and slow down at the end. The speed of these currents is approximately one order of magnitude lower than the maximum flow speed, which is found at the interface, and therefore comparable to the speed of the flow in the drop’s core.

Figure 12. Local volumetric helicity $\mathscr{H}=\boldsymbol{u}\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}\times \boldsymbol{u})$ for the azimuthal current emerging in the deformed drop with elliptical contact area. The arrow shows the direction of motion.

5.2. Irregular contact area

We now study the case wherein an evaporating spherical drop becomes deformed after the first contact-line jump, which is typically observed in the experiments at $t=0.7t_{tot}$ . In the range of temperatures investigated in the laboratory, the contact-line motion begins when a certain critical disequilibrium point is reached in the dynamics. The critical conditions change slightly with $\hat{T}_{w}$ but, in the range of temperatures examined in the laboratory $(40^{\circ }\leqslant \hat{T}_{w}\leqslant 70^{\circ })$ , these fall into the range $22^{\circ }\leqslant {\it\theta}_{c}\leqslant 26^{\circ }$ , which in terms of drop height becomes $0.37~\text{mm}\leqslant {\hat{H}}_{c}\leqslant 0.44~\text{mm}$ (figure 3), and $0.24\leqslant H_{c}\leqslant 0.28$ in dimensional and dimensionless forms, respectively. For the sake of simplicity, we round $H_{c}$ and consider that the drop depinning occurs when $H=0.3$ .

The experiments show that during the three-dimensional contact-area motion, a section of the triple line remains anchored to its original position (top-left quadrant in figure 2 d) while the rest suddenly moves to regain equilibrium by reducing the contact area. In our simulations, we generate the deformed drop in a similar way (see figure 13 a). Starting from an axisymmetric flow with $H=0.3$ and $R=1.2$ in three dimensions (dashed line in figure 13 a), the contact line is allowed to move freely everywhere except in the second quadrant, $90^{\circ }\leqslant {\it\psi}\leqslant 180^{\circ }$ , where it is pinned. Consequently, surface tension brings the MCL towards its new pinned position (solid line in figure 13 a) to satisfy interface area minimization. The new pinned position has been selected to approximate experimental observations, and the model implemented to allow contact-line motion is described in § 6. To maintain the focus of the investigation, this is not presented here. Note that the focus here is on elucidating the dynamics of evaporation from a deformed drop; thus modelling the complex motion of the contact line that brings about this deformation is beyond the scope of this work.

Figure 13. (a) Irregular contact area of the deformed drop found after the first contact jump (solid line) compared with that of its precursor spherical drop (dashed line) with $H=0.3$ and $R=1.2$ . (b) Contact-angle ${\it\theta}$ and horizontal-curvature $1/R_{xy}$ distribution along the contact line. This is given in terms of the azimuthal angle ${\it\psi}$ defined as shown in (a). Here, ${\it\theta}_{max}=64^{\circ }$ and ${\it\theta}_{min}=51^{\circ }$ . The key points are A $({\it\psi}_{A}=82^{\circ })$ , B $({\it\psi}_{B}=135^{\circ })$ , C $({\it\psi}_{C}=188^{\circ })$ and D $({\it\psi}_{D}=315^{\circ })$ .

Contrary to the previous drop geometry with elliptical contact area, this more irregular drop only presents one vertical plane of symmetry, namely that across points B and D in figure 13(a). The resulting ${\it\theta}$ distribution, shown in figure 13(b), is more intricate but follows similar principles. During the readjustment, ${\it\theta}$ increases from $28^{\circ }$ (spherical drop) to $51^{\circ }\leqslant {\it\theta}\leqslant 64^{\circ }$ , which is in excellent agreement with the range recorded experimentally, as shown in figure 3(a). Once again, ${\it\theta}$ is related to the contact-line curvature $1/R_{xy}$ as the points of minimum ${\it\theta}$ coincide fairly well with the locations where $1/R_{xy}$ is maximum, i.e. points A and C, and vice versa. Note that ${\it\theta}$ and $1/R_{xy}$ are perfectly in phase at B and D whereas there is a small mismatch for A and C. The reason lies in the drop symmetry. The angles ${\it\theta}_{B}$ and ${\it\theta}_{D}$ result from the equilibrium of interfacial forces acting on each side of the line BD, which are symmetric. In the case of A and C, however, the interface forces on each side of the line that connects these points are different due to the drop shape. Note that $R_{xy}$ corresponds to one of the principal interface radii of curvature for B and D but not necessarily for A and C.

The depinning of the drop occurs when it is sufficiently out of equilibrium, that is, when its interface surface area is very large in comparison with its volume. We find that for the same $V$ , the geometric readjustment leads to an interface whose area is significantly smaller (30.4 %) than the interface area of the spherical drop found right before the depinning event. In a similar manner, the evaporation rate $m(t)$ decreases by the same order of magnitude: $m(t)$ is 20.3 % lower after the contact-area jump. The reasons why there is a relative mismatch between the change of surface area and evaporation rate appear to be associated with the rapid transient nature of the flow. Careful analysis of the instantaneous $m(t)$ reveals that, during the transient stage, $m(t)$ in the irregular drop decreases at a rate $O(1)$ faster than that in the spherical drop used as reference. This indicates that the deviation between the change in area and evaporation rate (10.1 %) is being corrected with time and therefore is probably a short-term consequence of the complex transient flow. It should be noted that a complete correction is unlikely given that the contact-line motion also entails a significant redistribution of drop’s thermal field, along with a change in the strength of the Marangoni effect. Hence, a residual mismatch is expected but, once again, the surface area appears to be the leading factor in dictating changes in $m(t)$ .

Figure 14. Three-dimensional flow and temperature distribution in the irregular drop found after the first contact-line jump defined in figure 13. (a) 3D, (b) top, (c) side and (d) front views. The arrows illustrate the direction of motion.

During the contact-line jump, the interface motion severely deforms the preceding axisymmetric flow. This loses its quasi-steady condition, triggering a rapid transient stage until another equilibrium state is achieved for the new geometric configuration. The three-dimensional streamlines and $T$ distribution within the drop during this transitional phase are illustrated in figure 14. The Marangoni force develops a preferential direction, line BD in figure 13(a), along which warm liquid is transported from D towards B, and then recirculated downward and back along the periphery via the points A and C. Consequently, a pair of counter-rotating vortices emerge in the bulk flow: see figure 14. These are centred in the colder regions and are essentially parallel to the horizontal plane. The strength of these vortices diminishes with time as the interface temperature is homogenized by the central stream of warm fluid as well as by the heat being transported towards the cold vortices from the contact line due to both conduction and convection. Eventually the vortices disappear, along with the preferential direction of the Marangoni convection, and the drop reaches its new quasi-steady equilibrium wherein the apex is the coldest region and the interface flow is mainly radial towards the apex. Azimuthal components also emerge in this new quasi-steady equilibrium, but these are weaker than those reported for the drop with an elliptical contact area given that the interface geometry is less deformed. The interface temperature distribution during the intermediate transient stage is depicted in figure 15(a). Capturing these vortices in the laboratory is very difficult due to the relatively small range of the $T$ and because the transient stage after depinning for our chosen experiment lasts a very short period of time ( ${\sim}0.8~\text{s}$ ). Note that the drop volume at this point is quite small and thus any imbalance is quickly corrected by the substrate, whose dominance over the drop’s thermal field increases for decreasing $V$ . However, similar vortices are clearly observed during the initial transient stage of irregularly shaped colder drops placed on the heated substrate: see figure 15(b). By reducing the initial liquid temperature, the vortices emerge more vigorously and persist for longer times, which makes it easier to record them in the laboratory. The qualitative agreement between the numerical simulations (figure 15 a) and the experiments (figure 15 b) is excellent. Given the importance of these vortices in regulating the initial dynamics of evaporating drops, a follow-up investigation is being undertaken, the focus of which is on the detailed numerical and experimental study of these and similar three-dimensional features emerging in non-spherical drops with more diverse and complex contact areas.

Figure 15. Interface temperature showing the pair of counter-rotating vortices emerging during the rapid thermal adjustment in an irregular drop: (a) simulations, (b) IR image showing experimentally the emergence of counter-rotating vortices in evaporating drops with $\hat{T}_{0}=10\,^{\circ }\text{C}$ and $\hat{T}_{w}=70\,^{\circ }\text{C}$ .