Film deposition and dynamics of a self-propelled wetting droplet on a conical fibre

Abstract

We study the dynamic wetting of a self-propelled viscous droplet using the time-dependent lubrication equation on a conical-shaped substrate for different cone radii, cone angles and slip lengths. The droplet velocity is found to increase with the cone angle and the slip length, but decrease with the cone radius. We show that a film is formed at the receding part of the droplet, much like the classical Landau–Levich–Derjaguin film. The film thickness $h_f$ is found to decrease with the slip length $\lambda$. By using the approach of matching asymptotic profiles in the film region and the quasi-static droplet, we obtain the same film thickness as the results from the lubrication approach for all slip lengths. We identify two scaling laws for the asymptotic regimes: $h_fh''_o \sim Ca^{2/3}$ for $\lambda \ll h_f$ and $h_f h''^{3}_o\sim (Ca/\lambda )^2$ for $\lambda \gg h_f$; here, $1/h''_o$ is a characteristic length at the receding contact line and $Ca$ is the capillary number. We compare the position and the shape of the droplet predicted from our continuum theory with molecular dynamics simulations, which are in close agreement. Our results show that manipulating the droplet size, the cone angle and the slip length provides different schemes for guiding droplet motion and coating the substrate with a film.

1. Introduction

Coating a film onto a substrate as a liquid is forced to move along it is a technique used in painting and industrial applications such as lithography, which has been studied since the early twentieth century (Quéré 1999). Dip coating is one way to coat a plate as it is withdrawn from a liquid reservoir above a critical plate velocity (Snoeijer et al. 2006; Maleki et al. 2011; Gao et al. 2016). A mathematical model describing this film coating was developed in the seminal works by Landau & Levich (1942) and Derjaguin (1943). These theoretical works sparked a great interest in film coating, later adopted for a range of solid geometries, e.g. cylindrical fibres (White & Tallmadge 1966; Wilson 1988; De Ryck & Quéré 1996), coating by rollers (Taylor 1963; Wilson 1982) and coating the inner surface of a channel/tube (Bretherton 1961; Tabeling & Libchaber 1986). Some studies have focused on how other physical effects influence the film deposition such as gravity (Derjaguin 1943; Snoeijer et al. 2008), inertia (De Ryck & Quéré 1996; Orsini & Tricoli 2017), surfactants (Carroll & Lucassen 1973), particles on the interface (Colosqui, Morris & Stone 2013; Dixit & Homsy 2013a,b), van der Waals forces for deposited films of nanometric scales (Quéré, di Meglio & Brochard-Wyart 1989) as well as effects of substrate roughness (Krechetnikov & Homsy 2005) and confinement due to the reservoir (Kim & Nam 2017). Much of the extensive literature on the film deposition dynamics has been summarized in several review articles (Ruschak 1985; Quéré 1999; Weinstein & Ruschak 2004; Rio & Boulogne 2017).

The classical theory by Landau & Levich (1942) and Derjaguin (1943) gives a fundamental description of thin film coating, where the deposited film thickness is so thin that gravity can be neglected. The flow inside the film region is maintained by the balance of capillarity, characterized by the liquid/air surface tension coefficient $\gamma$, and the viscous forces, characterized by the liquid viscosity $\eta$. The film region is connected to a quasi-static liquid reservoir of a length scale that is much larger than the thickness of the deposited film. When a plate is withdrawn from a reservoir, this length is set by the capillary length $\ell _c\equiv (\gamma /\rho g)^{1/2}$, with $\rho$ the liquid density and $g$ the gravitational acceleration. By using the method of asymptotic matching, the thickness of the film $h_f$, denoted as the Landau–Levich–Derjaguin (LLD) film, is shown to have a universal scaling with respect to the plate velocity $U$ as $h_f/\ell _c\sim Ca^{2/3}$, where the capillary number $Ca\equiv \eta U/\gamma$ is the ratio between the viscous and the surface tension forces. Remarkably, this $Ca^{2/3}$ power law has been demonstrated to be a robust relation in many different systems when a fluid film is deposited. The only required change in the scaling relation is to replace the capillary length by the corresponding characteristic length of the system. For example, the film thickness is rescaled by the fibre radius for the case when a cylindrical fibre is withdrawn from a bath (White & Tallmadge 1966; James 1973; Wilson 1988); while in the case of a long bubble moving in a tube, also known as Bretherton's problem, the corresponding length is the tube radius (Bretherton 1961).

Common to fluid coating processes is that they often require an external driving force to displace the fluid. However, when a droplet with a size smaller than the capillary length comes in contact with a conical fibre, it moves spontaneously from the tip to the base of the cone due to capillarity (Lorenceau & Quéré 2004; Li & Thoroddsen 2013). In nature, this self-propelled mechanism has been exploited by plants (Liu et al. 2015) and animals (Zheng et al. 2010; Wang et al. 2015) to facilitate water transport at small scales. When a droplet is translating above a critical velocity, a layer of liquid film is expected to be deposited on the conical surface at the receding part of the droplet. It has been discovered recently on the trichome of the Sarrancenia that the deposited film provides a wetted surface, enabling later water droplets to be transported at a velocity several orders of magnitude larger than found in other plants (Chen et al. 2018). Despite the importance of understanding the film deposition and potential implications for biological evolution in plants and giving a path to very fast droplet transport, the film deposition has not been studied before on conical geometries. Previous fluid coating studies have assumed a no-slip condition at the fluid–solid boundary, i.e. no relative motion between the fluid and the solid boundary. Interestingly, slip lengths have been reported to be as large as a few micrometres for fluids such as polymer melts (Bäumchen, Fetzer & Jacobs 2009) and for superhydrophobic surfaces (Rothstein 2010). When the droplet size is decreased to a few micrometres or below, effects due to the fluid slip on solid surface may become significant, as demonstrated in dynamical fluidic systems (Lauga, Brenner & Stone 2007; Bocquet & Charlaix 2009) such as the dewetting of microdroplets (McGraw et al. 2016; Chan et al. 2017) and of liquid films (Fetzer et al. 2005). However, the influence of slip on the droplet dynamics and film deposition is not known for the directional droplet motion on a cone.

Although conical solid structures are common in nature and appear as a component in industrial processes, modelling of the droplet fluid flow on such geometries is lacking. In our recent article (Chan, Yang & Carlson 2020), we provide a physical picture of the spontaneous motion of the droplet based on the mismatch between the equilibrium contact angle and the apparent contact angles. This generates flow in the contact line regions and maintains the droplet motion. In this study, we implement the time-dependent lubrication equation developed in Chan et al. (2020) to investigate the evolution of the liquid–air interface of the capillary driven droplet motion on a smooth conical fibre. The properties of a deposited film generated by a self-propelled droplet are studied for small cone angles and for a wide range of slip lengths. Apart from the continuum approach, the simple geometry of a conical shape allows us to study fluid flow using molecular dynamics (MD) simulations. Results from the MD simulations will be used as a verification for the droplet shape predicted by the lubrication model. In fact, the approach of MD simulations has previously been implemented to study the wetting dynamics at the nanoscale (Nakamura et al. 2013), the slip condition at a contact line region (Qian, Wang & Sheng 2003), the frictional force on a sliding droplet (Koplik 2019) and the influence of physico-chemistry of the water/substrate interface on the droplet dynamics (Johansson, Carlson & Hess 2015).

2. Mathematical formulation

An axisymmetric viscous droplet with a volume $V$ is placed in contact with a wetted surface of a conical fibre with a cone angle $\alpha \ll 1$, see figure 1(a). We consider a fibre surface prewetted with a thin layer of the same fluid of thickness $\epsilon$. The prewetted layer can be deposited or interpreted as a microscopic precursor film for a perfectly wetting droplet, i.e. equilibrium contact angle $\theta _e=0^{\circ }$. The profile of the liquid–air interface is described by $h(r,t)$, the distance between the interface and the substrate, as a function of the distance from the vertex of the cone along its surface $r$ and time $t$. For droplets with a Bond number $Bo\equiv \rho g V^{2/3}/\gamma \ll 1$, gravitational effects can be ignored. We consider the Reynolds number $Re\equiv \rho U V^{1/3}/\eta \ll 1$ and the flow inside the droplet is described by the Stokes equations and the continuity equation.

Figure 1. $(a)$ A description of the system at study, where a droplet is moving across (from left to right) a conically shaped fibre (grey shaded region). Droplet profiles $h(r,t)$ are shown at three different times on the fibre with a cone angle $\alpha$. From $t=0$ to $t=0.89$, there is a fast relaxation of the droplet before it slowly spreads across the fibre. $(b)$ The centre of mass of the droplet $x_c$ plotted as a function of time. The diamond markers correspond to the two profiles shown in $(a)$, with a cone angle $\alpha =0.01\ \textrm {rad}$, a prewetted film thickness $\epsilon =10^{-3}$ and a slip length $\lambda =0$.

2.1. Lubrication approximation on a cone (LAC)

Consider the flow in the droplet as $\boldsymbol {u}(r,\theta )$, here $\theta$ is the polar angle measured from the axis of rotation. Supposing the polar angle of the free surface of the droplet is very small, the flow is primarily in the radial direction. By using these approximations, the Stokes equations reduce to the lubrication equations here given in spherical coordinates (Chan et al. 2020),

(2.1)\begin{gather} \frac{\partial p}{\partial r}=\frac{\eta}{r^2\theta}\frac{\partial}{\partial \theta}\left(\theta \frac{\partial u}{\partial \theta}\right), \end{gather}

(2.2)\begin{gather} \frac{\partial p}{\partial \theta}=0, \end{gather}

where $p$ is the pressure and $u$ is the radial velocity inside the droplet/film.

To describe the fluid flow (2.1) and (2.2) need to be accompanied by several boundary conditions. At the liquid–air surface, the tangential stress is zero as we neglect viscous effects in the air

(2.3)\begin{equation} \frac{\partial u}{\partial \theta}=0\quad \text{at}\ \theta=\alpha+\phi. \end{equation}

At the wetted substrate, the normal velocity is zero and we assume a radial velocity described by the Navier-slip condition (Lauga, Brenner & Stone 2008)

(2.4)\begin{equation} \frac{u}{\lambda}=\frac{1}{r}\frac{\partial u}{\partial \theta} \quad \text{at} \ \theta=\alpha, \end{equation}

where $\lambda$ is the slip length.

Solving (2.1) and (2.2) with the boundary conditions (2.3) and (2.4) gives the velocity, and by imposing mass conservation of the liquid we get

(2.5)\begin{align} \frac{\partial h}{\partial t}&+\frac{1}{r\alpha+ h}\frac{\partial}{\partial r} \left[ \frac{r^4\alpha^4}{2\eta}\frac{\partial p}{\partial r}\left\{\frac{1}{8}\left[3\left(1+\frac{h}{r\alpha}\right)^4-4\left(1+\frac{h}{r\alpha}\right)^2+1\right] \right.\right.\nonumber\\ &\quad -\left.\left.\frac{1}{2}\left(1+\frac{h}{r\alpha}\right)^4\ln\left(1+\frac{h}{r\alpha}\right)-\frac{\lambda h^2}{2r^3\alpha^3}\left(2+\frac{h}{r\alpha}\right)^2\right\} \right]=0. \end{align}

The pressure gradient inside the liquid is generated by the Laplace pressure $p=-\gamma \kappa$ where $\kappa$ is the curvature of the liquid–air interface, which for $\alpha \ll 1$ simplifies as

(2.6)\begin{equation} \kappa=\frac{h''}{\left(1+h^{'2}\right)^{3/2}}-\frac{1-\alpha h'}{(r\alpha +h )\left(1+h'^{2}\right)^{1/2}}, \end{equation}

with $()'\equiv {\partial ()}/{\partial r}$. The second term of the curvature is derived by using a rotation matrix with the cone angle $\alpha \ll 1$. We keep the $h'$ terms as the interface slope is not always small at the droplet scale. We will see in § 3.1 that the viscous effect is weak in the droplet and it quickly adopts a quasi-static shape at the leading order, which is determined by the uniform pressure condition i.e. $\kappa = \textrm {constant}$. Although the flow field computed from the lubrication equation is inaccurate at the droplet scale, the correct quasi-static shape determines the flows in the contact line regions where lubrication approximation does work. Hence, (2.5) is still valid for computing the evolution of the interface.

2.2. Finite element method

We solve a coupled system of equations consisting of (2.5) and the Laplace pressure equation $p=-\gamma \kappa$ numerically by using the finite element method. For the pressure equation, we split it into the two following equations:

(2.7)\begin{equation} p=-\gamma\left[\frac{q'}{\left(1+h^{'2}\right)^{3/2}}-\frac{1-\alpha h'}{(r\alpha +h )\left(1+h'^{2}\right)^{1/2}}\right] \end{equation}

and

(2.8)\begin{equation} q=h'. \end{equation}

The variables we solve for are $h(r,t)$, $p(r,t)$ and $q(r,t)$. These fields are discretized with linear elements and solved as a coupled equation set by using Newton's method in the FEniCS library (Logg, Mardal & Wells 2012). We use both an adaptive time stepping routine and an adaptive spatial discretization to refine the spatial resolution around the receding tail and the advancing front of the droplet with a resolution of $\Delta r=10^{-4} V^{1/3}$; here, $\Delta r$ is the difference of $r$ between two nodal points. The numerical simulations are initialized with the initial profile $h(r,t=0)=\epsilon +A[1-\tanh (r-r_i)^2]$ where $A$ determines the volume of the droplet and $r_i$ determines the initial position of the droplet's geometric centre. The simulations are insensitive to the initial droplet shape after a very short initial relaxation, see figure 1 and appendix B. Further, we impose the following boundary conditions at the boundary $\partial \varOmega$ of the numerical domain: $p(r=\partial \varOmega ,t)=p(r=\partial \varOmega ,t=0)$ and $\partial h(r=\partial \varOmega , t)/\partial r= 0$.

2.3. Molecular dynamics simulations

We can test our hydrodynamic model by means of a ‘numerical experiment’ – a classical MD simulation of a liquid drop placed on a solid cone, based on standard methods (Frenkel & Smit 2002). We consider a generic viscous liquid consisting of spherically symmetric atoms with a Lennard–Jones (LJ) interaction, bound into linear tetramer molecules by a FENE (finitely extensible nonlinear elastic) potential

(2.9a,b)\begin{equation} V_{LJ}(\mathcal{X}) = 4E \left[\left(\frac{\mathcal{X}}{\sigma} \right)^{-12} - \left(\frac{\mathcal{X}}{\sigma} \right)^{-6}\right] \quad V_{FENE}(\mathcal{X})=-\frac{1}{2} k_F \mathcal{X}_0^2 \ln\left(1-\frac{\mathcal{X}^2}{\mathcal{X}_0^2}\right), \end{equation}

where $\mathcal {X}$ is the separation between the centre of mass of two atoms. The LJ potential acts between all pairs of atoms within a cutoff distance $2.5\sigma$, and is shifted by a linear term so that the force vanishes at the cutoff. The FENE interaction (with parameters $k_F=30E/\sigma ^2$ and $\mathcal {X}_0=1.5\sigma$), following Grest & Kremer (1986) acts between adjacent atoms on the chain. The advantage of a molecular rather than a monatomic liquid is that the vapour pressure is very low and the liquid/vapour interface is relatively sharp and easy to visualize. The solid is a conical section of a regular lattice whose atoms are mobile but bound to their lattice sites by linear springs with stiffness $100E/\sigma ^2$. The simulations are conducted in a canonical/NVT ensemble, where the temperature is fixed at $0.8E/k_B$ using a Nosé–Hoover thermostat. This particular solid/liquid system has been used in a number of previous simulations (Busic, Koplik & Banavar 2003; Koplik et al. 2006; Koplik & Zhang 2013; Koplik & Maldarelli 2017), and its properties are well characterized. The liquid has bulk number density $0.857\sigma ^{-3}$, viscosity $5.18m/(\sigma \tau )$ and liquid–vapour surface tension $0.668E/\sigma ^2$, where $m$ is the mass of the liquid atoms and $\tau =\sigma (m/E)^{1/2}$ is the natural time scale based on the LJ parameters. Furthermore, the liquid is completely wetting: a drop placed on flat solid surface with the same density and interactions spreads completely into a thin film.

The simulation begins with all atoms on face-centred cubic (fcc) lattice sites, within a rectangular box of length $256.5\sigma$ and sides $171\sigma$, with repulsive confining walls in the long ($x$) direction and periodic boundary conditions on the sides. For the cone we select all atoms in an fcc lattice of number density $1.06\sigma ^{-3}$ within a radius $R_{co}(x)=3+x\tan \alpha$ (in the unit of $\sigma$) of the central axis, which runs in the $x$-direction through the centre of the box, and where $\alpha =0.1\ \textrm {rad}$. The resulting solid has 74,362 atoms. The liquid initially occupies a disc-shaped region near the left edge, $x_1<x<x_2$ and $R_{co}(x)<\mathcal {R}<R_{out}$ consisting of all atoms outside the cone but inside an outer radius $R_{out}$. We have studied two cases (A, B) where the remaining cone surface is initially dry and one case (C) where there is also a liquid (prewetting) film of thickness $4\sigma$. The parameters for the various cases and the number of fluid atoms $N$ is given in table 1. The simulation temperature starts at a low value in the solid phase and increases linearly to the final value, $0.2\to 0.8E/k_B$ over $250\tau$, to prevent the liquid atoms from leaving the cone, and subsequently the drop is allowed to evolve freely at the final temperature.

Table 1. Initial geometry of the drops.

3. Results and discussions

3.1. Numerical solutions of the LAC

We first present numerical simulation results for the typical evolution of the droplet profile obtained from the lubrication approximation on a cone. In the following, all lengths are rescaled by $V^{1/3}$ and time is rescaled by $V^{1/3}\eta /\gamma$. The dimensionless parameters are the cone angle $\alpha$, the thickness of the prewetted layer $\epsilon$ and the slip length $\lambda$.

A typical dynamical process is shown in figure 1(a). First the droplet relaxes from an initial shape ($t=0$) to a quasi-static shape ($t=0.89$) in a short time. At $t=0$, the initial shape gives a non-uniform curvature and hence a non-uniform pressure inside the whole droplet region. The pressure gradient generates flow inside the droplet. At $t=0.89$, a nearly uniform pressure distribution is achieved in the bulk of the droplet, but a large pressure gradient is created at the two edges of the droplet, which are commonly referred to as the ‘contact line regions’, see figure 2(b) for the pressure distribution at $t=256$. A concentration of stresses at the contact line would be expected (Huh & Scriven 1971). After the quick initial relaxation, the droplet starts to propagate toward the thicker part of the cone. The position of the droplet is described by the centre of mass of the droplet, for $\alpha \ll 1$, defined as

(3.1)\begin{equation} x_c={\rm \pi}\int^{r_a}_{r_r}h(h+2\alpha r) r\,\textrm{d}r, \end{equation}

which is plotted as a function of time in figure 1(b). Here $r_r$ is the apparent receding and $r_a$ is the apparent advancing contact line positions, which are defined in appendix A. The numerical simulations suggest that the droplet adopts a quasi-static shape during the directional spreading. For example, the profile for $t=256$ is plotted as the black solid line in figure 2(a), and the pressure distribution is shown in figure 2(b). Given the uniform pressure/curvature condition, one can solve for a static droplet profile, see the details of the computation in § 3.5. The static profile obtained in this way is plotted as the red dashed curve in figure 2(a) for the same droplet position as obtained from the LAC at $t=256$. The agreement between the two profiles again demonstrates that the droplet profile is quasi-static on the droplet scale. As the droplet shape appears more round than flat, it may affect the validity of the lubrication approximation in the contact line regions where the viscous effects are significant. We zoom into the advancing and receding contact line regions and compute the pressure gradient $p'$ rescaled by its maximum magnitude and the interfacial slope $h'$, which are shown in figures 2(c) and 2(d). We observe that when approaching the contact line regions, the pressure gradient $p'$ increases from almost zero in the bulk of the droplet, and along with it the interface slope decreases. The maximum magnitude of the pressure gradient in both the receding and the advancing regions corresponds to an interfacial slope of magnitude $<0.1\ \textrm {rad}$ ($5.7^{\circ }$) which presumably fulfils the small slope assumption.

Figure 2. $(a)$ Solid line: droplet shape on a conical fibre (grey shaded region) with a cone angle $\alpha =0.01\ \textrm {rad}$ at $t=256$ obtained from a numerical solution of the LAC with $\epsilon =10^{-3}$ and $\lambda =0$. Red dashed line: static droplet shape obtained from solving the uniform curvature condition $\kappa = \textrm {constant}$. $(b)$ Solid line: the Laplace pressure $p$ as a function $x$ obtained from LAC. Red dashed line: the Laplace pressure of a static droplet, where the domain of the static droplet is between $x=3.36$ and $x=4.83$. In (c,d) the pressure gradient is rescaled by its maximum value, denoted as $\bar {p}'$ and the interface slope $h'$ in the receding region in $(c)$ and the advancing region in $(d)$.

As shown in figures 2(c) and 2(d), the pressure gradients at the receding and the advancing contact line regions are large. By zooming into these regions of the droplet (see figure 3) at early times, we observe large interface curvatures, consistent with what one would expect from the results of a large pressure gradient. We want to highlight that as the droplet starts to move across the cone, a film is formed at the receding region. Since the later self-propelling state is independent of the initial conditions, the droplet properties such as the deposited film thickness and the droplet velocity are a function of the droplet position on the cone.

Figure 3. $(a)$ Interface dynamics at the receding region of the droplet, showing the formation of a deposited film from $t>2.18$. $(b)$ The interface dynamics at the advancing region of the droplet. The far field conditions of the profiles at both ends match to a constant prewetted fluid layer of thickness $h=\epsilon = 10^{-3}$ and we have $\alpha =0.01\ \textrm {rad}$ and $\lambda =0$.

3.2. Comparison of the numerical solutions of the LAC with molecular dynamics simulations

We compare the results for the case of cone angle $\alpha =0.1\ \textrm {rad}$ from the numerical solutions of the LAC and the MD simulations. The comparison serves also as a verification of our lubrication model and can help reveal nanoscopic physical effects. We first compare the centre of mass of the droplet $x_c$ as a function of time in figure 4(a). For the LAC, we have used two different values of prewetted layer thickness, i.e. $\epsilon =10^{-4}$ and $10^{-3}$, to highlight their weak influence on the results. We have three cases for the MD simulations. We see that droplets on the dry surface (case A and case B) move slower than the droplet on the wet surface (case C, with the rescaled thickness of the prewetted layer = 0.065). This is consistent with the expectation that the wetted layer reduces the frictional force between the droplet and the substrate. When comparing with the LAC results, we find that the results from the MD simulation (for both the dry surface and wet surface cases) have a larger non-dimensional velocity. When we multiply the time scale in MD by a prefactor of 2.8 for case C (with a wetted layer), we effectively shift the data from MD horizontally to the right and obtain the results represented by red circles, which makes the two models give the same results. The difference in time scale is equivalent to the difference in drop velocity, and this is related to the strength of the ‘driving force’. In the MD calculation the liquid/solid interaction strength controls the speed of the drop (as well as the slip length), but this parameter does not appear directly in the continuum calculations. Indeed, in continuum calculations of wetting usually the contact angle (for partially wetting cases) and the van der Waals interaction (for the completely wetting case) controls the motion. Here, the prewetting film thickness and the assumption that it is asymptotically constant indirectly incorporate this information. There is no obvious way to match this aspect of the molecular and continuum calculations in advance and it is no surprise that they differ. Furthermore, the presence of slip in the simulations is rather unclear because the translation velocity of the drop is much smaller than the thermal velocity of the atoms, by a factor of 10$^{-3}$ or less, and it is not possible to resolve the flow field inside the drop. However, simulations of the same liquid in shear flow along a planar solid of the same structure as the cone, under otherwise identical conditions, have a velocity field which extrapolates to zero roughly halfway between the innermost liquid and outermost solid atoms. If one (naturally) identifies the latter point as the liquid/solid boundary then the slip length is at most a small fraction of an atomic diameter, which is essentially zero. The shapes of the droplet are shown in figure 4(b–d), where the same droplet shapes are predicted by the LAC and the MD when comparing for the same centre of mass of the droplet. However, no deposited film is observed for all cases in the MD simulations.

Figure 4. $(a)$ The centre of mass of the droplet $x_c$ as a function of time $t$ obtained from the MD simulations (symbols) for three different cases and the numerical solutions of the LAC (lines) for $\lambda =0$ and $\alpha =0.1\ \textrm {rad}$. The red circles are results obtained by multiplying the time in case C with a prefactor 2.8. (b–d) Comparison between droplet profiles obtained from the LAC and the MD simulations. For each comparison, the profiles are chosen such that $x_c$ is the same. Both (b,c) are for the case C of the MD simulations (wetted substrate) but at two different droplet positions; (d) is for case A (dry substrate) of the MD simulations. Red dashed curves: profiles from LAC. Green dots: liquid molecules of the droplet. Pink colour: the cone substrate.

The effects of the thin prewetted layer are illustrated when comparing a drop advancing on a cone at the same centre of mass position ($x_c=3.45$) for the wetted case in figure 4(c) with the dry case in figure 4(d). The lighter colouring of the liquid region as compared to figure 4(b,c) reflects the fact that there are fewer liquid molecules present in the dry case. For the dry case, the advancing meniscus of the drop is irregular at molecular scales, corresponding to individual molecules hopping to attractive sites on the surface, which is also the case for wetting drops advancing on a dry flat surface (D'Ortona et al. 1996). The receding meniscus region is an uneven film as well, zero to two molecules in thickness, and this behaviour is also present in the prewetted case. The absence of a continuous trailing film for these drops is surprising because one would expect a completely wetting liquid to remain in contact with a solid unless removed by an external force, and the lubrication calculations in this paper incorporate this assumption. One possible explanation is the finite size of the simulated droplets, which may not have enough molecules to exhibit all features of continuum behaviour. A second, more specific explanation involves the curvature of the surface. Liquid adjacent to a flat surface tends to form pronounced layers and, at least for a crystalline solid, there is an ordered structure in each layer because the molecules favour positions in register with the lattice. High curvature disrupts the usual lattice structure and could thereby weaken the liquid-solid attraction. In this vein, it is known that solid curvature has a significant effect on slip lengths, which are controlled by the same interaction (Chen, Zhang & Koplik 2014; Guo, Chen & Robbins 2016).

3.3. Droplet velocity

We consider here only cases in which the prewetted layer is much thinner than the droplet, i.e. $\epsilon \in [10^{-4}, 10^{-3}]$, but we vary the slip length across a wide range $\lambda \in [0, 20]$. When the prewetted layer is thick, e.g. $\epsilon >10^{-2}$, it becomes unstable quickly due to Rayleigh–Plateau instability (Eggers & Villermaux 2008). To investigate the droplet dynamics, we define the capillary number as the dimensionless velocity of the droplet, namely

(3.2)\begin{equation} Ca\equiv \frac{\textrm{d}x_c}{\textrm{d}t}. \end{equation}

3.3.1. Dependence of droplet velocity on the thickness of the prewetted layer and the slip length

We start by looking at cases when both the prewetted layer and the slip length are small, i.e. $\epsilon \ll 1$ and $\lambda \ll 1$. In models of dynamical wetting (Bonn et al. 2009; Snoeijer & Andreotti 2013), these microscopic lengths act as a cutoff length scale for moving contact line singularity and the length scales appear in a logarithmic term of the viscous dissipation. In Chan et al. (2020), by using asymptotic matching, it is shown that the capillary number scales as $Ca= \theta _a^3/9\ln (c/\lambda )$, where $\theta _a$ is the advancing apparent contact angle of the corresponding static droplet and $c$ is a fitting parameter. We here propose a similar relation but include the prewetted layer thickness $\epsilon$ as

(3.3)\begin{equation} Ca=\frac{\theta_a^3}{9\ln(c/[\lambda+\epsilon])}. \end{equation}

With an adjustment of the fitting parameter $c$, this relation describes well the results from the LAC as shown in figure 5, but gradually becomes invalid when $\lambda$ is no longer small.

Figure 5. Lines: the capillary number $Ca$ as a function of the centre of mass of the droplet $x_c$ obtained from LAC. Symbols: the relation given by (3.3) with different values of $c$. The cone angle is $\alpha =0.01$.

When exploring a wider range of slip length $\lambda$, one expects a change in the flow profile inside the droplet, from a Poiseuille flow for the case of no-slip to a plug flow as we approach free slip (Münch, Wagner & Witelski 2005). Figure 6(a) shows $Ca$ as a function of $x_c$ for different slip lengths, where droplets move faster for larger slip lengths as the viscous dissipation is decreased. At large slip lengths, it is expected that the term with the slip length in the governing equation (2.5) dominates over the other terms, thus $\lambda$ can be scaled out from the equation by defining $\bar {t}\equiv \lambda t$. This implies $Ca$ scales linearly with $\lambda$. We plot $Ca$ rescaled by $\lambda$ in figure 6(b), and find that $Ca/\lambda$ collapse onto a single curve for $\lambda >1$, consistent with our expectation. An alternative derivation can also be made based on a balance between the rate of change of capillary energy and the viscous dissipation. The viscous stress scales as $\sim \eta U/\lambda$, giving a bulk dissipation $\sim \eta V U^2/\lambda ^2$, which is much smaller than the dissipation due to friction at the substrate $\sim \eta A_w U^2/\lambda$; here, $A_w$ is the wetted area. By balancing the dominant viscous dissipation with the rate of change of the surface energy $\partial (\gamma A_w)/\partial t \sim \gamma x_c\alpha U$ gives $U\sim \lambda$.

Figure 6. $(a)$ The capillary number $Ca$ plotted as a function of the centre of mass of the droplet $x_c$ for different slip lengths. $(b)$ The rescaled $Ca$ by $\lambda$ as a function of $x_c$. Parameters: $\epsilon =10^{-3}$ and $\alpha =0.01\ \textrm {rad}$. We rescale the data in $(a)$ as $Ca/\lambda$, which collapses the data onto a single curve for $\lambda >1$.

3.3.2. Dependence of droplet velocity on the cone angle and the droplet position

As there is no directional spreading when $\alpha =0$, it is natural to expect that a droplet moves faster at larger cone angles. This is true when comparing $Ca$ at the same cone radius, as shown in figure 7(a) in which $Ca$ is plotted as a function of $R_c\equiv x_c\tan \alpha$. The results from the LAC agree nicely with the matching results of (3.3) for the three different values of $\alpha$ using the same value of $c=2$. Remarkably, the agreement is good even when the apparent contact angle is as large as $\theta _a\approx 1\ \textrm {rad}$, for example when $R_c\approx 0.06$ and $\alpha =0.03\ \textrm {rad}$. Another feature we observe is that $Ca$ decreases when the droplet is at a position of larger cone radius for a fixed cone angle, namely the droplet slows down when moving to the thicker part of the cone. When plotting $Ca$ rescaled by $\alpha$ in figure 7(b), the results for the three different cone angles nearly collapse onto a single curve.

Figure 7. Lines: results from LAC. Symbols: the relation given by (3.3) with $c=2$. The droplet capillary number $Ca$ plotted as a function of $R_c\equiv x_c\tan \alpha$ in $(a)$, and as a function of the centre of mass of the droplet $x_c$ in $(c)$, whereas in $(b)$ $Ca$ is rescaled by $\alpha$ plotted as a function of $R_c$. Parameters: $\epsilon =10^{-3}$ and $\lambda =0$.

We have shown that the cone angle and the cone radius give opposite effects to the droplet velocity. It might be interesting to see how the droplet velocity depends on the distance from the tip of the cone, particularly the length of a fibre can be a more important parameter for certain functionality. In figure 7(c), we show $Ca$ as a function of the droplet's centre of mass $x_c$. Remarkably, a non-monotonic behaviour is observed. In the limit of large distances from the tip, droplets on cones with smaller cone angles move faster when comparing at the same $x_c$. When decreasing $x_c$, there are changes of relative strength of $Ca$. For example, at $x_c=2.3$, $Ca$ for $\alpha =0.03\ \textrm {rad}$ is even higher than that for $\alpha =0.01\ \textrm {rad}$. The reason for the non-monotonic behaviour is that two factors are playing roles when comparing at the same $x_c$, namely $\alpha$ and $R_c$. The influence of the cone radius $R_c$ is dominant over the cone angle effect when $x_c$ is large, thus droplets move faster at smaller $\alpha$. The cone angle effect becomes more important when $x_c$ is small. Our results demonstrate that a sensitive control of the geometrical factors such as $\alpha$ and $R_c$ is necessary for optimal droplet transport on these structures.

3.4. Film deposition

A film is formed at the receding region of the droplet when the droplet moves to the thicker part of the cone, as already shown in figure 3(a). We refer to the region that connects the prewetted layer and the deposited film as the film edge region. One can observe from figure 3(a) that the film edge region (around $r=1$) propagates much slower than the motion of the droplet. Hence a long deposited film is generated and the film profile is found to remain steady within the simulated time. However, the film would eventually become unstable due to the Rayleigh–Plateau instability, but the time scale for the growth of the disturbance is here greater than the time for the droplet to spread across the cone. For films of nanometric thickness, they can be stabilized by intermolecular forces (Quéré, di Meglio & Brochard-Wyart 1990).

3.4.1. Dependence of the deposited film on the cone angle

We first consider cases of no slip ($\lambda =0$). The profiles of the deposited films are shown in figure 8(a) for $\alpha =[0.01$, $0.03$, $0.05]\ \textrm {rad}$. It is found that the film thickness increases with both $\alpha$ and $r$. It is also important to understand the influence of the cone angle on the film thickness when comparing at the same cone radius. We hence plot in figure 8(b) the profiles of the films as a function of $R\equiv r\sin \alpha$. The film is thicker for larger cone angles. As will be explained in § 3.5, this is mainly due to the larger capillary number for larger cone angles.

Figure 8. $(a)$ The profiles of the deposited film $h$ plotted as a function of the distance from the tip of the cone $r$ for $\alpha =[0.01,0.03,0.05]\ \textrm {rad}$. $(b)$ The profiles of the deposited film $h$ plotted as a function of the cone radius $R$. Parameters: $\epsilon =10^{-3}$ and $\lambda =0$.

3.4.2. Dependence of the deposited film on the slip length

As droplets move faster at larger slip lengths, one may expect that a thicker film is deposited according to the LLD model. However, our analysis shows the opposite results (figure 9) with $\alpha =0.01\ \textrm {rad}$ and $\epsilon =10^{-3}$. The film thickness decreases with the slip length. We find two asymptotic film profiles. One is for the limit of small slip length ($\lambda <10^{-4}$). Another one is for the limit of large slip length ($\lambda >1$), which can be understood by the argument that the slip length is absorbed into the time variable $\bar {t}$ as explained in § 3.3. Hence, the droplet profile and the deposited film thickness become independent of the slip length. A dramatic change of film thickness appears for slip length in between $10^{-4}$ and $1$. The difference of film thickness between these two limits is of two orders of magnitude. Our results show that when the droplet size is too small, film deposition is not possible as the film thickness computed from our model can be of sub-molecular size. This is particularly relevant for large slip lengths, for which the deposited film is much thinner and the large slip regime can be realized usually for droplet size of micrometres or below.

Figure 9. The profiles of the deposited film $h$ plotted as a function of the distance from the tip of the cone $r$ for different slip lengths. Parameters: $\epsilon =10^{-3}$ and $\alpha =0.01\ \textrm {rad}$.

3.5. Asymptotic matching

Although we have shown the interface profiles $h(r)$ of the deposited films (figures 8 and 9), it is not clear yet how a particular film thickness is related to the capillary number. For droplets spreading on a cone, the motion of the droplet is self-propelled, and the capillary number is a function of the droplet position $x_c$ (or time $t$). Nevertheless, at each moment in time, the droplet deposits a portion of film with a particular thickness $h_f=h_f(x_c)$. Hence, we can link a particular film thickness to the corresponding $Ca$ at each droplet position on the cone. The procedure of determining $h_f$ is given in appendix A. A natural way to rescale the film thickness $h_f$ is by using the corresponding cone radius $R_f\equiv r_r\sin \alpha$ where the film is deposited. The rescaled $h_f/R_f$ is plotted as a function of $Ca$ for a wide range of slip length $\lambda$ in figure 10(a) with log–log axes. First, the $2/3$ scaling is not observed for any cases, even for the no-slip case which we would expect from the LLD model. Second, the local slope (in log scales) decreases with the slip length, and becomes negative when the slip length $>10^{-2}$.

Figure 10. $(a)$ The rescaled film thickness $h_f/R_f$ as a function of the capillary number $Ca$ for different slip lengths. Lines: results from the lubrication approach (LAC). Parameters: $\epsilon =10^{-3}$ and $\alpha =0.01\ \textrm {rad}$. Squares: results from the asymptotic matching (AM). The grey dashed line indicates a scaling $Ca^{2/3}$. $(b)$ Solid line: the second derivative at the apparent receding contact line of a static droplet $h''_o\equiv h''_s(r=r_{r})$ as a function of the cone radius of a deposited film $R_f$. The dashed line represents the asymptotic relation $h''_o=1/R_f$ when $R_f\ll 1$. $(c)$ The value of $a_2$ obtained from the two-dimensional lubrication equation for the film region as a function of the rescaled slip length $\bar {\lambda }\equiv \lambda /h_f$.

To understand better the numerical solutions of the LAC, we revisit the approach of asymptotic matching. We consider two regions of the liquid–air interface profile: the film region and the static droplet region, which are described by two different force balance equations. We then match the asymptotic profiles of these two regions to determine the deposited film thickness.

As the film thickness is much smaller than the cone radius, we propose that the profile in the film region locally is described by a steady solution $h=h_{2d}(x)$ of the two-dimensional lubrication equation. In the droplet frame, translating with a non-dimensional velocity $Ca$, the rescaled liquid–air interfacial profile $H(\xi )=h_{2d}/h_f$, here $\xi =xCa^{1/3}/h_f$, follows (Snoeijer et al. 2008)

(3.4)\begin{equation} \frac{\partial^3 H}{\partial \xi^3}=\frac{3}{H\left(H+3\bar{\lambda}\right)}\left(1-\frac{1}{H}\right), \end{equation}

where $\bar {\lambda }\equiv \lambda /h_f$. We impose a flat film boundary condition $H(\xi \rightarrow -\infty )=1$, and hence close to the flat film, we can write $H=1+\delta \exp [3^{1/3}\xi /(1+3\bar {\lambda })^{1/3}]$, with $\delta \ll 1$ (Oron, Davis & Bankoff 1997). The value of $\delta$ is arbitrary due to the translational invariance of (3.4). Here, we set $\delta =5 \times 10^{-7}$ when $\xi =0$. When $\xi \rightarrow \infty$, the profile of the film has to match to the droplet shape at the receding region, thus $H$ tends to $\infty$, the asymptotic solution of (3.4) is described by $H=a_1\xi +a_2\xi ^2$. The value of $a_1$ and $a_2$ are determined by the numerical solution of (3.4). A comparison between the similarity profile $H(\xi )$ and the rescaled profiles from LAC in the region connecting the flat film and the droplet is given in figure 11 for three different slip lengths: $\lambda =0$, $10^{-2}$ and $7$. The profiles from LAC are shifted manually by $r_o$ so that they match the best with the solution of (3.4). The value of $r_o$ is close to $r_r$ (difference within $2\,\%$). Note that for non-zero slip lengths, as $h_f$ varies with droplet positions, $\bar {\lambda }$ has different values at different droplet positions even though $\lambda$ is the same. We can see in figure 11 that the similarity profiles describe well the profiles from the LAC particularly at the region closer to the flat film. Away from the flat film, the profiles from the LAC bend to match the droplet shape.

Figure 11. (a,c,e) The interfacial profile $h$ as a function of the shifted radial coordinate $r-r_o$ obtained from LAC in the region connecting the flat film and the droplet for different slip lengths $\lambda$ and droplet positions (characterized by $r_r$ in b,d,f). The parameters are the cone angle $\alpha =0.01\ \textrm {rad}$ and the prewetted layer thickness $\epsilon =0.001$. (b,d,f) The rescaled profiles $h/h_f$ as a function of $(r-r_o)Ca^{1/3}/h_f$. The solid lines are the numerical solution of (3.4) with values of $\bar {\lambda }\equiv \lambda /h_f$ computed by the corresponding values of $\lambda$ and $h_f$ from the LAC.

In the static droplet region, the profile $h_s(r)$ is determined by the static equation of uniform curvature $\kappa _s$ obtained by substituting $h_s(r)=h(r)$ into (2.6), with a magnitude of $\kappa _s$ that depends on the droplet position on the cone. The problem is closed by including the boundary conditions $h_s(r=r_{r})=0$ and $h_s'(r=r_{r })=0$ at the substrate.

Now we are in a position of matching the two asymptotic profiles in the overlapping region. As we already impose the condition $h_s'(r=r_{r})=0$ for the droplet region. A natural matching condition is equating the second derivatives of the asymptotic profiles. In the film region,

(3.5)\begin{equation} h''_{2d}(\xi\rightarrow \infty)=2a_2 Ca^{2/3}/h_f. \end{equation}

Matching $h''_{2d}(\xi \rightarrow \infty )$ to the second derivative $h''_o\equiv h''_s(r=r_{r})$ in the static droplet region gives the film thickness

(3.6)\begin{equation} h_f=\frac{2a_2(\bar{\lambda})}{h''_o(R_f)}Ca^{2/3}. \end{equation}

Importantly, as all the cone angles are small, $h''_o$ is independent of $\alpha$ but only a function of the cone radius $R_f$, which is plotted in figure 10(b) in log scales. Note also that $a_2$ is determined from (3.4) and a function of $\bar {\lambda }\equiv \lambda /h_f$, which is plotted in figure 10(c). We note that $h''_o=1/R_f+\kappa _s$, and thus $h''_o\approx 1/R_f$ when $R_f\ll 1$, which is represented by the dashed line in figure 10(b). With the computed values of $h''_o$ and $a_2$, and using the values of $Ca$ for each droplet position obtained from LAC, we plot $h_f/R_f$ computed from (3.6) in figure 10(a) as square markers. Remarkably, the results from the asymptotic matching agree with the numerical results for all slip lengths. Provided the excellent agreement between the two approaches, we can understand our results in terms of the flow inside the film and the geometry of the droplet, and hence provide a better picture of the physical mechanism of film deposition by a droplet moving on a cone.

We first look at the no-slip case. When $\bar {\lambda }=0$, $a_2=0.669$, which is the same value as obtained from previous studies (Rio & Boulogne 2017). The description of the film region is the same as, for example, the dip-coating cases. Then why is the $2/3$ scaling not obtained when plotting $h_f/R_f$ as a function of $Ca$? One important aspect in our problem is that the droplets have a finite size. Hence, there are two length scales: the cone radius and the droplet radius. In terms of rescaled quantities, this means that the second derivative $h''_o$ is not a linear function of $1/R_f$, except when $R_f\ll 1$, which is already demonstrated in figure 10(b). When we rescale the results of $h_f$ obtained from the LAC by $1/h''_o$ and plot it as a function of $Ca$ in figure 12(a) for three different cone angles. The scaling $Ca^{2/3}$ is recovered and agrees well with the prediction from asymptotic matching especially for smaller cone angles.

Figure 12. $(a)$ The film thickness $h_f$ rescaled by $1/h''_o$ as a function of $Ca$ for $\lambda =0$. Symbols are results from LAC for three different cone angles. The prewetted layer thickness $\epsilon =10^{-3}$. The dashed line is the result from AM: $h_fh''_o=1.34 Ca^{2/3}$. $(b)$ The film thickness $h_f$ rescaled by $1/h''^{3}_o$ as a function of $Ca/\lambda$ for three different slip lengths $\lambda \gg h_f$. Symbols are results from LAC. The dashed line is the result from AM: $h_fh''^{3}_o=3.667 (Ca/\lambda )^{2}$.

Next we look at the slip dependence. From figure 10(c), we see $a_2$ is independent of $\bar {\lambda }$ when $\bar {\lambda }\ll 1$. Hence for $\lambda \ll h_f$, the film thickness becomes independent of $\lambda$. As the typical order of magnitude of $h_f$ for a no-slip case is $10^{-3}$, $h_f$ starts to depend on $\lambda$ significantly when $\lambda >10^{-3}$. This is consistent with our numerical results. For $\bar {\lambda }\gg 1$, we find that $a_2=0.771\bar {\lambda }^{-2/3}$. Substituting this expression of $a_2$ into (3.6), we obtain

(3.7)\begin{equation} h_f=\frac{3.667}{h''^{3}_o}\left(\frac{Ca}{\lambda}\right)^2. \end{equation}

This expression is in perfect agreement with our numerical results from the LAC for $\lambda \gg h_f$ in figure 12(b). For droplets moving on a conical fibre, we show already that $Ca \sim \lambda$ when $\lambda \gg 1$, (3.7) then suggests $h_f$ is independent of the slip length as shown in figure 9 for $\lambda \gg 1$.

One may expect that the droplet profile does not maintain a quasi-static shape when the slip length is not small due to significant viscous effects in the entire droplet. However, the excellent agreement between the results from the lubrication equation on a cone and the approach of asymptotic matching suggests that the quasi-static assumption is still valid. The reason might be the large length separation between the deposited film and the droplet height maintaining a very large difference in time scales, as one can observe from the mobility term which scales as $\sim \lambda h^2$. Thus there is sufficient time for the droplet to relax to a quasi-static shape when the apparent contact lines move. For the large slip length regime, elongational flow has been proposed to appear and dominate the viscous dissipation (Münch et al. 2005), which has been observed for dewetting droplets (McGraw et al. 2016; Chan et al. 2017). The large slip regime in our model, assuming Poiseuille flow as the dominating flow structure, is considered as the intermediate slip regime in the analysis of Münch et al. (2005). For a translating droplet, the effect of elongational flow is unclear, which requires additional experimental and theoretical studies.

For partially wetting surfaces, the droplet dynamics is the same as for the perfectly wetting cases if the prewetted layer is thick enough so that van der Waals forces between the liquid–air and the solid–liquid interfaces can be neglected. A partially wetting droplet moving on a cone without a prewetted layer will also deposit a film if it moves with a velocity above a critical value. The properties of the film are expected to be similar to the wetting cases, namely following the asymptotic relation (3.6), as long as the droplet maintains an axisymmetric shape.

4. Conclusions

The directional spreading of a viscous droplet on a conical fibre due to capillarity is investigated for small cone angles and for a wide range of slip lengths by using the lubrication equation on a cone. The droplet velocity increases with the cone angle and the slip length, but decreases as the cone radius becomes larger. At the receding part of the droplet, a film is deposited on the cone surface while the droplet is moving. When comparing with our MD simulations, we find that the droplet shapes obtained from these two approaches are the same. The velocity also shows a similar trend. However, no deposited film is observed in the MD simulations, which might be due to the nanoscopic size of the droplet.

The thickness of the deposited film observed in the LAC decreases from $h_f\approx 10^{-3}$ for the no-slip case ($\lambda =0$) to $h_f\approx 10^{-5}$ for $\lambda > 1$. We show that the film thickness obtained from the lubrication model can be understood by a similar approach of asymptotic matching used in the LLD model. For the no-slip limit, the standard $Ca^{2/3}$ scaling is recovered only when the length scale is given by $1/h''_o$ in the re-scaling. In the limit of $\lambda \gg h_f$, we find another asymptotic regime in which the film thickness scales as $h_f h''^{3}_o\sim (Ca/\lambda )^2$. For the problem we study here, the cross-over of these two regimes occurs at $\lambda \approx 10^{-4}-10^{-1}$. Our results show that manipulating the droplet size, the cone angle and the slip length provides different schemes for guiding droplet motion and coating the substrate with a film.