1 Introduction
The aggregation of many wet hairs into a series of clumps is familiar from everyday examples including wet paint brushes or eyelashes wetted by tears (see Bico et al. Reference Bico, Roman, Moulin and Boudaoud2004; Kim & Mahadevan Reference Kim and Mahadevan2006, for example). In these scenarios, the surface tension of a liquid acts to minimize the liquid surface area but is resisted by the bending stiffness of the hairs involved. The result is elastocapillary aggregation with a number of clumps of hairs, separated by larger gaps. At a microscopic scale, a similar phenomenon is observed in the manufacture of microelectromechanical systems (MEMS) where long thin elements are etched (using photolithography) and then rinsed (Tanaka, Morigami & Atoda Reference Tanaka, Morigami and Atoda1993). In the rinse step, the elastic elements are vulnerable to the formation of clumps that can, if the clumps are large enough, lead to fracture (see figure 1 a).
While, in many applications, elastocapillary aggregation is an undesirable feature of the process and is best avoided, in other situations it is exploited to form controlled patterns (see Chakrapani et al. Reference Chakrapani, Wei, Carrillo, Ajayan and Kane2004; Pokroy et al. Reference Pokroy, Kang, Mahadevan and Aizenberg2009, for example). Indeed, the clumps of carbon nanotubes that form when a nanotube ‘forest’ is wetted and the liquid evaporated (see figure 1 b) have sparked considerable interest as a design strategy at microscopic scales (see de Volder & Hart Reference de Volder and Hart2013; de Volder et al. Reference de Volder, Tawfick, Baughman and Hart2013, for reviews). Despite the importance of other forces (including van der Waals and electrostatic forces) at these small scales, it is clear that capillary forces from the liquid play a vital role: aggregation only happens if the forest has been wetted, and varying the surface tension coefficient changes the properties of the pattern (Tanaka et al. Reference Tanaka, Morigami and Atoda1993; de Volder & Hart Reference de Volder and Hart2013).
Figure 1. Elastocapillary aggregation driven by evaporation of a volatile liquid. (a) Scanning electron micrograph of a pattern created to be part of a MEMS device but damaged by surface tension forces in the rinse step. (Image reproduced from Tanaka et al. (Reference Tanaka, Morigami and Atoda1993), Copyright (1993) The Japan Society of Applied Physics.) (b) The formation of cellular foams from the elastocapillary collapse of a forest of carbon nanotubes. (Image reproduced from Chakrapani et al. (Reference Chakrapani, Wei, Carrillo, Ajayan and Kane2004), Copyright (2004) the National Academy of Sciences.) (c) The simplified model considered here in which rigid blocks, connected to their initial position by linear springs, replace flexible beams. Lubrication theory in the intervening liquid gaps leads to a second-order differential equation for the pressure within each gap,
$p_{j}(x,t)$
, which is solved subject to a no-flux condition at
$x=0$
and imposed capillary pressure at the meniscus,
$x=x_{j}(t)$
(see appendix A).
Most previous theoretical approaches aim to understand the pattern formation in these systems by focusing on energy minimization arguments (see Bico et al. Reference Bico, Roman, Moulin and Boudaoud2004; Py et al. Reference Py, Bastien, Bico, Roman and Boudaoud2007; de Volder et al. Reference de Volder, Park, Tawfick, Vidaud and Hart2011; Duprat et al. Reference Duprat, Protiére, Beebe and Stone2012; Jung, Clanet & Bush Reference Jung, Clanet and Bush2014, for example). This approach neglects the dynamic manner in which the pattern forms and so does not address any features of the pattern that might be controlled dynamically. Boudaoud, Bico & Roman (Reference Boudaoud, Bico and Roman2007) developed a toy ballistic aggregation model that gives predictions for the distribution of cluster sizes. However, only recently have detailed dynamic models been developed to capture the interaction between fluid flow and elasticity that occurs in this problem. Several early authors developed models that couple the elastic deflection with the fluid flow between a pair of elastic beams (Aristoff, Duprat & Stone Reference Aristoff, Duprat and Stone2011; Duprat, Aristoff & Stone Reference Duprat, Aristoff and Stone2011; Taroni & Vella Reference Taroni and Vella2012). In particular, Taroni & Vella (Reference Taroni and Vella2012) showed via their detailed model that even the simple system consisting of a fixed volume of liquid confined between two elastic beams generally exhibits multiple stable equilibria. Taroni & Vella (Reference Taroni and Vella2012) found that which of these equilibria is attained dynamically depends on which is ‘closer’ to the initial condition, rather than which configuration has the lower energy. From this observation, one might expect the final state exhibited by the many-body problem to also be sensitive to the history of deformation, including the initial condition.
Recently, there has been a focus on understanding dynamic aggregation in many-body systems, though to make progress it has been necessary to simplify the elastic problem. Gat & Gharib (Reference Gat and Gharib2013) and Singh, Lister & Vella (Reference Singh, Lister and Vella2014) did this by considering arrays of spring–block elements while Wei & Mahadevan (Reference Wei and Mahadevan2014), Wei et al. (Reference Wei, Schneider, Kim, Kim, Aizenberg and Mahadevan2015) considered rigid rods that rotate, resisted by a torsional spring. These models are sufficiently simple that they allow for the development of analytical understanding of the underlying problem while also being able to replicate some of the more complicated behaviour observed in experiments and numerical simulations. For example, Singh et al. (Reference Singh, Lister and Vella2014) showed that multiple equilibria exist, as in the full problem (Taroni & Vella Reference Taroni and Vella2012), and that the distribution of cluster size predicted numerically is similar to that observed experimentally by Boudaoud et al. (Reference Boudaoud, Bico and Roman2007). Similarly, Munro (Reference Munro2014) demonstrated novel scaling laws for meniscus motion in elastocapillary rise that resemble those found experimentally (Duprat et al. Reference Duprat, Aristoff and Stone2011).
Linear stability analysis of the simplified models of dynamic aggregation have all revealed that pairwise aggregation is the mode of aggregation with the largest growth rate (Gat & Gharib Reference Gat and Gharib2013; Singh et al. Reference Singh, Lister and Vella2014; Wei et al. Reference Wei, Schneider, Kim, Kim, Aizenberg and Mahadevan2015). However, simulations and experiments show that larger clusters form. This apparent discrepancy is believed to be resolved either as the result of iterated clustering (in which quadruples form from a pair of pairs etc., as suggested by Gat & Gharib (Reference Gat and Gharib2013) and Wei et al. (Reference Wei, Schneider, Kim, Kim, Aizenberg and Mahadevan2015)) or by the propagation of a disturbance from a localized perturbation, which leaves behind clusters with a well-defined size (as suggested by Singh et al. (Reference Singh, Lister and Vella2014)).
The models of elastocapillary aggregation discussed above have common features, most notably that the amount of liquid is conserved or is decreased quasi-statically. However, in many applications the liquid is volatile and therefore only wets the elastic elements transiently. Once the liquid has evaporated, elements can only ‘stick’ together if they were brought sufficiently close together while wet that dry adhesive forces (e.g. van der Waals forces) come into play. A key question is then: how close to each other do the elastic objects come during evaporation? In this paper we study a reduced model of this dynamic aggregation accounting for the effect of evaporation.
The fluid mechanics of evaporation has seen a great deal of interest in recent years, with particular emphasis being placed on the evaporation of droplets (Cazabat & Guena Reference Cazabat and Guena2010). Much of this work was motivated by the observation that the dark ring left when a drop of coffee evaporates, the so-called ‘coffee-ring effect’, is caused by the pinning of the contact line and a singularity in the evaporative flux close to the contact line (Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten1997; Popov Reference Popov2005). In other cases, the contact line is not pinned so that the evolution of both the radius and contact angle of the droplet need to be determined dynamically (McHale et al. Reference Mchale, Aqil, Shirtcliffe, Newton and Erbil2005; Stauber et al. Reference Stauber, Wilson, Duffy and Sefiane2014). A variety of different evaporation models are available depending on, among other factors, whether the rate of evaporation is limited by diffusion in the vapour phase (Murisic & Kondic Reference Murisic and Kondic2011; Ledesma-Aguilar, Vella & Yeomans Reference Ledesma-Aguilar, Vella and Yeomans2014), conduction of heat from a solid boundary (Dunn et al. Reference Dunn, Wilson, Duffy, David and Sefiane2009) or convection above the drop (Kelly-Zion et al. Reference Kelly-Zion, Pursell, Hasbamrer, Cardozo, Gaughan and Nickels2013). However, even very simple models of evaporation give rise to extremely intricate fluid mechanical problems whose solution can be difficult to distinguish from experimental observations (Oliver et al. Reference Oliver, Whiteley, Saxton, Vella, Zubkov and King2015).
The coffee-ring effect itself may be viewed as the relatively dilute limit of the drying of a colloidal suspension (Routh Reference Routh2013). However, at higher concentrations, the densely packed layer of particles left after evaporation cracks just as mud cracks when it dries (Lecocq & Vandewalle Reference Lecocq and Vandewalle2002; Dufresne et al. Reference Dufresne, Stark, Greenblatt, Cheng, Hutchinson, Mahadevan and Weitz2006). Cracks may only propagate in these systems when the elastic energy that is released by fracture is sufficient to pay the associated surface energy penalty, much like a crack in a more classical elastic medium (Dufresne et al. Reference Dufresne, Stark, Greenblatt, Cheng, Hutchinson, Mahadevan and Weitz2006). This problem therefore contains similar ingredients to elastocapillary aggregation, namely elasticity, capillary-induced flow and evaporation.
In this paper, we consider the elastocapillary aggregation problem accounting for the motion of individual elastic elements within the system, their displacement relative to their initial position, and the fluid flow between elements that controls the dynamics of the process. We shall focus on the effect that the rate of evaporation has on the pattern that is ultimately formed.
2 Governing equations
To mimic the elasticity of the beams typical in many-body elastocapillary aggregation, we consider an array of rigid blocks, each connected to their initial position by a linear spring of stiffness
$k$
(see figure 1
c). The gap between the blocks is considered to be filled with liquid of initial volume
$V_{0}$
(per unit depth into the page); the initial thickness of each gap is
$w$
.
The surface tension of the liquid,
${\it\gamma}$
, causes a capillary suction beneath each meniscus, which acts to bring the blocks together. This aggregation is resisted by the springs, mimicking the balance between capillarity and elasticity common to elastocapillary aggregation in many systems.
To model the dynamics of aggregation, we assume that the gap is relatively thin,
$V_{0}/w^{2}\gg 1$
. In this case, lubrication theory can be used to determine the pressure profile
$p_{j}(x)$
within each gap, subject to the boundary conditions that there be no flux through the base (
$x=0$
) and that the pressure is equal to the capillary suction at the meniscus,
$x=x_{j}$
. More details of this calculation are given in appendix A; the key result is that the pressure profile within each gap,
$p_{j}(x)$
, is determined in terms of the values of
$x_{j}$
,
$h_{j}$
and the rate at which the gap is growing/shrinking,
$\text{d}h_{j}/\text{d}t$
. The pressure profile
$p_{j}(x)$
can then be integrated along the wetted length of each side of the gap to give the force on the two neighbouring blocks exerted by the liquid in the
$j\text{th}$
gap:
$$\begin{eqnarray}F_{j}=\frac{x_{j}^{3}}{h_{j}^{3}}\frac{\text{d}h_{j}}{\text{d}t}+\frac{x_{j}}{h_{j}}.\end{eqnarray}$$
Here the gap thickness and position of the meniscus have been non-dimensionalized by
$w$
and
$V_{0}/w$
respectively, force has been non-dimensionalized by
$2{\it\gamma}V_{0}/w^{2}$
and time by the capillary–viscous time
$$\begin{eqnarray}t_{cv}=\frac{2{\it\mu}V_{0}^{2}}{w^{3}{\it\gamma}}.\end{eqnarray}$$
In (2.1) the first term represents the lubrication force provided by the liquid in the gap (which resists the motion if
$\text{d}h_{j}/\text{d}t<0$
) while the second term represents the net attractive force that results from capillary suction.
Now, each block has two forces acting on it (from the liquid gaps on either side) and the difference of these must be balanced by the net spring force on the block. Subtracting these net forces for the
$j\text{th}$
and
$(j-1)\text{th}$
gaps, we find that
$$\begin{eqnarray}2K(h_{j}-1)=F_{j+1}-2F_{j}+F_{j-1}.\end{eqnarray}$$
Here the dimensionless spring stiffness
$$\begin{eqnarray}K=\frac{kw^{3}}{4{\it\gamma}V_{0}}.\end{eqnarray}$$
To close the problem we must also give an equation for the evolution of the position of the meniscus,
$x_{j}$
, within the
$j\text{th}$
gap. Previous models have assumed that the total amount of liquid within each gap is conserved, that is
$\text{d}(x_{j}h_{j})/\text{d}t=0$
(Gat & Gharib Reference Gat and Gharib2013; Singh et al.
Reference Singh, Lister and Vella2014). The dynamics of evaporation are, in general, complicated depending on the presence of convection in the atmosphere, the curvature of the meniscus and evaporation-induced Marangoni effects among others. Furthermore, in the two-dimensional geometry of interest here, the classic approach of considering evaporation limited by steady-state diffusion is ill-posed (essentially because of the logarithmic term in the solution of the two-dimensional Laplace equation). For simplicity, therefore, we assume that evaporation occurs at a (constant) dimensional rate
$e$
per unit surface area. It follows that the volume of liquid within each liquid-filled gap evolves according to
$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}(x_{j}h_{j})=-Eh_{j},\end{eqnarray}$$
where the dimensionless evaporation rate
$E$
is defined by
$$\begin{eqnarray}E=\frac{ewt_{cv}}{V_{0}}.\end{eqnarray}$$
Our non-dimensionalization of the problem ensures that a uniform initial condition may be written as
$x_{j}(0)=h_{j}(0)=1$
for each liquid gap, i.e. for each
$j$
. The problem therefore has only two dimensionless parameters: the spring stiffness,
$K$
, and the evaporation rate,
$E$
. Ultimately, we shall be interested in how the final pattern, particularly the maximum number of blocks in any cluster, depends on these two parameters for systems consisting of a large number of spring–block elements. However, to gain some understanding of the behaviour of the system in different regions of
$(E,K)$
parameter space, we first consider the simpler problem of a pair of spring–block elements with a single fluid-filled gap between them.
3 An isolated pair: the two-body problem
For a pair of spring–block elements, there is only a single liquid gap,
$j=1$
, to consider. In this case
$F_{0}=F_{2}=0$
so that (2.3) simplifies considerably. Denoting the (only) gap thickness by
$h(t)$
and the meniscus position by
$x(t)$
we find that
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{x^{3}}{h^{3}}\frac{\text{d}h}{\text{d}t}=-\frac{x}{h}+K(1-h), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}x}{\text{d}t}=-E-\frac{\text{d}h}{\text{d}t}\frac{x}{h}, & \displaystyle\end{eqnarray}$$
with initial conditions
$x(0)=h(0)=1$
.
3.1 Phase-plane analysis
The system of equations (3.1)–(3.2) is autonomous so it is natural to examine the behaviour of its phase plane. To facilitate this, we introduce a rescaled gap aspect ratio,
${\it\xi}(t)=x(t)/[Kh(t)]$
and a time-like coordinate
$s$
such that
$\text{d}/\text{d}s=K^{2}h{\it\xi}^{3}\,\text{d}/\text{d}t$
with
$s(t=0)=0$
. The system (3.1)–(3.2) can then be written
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}h}{\text{d}s}=h(1-h-{\it\xi}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}{\it\xi}}{\text{d}s}=-{\it\xi}[KE{\it\xi}^{2}+2(1-h-{\it\xi})]. & \displaystyle\end{eqnarray}$$
From (3.3)–(3.4) we see that this (slightly opaque) transformation has eliminated individual occurrences of the dimensionless parameters
$K$
and
$E$
; instead only their product
$KE$
appears. The dynamics does still depend on the separate values of
$K$
and
$E$
, however, since the initial conditions are now
${\it\xi}(0)=1/K$
,
$h(0)=1$
. Nevertheless, to understand the phase plane, this rescaling, and the associated reduction of the number of effective parameters, is extremely useful.
It is clear from (3.3)–(3.4) that the (non-trivial) nullclines in the
$({\it\xi},h)$
plane are
$h=1-{\it\xi}$
and
$h=1-{\it\xi}+KE{\it\xi}^{2}/2$
. Combining these with the trivial nullclines
${\it\xi}=0,h=0$
, we find that if
$KE<1/2$
, then the dynamical system (3.3)–(3.4) has two stable fixed points:
${\it\xi}=0,h=1$
(corresponding to two dry blocks in their undeformed, separated configuration) and
${\it\xi}={\it\xi}_{\ast }\equiv [1+(1-2KE)^{1/2}]/KE$
,
$h=0$
(corresponding to the two blocks coming into contact as the liquid between them evaporates). We refer to the fixed point
$({\it\xi}_{\ast },0)$
as the ‘stuck’ configuration since we expect short-ranged forces to maintain contact even after all of the liquid has evaporated. There are also two saddle points for
$KE<1/2$
: one at
$({\it\xi}_{s},0)$
with
${\it\xi}_{s}=[1-(1-2KE)^{1/2}]/KE$
and another at
$(0,0)$
.
Figure 2. Phase planes and numerically computed trajectories showing the behaviour of the system (3.3)–(3.4) for
$KE=0.49$
(a,b) and
$KE=1$
(c,d). In (a,c) the phase planes show various trajectories (solid curves), corresponding to changing the value of
$K$
with
$KE$
fixed, with arrows indicating the direction of increasing time-like coordinate
$s$
(defined such that
$\text{d}/\text{d}s=K^{2}h{\it\xi}^{3}\text{d}/\text{d}t$
). The nullclines
$h=1-{\it\xi}$
and
$h=1-{\it\xi}+KE{\it\xi}^{2}/2$
are shown by the dotted straight line and dashed curve, respectively. For
$KE\leqslant 1/2$
there are stable equilibria at
$(0,1)$
(separated blocks, indicated by an open circle) and
$({\it\xi}_{\ast },0)$
(stuck blocks, indicated by a filled circle), and two saddle points at
$(0,0)$
(indicated by
$+$
) and
$({\it\xi}_{s},0)$
(indicated by
$\times$
). For
$KE>1/2$
there is only the equilibrium point at
$(0,1)$
(open circle) and the saddle at
$(0,0)$
(
$+$
). (b) Trajectories of
$x(t)$
(solid curves) and
$h(t)$
(dashed curves) for
$K=5$
(labelled A) and
$K=0.5$
(labelled B); in each case
$KE=0.49$
and we observe sticking only for sufficiently small
$K$
. (d) Trajectories of
$x(t)$
(solid curves) and
$h(t)$
(dashed curves) for
$K=5$
(labelled C) and
$K=0.5$
(labelled D); in each case
$KE=1$
and sticking is never observed.
A typical phase plane for
$KE<1/2$
is shown in figure 2(a) (with
$KE=0.49$
); the evolution of the physical variables for two example trajectories are shown in figure 2(b). The phase plane, figure 2(a), shows that for
${\it\xi}(0)=1/K$
sufficiently small, i.e. for sufficiently stiff springs, the system ultimately tends to the separated fixed point,
$(0,1)$
. However, for
${\it\xi}(0)=1/K$
larger than some threshold, i.e. for sufficiently weak springs, the system instead tends to the stuck fixed point,
$({\it\xi}_{\ast },0)$
. This is illustrated by the evolution of the physical variables shown in figure 2(b). The critical value of
$K$
at which the transition from separated to sticking occurs,
$K_{\mathit{crit}}(KE)$
, can be determined numerically from the stable manifold of the saddle point at
$({\it\xi}_{s},0)$
(see appendix B). For the case
$KE=1/2$
we find that the system tends to the stuck fixed point if
$K<3.30$
(i.e.
$E>0.15$
) but otherwise becomes separated. In the limit
$KE\ll 1$
we find that
$K_{\mathit{crit}}=27/4-O(E^{2/3})$
, which agrees with the result of Singh et al. (Reference Singh, Lister and Vella2014) in the absence of evaporation.
For
$KE>1/2$
,
${\it\xi}_{\ast }$
is complex so that the only stable fixed point is at
${\it\xi}=0,h=1$
; there is also a single saddle point at
$(0,0)$
. Physically, we deduce that, for sufficiently fast evaporation or sufficiently strong springs, the system must always end up with the blocks separated and all of the liquid evaporated: the blocks are ‘unstuck’. An example of the trajectories of the system, plotted in
$({\it\xi},h)$
space, is shown in figure 2(c); as expected, this shows that all trajectories (i.e. regardless of the value of
$K$
) ultimately reach the unstuck state. This is confirmed by the physical evolution along two example trajectories, as shown in figure 2(d). Note from the case
$K=0.5$
,
$E=2$
(labelled ‘D’ in figure 2
d) that once the blocks start to separate they may do so very rapidly: both evaporation and the widening gap reduce the hydraulic resistance to motion while the total capillary suction in the gap is also dramatically decreased by a decrease in
$x$
.
In summary, our analysis of the phase plane shows that for
$KE>1/2$
the blocks must ultimately tend to the separated state, while for
$KE\leqslant 1/2$
the blocks stick if
$K\leqslant K_{\mathit{crit}}(KE)$
and otherwise separate. This behaviour is summarized in figure 3(a), which shows the regions of
$(K,E)$
parameter space for which the blocks are ultimately stuck or separated. However, we note that in either case the system tends to a state in which all of the liquid has evaporated,
$x=0$
, whether that be with
$h=1$
(the separated state) or with
$h=0$
,
$x/h=K{\it\xi}_{\ast }$
(the stuck state).
Figure 3. Regime diagrams showing the behaviour of the system as functions of the spring stiffness
$K$
and evaporation rate
$E$
. (a) The regime diagram for the two-body problem in which only stuck and separated states exist. The solid curve shows the numerically determined boundary between stuck and separated while the results with
$E=0$
,
$K_{\mathit{crit}}=27/4$
(dotted vertical line), and
$K_{\mathit{crit}}=1/2E$
(dashed line), are also shown. (For
$K\lesssim 3.30$
, the boundary between stuck and separated states is
$KE=1/2$
, see the text.) (b) Results of numerical simulations with
$N=2000$
blocks with randomized initial perturbations showing the largest number of blocks in a cluster,
$N_{\mathit{max}}$
, at various
$(K,E)$
. Here green shows situations in which no blocks aggregate ultimately while red and blue are used to distinguish results at neighbouring points of
$(K,E)$
space. The dashed line again shows the result
$K_{\mathit{crit}}=1/2E$
as a hint that the two-body problem is relevant to understanding the many-body problem.
3.2 Contact occurs in finite time or not at all
In the limit of no evaporation,
$E=0$
, Singh et al. (Reference Singh, Lister and Vella2014) showed that if the springs are sufficiently weak for a pair to stick, then the gap thickness
$h$
decreases with time according to
$h\sim (3t)^{-1/3}$
for
$t\gg 1$
. In this idealized system, therefore, the blocks do not actually contact one another (though in reality van der Waals attraction or the roughness of the blocks would ensure that contact does eventually occur).
The situation with even a small amount of evaporation,
$E>0$
, is qualitatively different in that contact occurs within a finite time. To see this we let
$h\ll 1$
,
$x/h\approx K{\it\xi}_{\ast }$
in (3.1), yielding
$$\begin{eqnarray}\frac{\text{d}h}{\text{d}t}\approx \frac{1-{\it\xi}_{\ast }}{K^{2}{\it\xi}_{\ast }^{3}}=-\frac{E}{2K{\it\xi}_{\ast }}.\end{eqnarray}$$
Crucially, this is constant so that the gap thickness decreases at a constant rate and, hence, the gap must close and the blocks touch.
Note that in the limit of small evaporation,
$E\ll 1$
, and with
$0<K<K_{\mathit{crit}}(KE)$
so that sticking occurs, we have
${\it\xi}_{\ast }\approx 2/KE$
leading to
$$\begin{eqnarray}\frac{\text{d}h}{\text{d}t}\approx -\frac{E^{2}}{4}.\end{eqnarray}$$
We therefore expect the time at which contact occurs to diverge as
$E\rightarrow 0$
, consistent with the previous observation that contact does not occur in the absence of evaporation (Singh et al.
Reference Singh, Lister and Vella2014). This expectation may also be confirmed for the case of no springs,
$K=0$
, where a simple calculation shows that contact occurs at
$t=t_{\ast }$
with
$$\begin{eqnarray}t_{\ast }=\frac{1}{E}+\frac{2}{E^{3/2}(2-E)^{1/2}}\tan ^{-1}\left[\left(\frac{2}{E}-1\right)^{1/2}\right]\sim \frac{{\rm\pi}}{\sqrt{2}}E^{-3/2}\quad \text{as}~E\rightarrow 0.\end{eqnarray}$$
4 The many-body problem
4.1 Numerical results
Our primary interest is in the dynamics of aggregation for a large number of spring–block elements. To gain some insight into the behaviour of these systems, we solved the system of differential equations (2.3)–(2.5) with up to
$N_{T}=2000$
elements. We used symmetry boundary conditions at the edge of the system, i.e.
$F_{-1}=F_{1}$
and
$F_{N_{T}+1}=F_{N_{T}-1}$
(Singh et al.
Reference Singh, Lister and Vella2014), though similar results were obtained with free edges,
$F_{0}=F_{N_{T}}=0$
. Using symmetry boundary conditions has the advantage that the uniform state
$h_{j}=1$
is an equilibrium. The initial condition was then taken to be a small random perturbation to this state, i.e.
$$\begin{eqnarray}h_{j}(0)=1+{\it\epsilon}\mathscr{R}_{j}\end{eqnarray}$$
with
$\mathscr{R}_{j}$
a random number drawn from
$N(0,1)$
. For the results presented here, we took
${\it\epsilon}=10^{-2}$
. The initial meniscus position was
$x_{j}(0)=1/h_{j}(0)$
so that each gap contained the same volume of liquid initially.
Two adjustments to the codes used by Singh et al. (Reference Singh, Lister and Vella2014) were made to account for the ultimate disappearance of liquid inherent in evaporation. First, the liquid level
$x_{j}(t)$
was prevented from becoming negative by enforcing
$x_{j}(t)\geqslant 10^{-10}$
for all times. Second, a short range van der Waals interaction was introduced to ensure that once elements come sufficiently close to contact they remain close, but do not pass through one another; this mimics the sticking upon contact that is observed in MEMS devices. Simulations were run until
$\sum _{j}|{\dot{h}}_{j}|<10^{-10}$
, so that each gap has reached an equilibrium to within numerical error.
Figure 4. Space–time diagrams showing the position of
${\approx}100$
blocks within a larger simulation of the elastocapillary aggregation of
$N_{T}=500$
blocks. The spring stiffness
$K=0.1$
in each case and the evaporation rate
$E$
is varied: (a)
$E=0$
; (b)
$E=0.1$
; (c)
$E=1$
; (d)
$E=3$
.
Figure 4 shows space–time diagrams for the evolution of around 100 block positions during aggregation. Results are shown for a range of values of the dimensionless evaporation rate
$E$
but with
$K=0.1$
fixed. Figure 4(a) (in the absence of evaporation) shows that clusters first form and then shrink down as time increases. Introducing a small amount of evaporation (figure 4
b), we see that the typical size of clusters is broadly similar but that the clusters seem to shrink to contact within a finite time. We also observe that there seem to be more scenarios in which clusters start to form but ultimately split: what we term a ‘cluster bifurcation’. (These cluster bifurcations are also associated with a ‘kink’ in the space–time diagram because the force balance is radically altered as the cluster splits: the new sub-clusters shift to positions where all of the spring forces are reduced.) Increasing the evaporation rate still further (figure 4
c,d), we see, that the clusters that form shrink to contact more quickly than with lower evaporation rates. We also see that there are many more cluster bifurcations and that the final clusters tend to be smaller than those observed with smaller evaporation rates.
Figure 5. The statistical properties of the cluster size
$N$
for a range of evaporation rates
$E$
, with
$K=0.0137$
fixed. Inset: the mean final cluster size,
$\langle N\rangle$
, as estimated, for each value of
$E$
, from 30 numerical experiments with initial conditions
$h_{j}=1+{\it\epsilon}\mathscr{R}_{j}$
, where
${\it\epsilon}=10^{-2}$
and
$\mathscr{R}_{j}$
are random numbers drawn from
$N(0,1)$
. The vertical error bars represent the standard error. The mean cluster size with no evaporation (
$E=0$
) is shown by the horizontal dashed line. Main figure: the probability distribution function for the cluster size as found in the numerical experiments. Different evaporation rates are indicated as follows:
$E=0$
(▸),
$E=10^{-2}$
(◂),
$E=0.5$
(▴),
$E=1$
(▾),
$E=5$
(▪) and
$E=30$
(●). Horizontal error bars show (for large
$E$
and small
$N$
) the discrete nature of the data, while vertical error bars show the standard error. The solid curve shows a normal distribution with mean
$1$
and standard deviation 0.3, as suggested previously (Singh et al.
Reference Singh, Lister and Vella2014); the dashed curve shows the distribution found by Boudaoud et al. (Reference Boudaoud, Bico and Roman2007) from a mean-field aggregation model with a single parameter fitted to match experiments.
Another perspective on the effect of evaporation is provided by examining the statistics of the cluster size
$N$
as
$E$
varies with
$K$
fixed. The inset of figure 5 shows that the mean cluster size
$\langle N\rangle$
is a decreasing function of the evaporation rate, as should be expected from the earlier observations of figure 4. The main portion of figure 5 shows the probability distribution function for a range of evaporation rates,
$0\leqslant E\leqslant 30$
, calculated by averaging the results of 30 runs at each evaporation rate. To allow for a comparison of the shape of the distribution as the evaporation rate varies, we have normalized the
$x$
-axis in figure 5 by
$\langle N\rangle$
. In these normalized coordinates, the typical (normalized) width of the distribution does not change noticeably as the evaporation rate increases (i.e. the width of the true distribution scales with the mean
$\langle N\rangle$
and, hence, does become narrower as
$E$
increases and
$\langle N\rangle$
decreases). Surprisingly the behaviour in the tails of the distribution seems to follow quite closely that found in the absence of evaporation (this is true even at
$E=30$
where only pairs of blocks or single blocks are observed). However, we note that the peak around
$N/\langle N\rangle \approx 0.8$
appears to become more pronounced as the evaporation rate
$E$
increases up to
$E=1$
; at the largest values of
$E$
the discrete nature of the distribution and the decrease in
$\langle N\rangle$
blunts the peak.
In summary, the simulation results shown in figures 4 and 5 suggest that the presence of evaporation modifies the dynamics of aggregation in three important ways: (i) any clusters that form shrink down more quickly than would be the case without evaporation, (ii) clusters are more prone to split up, or bifurcate, when evaporation is included and (iii) the final state generally consists of smaller clusters. The first of these observations recalls what was seen in the two-body problem where contact occurs in finite time with individual blocks moving at a constant speed, see (3.6). We shall see that the second and third observations are intimately related.
To better understand the third observation, numerical simulations with a range of values of
$K$
and
$E$
were performed and the maximum number of blocks in any cluster noted, see figure 3(b). This plot confirms the trend observed in (iii). In particular, we note that for sufficiently large evaporation rates
$E$
, no clusters are formed (the green points marked ‘1’ in figure 3
b). Comparing the transition between no clustering and clustering observed with many blocks (figure 3
b) to the transition between separated and stuck states observed in the two-body problem (figure 3
a) we observe a strong similarity; in particular, the line
$KE=1/2$
seems to be important in both transitions. Away from this line, the maximum number of elements in a cluster varies as the evaporation rate
$E$
varies (at fixed
$K$
). To try to understand these results, we begin by considering what happens in the limit of no springs.
4.2 No springs:
$K=0$
With no springs,
$K=0$
, the problem simplifies considerably. To see this, we note that (2.1) then gives
$$\begin{eqnarray}F_{j+1}-2F_{j}+F_{j-1}=0.\end{eqnarray}$$
If the system is unforced at its edges,
$F_{0}=F_{N+1}=0$
, then this may be inverted to give
$$\begin{eqnarray}F_{j}=\frac{x_{j}^{3}}{h_{j}^{3}}\frac{\text{d}h_{j}}{\text{d}t}+\frac{x_{j}}{h_{j}}=0\end{eqnarray}$$
for
$j=1,\ldots ,N$
. The motion within each gap therefore decouples from that in every other gap and we recover the two-block problem (3.1)–(3.2) with
$K=0$
. We conclude that in this case all of the blocks will aggregate with contact occurring at
$t=t_{\ast }$
with
$t_{\ast }$
given by (3.7).
4.3 Finite strength springs
$K>0$
Analytical progress is possible in the limit of no springs because the evolution of the liquid within each gap becomes decoupled from that in every other gap. However, with non-zero spring stiffness, the force balance (2.3) must hold; this couples the evolution within the gaps to one another.
4.3.1 Analytical results close to contact
To make analytical progress, we consider the limit of small gap thickness,
$h_{j}\ll 1$
, so that a cluster of blocks are already close to contact. In this limit additional displacements are small and the force arising from the springs is approximately constant, simplifying the problem considerably. In particular, when the blocks are close to contact,
$h_{j}\ll 1$
, the force-balance equation (2.3) may be approximated by
$$\begin{eqnarray}-2K=F_{j+1}-2F_{j}+F_{j-1}.\end{eqnarray}$$
Physically, this equation expresses the fact that in the small-gap limit, the net hydrodynamic force on a block,
$F_{j+1}-F_{j}$
, exceeds that on its neighbour by
$2K$
since it is approximately one unit further from its equilibrium position.
Given two particular
$F_{j}$
somewhere within the system, (4.4) may readily be inverted to give an expression for general
$F_{j}$
. Motivated by the repeated cluster-bifurcation events seen in figure 4(c,d), we consider whether a cluster of
$N$
spring–block elements (enclosing
$N-1$
liquid-filled gaps) that is close to complete collapse will bifurcate. Labelling the liquid gaps from one side of the cluster, we can approximate the hydrodynamic force at its edges by zero, i.e.
$F_{0}=F_{N}=0$
, since the cluster will be comparatively well separated from any neighbouring elements. Inverting (4.4) subject to these edge conditions, we find that
$$\begin{eqnarray}F_{j}=Kj(N-j).\end{eqnarray}$$
Thus the hydrodynamic force
$F_{j}$
in a gap is given by a constant value
$\mathscr{K}(j;N)\equiv Kj(N-j)$
that depends only on the label
$j$
of that gap within a particular cluster, and the number of blocks
$N$
in the cluster. (The detailed gap thicknesses do not appear because the displacements of blocks from their equilibrium positions, and the consequent spring forces, are nearly constant when the blocks are close to contact.)
Now, recalling from (2.1) that the hydrodynamic force
$F_{j}$
is the sum of lubrication and capillary components, we can write (4.5) together with the equation describing mass conservation as the system
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{x_{j}^{3}}{h_{j}^{3}}\frac{\text{d}h_{j}}{\text{d}t}=-\frac{x_{j}}{h_{j}}+\mathscr{K}(j;N) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}x_{j}}{\text{d}t}=-E-\frac{\text{d}h_{j}}{\text{d}t}\frac{x_{j}}{h_{j}}. & \displaystyle\end{eqnarray}$$
These equations are equivalent to those describing the evolution of a single gap, (3.1)–(3.2), in the limit of
$h\ll 1$
, but with an effective spring stiffness
$\mathscr{K}(j;N)$
. Hence, the many-body nature of the cluster problem enters only through the effective spring stiffness
$\mathscr{K}(j;N)$
, which varies through the cluster from relatively soft springs at the edge to relatively stiff springs in the centre. Within the cluster the presence of nearby elements increases the spring force acting to open an individual gap. This is reminiscent of the ‘collaborative stiffening’ that occurs when elastic plates are stuck together by surface tension: clusters interact just as individual elements do, albeit with a larger effective stiffness that comes from having to bend many plates (Bico et al.
Reference Bico, Roman, Moulin and Boudaoud2004). Collaborative stiffening here corresponds to a multi-block cluster whose springs have total stiffness
${\sim}NK$
and typical displacement proportional to
$N$
. Hence, the typical force in the cluster
${\sim}N^{2}K$
. In the two-body problem the displacement close to contact is unity and so, to achieve the same force, the effective stiffness satisfies
$\mathscr{K}\sim N^{2}K$
, as seen in (4.5) and (4.6).
In § 3.1 we showed for the two-body problem that if
$K<K_{\mathit{crit}}(E)$
the gap thickness will shrink to zero, and the blocks stick; otherwise the blocks will separate and return to their original positions. From the similarity of (4.6)–(4.7) to (3.1)–(3.2) at small gap thicknesses, it seems reasonable to expect that, for a cluster containing
$N$
blocks to collapse completely to a final stuck cluster, the effective stiffness
$\mathscr{K}(j;N)$
must be less than
$K_{\mathit{crit}}(E)$
for each of the gaps within that cluster, i.e.
$$\begin{eqnarray}\max _{j}\{\mathscr{K}(j;N)\}<K_{\mathit{crit}}(E).\end{eqnarray}$$
(We assume that any differences between the systems at large gap thicknesses have no significant effect on
$K_{\mathit{crit}}(E)$
.)
Noting that
$\max _{j}\{\mathscr{K}(j;N)\}=KN^{2}/4$
, we can then restate the criterion in (4.8) as collapse of an
$N$
-block cluster (for
$N\geqslant 2$
) can only occur if
$$\begin{eqnarray}N<N_{\mathit{max}}=\left(\frac{4K_{\mathit{crit}}(E)}{K}\right)^{1/2}.\end{eqnarray}$$
Thus, if
$N_{\mathit{max}}<2$
, which occurs at sufficiently large spring stiffnesses, all blocks end up separate (in a trivial
$N=1$
cluster). More generally, equation (4.9) is a prediction for the maximum number of blocks in a cluster as a function of the evaporation rate
$E$
and spring stiffness
$K$
. The only assumption used in its derivation was that whether a cluster collapses, and the blocks stick, or splits up should be determined by the behaviour of the system with small gap thicknesses,
$h_{j}\ll 1$
.
Figure 6. Rescaling of the data in figure 3(b) to test the theoretical prediction (4.9) for the maximum number of blocks in a cluster. (a) The maximum number of blocks in any cluster,
$N_{\mathit{max}}$
, as a function of
$K_{\mathit{crit}}(E)/K$
(points); equation (4.9) suggests that, provided
$N_{\mathit{max}}>1$
, this should follow the behaviour
$y=2x^{1/2}$
(solid line). (b) An alternative rescaling showing the dependence on the evaporation rate,
$E$
, directly for data with
$N_{\mathit{max}}\geqslant 2$
. Here, the solid curve shows
$2K_{\mathit{crit}}(E)^{1/2}$
, the dashed horizontal line shows
$N_{\mathit{max}}K^{1/2}=2^{7/2}{\rm\pi}/9\approx 3.95$
(predicted by Singh et al. (Reference Singh, Lister and Vella2014), in the limit
$E=0$
) and the dotted line shows
$N_{\mathit{max}}K^{1/2}=(2/E)^{1/2}$
, valid for
$E\gtrsim 0.15$
. In both panels, colours and symbols correspond to results obtained at different evaporation rates:
$E=10^{-3}$
(red squares),
$E\approx 10^{-2}$
(green triangles),
$E\approx 0.1$
(blue circles),
$E\approx 1$
(black inverted triangles),
$E=5$
(magenta ‘
$\times$
’),
$E=10$
(cyan ‘
$+$
’) and
$E=10^{2}$
(yellow side triangles).
4.3.2 Reassessment of numerical results
Figure 6 replots data for the maximum number of blocks in any cluster from figure 3(b), at seven different evaporation rates
$E$
and with a variety of spring stiffnesses
$K$
(i.e. the values of
$N_{\mathit{max}}$
along seven different horizontal lines in figure 3
b). In figure 6, this data has been rescaled based on the analytical prediction (4.9): figure 6(a) highlights how the raw values of
$N_{\mathit{max}}$
vary with
$K_{\mathit{crit}}(E)/K$
while figure 6(b) highlights the dependence on
$E$
.
We note first that both rescalings in figure 6 lead to reasonable collapses of the data from a two-dimensional parameter space (figure 3
b) onto a single master curve. Since both of these rescalings follow naturally from (4.9), we conclude that the reduction from the many-body problem to the two-body problem yields useful insight. However, we also note that for small evaporation rates,
$E\ll 1$
, the maximum cluster size
$N_{\mathit{max}}$
lies somewhat below that predicted on the basis of (4.9). Instead, it seems that for
$E\ll 1$
the results are closer to the prediction from the no-evaporation problem,
$E=0$
, (Singh et al.
Reference Singh, Lister and Vella2014) where
$$\begin{eqnarray}N_{\mathit{max}}K^{1/2}=\frac{2^{7/2}{\rm\pi}}{9}\approx 3.95.\end{eqnarray}$$
Equation (4.10) was obtained from an analysis of front propagation in an unstable medium, which was found to leave behind clusters of this well-defined size. Though larger clusters would, in principle, be able to collapse and not bifurcate, they are apparently not formed by the initial front-propagation mechanism. For all values of
$E$
, we find that clusters do not reach the maximum size predicted by the linear instability of a nearly uniform state without evaporation, which in this notation reads
$N_{\mathit{max}}K^{1/2}=2{\rm\pi}$
(Singh et al.
Reference Singh, Lister and Vella2014).
5 Discussion
We have studied the effect of evaporation on a simple model of elastocapillary aggregation, which revealed two important differences from the same model in the absence of evaporation. First, evaporation acts to bring the elements together more quickly and leads to contact within a finite time. Second, the maximum number of elements in a final cluster decreases with increasing evaporation rate.
The first of these observations is important for applications since it suggests that evaporation can drive objects together more quickly than would be the case without evaporation, thereby increasing the risk of, for example, MEMS stiction (Tanaka et al. Reference Tanaka, Morigami and Atoda1993). However, the practical consequences may be (partially) mitigated by the second observation: if high evaporation rates cause fewer elements to come into contact, then those elements that do come into contact are less likely to be deformed enough to fracture (which is one of the main manufacturing concerns, see figure 1 a).
Our analysis of the mathematical model allowed us to reduce the problem of a cluster that is close to collapsed (i.e. small gap thicknesses) to a two-body problem with an effective spring stiffness. If a cluster starts to form then this effective stiffness is smallest at the edges and largest at the centre of the cluster; thus, if the cluster breaks up (because it is too large), then we expect a void to form in the centre. This is similar to observations of microscopic pillars made up of many carbon nanotubes: if the diameter of the initial pillar (cluster) is small enough then the entire pillar will collapse during evaporation; if not, then voids form in the centre (de Volder & Hart Reference de Volder and Hart2013).
We have shown that the number of elements in a cluster is controlled by the dynamics of evaporation. As such, the final pattern observed depends on the dynamics of pattern formation and not simply which pattern is the global energy minimizer, as is often assumed. However, we have also seen that the effect of evaporation only becomes significant once
$E\gtrsim 0.1$
, i.e. when the evaporation dynamics occur on a time scale comparable to, or faster than, that of the lubrication-type dynamics itself. From (2.6) the dimensionless evaporation rate
$E=2\mathit{Ca}_{e}(x_{0}/w)$
, where
$\mathit{Ca}_{e}={\it\mu}e/{\it\gamma}$
is an evaporative capillary number. Written in this way, we see that
$E$
is the product of a property of the liquid filling the gaps (
$\mathit{Ca}_{e}$
) and the initial aspect ratio of the liquid gaps (
$x_{0}/w$
). For the evaporation of pure, still water at room temperature and ambient humidity,
$e\approx 60~\text{nm s}^{-1}$
(Li et al.
Reference Li, Cabane, Sztucki, Gummel and Goehring2012; Routh Reference Routh2013),
${\it\mu}=10^{-3}~\text{Pa s}$
,
${\it\gamma}=72~\text{mN m}^{-1}$
so that
$\mathit{Ca}_{e}\approx 10^{-9}$
. Other liquids evaporate more quickly and hence may have a larger value of
$\mathit{Ca}_{e}$
. For example, HFE-7100 (3M) has an evaporation rate
$e\approx 20~{\rm\mu}\text{m s}^{-1}$
(see Lyulin & Kabov Reference Lyulin and Kabov2013; Machrafi et al.
Reference Machrafi, Sadoun, Rednikov, Dehaeck, Dauby and Colinet2013, for example),
${\it\gamma}=13.6~\text{mN m}^{-1}$
and
${\it\mu}=6.1\times 10^{-4}~\text{Pa s}$
so that
$\mathit{Ca}_{e}\approx 10^{-5}$
. Such small values of the evaporative capillary number mean that to obtain
$E\gtrsim 0.1$
would require
$x_{0}/w\gtrsim 10^{4}$
. While this is at least consistent with the requirement
$x_{0}/w\gg 1$
for lubrication theory to be valid, it is unlikely ever to be the case for the micrometre-scale plates and pillars often considered experimentally (see Wei et al.
Reference Wei, Schneider, Kim, Kim, Aizenberg and Mahadevan2015, for example). However, such large aspect ratios are possible for nanotube forests where the typical length of a nanotube is on the order of
$100~{\rm\mu}\text{m}$
but the radius, and hence the typical separation, is
$O(10~\text{nm})$
(Chakrapani et al.
Reference Chakrapani, Wei, Carrillo, Ajayan and Kane2004).
The geometrical enhancement of the evaporation rate suggests that our prediction that clusters typically contain fewer elements as the evaporation rate increases may be observable in sufficiently dense ‘forests’. Indeed, from their experiments on the formation of cells in carbon nanotube forests (see figure 1 b), Chakrapani et al. (Reference Chakrapani, Wei, Carrillo, Ajayan and Kane2004) report that ‘[a faster rate of evaporation] favors crack growth and results in a decrease in the average cell width’ (i.e. fewer elements form in the clusters between cracks). However, they do not present quantitative data to support this comment.
We hope that the theoretical results presented here may motivate experimental efforts to better understand how evaporation rate influences pattern formation in these systems. At the same time, there are a number of improvements that should be made to the model developed here. For example, experiments are bound to take place in a three-dimensional setting where the force from surface tension will require a more detailed calculation (Wei et al. Reference Wei, Schneider, Kim, Kim, Aizenberg and Mahadevan2015) and the resistive forces from the lubricating flow may be reduced since liquid can flow horizontally, as well as vertically.
Even in situations where the delicate balance between evaporation and lubrication dynamics discussed here is not appropriate, other, slower timescales exist with which evaporation can more readily compete. For example the viscoelastic relaxation time of the elements themselves may be important (Wei et al. Reference Wei, Schneider, Kim, Kim, Aizenberg and Mahadevan2015). We expect that the techniques developed in this paper may be applied, with suitable modifications, to understand the interplay between evaporation and dynamics in such systems.
Acknowledgements
K.S. and D.V. wish to acknowledge the support of the King Abdullah University of Science and Technology (KAUST; Award No. KUK-C1-013-04), and the John Fell Oxford University Press (OUP) Research Fund.
Appendix A. Derivation of model
In this appendix, we outline the derivation of the expression for the force exerted by the
$j$
th liquid gap on each of its neighbouring blocks (2.1). More details are given by Singh et al. (Reference Singh, Lister and Vella2014). The well-known lubrication equation, together with the local conservation of mass leads to a relationship between the pressure profile within the gap,
$p_{j}(x)$
, and the rate at which the gap thickness,
$h_{j}$
, is changing with time. In particular,
$$\begin{eqnarray}\frac{\text{d}h_{j}}{\text{d}t}=\frac{h_{j}^{3}}{12{\it\mu}}\frac{\partial ^{2}p_{j}}{\partial x^{2}}.\end{eqnarray}$$
This equation for
$p_{j}$
is to be solved subject to the boundary conditions of no flux through the base, i.e.
$\partial p_{j}/\partial x|_{x=0}=0$
, and that the pressure just beneath the meniscus is equal to the capillary pressure, i.e.
$p_{j}(x_{m}^{-})=-2{\it\gamma}/h_{j}$
. We find that
$$\begin{eqnarray}p_{j}(x)=\frac{6{\it\mu}}{h_{j}^{3}}\frac{\text{d}h_{j}}{\text{d}t}(x^{2}-x_{j}^{2})-\frac{2{\it\gamma}}{h_{j}}.\end{eqnarray}$$
The force on each of the blocks that bound the
$j$
th gap is then given by
$$\begin{eqnarray}F_{j}=-\int _{0}^{x_{j}}p_{j}\,\text{d}x=\frac{4{\it\mu}x_{j}^{3}}{h_{j}^{3}}\frac{\text{d}h_{j}}{\text{d}t}+\frac{2{\it\gamma}x_{j}}{h_{j}},\end{eqnarray}$$
which is the dimensional version of (2.1).
Appendix B. Phase-plane analysis to determine
$K_{\mathit{crit}}(E)$
For
$0<KE<1/2$
the phase plane of (3.3)–(3.4) is topologically equivalent to figure 2(a): trajectories generically tend either to the stable node
$(0,1)$
(separated blocks) or to the stable node
$({\it\xi}_{\ast },0)$
(stuck blocks). The dividing trajectory between these outcomes is the stable manifold of the saddle point
$({\it\xi}_{s},0)$
, which is given by
${\it\xi}\sim {\it\xi}_{s}+2h{\it\xi}_{s}/({\it\xi}_{s}-3)+O(h^{2})$
as
$h\rightarrow 0$
. Starting with this limit, we can integrate back along the stable manifold to
$h=1$
, which defines a point
$(K_{\mathit{crit}}^{-1},1)$
. Then if
$K\leqslant K_{\mathit{crit}}$
an initial condition
$(K^{-1},1)$
leads to a stuck state.
As
$KE{\nearrow}1/2$
the points
$({\it\xi}_{s},0)$
and
$({\it\xi}_{\ast },0)$
merge into a saddle-node at
$(2,0)$
and
$K_{\mathit{crit}}=3.30+O(1-2KE)^{1/2}$
. For
$KE>1/2$
all trajectories lead to the separated state.
At
$KE=0$
the (non-trivial) nullclines of (3.3)–(3.4) coincide as a line of fixed points
$h+{\it\xi}=1$
, corresponding to a mechanical equilibrium
$K(1-h)=x/h$
between the spring and capillary forces in the absence of evaporation. Away from this equilibrium, the trajectories satisfy
$h^{2}{\it\xi}=\text{const.}$
, corresponding to mass conservation
$hx=\text{const}$
.
To consider the useful limit of slow but non-zero evaporation, we let
${\it\epsilon}=KE\ll 1$
and note that the stable manifold of the saddle point,
${\it\xi}=1+{\it\epsilon}/2-h(1+3{\it\epsilon}/4)+O({\it\epsilon}^{2},h^{2})$
, starts almost parallel to the two nearly coincident nullclines. This motivates setting
${\it\eta}=(h+{\it\xi}-1)/{\it\epsilon}$
so that (3.3)–(3.4) become
$$\begin{eqnarray}\frac{\text{d}{\it\xi}}{\text{d}s}={\it\epsilon}{\it\xi}[2{\it\eta}-{\it\xi}^{2}],\quad \frac{\text{d}{\it\eta}}{\text{d}s}=(3{\it\xi}-1){\it\eta}-{\it\xi}^{3}-{\it\epsilon}{\it\eta}^{2}.\end{eqnarray}$$
Since
${\it\epsilon}\ll 1$
, we expect the evolution of
${\it\xi}$
to be slow and hence, integrating backwards,
${\it\eta}$
collapses onto a ‘slow’ manifold given by
$\text{d}{\it\eta}/\text{d}s=O({\it\epsilon})$
, or
$$\begin{eqnarray}{\it\eta}\sim \frac{{\it\xi}^{3}}{3{\it\xi}-1},\end{eqnarray}$$
over the range
$1/3<{\it\xi}\leqslant 1$
. As
${\it\xi}$
approaches
$1/3$
,
${\it\eta}$
diverges and we need to rescale the equations to understand how the stable manifold leaves the
$O({\it\epsilon})$
neighbourhood of
$h+{\it\xi}=1$
.
A suitable rescaling is to set
$x=3(3{\it\xi}-1)/2{\it\epsilon}^{1/3}$
,
$y=9{\it\eta}{\it\epsilon}^{1/3}=9({\it\xi}+h-1)/{\it\epsilon}^{2/3}$
. At leading order (B 1) becomes
$$\begin{eqnarray}\frac{3}{{\it\epsilon}^{1/3}}\frac{\text{d}x}{\text{d}s}=y,\quad \frac{3}{{\it\epsilon}^{1/3}}\frac{\text{d}y}{\text{d}s}=2xy-1,\quad \text{and thus}\quad \frac{\text{d}y}{\text{d}x}=2x-\frac{1}{y}.\end{eqnarray}$$
For the solution to (B 3) to match onto (B 2), we require
$y\sim 1/2x$
as
$x\rightarrow \infty$
. This condition specifies a unique solution to (B 3), given implicitly by
$x=-\text{Ai}^{\prime }(z)/\text{Ai}(z)$
, where Ai denotes the Airy function and
$z=x^{2}-y$
. As
$x\rightarrow -\infty$
,
$y\sim x^{2}-z_{0}$
, where
$z_{0}=-2.338$
is the first zero of
$\text{Ai}(z)$
. This asymptotic behaviour matches onto the trajectory
$h^{2}{\it\xi}=(4/27)(1-{\it\epsilon}^{2/3}z_{0}/3)+O({\it\epsilon})$
, which continues to
$(K_{\mathit{crit}}^{-1},1)$
. Hence,
$K_{\mathit{crit}}\sim (27/4)(1+{\it\epsilon}^{2/3}z_{0}/3)\sim (27/4)-18.8E^{2/3}$
. The limiting value
$K_{\mathit{crit}}=27/4$
for
$E=0$
was obtained by Singh et al. (Reference Singh, Lister and Vella2014).
