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 Reference Diddens2017). 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 Reference Christy, Hamamoto and Sefiane2011; Bennacer & Sefiane Reference Bennacer and Sefiane2014; Diddens et al. Reference Diddens, Tan, Lv, Versluis, Kuerten, Zhang and Lohse2017; Lopez de la Cruz et al. Reference Lopez de la Cruz, Diddens, Zhang and Lohse2021) or condensation of water onto a water–glycerol droplet (Shin, Jacobi & Stone Reference Shin, Jacobi and Stone2016; Diddens Reference Diddens2017). A more detailed explanation of such instabilities can be found in Diddens (Reference Diddens2017).

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. Reference Salmon, Doumenc and Guerrier2017). 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.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 Reference Popov2005; Stauber et al. Reference Stauber, Wilson, Duffy and Sefiane2014):

(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 Reference Popov2005; Stauber et al. Reference Stauber, Wilson, Duffy and Sefiane2014). 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 (Reference Heil and Hazel2006), 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. (Reference D'Errico, Ortona, Capuano and Vitagliano2004), while the mass density $\rho$ was fitted according to the data of Takamura, Fischer & Morrow (Reference Takamura, Fischer and Morrow2012). The activity coefficient of water was calculated by AIOMFAC (Zuend et al. Reference Zuend, Marcolli, Booth, Lienhard, Soonsin, Krieger, Topping, McFiggans, Peter and Seinfeld2011).

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.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. Reference Incropera, Dewitt, Bergman and Lavine2007)

(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 Reference Crank1975)

(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$.