Evaporation of binary liquids from a capillary tube

Abstract

Evaporation of multi-component liquid mixtures in confined geometries, such as capillaries, is crucial in applications such as microfluidics, two-phase cooling devices and inkjet printing. Predicting the behaviour of such systems becomes challenging because evaporation triggers complex spatio-temporal changes in the composition of the mixture. These changes in composition, in turn, affect evaporation. In the present work, we study the evaporation of aqueous glycerol solutions contained as a liquid column in a capillary tube. Experiments and direct numerical simulations show three evaporation regimes characterised by different temporal evolutions of the normalised mass transfer rate (or Sherwood number $Sh$), namely $Sh (\tilde{t} ) = 1$, $Sh \sim 1/\sqrt {\tilde{t} }$ and $Sh \sim \exp (-\tilde{t} )$, where $\tilde {t}$ is a normalised time. We present a simplistic analytical model that shows that the evaporation dynamics can be expressed by the classical relation $Sh = \exp ( \tilde{t} )\,\mathrm {erfc} ( \sqrt {\tilde{t} })$. For small and medium $\tilde{t}$, this expression results in the first and second of the three observed scaling regimes, respectively. This analytical model is formulated in the limit of pure diffusion and when the penetration depth $\delta (t)$ of the diffusion front is much smaller than the length $L(t)$ of the liquid column. When $\delta \approx L$, finite-length effects lead to $Sh \sim \exp (-\tilde{t} )$, i.e. the third regime. Finally, we extend our analytical model to incorporate the effect of advection and determine the conditions under which this effect is important. Our results provide fundamental insights into the physics of selective evaporation from a multi-component liquid column.

1. Introduction

Evaporation of multi-component volatile liquids into a gaseous phase is ubiquitous in both biological and industrial processes (Erbil 2012; Lohse & Zhang 2020; Bourouiba 2021; Morris et al. 2021; Lohse 2022). A multi-component liquid can consist of multiple solvents, surfactants, polymers, colloids and salts. Evaporation from such systems leads to a plethora of phenomena, such as instabilities (Li et al. 2020), phase separation (Tan et al. 2016), altering of deposition patterns (Palacios et al. 2012), crystallisation (Shahidzadeh-Bonn et al. 2008), stratification (Hooiveld et al. 2023) and evaporation-driven flows (Deegan et al. 1997; Hu & Larson 2006). In addition to the composition of the liquid, geometrical confinement also affects evaporation significantly. The geometrical confinement can be in the form of a droplet (Picknett & Bexon 1977; Lohse & Zhang 2020), a liquid film (Okazaki et al. 1974), a porous membrane (Prat 2002), a shallow pit (D'Ambrosio et al. 2021) or a capillary (Chauvet et al. 2009).

Understanding the evaporation of liquids from capillaries is crucial for applications such as microfluidics (Zimmermann et al. 2005; Bacchin, Leng & Salmon 2022), inkjet printing (Lohse 2022; Rump et al. 2023), heat pipes (Chen et al. 2016), chromatography (Kamp et al. 2005) and the measurement of material properties (Roger, Sparr & Wennerström 2018; Merhi et al. 2022; Nguyen, Bouchaudy & Salmon 2022). Capillaries are also considered to be idealised systems for modelling porous structures (Yiotis et al. 2007; Chauvet et al. 2009; Chen et al. 2022; Le Dizès Castell et al. 2023), the transport of water across skin (Sparr & Wennerström 2000; Roger et al. 2016) and film drying (Guerrier et al. 1998; Salmon, Doumenc & Guerrier 2017).

Capillary evaporation is determined mainly by the behaviour of the volatile liquid meniscus. There have been several experimental and numerical studies to determine the evaporation from a liquid meniscus. These studies describe the evaporation rate (Wayner & Coccio 1971), heat transfer coefficients (Wayner, Kao & LaCroix 1976; Park & Lee 2003; Dhavaleswarapu, Murthy & Garimella 2012; Zhou et al. 2018), shape of the meniscus (Potash & Wayner 1972; Swanson & Herdt 1992), capillary flows that replenish the evaporated liquid (Potash & Wayner 1972; Ransohoff & Radke 1988; Park & Lee 2003) and additional flows that might be driven by evaporation-induced surface tension gradients (Schmidt & Chung 1992; Buffone & Sefiane 2003; Dhavaleswarapu et al. 2007; Cecere, Buffone & Savino 2014) or buoyancy (Dhavaleswarapu et al. 2007; Buffone 2019).

Broadly speaking, the evaporation of a single-component liquid from a capillary can be divided into two main classes depending on the location of the liquid–air meniscus with respect to the open end of the capillary (henceforth referred to as its ‘mouth’; see figure 1a). In the first class of problems, the liquid–air interface is far away from the mouth of the capillary. In such a configuration, the evaporation rate is limited by the transport of vapour from the liquid–air interface to the mouth of the capillary tube, and it varies approximately as the inverse square root of time (Stefan 1873, 1889).

Figure 1. Schematic of the different configurations of evaporation from a capillary: (a) first class of problems, where the liquid meniscus is away from the mouth of the capillary $(l\gg 0)$; (b,c) second class of problems, where the liquid meniscus is close to the mouth of the capillary $(l \approx 0)$, with (b) $\theta >90^\circ$ and (c) $\theta < 90^\circ$. Red arrows represent the evaporative flux.

In the second class of problems, the liquid meniscus is at (or relatively close to) the mouth of the capillary (Chauvet et al. 2009; Gazzola, Franchi Scarselli & Guerrieri 2009). Within this second class of problems, if the contact angle inside the liquid between the liquid–gas interface and the capillary wall is $\theta \geq 90^\circ$, then one can realise immediately the resemblance to a sessile droplet (figure 1b). In such a case, one can use the Popov model (Popov 2005; Li et al. 2019) to predict the evaporation rate. For $130^\circ > \theta \geq 90^\circ$ (which is equivalent to a contact angle of $40^\circ > \theta _{{drop}} \geq 0^\circ$ in the case of a sessile droplet), the evaporation rates will be practically independent of the contact angle (Sobac & Brutin 2011) and depend mainly on the base radius, ambient humidity and properties of the liquid.

However, when $\theta \leq 90^\circ$, the droplet model for evaporation is no longer applicable (figure 1c). The evaporation under such conditions shows two distinct regimes – a ‘constant rate period’ and a ‘falling rate period’ (Chauvet et al. 2009; Keita et al. 2014, 2016) – similar to that observed during the drying of a porous medium (Coussot 2000). During the constant rate period, evaporation is still determined mainly by the ambient conditions. To replenish the liquid lost by evaporation, upstream liquid is driven by capillary pressure to the mouth of the capillary through thin liquid films (Ransohoff & Radke 1988; Chauvet, Duru & Prat 2010). However, even during the so-called constant rate period, the evaporation rate decreases slightly (Coussot 2000; Chauvet et al. 2009). For square capillaries, this slight decrease is due to the thinning of liquid films at the mouth of the capillary (Chauvet et al. 2009). As evaporation proceeds further, the depinning of the liquid films from the mouth of the capillary leads to the falling rate period (Chauvet et al. 2009). When the liquid meniscus is sufficiently far away from the mouth of the capillary, the situation reverts to the first class of problems (figure 1a).

In addition to the studies on evaporation of single-component liquids from capillaries, there have also been several studies on evaporation in multi-component systems, such as evaporation of binary liquid mixtures (Duursma, Sefiane & Clarke 2008; Cecere et al. 2014; Salmon et al. 2017; Zhou et al. 2018), surfactant solutions (He et al. 2015; Roger et al. 2016, 2018), salt solutions (Camassel et al. 2005; Naillon et al. 2015) and colloidal dispersions (Abkarian, Nunes & Stone 2004; Kamp et al. 2005; Sarkar & Tirumkudulu 2009; Wang et al. 2020; Roger & Crassous 2021) in capillaries. In all these scenarios, evaporation leads to spatio-temporal variations in the composition of the mixture. Nevertheless, it is generally possible to model the capillary evaporation of a multi-component system as a one-dimensional transport problem. It is unsurprising, but perhaps interesting to note, that in the absence of instabilities (de Gennes 2001), the one-dimensional model of evaporation of polymeric liquid films (Okazaki et al. 1974; Guerrier et al. 1998; Saure, Wagner & Schlünder 1998; Okuzono, Ozawa & Doi 2006) is mathematically equivalent to the evaporation of binary solutions from capillaries (Salmon et al. 2017).

In multi-component systems, especially binary systems, the evaporation rate can also show a constant rate period followed by a falling rate period (Okazaki et al. 1974; Salmon et al. 2017; Huisman et al. 2023), similar to pure liquids. However, in multi-component systems, the different evaporation regimes are determined additionally by changes in the composition of the mixture. Hence a complete evaporation model must include the spatio-temporal variations in the composition and properties of the mixture. Recently, Salmon et al. (2017) showed how a steep decrease in thermodynamic activity and the diffusion coefficient of water at high solute concentrations can lead to evaporation being almost independent of ambient humidity for certain molecularly complex fluids. These authors modelled the variable diffusion coefficient as a piecewise-constant function, but bypassed the necessity of using an analytical expression for the thermodynamic activity of water. Moreover, Salmon et al. (2017) considered the parameter range where the medium can be approximated as semi-infinite.

In this work, we study the evaporation of binary liquids in capillaries with experiments, direct numerical simulations and analytical modelling. We perform experiments for the evaporation of water–glycerol mixtures in a cylindrical capillary tube under controlled humidity conditions. Our direct numerical simulations show excellent agreement with the experiments. Further, to unravel the physics of the evaporation dynamics, we develop a one-dimensional analytical model. We introduce a linear approximation for the thermodynamic activity of water as a function of its weight fraction. We also identify the condition under which the semi-infinite assumption breaks down and accordingly take it into account in our modelling. Finally, we discuss the explicit role of the advective mass transport compared with the diffusive mass transfer for our particular system. We show that our model predicts accurately the relevant scaling laws observed in the experiments and the numerical simulations.

The paper is organised as follows. In § 2, we describe the experimental set-up and observations. The governing equations of our system and the numerical method are described in § 3. In § 4, we provide three simplified analytical models of the problem, each with an added level of complexity over its predecessor, and compare their predictions with the results obtained from the experiments and the simulations. The paper culminates in a summary of the results and an outlook in § 5.

2. Experiments

2.1. Experimental set-up

Aqueous solutions of different mass fractions of glycerol (Sigma-Aldrich) were used as the probe liquids in the present experiments. The use of glycerol has the following advantages. First, since glycerol has a very low vapour pressure, it is practically non-volatile under the current experimental conditions. Thus we need only account for the evaporation of water. Further, since the more volatile liquid (water) has higher surface tension, its evaporation should not lead to any flow instabilities close to the interface (Diddens 2017). Such instabilities occur whenever the mass transfer (evaporation or condensation) leads to an increase in surface tension, such as for evaporating water–ethanol mixtures (Christy, Hamamoto & Sefiane 2011; Bennacer & Sefiane 2014; Diddens et al. 2017; Lopez de la Cruz et al. 2021) or condensation of water onto a water–glycerol droplet (Shin, Jacobi & Stone 2016; Diddens 2017). A more detailed explanation of such instabilities can be found in Diddens (2017).

In the present experiments, we study the evaporation dynamics of aqueous solutions of glycerol in a thin cylindrical capillary tube (figure 2; inner diameter 1 mm, outer diameter 1.2 mm, length 100 mm; Round Boro Tubing, CM Scientific). The liquid column inside the capillary had an initial height of $19 \pm 2\,{\rm mm}$. The initial weight fraction of water, $w_i$, in the water–glycerol mixture was varied as 0.2, 0.6, 0.9 and 1.0, to cover a wide range of initial compositions. The lower end of the capillary was placed inside an in-house-developed, optically transparent, humidity-controlled chamber at room temperature. The humidity and temperature inside the chamber were monitored using a temperature–humidity sensor (HIH6121, Honeywell). The relative humidity $(H_r)$ in the chamber was maintained at $H_r = 10 \pm 5\,\%$. The upper end of the capillary tube was exposed to room humidity (${>}50\,\%$). The evaporation or condensation of water at the upper meniscus was negligible compared to the evaporation from the bottom. This is because of the relatively large distance between the liquid's upper meniscus and the capillary's upper end (see Appendix A for a detailed discussion).

Figure 2. (a) Schematic of the experimental set-up. (b) Geometry used for numerical simulation and analytical modelling. (c) Left: schematic of the relative length of the liquid column with respect to the length of the capillary. Right: time-lapsed experimental snapshots of water–glycerol mixtures in the capillaries for an initial weight fraction of water $w_{i} = 0.9$. The top interface of the liquid mixture keeps moving downwards due to the evaporation of water from the bottom of the capillary. Red arrows in (a,b) represent the evaporative flux.

The contact line of the lower meniscus remained pinned at the lower mouth of the capillary. Thus the loss of water by evaporation from the lower mouth of the capillary leads to a decrease in the length $L$ of the liquid column (figure 2c). To study this evaporation process quantitatively, time-lapsed images of the liquid column were captured using a DSLR camera (D750, Nikon) equipped with either a long-distance microscope (Navitar 12$\times$) or a macro lens (50 mm DG Macro D, Sigma), while the capillary tube was back-illuminated with a cold LED light source (Thorlabs). For pure water, the velocity $v_y$ of the upper interface,

(2.1)\begin{equation} v_y = \frac{\mathrm{d} y}{\mathrm{d}t}= \frac{1}{\rho_w}\,\frac{\mathrm{d}M^{\prime \prime}}{\mathrm{d}t} , \end{equation}

is a direct measure of the evaporation rate $\mathrm {d}M^{\prime \prime }/\mathrm {d}t$ of water, where $y$ is the displacement of the top interface, $M^{\prime \prime }$ is the mass per unit area and $\rho _w$ is the density of water.

At a later time, the contact line of the lower liquid meniscus eventually depins. At this point, we stop the measurements because $v_y$ is thereafter no longer a correct measure of the evaporation rate. Additionally, as the lower interface moves inwards into the capillary tube after depinning, the evaporation boundary condition at the lower interface also changes (see § 3.1 for details of the boundary conditions).

2.2. Experimental results

The discrete data points in figure 3(a) denote the experimentally obtained vertical displacement $y$ of the upper interface with time $t$ for different initial weight fractions of water, $w_i$. Only every fifth data point is plotted in figure 3(a) to avoid overcrowding of the plot. The markers represent the mean of at least three independent experimental realisations, while the error bars (denoting $\pm$ one standard deviation) reflect the uncertainty due to experimental variabilities and the resolution of the imaging. The velocity $v_y = \mathrm {d} y/ \mathrm {d}t$ of the top interface, calculated from $y$, is plotted as discrete data points in figure 3(b).

Figure 3. (a) Displacement $y$ and (b) velocity $v_y$ of the top interface of the aqueous solutions of glycerol observed in experiments (discrete data points) and numerical simulations (continuous lines), for initial weight fractions of water $w_{i} = 0.2, 0.6, 0.9\text { and }1.0$.

Obviously, the motion of the top interface depends on the initial composition of the liquid mixture (see figure 3 and supplementary movie 1 available at https://doi.org/10.1017/jfm.2024.122). For pure water ($w_i=1$), the interface moves at an almost constant velocity ($y=v_y t$ and ${v_y(t)} \approx 10^{-3}\,{\rm mm}\,{\rm s}^{-1}$; see figure 3). For $w_i = 0.2$ and $0.6$, the experiments suggest $y \sim t^{1/2}$ and $v_y \sim t^{-1/2}$ scaling relations. However, for $w_i=0.9$, the experiments show three different regimes: ${v_y(t)} \approx 0.8 \times 10^{-3}\,{\rm mm}\,{\rm s}^{-1}$, similar to $w_i = 1$, in approximately the first $800$ s; $v_y \sim t^{-1/2}$, similar to $w_i=0.2$ and $0.6$, at intermediate times ($t \lesssim 5 \times 10^4$ s); and a decrease in velocity (which is steeper than $v_y \sim t^{-1/2}$) at long times.

The addition of glycerol reduces primarily the local concentration of water at the lower interface, which in turn leads to a reduced evaporation rate and a decrease in the velocity of the upper interface. Thus we get the so-called constant rate period and the so-called falling rate period as described in the literature (Salmon et al. 2017). Our experimental results, however, raise two important questions. First, why does the constant rate period not appear for $w_i=0.2$ and $0.6$ in the experiments? Second, why does the falling rate period show two sub-regimes for $w_i=0.9$? To answer these questions, and to understand how the evaporation dynamics of a water–glycerol mixture in a capillary changes with time and the initial composition of the mixture, we develop a theoretical model in the next section.

3. Problem formulation

3.1. Theoretical model

We model the column of water–glycerol mixture as a one-dimensional isothermal system with length $L(t)$ (figure 2b). The lower end of the capillary is located at $z=0$. We can express the liquid composition using the weight fraction of water $w(z,t)$; the weight fraction of glycerol then is simply $1- w(z,t)$. The spatio-temporal variations in the concentration of water in the liquid column can be determined by solving the one-dimensional continuity and advection–diffusion equations:

(3.1)\begin{equation} \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial z} (\rho u) = 0 \end{equation}

and

(3.2)\begin{equation} \frac{\partial}{\partial t} (\rho w) + u\,\frac{\partial}{\partial z} (\rho w) = \frac{\partial}{\partial z} \left( \rho D\,\frac{\partial w}{\partial z} \right) , \end{equation}

where $\rho (w)$ is the local density of the mixture, $D(w)$ is the diffusion coefficient of the water–glycerol mixture and $u$ is the fluid velocity along the $z$-direction.

Initially (at $t=0$), the composition in the liquid column is uniform and equal to $w_i$. Since the evaporation of water from the upper interface is negligible (see Appendix A), we model the upper interface as non-evaporative. The lower interface of the liquid column is exposed to the ambient constant relative humidity $H_{r}$ and loses water by evaporation (per unit area) at rate $\mathrm {d}M^{\prime \prime }/\mathrm {d}t$. The water lost due to evaporation is replenished by the diffusive and advective transport of water from the bulk liquid inside the capillary. These considerations lead to one initial and two boundary conditions, as follows:

(3.3)\begin{equation} \left. \begin{gathered} \displaystyle w = w_i\quad \text{at}\ t=0 ,\\ \displaystyle\frac{\partial w}{\partial z} = 0\quad \text{at}\ z = L(t) , \\ \displaystyle -\rho u w + \rho D\,\frac{\partial w}{\partial z} = \frac{\mathrm{d}M^{\prime \prime}}{\mathrm{d}t} \quad \text{at}\ z=0 . \end{gathered} \right\} \end{equation}

The evaporation of water from the bottom interface can be approximated using a quasi-steady diffusion-limited model of evaporation of a pinned sessile droplet having zero contact angle (Popov 2005; Stauber et al. 2014):

(3.4)\begin{equation} \frac{\mathrm{d}M^{\prime \prime}}{\mathrm{d}t} = {\rm \pi}D_{v,a}\,R (c_{w,s}-c_{w,\infty})\, f(\theta_{{drop}})\,\frac{1}{{\rm \pi} R^2} , \end{equation}

where $D_{v,a}$ is the diffusion coefficient of water vapour in air, $R$ is the inner radius of the capillary tube, $c_{w,s}$ is the concentration (vapour mass per volume) of water vapour in the air at the lower interface, $c_{w,\infty }=H_r c_{w,s}^o$ is the concentration of water vapour in the ambient air (far away from the capillary), $c_{w,s}^o$ is the saturation concentration of water vapour at the surface of pure water, $\theta _{{drop}}$ is the angle between the horizontal plane and the liquid–air interface at the contact line, and $f(\theta _{{drop}})$ is a known function of $\theta _{{drop}}$. For $\theta _{{drop}}=0$ as in our model here, $f(\theta _{{drop}}) = 4/{\rm \pi}$ (Popov 2005; Stauber et al. 2014). Thus

(3.5)\begin{equation} \frac{\mathrm{d}M^{\prime \prime}}{\mathrm{d}t} =\frac{4 D_{v,a}}{{\rm \pi} R}\,(c_{w,s}-c_{w,\infty}) . \end{equation}

Here, $c_{w,s}$ depends on the composition of the liquid mixture close to the interface (at $z=0^+$) and is given by Raoult's law as

(3.6)\begin{equation} c_{w,s} = a c_{w,s}^o = x_{o} \psi_{o} c_{w,s}^o , \end{equation}

where $a(w)$ is the thermodynamic activity of water, $x_{o}$ is the mole fraction of water in the liquid at $z=0^+$ and $\psi _o(x_o)$ is the activity coefficient of water corresponding to $x_o$. As long as $c_{w,s}>c_{w,\infty }$, water will evaporate from the lower interface.

3.2. Numerical solution

We solve numerically the equations pertaining to the aforementioned one-dimensional theoretical model using finite element method simulations with initial length $L=20\,\textrm {mm}$ to obtain the evaporation rates and the spatio-temporal distribution of the concentration of water, $w$, in the capillary. To that end, a line mesh initially consisting of 100 second-order Lagrangian elements is created to cover the initial length $L$. The motion of the top interface is realised by moving the mesh nodes along with the top interface, i.e. by an arbitrary Lagrangian–Eulerian (ALE) method with a Laplace-smoothed mesh. The mesh displacement at the top interface $z=L(t)$, i.e. $\dot {L}(t)=u(L,t)$, is enforced by a Lagrange multiplier acting on the node position at the top. Likewise, the evaporation at the bottom is considered, but here the Lagrange multiplier is acting on the velocity $u$ at $z=0$, whereas the bottom node remains fixed at $z=0$.

The implementation of (3.2) and (3.1) along with the boundary conditions (3.3) and (3.5) is achieved by the conventional weak formulations of advection–diffusion equations, including the ALE corrections for the time derivatives. A posteriori spatial adaptivity based on the jumps in the slopes of $w$ across the elements is considered. Also, a posteriori temporal adaptivity is considered by calculating the difference between the freshly calculated value of each field at each point and its prediction. For the prediction, values from the previous time steps are extrapolated to the current time. If the difference is large, then this means that the system changes excessively during a time step. In that case, the current time step is rejected and calculations are done again with a smaller time step. The implementation has been done with the finite element library oomph-lib by Heil & Hazel (2006), which solves monolithically the coupled equations with a backward differentiation formula of second order for the temporal integration.

The variation of the diffusion coefficient $D$ as a function of the local composition $w$ is considered in the simulations based on the experimental data of D'Errico et al. (2004), while the mass density $\rho$ was fitted according to the data of Takamura, Fischer & Morrow (2012). The activity coefficient of water was calculated by AIOMFAC (Zuend et al. 2011).

The results of the direct numerical simulations are shown by the continuous lines in figure 3. Excellent agreement between the experiments and the numerical simulations is observed. In particular, figure 3(b) shows that the simulations can reproduce the experimentally observed $v_y \sim t^{-1/2}$ scaling for $w_i=0.2$ and $0.6$, and all three velocity scalings for $w_i=0.9$. Interestingly, the simulations also show that for very early times, $v_y$ is almost constant for $w_i=0.2$ and $0.6$ as well. However, we cannot access these time scales in experiments because of limitations arising from the lack of spatio-temporal resolution. Overall, the quantitative match between experiments and numerical simulations shows that our theoretical model incorporates all the relevant physics of the problem. In figure 10 of Appendix B, we also show the spatial variation in the axial concentration profiles for the different initial concentrations $w_i$. In the next section, we will use additional simplifying assumptions to formulate a simplistic model that captures the essential physics of the system and recovers the various evaporation regimes.

4. Analytical model

We present here a simplified description of the problem with the purpose of elucidating the physical mechanisms behind the different regimes observed in the experiments and the direct numerical simulations. We introduce some assumptions that will allow us to treat the resulting problem analytically. As we will show below, despite these simplifications, the quantitative comparison between the model and the experiments and simulations is reasonably good. Our one-dimensional analytical model relies on the following assumptions.

1. (i) Constant properties: we assume that the properties of the water–glycerol mixture, namely the density $\rho$ and diffusion coefficient $D$, are constant and equal to the values corresponding to the initial composition. These properties can be obtained from figure 12 in Appendix D by setting $w = w_i$. Setting the density to a constant in the continuity equation (3.1) yields that $u$ is independent of $z$ and depends only on $t$. Hence

   (4.1)\begin{equation} u(z,t) ={-} v_y(t) . \end{equation}

2. (ii) Linearisation of the water vapour concentration difference: the concentration of water vapour at the liquid–gas interface depends on the concentration of water at $z=0^+$ in (3.6). To solve the model analytically, we linearise the expression in (3.5) for the difference in concentration of the water vapour between $c_{w,s}(w_i)$ and $c_{w,s}(w_{eq})$ in terms of $w$:

   (4.2)\begin{equation} c_{w,s} - c_{w,\infty} = c_{w,s}^0\,\frac{x_i \psi_i - H_r}{w_i - w_{eq}}\,(w \vert_{z=0} - w_{eq}), \end{equation}
   where $x_i$, $\psi _i$ and $w_i$ are the initial mole fraction, activity coefficient and weight fraction of water in the liquid mixture, respectively, and $w_{eq}$ is the weight fraction of water at equilibrium, i.e. when $c_{w,s}$ becomes equal to $c_{w,\infty }$ and evaporation stops. In (4.2), we have effectively linearised $c_{w,s}-c_{w,\infty}$ in terms of $w$, between the initial ($w=w_i$) and final ($w=w_{eq}$) concentrations of water (see detailed derivation in Appendix C). Combining (3.5) and (4.2), we get
   (4.3)\begin{equation} \frac{\mathrm{d}M^{\prime \prime}}{\mathrm{d}t} = h^{{\ast}} (w \vert_{z=0} - w_{eq}) , \end{equation}
   where $h^{\ast }$, defined as
   (4.4)\begin{equation} {h^{{\ast}} = \frac{4D_{v,a} c_{w,s}^0 }{{\rm \pi} R\,\Delta w_i}\,(x_i \psi_i - H_r) ,} \end{equation}
   is a modified mass transfer coefficient, with $\Delta w_i = w_i - w_{eq}$. We put an asterisk in $h^{\ast }$ to indicate that its units (kg m$^{-2}$ s) are different from those of the conventional mass transfer coefficient $h$ (Incropera et al. 2007), which is related to $h^*$ through $h=h^*/\rho$ (unit m s$^{-1}$).

3. (iii) Velocity of meniscus: assuming that the volume of the water–glycerol mixture is the sum of the glycerol and water volumes, the velocity at which the length of the liquid column recedes is given by

   (4.5)\begin{equation} v_y ={-}\frac{\mathrm{d}}{\mathrm{d}t}\int_0^L \frac{\rho w}{\rho_w}\,\mathrm{d}z, \end{equation}
   i.e. the negative of the time derivative of the volume occupied by the water, since the glycerol volume is constant. Since $\mathrm {d}M^{\prime \prime }/\mathrm {d}t$ denotes the rate at which the water mass per unit cross-section of capillary is lost, we get
   (4.6)\begin{equation} v_y = \frac{1}{\rho_w}\,\frac{\mathrm{d}M^{\prime \prime}}{\mathrm{d}t} . \end{equation}
   Note that this equation is identical to the exact expression (2.1) obtained for the case where the liquid contains just water. This is a direct consequence of the nearly ideal character of the water–glycerol mixtures.

We define a diffusive length scale $l_D$ (as also done by Salmon et al. 2017) based on the mass transfer coefficient by considering a balance between evaporation and diffusion at the lower interface ($D\,\Delta w / l_D \sim h\,\Delta w$; (3.3) and (4.3)), yielding

(4.7)\begin{equation} l_D = \frac{D}{h}. \end{equation}

Based on the aforementioned considerations, we non-dimensionalise the variables as follows:

(4.8a–d)\begin{align} {\tilde{w} = \frac{w_i - w}{w_i - w_{eq}} = \frac{w_i-w}{\Delta w_i}}, \quad \tilde{t} = \frac{t}{l_D^2/D} = \frac{h^2 \, t}{D}, \quad \tilde{z} = \frac{z}{l_D} = \frac{h \, z}{ D}, \quad \tilde{L} = \frac{L}{l_D} = \frac{h \, L}{D}. \end{align}

Further, since the velocity of the interface is in fact a proxy for the mass transfer (evaporation) rate of water, as follows from mass conservation, we can describe the system in terms of the Sherwood number, $Sh$. This parameter denotes the non-dimensional velocity or non-dimensional mass transfer rate:

(4.9)\begin{equation} Sh = \frac{\rho_w v_{y}}{h^{{\ast}}\,{\Delta w_i}} . \end{equation}

Substituting (4.1)–(4.9) into the governing differential equations (3.1) and (3.2) of mass transport in the capillary, the initial and boundary conditions (3.3), and the equation governing evaporation of water into air (3.5), we get the following system of equations and initial and boundary conditions:

(4.10) \begin{gather} \frac{\partial \tilde{w}}{\partial \tilde{t}} - {\frac{\rho\,\Delta w_i\,Sh}{\rho_w}}\,\frac{\partial \tilde{w}}{\partial \tilde{z}} = \frac{\partial^2 \tilde{w}}{\partial \tilde{z}^2} , \end{gather}

(4.11)\begin{gather} \left. \begin{gathered} \displaystyle\tilde{w} = 0 \quad \text{at}\ \tilde{t}=0,\\ \displaystyle\frac{\partial \tilde{w}}{\partial \tilde{z}} = 0\quad \text{at}\ \tilde{z} = \tilde{L}, \\ {\displaystyle \frac{\partial \tilde{w}}{\partial \tilde{z}} - \left(1 - \frac{\rho\,\Delta w_i\,Sh}{\rho_w}\right) \tilde{w} + \left( 1 - \frac{\rho w_i\,Sh}{\rho_w}\right) = 0} \quad \text{at}\ \tilde{z}=0 , \end{gathered} \right\} \end{gather}

(4.12)\begin{gather} {Sh = 1 - \tilde{w}(\tilde{z}=0, \tilde{t}) . } \end{gather}

Finally, for the case of pure water, we do not need to solve the model for the binary mixture, since there is no change in the concentration. The velocity of the interface can then be obtained directly by substituting $w_i=1$ in (4.3):

(4.13)\begin{equation} v_y = \frac{4D_{v,a}}{{\rm \pi} R \rho_w}\,c_{w,s}^0 (1 - H_r ), \end{equation}

which corresponds to $Sh=1$.

4.1. Semi-infinite transient diffusion model

As the zeroth-order simplification, we assume that advection is negligibly small. Further, we approximate the liquid column as a semi-infinite medium. Thus the governing equation and the initial and boundary conditions ((4.10) and (4.11)) can be reduced to

(4.14)\begin{equation} \frac{\partial \tilde{w}}{\partial \tilde{t}} = \frac{\partial^2 \tilde{w}}{\partial \tilde{z}^2} , \end{equation}

with

(4.15)\begin{equation} \left. \begin{gathered} \displaystyle\tilde{w} = 0\quad \text{at}\ \tilde{t}=0,\\ \displaystyle\frac{\partial \tilde{w}}{\partial \tilde{z}} = 0\quad \text{at}\ \tilde{z} \rightarrow \infty, \\ \displaystyle \frac{\partial \tilde{w}}{\partial \tilde{z}} = \tilde{w} - 1\quad \text{at}\ \tilde{z}=0 . \end{gathered} \right\} \end{equation}

This is a classical transient diffusion problem with mixed boundary conditions, whose solution is given by (Incropera et al. 2007)

(4.16)\begin{equation} \tilde{w} \left( \tilde{z},\tilde{t} \right) = \mathrm{erfc}\left( \frac{\tilde{z}}{2\sqrt{\tilde{t}}} \right) -\exp \left(\tilde{t} + \tilde{z}\right)\mathrm{erfc}\left(\frac{\tilde{z}}{2\sqrt{\tilde{t}}} + \sqrt{\tilde{t}}\right) . \end{equation}

The velocity of the top interface can now be evaluated from the evaporation rate of water using (4.9), (4.12) and (4.16) to yield

(4.17)\begin{equation} v_y = \frac{h^{{\ast}}\,{\Delta w_i}}{\rho_w} \exp(\tilde{t})\,\mathrm{erfc}\left(\sqrt{\tilde{t}}\right) , \end{equation}

or

(4.18)\begin{equation} Sh = \exp ( \tilde{t} )\,\mathrm{erfc}\left( \sqrt{\tilde{t}}\right) . \end{equation}

For $\tilde{t} \rightarrow 0$, $Sh \rightarrow 1$, whereas for $\tilde{t} \gg 1$, $Sh \approx 1/\sqrt {{\rm \pi} \tilde{t} }$.

The predictions of the semi-infinite transient diffusion model (4.18) are compared in figure 4 with the experimental measurements (discrete data points) and the numerical simulations (continuous lines). It can be observed that (4.18) predicts correctly the early-time limit $Sh=1$ for all cases, and $Sh \sim 1/\sqrt {\tilde{t} }$ at intermediate times for $w_i=0.2$ and $0.6$. Moreover, figure 4 also shows that the predicted value of $Sh$ from the model agrees reasonably well with the experiments and the simulations for $w_i=0.2$ and $0.6$. However, (4.18) severely underpredicts $Sh$ for $w_i=0.9$ until $\tilde{t} \approx 80$. Furthermore, (4.18) also fails to capture the steep decay in $Sh$ seen in the simulations for $w_i=0.6$ and $0.9$ at long times (figures 4b,c). This calls for a careful re-examination of the assumptions made in the model.

Figure 4. (a–d) Comparisons of the normalised evaporation rates (Sherwood number $Sh$) against normalised time ($\tilde{t}$) obtained from the experiments (discrete data points) and the numerical simulations (continuous lines) with theoretical values (dashed lines) obtained using the semi-infinite transient diffusion model (4.18), for different initial weight fractions of water $w_i=0.2,0.6,0.9,1.0$, respectively. Notice that in the $w_i = 1.0$ case in (d), the theoretical curve corresponds to $Sh = 1$, not that given by (4.18). Note also that all the theoretical curves for $w_i<1.0$ are the same, as (4.18) does not take into account $w_i$. Insets highlight the comparisons between the experiments, the numerical simulations and the simplified analytical modelling.

4.2. Semi-infinite transient diffusion model with advection – asymptotic solution for $\tilde{t} \gg 1$ and approximate solution for all times

In order to improve the agreement of our model with the experiments and full numerical simulations, we turn our attention to the effect of advection. Indeed, advection is expected to become important at short times and for values of $w_i$ close to unity. To make this clear, we define the Péclet number of the problem as the coefficient of the advection term in (4.10), namely $\Delta w_i\,Sh\,\rho / \rho _w = \rho v_y / h^{\ast } = (w \vert _{z=0} - w_{eq}) \rho /\rho _w$. Notice that in our problem, the Péclet number is nothing more than the Sherwood number with a coefficient that modulates the importance of the initial water concentration. However, we find it useful to work with this quantity in order to evaluate the effect of advection. This Péclet number is plotted in figure 5, where we can see that it becomes of order unity during the first stage of the evaporation process for $w_i \gtrsim 0.6$.

Figure 5. Plots showing the variation of the Péclet number $\rho \,\Delta w_i\,Sh/ \rho _w$ with the normalised time $\tilde{t}$ for three different initial water concentrations, $w_{i} = 0.2$, 0.6 and 0.9.

To predict quantitatively the evaporation dynamics including the effect of advection, we propose an improvement over our simplistic analytical model. We will first develop an asymptotic solution, valid formally in the limit $\tilde{t} \gg 1$, and then present an approximate expression for $Sh$ that converges uniformly to the asymptotic solution for $\tilde{t} \gg 1$ and to $Sh \approx 1$ for $\tilde{t} \ll 1$.

We start from the problem formulated in (4.10) and (4.11), with $\tilde {L}\rightarrow \infty$. In order to make the problem tractable, we assume that the mixture's density is constant and equal to the initial one, $\rho = \rho (w_i)$. The transport problem (4.10) with advection becomes

(4.19)\begin{equation} {\frac{\partial\tilde{w}}{\partial\tilde{t}} - \Delta w_i\,Sh\,\frac{\rho}{\rho_w}\,\frac{\partial\tilde{w}}{\partial\tilde{z}} = \frac{\partial^2\tilde{w}}{\partial\tilde{z}^2},} \end{equation}

with one initial and two boundary conditions (4.11), namely

(4.20)\begin{equation} \left. \begin{gathered} \displaystyle {\tilde{w} = 0} \quad \text{at}\ {\tilde{t} = 0 ,}\\ \displaystyle {\frac{\partial\tilde{w}}{\partial\tilde{z}} = 0} \quad \text{at}\ {\tilde{z} \rightarrow \infty,}\\ \displaystyle {\frac{\partial\tilde{w}}{\partial\tilde{z}} - \left(1 - \Delta w_i\,\frac{\rho}{\rho_w}\, Sh\right)\tilde{w} + 1 - \frac{\rho}{\rho_w}\,w_i\,Sh = 0} \quad \text{at}\ {\tilde{z} = 0,} \end{gathered} \right\} \end{equation}

and with the Sherwood number $Sh$ (the water evaporation rate at the interface) given by (4.12),

(4.21)\begin{equation} {Sh = 1 - \tilde{w}(\tilde{z}=0, \tilde{t}).} \end{equation}

4.2.1. Asymptotic solution

We will now look for solutions at long times. Inspired by both experiments and numerical simulations, we investigate solutions where

(4.22)\begin{equation} {Sh = A\tilde{t}^{-\lambda}}, \end{equation}

with $\lambda > 0$. Further, since at long times we expect $\tilde{w}$ to approach unity everywhere, we propose the following change of variables:

(4.23)\begin{equation} {W = 1 - \tilde{w} - Sh .} \end{equation}

Thus $W \rightarrow 0$ as $\tilde{w} \rightarrow 1$ and $Sh \rightarrow 0$, i.e. at long times. Reformulating the problem in terms of the new variable and neglecting terms of $\mathcal {O}(Sh^2)$, (4.19)–(4.21) can be rewritten as

(4.24)$$\begin{gather} {\frac{\mathrm{d}\,Sh}{\mathrm{d}\tilde{t}} + \frac{\partial W}{\partial\tilde{t}} - \Delta w_i\,\frac{\rho}{\rho_w}\,Sh\,\frac{\partial W}{\partial\tilde{z}} = \frac{\partial^2 W}{\partial\tilde{z}^2} ,} \end{gather}$$

(4.25)$$\begin{gather}{W = 1 - Sh} \quad \mathrm{{at}}\ {\tilde{t} = 0 } , \end{gather}$$

(4.26)$$\begin{gather}{\frac{\partial W}{\partial\tilde{z}} = 0} \quad \mathrm{{at}}\ {\tilde{z} \rightarrow \infty ,} \end{gather}$$

(4.27)$$\begin{gather}{-\frac{\partial W}{\partial\tilde{z}} + Sh \left(1 + w_{eq}\,\frac{\rho}{\rho_w}\right)= 0} \quad \mathrm{{at}}\ {\tilde{z} = 0} , \end{gather}$$

(4.28)$$\begin{gather}{W(\tilde{z} = 0, t) = 0 .} \end{gather}$$

We seek self-similar solutions of the type $W(\tilde{z}, \tilde{t} ) = F(\eta )$, with $\eta = \tilde{z} / \tilde{t} ^\lambda$, in the limit $\tilde{t} \gg 1$. For such a self-similar solution to exist, the exponent $\lambda$ of $\tilde {t}$ in the self-similar variable $\eta$ must be the same as that in the definition of $Sh$ by virtue of (4.27). Indeed, $\partial W/\partial \tilde{z} = F^{\prime } \tilde{t} ^{-\lambda }$. So if the time exponent in $Sh$ and $\eta$ were different, then (4.27) could not be made self-similar.

Substituting the proposed ansatz into (4.24), we get

(4.29)\begin{equation} {-\lambda A \tilde{t}^{-\lambda-1} - \lambda \eta \tilde{t}^{{-}1}F^{\prime} - \Delta w_i\,\frac{\rho}{\rho_w}\,A \tilde{t}^{{-}2\lambda} F^{\prime} = \tilde{t}^{{-}2\lambda}F^{\prime \prime} .} \end{equation}

For $\tilde{t} \gg 1$, the first term is always negligible compared to the second one, since $\lambda > 0$. Then, balancing the second term with the last two terms, we obtain that the equation becomes self-similar if $\lambda = 1/2$, as expected from the experiments and the simulations, resulting in

(4.30)\begin{equation} {-\left(\frac{1}{2}\eta + A\,\Delta w_i\,\frac{\rho}{\rho_w}\right) F^{\prime} = F^{\prime \prime} .} \end{equation}

This differential equation can be solved with the boundary conditions $F(0) = 0$ (from (4.28)) and $F^{\prime }(0) = A$ (from (4.27)) to yield

(4.31)\begin{align} {F(\eta)} &= {-\sqrt{\rm \pi}A\left(1 + w_{eq}\,\frac{\rho}{\rho_w}\right)\exp\left(A^2\,\Delta w_i^2 \left(\frac{\rho}{\rho_w}\right)^2\right) }\nonumber\\ &\quad \times {\left(\mathrm{erf}\left(A\,\Delta w_i\,\frac{\rho}{\rho_w}\right)-\mathrm{erf}\left(A\,\Delta w_i\,\frac{\rho}{\rho_w} + \frac{\eta}{2}\right)\right) .} \end{align}

Finally, an additional condition is needed to determine the value of $A$. This condition stems from the behaviour of $W$ far away from the evaporation boundary $\tilde{z} = 0$, where $W \rightarrow 1$ (from (4.25) while neglecting $Sh \ll 1$ against unity). Introducing this condition into (4.31), we finally get

(4.32)\begin{equation} {\sqrt{\rm \pi}A \left(1 + w_{eq}\,\frac{\rho}{\rho_w}\right) \exp\left(A^2\,\Delta w_i^2 \left(\frac{\rho}{\rho_w}\right)^2\right) \mathrm{erfc}\left(A\,\Delta w_i\, \frac{\rho}{\rho_w}\right) - 1 = 0 .} \end{equation}

The value of $A$ can be evaluated numerically. Plotting $A$ against $w_i$, we see that $A$ grows monotonically with the initial water concentration, recovering the diffusion-driven asymptotic solution $A = 1/\sqrt {{\rm \pi} }$ for $w_i \ll 1$ (if the additional assumption that $w_{eq} = 0$ is made; figure 6). This means that the higher the initial concentration of water $w_i$, the greater the advective enhancement of mass transfer compared to pure diffusion.

Figure 6. Coefficient $A$ in $Sh = A \tilde{t} ^{-1/2}$, obtained by solving (4.32) numerically, as a function of the initial water concentration. At low initial water concentrations, and under the assumption that $w_{eq} \approx 0$, we recover the pure diffusion regime, $A = 1/\sqrt {{\rm \pi} }$ (blue dashed line). For large water concentrations, $A$ grows unbounded, indicating that the regime $Sh \sim \tilde{t} ^{-1/2}$ is never reached, as $Sh = 1$ in this limit.

4.2.2. Uniform approximation

For practical applications, it is desirable to have an approximate expression that converges to the asymptotic solutions in the limits $t \ll 1$ ($Sh \approx 1$) and $t \gg 1$ ($Sh \approx A\tilde{t} ^{-1/2}$). To this end, we notice that the exact solution for the problem without advection,

(4.33)\begin{equation} Sh = \exp\left(\tilde{t}\right)\mathrm{erfc}\left(\sqrt{\tilde{t}}\right), \end{equation}

captures the behaviour of the numerical solution with advection, except that the prefactor of the equation $Sh \sim \tilde{t} ^{-1/2}$ in (4.33) at $\tilde{t} \gg 1$ is $1/\sqrt {{\rm \pi} }$, instead of $A$ (as in (4.22)). Thus we propose the following expression to approximate the full numerical solution uniformly at all times:

(4.34)\begin{equation} Sh = \exp\left(\frac{\tilde{t}}{{\rm \pi} A^2}\right)\mathrm{erfc}\left(\sqrt{\frac{\tilde{t}}{{\rm \pi} A^2}}\right). \end{equation}

We observe that this approximation reproduces much more successfully the numerical simulations than the solution without advection, as shown in figure 7.

Figure 7. Variation of the normalised evaporation rates (Sherwood number $Sh$) with normalised time ($\tilde{t}$) obtained from the experiments (discrete data points), the numerical simulations (continuous lines), the asymptotic solution for long times (dash–dotted lines) and the approximate expression ((4.34), dashed lines), for different initial weight fractions of water $w_{i}$: (a) $w_{i} = 0.2$, (b) $w_{i} = 0.6$, (c) $w_{i} = 0.9$ and (d) $w_{i} = 1.0$. Insets highlight the comparisons between the experiments, the numerical simulations and the simplified analytical modelling.

To conclude this subsection, we point out that another uniform approximation for computing the Sherwood number $Sh$, similar to (4.34), can be obtained by solving the simplified problem (4.19)–(4.20) analytically considering advection as being quasi-constant. For completeness, this solution is described in Appendix E. The fact that treating advection as quasi-constant yields an expression that works fairly well is interesting, as it points out that advection is important mostly when it is nearly constant, that is, at short times (see figure 5). This makes sense, since it is at this stage when the Péclet number reaches the largest value. This idea could be useful in pursuing further analytical approaches for similar problems.

4.3. Transient diffusion model with finite-length effects

The models described in §§ 4.1 and 4.2 can explain faithfully both $Sh=1$ and $Sh \sim \tilde{t} ^{-1/2}$ behaviours seen in the experiments and the simulations. However, we are yet to explain the sharp deviation from $Sh \sim \tilde{t} ^{-1/2}$ seen for very late times in the simulations at $w_i = 0.6$ and $0.9$ (figures 7b,c). To this end, we turn our attention to the semi-infinite assumption.

For the semi-infinite assumption to hold, the penetration depth of the diffusion front $\delta (t) = \sqrt {D t}$ should be much smaller than the length of the liquid column $L(t)$. We plot the variation of $\delta /L\ (= \sqrt {D t}/L)$ with $\tilde{t}$ as obtained from the experiments and the numerical simulations in figure 8 for different $w_i$. It can be observed that for $w_i=0.9$, $\delta /L=1$ at $\tilde{t} \approx 60$, which agrees approximately with the time when the slope of $Sh( \tilde{t} )$ starts to deviate from $Sh \sim \tilde{t} ^{-1/2}$ in figure 7(c). The same holds true for $w_i=0.6$ at $\tilde{t} \approx 10^3$ (figure 7b). Thus we conclude that although the semi-infinite assumption holds at early times, finite-length effects should be included at later times for $w_i=0.6$ and $0.9$ in order to capture accurately the physics of the problem.

Figure 8. Variation of the penetration depth $\delta$ (normalised by the length $L$ of the liquid column) with the normalised time $\tilde{t}$. Discrete data points are based on the experiments, and continuous lines are based on the numerical simulations. The dotted line at $\delta /L=1$ indicates when the penetration depth is equal to the size of the liquid column. The penetration depth is defined as $\delta (t) = \sqrt {D t}$.

We show in this subsection that this is an effect of the finite length of the capillary, which becomes relevant at long times. To include the effects of the finite length, we use the original boundary conditions of (3.3).

We will first check how the three terms in the governing equation (4.10), reproduced here for convenience, compare with each other for the late times when the finite-length effects start playing a role:

(4.35)\begin{equation} \frac{\partial\tilde{w}}{\partial\tilde{t}} - \frac{\rho\,\Delta w_i}{\rho_w}\,Sh\,\frac{\partial\tilde{w}}{\partial\tilde{z}} = \frac{\partial^2\tilde{w}}{\partial\tilde{z}^2} . \end{equation}

When the boundary layer becomes of the order of the length of the domain, $\delta \sim L$, the first and second spatial derivatives in (4.35) scale as $\Delta \tilde{w} /\tilde{L}$ and $\Delta \tilde{w} /\tilde{L} ^2$, respectively. For $w_i < 1$, the volume occupied by water at long times is small compared to that occupied by glycerol, so $\tilde{L}$ tends asymptotically to a constant as the total liquid volume changes very slowly. This means that the orders of magnitude of the first and second spatial derivatives differ by a constant factor $\tilde{L}$. In this situation, since the prefactor of the advective term goes to zero as $Sh \rightarrow 0$, the advective term is going to be a factor of $Sh$ smaller than the diffusive one and can be neglected. Thus we formulate a quasi-constant-length transient diffusion model valid for $Sh \ll 1$.

The governing equations and the initial and boundary conditions ((4.10)–(4.12)) can now be written as

(4.36)\begin{gather} \frac{\partial \tilde{w}}{\partial \tilde{t}} = \frac{\partial^2 \tilde{w}}{\partial \tilde{z}^2} , \end{gather}

(4.37) \begin{gather} \left. \begin{gathered} \displaystyle\tilde{w} = 0\quad \text{at}\ \tilde{t}=0,\\ \displaystyle\frac{\partial \tilde{w}}{\partial \tilde{z}} = 0\quad \text{at}\ \tilde{z} = \tilde{L} , \\ \displaystyle \frac{\partial \tilde{w}}{\partial \tilde{z}} = \tilde{w}-1\quad \text{at}\ \tilde{z}=0 , \end{gathered} \right\} \end{gather}

(4.38)\begin{gather} {Sh = 1 - \tilde{w}(\tilde{z}=0, \tilde{t}) .} \end{gather}

The solution of this system of equations is given by (Crank 1975)

(4.39)\begin{gather} \tilde{w}(\tilde{z},\tilde{t}, \tilde{L}) = 1 - \sum_{n=1}^{\infty} \frac{2 \tilde{L} \cos \left( \lambda_n \left( 1 - \tilde{z} / \tilde{L} \right)\right) \exp \left(-\lambda_n^2 \tilde{t} / \tilde{L}^2\right)}{\left( \lambda_n^2 + \tilde{L}^2 + \tilde{L} \right) \cos(\lambda_n)}, \end{gather}

(4.40)\begin{gather} { Sh = 1 - \tilde{w}(\tilde{z}=0, \tilde{t}) = \sum_{n=1}^\infty \frac {2\tilde{L} \exp \left( -\lambda_n^2 \tilde{t} / \tilde{L}^2 \right)}{\lambda_n^2 + \tilde{L}^2 + \tilde{L}} ,} \end{gather}

where $\lambda _n$ is the $n{\textrm {th}}$ root of the equation $\lambda \tan (\lambda ) = \tilde{L}$. Combining (4.9) and (4.40), the velocity of the top interface can now be written as

(4.41)\begin{equation} \frac{\mathrm{d}L}{\mathrm{d}t} = v_y = \frac{h^{{\ast}}\,\Delta w_i}{\rho_w} \sum_{n=1}^\infty \frac{2\tilde{L} \exp \left( -\lambda_n^2 \tilde{t} / \tilde{L}^2 \right)}{\lambda_n^2 + \tilde{L}^2 + \tilde{L}} . \end{equation}

Equation (4.41) is integrated in time using the built-in Matlab function ode45 (which implements an adaptive-time-step Runge–Kutta algorithm of fourth order) to obtain $L(t)$ and subsequently the variation of $Sh$ with $\tilde{t}$ (figure 9c). For $\tilde{t} \gg 1$, $v_y$ can be approximated by the first term in the series:

(4.42)\begin{equation} v_y = \frac{h^{{\ast}}\,\Delta w_i}{\rho_w}\,\frac{2\tilde{L} \, \exp \left( -\beta_1^2 \tilde{t} / \tilde{L}^2 \right)}{\beta_1^2 + \tilde{L}^2 + \tilde{L}} . \end{equation}

When the changes in the length $L$ of the liquid column are slow, the velocity decreases exponentially with time (as per (4.42)), thus correctly predicting the late-time $Sh (\tilde{t} )$ behaviour of $w_i=0.9$ in the simulations (figure 9c). Note that observing this regime experimentally is complicated from the practical point of view, due to the very small velocities associated with it. Although we can observe it easily in the numerical simulations, we are able to see the beginning of this regime in the experiments only for the most favourable case, $w_i = 0.9$, for which this regime appears at velocities of the order of less than a micron per second (see figure 3). In conclusion, this model captures successfully the deviation from $Sh\sim \tilde{t} ^{1/2}$ due to the finite-length effects.

Figure 9. Plots comparing the normalised evaporation rates (Sherwood number $Sh$) versus normalised time ($\tilde{t}$) obtained from the experiments (discrete data points) and the numerical simulations (continuous lines) with the theoretical values obtained using (a) the semi-infinite transient diffusion model (dashed line), (b) the semi-infinite uniform approximation model (dotted line) and (c) the transient diffusion model with finite-length effects (dash–dotted line) for an initial weight fraction of water $w_i = 0.9$.

5. Conclusions and outlook

In this work, we have studied the evaporation of aqueous glycerol solutions in cylindrical capillaries having a circular cross-section. We characterised the drying behaviour in terms of a normalised mass transfer rate (Sherwood number $Sh$) and a normalised time $\tilde{t}$. Our experiments demonstrate quantitatively how the addition of glycerol reduces the evaporation rate of water. The corresponding direct numerical simulations indicate that modelling the system as a one-dimensional advection–diffusion mass transfer problem with composition-dependent properties can reproduce quantitatively the evaporation behaviour observed in the experiments. The evaporation of water shows three main regimes: (i) $Sh = 1$, (ii) $Sh \sim 1/\sqrt {\tilde{t} }$ and (iii) $Sh \sim \exp (-\tilde{t} )$. We describe the physical origins of these regimes using a one-dimensional simplistic analytical model with constant material properties and a linearised composition-dependent activity of water.

Modelling the system as a problem of pure diffusion in a semi-infinite medium reproduces $Sh = 1$ and $Sh \sim 1/\sqrt {\tilde{t} }$ as the early-time and late-time behaviours, respectively. In the early time limit, $Sh=1$ corresponds to a rapid replenishment of water at the evaporating interface, leading to a constant evaporation rate and constant interfacial concentration of water. In the late-time regime, replenishment of the interfacial concentration of water is limited by diffusion, leading to the classical diffusion-like $Sh \sim 1/\sqrt {\tilde{t} }$ behaviour.

However, we also show that even if the pure diffusion model captures the scaling relations of $Sh(\tilde{t} )$ correctly, advective replenishment of water needs to be considered for more precise prediction of the evaporation rates. A Péclet number defined as $\Delta w_i\,Sh\,\rho / \rho _w = \rho v_y / h^{\ast } = (w \vert _{z=0} - w_{eq}) \rho /\rho _w$ (the coefficient of the advective term in (4.10)) is the relevant parameter dictating the importance of advection. Thus advection is small when the interfacial concentration of water $w \vert _{z=0}$ is close to the equilibrium concentration $w_{{eq}}$, and high otherwise.

Finally, we show that as the diffusive penetration depth $\delta = \sqrt {Dt}$ increases and the length $L(t)$ of the system decreases, they can become comparable in magnitude. In such a scenario, the semi-infinite approximations hold only as long as the late-time behaviour is modified to $Sh \sim \exp (-\tilde{t} )$. This change in the evaporation regime essentially reflects the effect of the finite size of the liquid column.

Even though a model with constant material properties was used to describe the evaporation of a binary mixture, the spatio-temporal changes in properties such as the density and the diffusion coefficient affect the precise quantitative prediction of evaporation rates. The direct numerical simulations are devoid of these deficiencies. However, the simplified analytical models provide valuable insight into the essential physics of the system and can provide predictions for more complex liquid mixtures.

We also note that for predicting the evaporation rate, we have used the expression corresponding to that of a thin droplet on a substrate ($\theta _{drop}=0^\circ$), giving excellent predictions that match the experimental observations. However, one can include corrections to the mass transfer coefficient $h^{\ast }$ to account for the difference between the air–liquid interface of a droplet and the air–liquid interface at the mouth of a capillary tube (Li et al. 2019). Moreover, during evaporation, the shape of the lower interface changes ($\theta <90^\circ$) until the meniscus depins. This change in the shape of the interface may also require a small correction to $h^{\ast }$, and might explain the very small decrease in $Sh$ seen in the experiments on pure water (figure 3b). However, further studies along the lines of D'Ambrosio et al. (2021) are required to confirm this hypothesis.

Finally, all measurements were limited to times during which the lower meniscus was pinned at the mouth of the capillary. When the lower meniscus depins and propagates into the capillary, the rate-limiting step of evaporation would change from three-dimensional vapour diffusion to one-dimensional vapour diffusion. In the case of a single-component liquid evaporating from a square capillary, there have been efforts to predict the time when the lower meniscus depins (Chauvet et al. 2010). Further studies are required to predict the same for multi-component liquids and capillaries of various geometries and configurations (e.g. inclination with respect to gravity). Nonetheless, in this study, we have shown that these kinds of phenomena may not be essential to predict, with a reasonable degree of accuracy, the mass transfer rate in a relatively complex system such as the one that we consider here.

The aforementioned results can also be applied directly to predict the evaporation of multi-component liquids from porous structures, which can be modelled as bundles of thin capillaries. Evaporation from capillaries can also help us to understand the evaporative behaviour of biological fluids (such as blood, saliva or liquids in respiratory droplets; Merhi et al. 2022; Seyfert et al. 2022) or novel liquid mixtures (for applications such as evaporative cooling and spray drying). Finally, studying evaporation from capillaries can also be useful in the case of inkjet printing, where the evaporation from the tip of the printing nozzle can lead to changes in the composition of the ink (Rump et al. 2023). Our model can provide insight into the changes in the composition at the nozzle tip for a given time scale and assist in choosing carefully the properties of the ink.