Sliding drops - ensemble statistics from single drop bifurcations

TL;DR

This paper investigates the complex dynamics of sliding drops on inclined surfaces, revealing how merging and splitting processes influence the steady-state size distribution through numerical simulations and a statistical model.

Contribution

It introduces a combined approach of numerical simulations and bifurcation analysis to understand ensemble drop dynamics and develops a Smoluchowski-type model for these processes.

Findings

Drop ensembles exhibit a balance of coalescence and breakup.

The size distribution reaches a stationary state.

The statistical model agrees well with direct simulations.

Abstract

Ensembles of interacting drops that slide down an inclined plate show a dramatically different coarsening behavior as compared to drops on a horizontal plate: As drops of different size slide at different velocities, frequent collisions result in fast coalescence. However, above a certain size individual sliding drops are unstable and break up into smaller drops. Therefore, the long-time dynamics of a large drop ensemble is governed by a balance of merging and splitting. We employ a long-wave film height evolution equation and determine the dynamics of the drop size distribution towards a stationary state from direct numerical simulations on large domains. The main features of the distribution are then related to the bifurcation diagram of individual drops obtained by numerical path continuation. The gained knowledge allows us to develop a Smoluchowski-type statistical model for the ensemble dynamics that well compares to full direct simulations.

Equations

$$\partial_{t}h=-\nabla\cdot\left[\frac{h^{3}}{3}\nabla\bigr{(}\Delta h+\Pi(h)\bigr{)}+\frac{h^{3}}{3}G\begin{pmatrix}\alpha\\ 0\end{pmatrix}\right].$$

$$U(V)=a_0\alpha V^{{\beta }_0},$$

$$V_{\mathrm{cr}}=a_1{\alpha }^{{\beta }_1},$$

$${\partial }_{t}f(V,t)=-\underset{\mathrm{loss}\mathrm{due}\mathrm{to}\mathrm{coalescence}}{\underbrace{\underset{0}{\overset{\infty }{\int }}K(V,\tilde{V})f(V)f(\tilde{V})\mathrm{d}\tilde{V}}}$$

$$+\underset{\mathrm{gain}\mathrm{due}\mathrm{to}\mathrm{coalescence}}{\underbrace{\underset{0}{\overset{V}{\int }}\frac{1}{2}K(\tilde{V},V-\tilde{V})f(\tilde{V})f(V-\tilde{V})\mathrm{d}\tilde{V}}}$$

$$-\underset{\mathrm{loss}\mathrm{due}\mathrm{to}\mathrm{fragmentation}}{\underbrace{\underset{0}{\overset{\infty }{\int }}\frac{1}{2}J(V,\tilde{V})f(V)\mathrm{d}\tilde{V}}}+\underset{\mathrm{gain}\mathrm{due}\mathrm{to}\mathrm{fragmentation}}{\underbrace{\underset{0}{\overset{\infty }{\int }}J(\tilde{V},V)f(\tilde{V})\mathrm{d}\tilde{V}}},$$

$$\mathrm{where}K(V_a,V_b)=\frac{2k_1}{L}\left|U(V_a)-U(V_b)\right|,$$

$$J(V_a,V_b)=j\sigma (V_a,V_{\mathrm{cr}})\Theta \left(V_a-V_b\right).$$

Full text

Introduction

The coarsening of small-scale structures as, for instance, clusters, crystals, drops or quantum dots, into larger ones is a fundamental and widely investigated physical process common in nature and technology Bray (1994); Nepomnyashchy (2015). Amongst the first investigations of coarsening dynamics is Ostwald’s work on the growth of larger crystals or particles in solution at the expense of smaller ones Ostwald (1896). Often, these processes of Ostwald ripening can be described by scaling laws for the time evolution of typical length scales. The power law scaling for cluster growth was explained by Lifshitz and Slyozov Lifshitz and Slyozov (1961) and independently by Wagner Wagner (1961). In their derivation, the dynamics of the individual objects is related to the dynamics of the entire ensemble.

A particular soft matter example is an ensemble of liquid drops on a solid substrate which naturally exhibits coarsening. The statistical description of condensing and coarsening drops on horizontal substrates, i.e., their evolution towards equilibrium, was addressed by Meakin and coworkers in terms of particle based statistical models Meakin (1992) and also by Smoluchowski-type integro-differential equations for volume distribution functions Smoluchowski (1916). They also consider the case of inclined substrates, where the drops are initially pinned and then depin at a critical volume in an avalanche process. The coarsening and migration of liquid drops on horizontal substrates was also addressed in detail for the one- (1D) Gratton and Witelski (2009); Kitavtsev and Wagner (2010); Kitavtsev (2014); Glasner and Witelski (2003) and two-dimensional (2D) case Pismen and Pomeau (2004); Otto et al. (2006) employing a lubrication or long-wave model Oron et al. (1997); Thiele (2007). The relation to Ostwald ripening is also discussed in Glasner et al. (2009).

Here, we analyze the coarsening dynamics of liquid drops that due to gravitation slide down a plate of fixed inclination. The lateral motion of the drops with respect to each other depends strongly on differences in drop size resulting in a fast relative transport that facilitates coarsening, i.e., the coalescence of smaller drops into larger ones. At the same time, large drops above a certain critical size are unstable with respect to break-up into smaller ones, due to the so-called pearling instability. Similarly, drops of fixed size are unstable above a critical substrate inclination Podgorski et al. (2001); Engelnkemper et al. (2016). We investigate the interplay of the accelerated coarsening and the pearling instability and elucidate the resulting statistical properties of large ensembles of sliding drops. To this end, we employ a long-wave film height evolution equation and conduct large-scale direct numerical simulations (DNS) of sliding drop ensembles to extract the dynamics of statistical measures like the drop size distribution. Next, the resulting stationary distribution of the ensemble is related to the bifurcation diagram and stability properties of individual drops obtained by numerical path continuation techniques Dijkstra et al. (2014). Finally, we merge the numerically obtained single-drop information including several scaling laws and develop an augmented Smoluchowski coagulation equation as simple statistical model that describes the dynamics of the drop size distribution.

Modelling and Numerical Implementation

A non-dimensional long-wave equation is used to model the time evolution of the height profile $h(x,y,t)$ that describes drops of a simple liquid on a partially wetting substrate, cf. Engelnkemper et al. (2016):

[TABLE]

The model accounts for the surface tension of the liquid via a Laplace pressure, substrate-liquid interactions such as wettability via a Derjaguin (or disjoining) pressure $\Pi(h)=-\frac{1}{h^{3}}+\frac{1}{h^{6}}$ and for the lateral driving where $G$ and $\alpha$ are a non-dimensional gravity parameter and the inclination angle of the substrate, respectively 111Starting from a dimensional form $\tilde{\Pi}(\tilde{h})=-A/\tilde{h}^{3}+B/\tilde{h}^{6}$ for the disjoining pressure, we employ scales $h_{\mathrm{eq}}=(B/A)^{1/3}$ for height, $l_{0}=\sqrt{3}h_{\mathrm{eq}}/\sqrt{5}\theta_{\mathrm{eq}}$ for lateral lengths, and $t_{0}=9\eta h_{\mathrm{eq}}/25\gamma\theta_{\mathrm{eq}}^{4}$ for time. Then $\theta_{\mathrm{eq}}=\sqrt{3A/5\gamma h_{\mathrm{eq}}^{2}}$ is the equilibrium contact angle and $G=3\rho gh^{2}_{\mathrm{eq}}/5\gamma\theta_{\mathrm{eq}}^{2}$ is the gravitation number. Here, we use $G=10^{-3}$.. The employed Derjaguin pressure Pismen (2001) results in the presence of a thin adsorption layer in the whole domain on which the drops slide. DNS of this model are conducted on a large spatial domain $\Omega=[0,4000]\times[0,4000]$ with periodic boundary conditions using a finite-element method on a quadratic mesh with bilinear ansatz functions and a 2nd-order implicit Runge-Kutta scheme for time-stepping, implemented using the DUNE PDELab framework Bastian et al. (2010, 2008a, 2008b) (for more numerical details see Engelnkemper et al. (2016)).

Properties of the Drop Ensembles

Figure 1 presents simulation snapshots at different times $t$ that contrast coarsening on a horizontal ($\alpha=0$, top row) and an inclined ($\alpha=1.0$, bottom row) substrate. Up to $t\approx 0.5\cdot 10^{4}$, the coarsening proceeds very similarly in both cases, however, at non-zero inclination the later stages are dominated by a faster coarsening process that results in larger drop sizes. This continues until a certain time ($t_{\mathrm{c}}\approx 6.5\cdot 10^{4}$ at $\alpha=1.0$), after which the typical drop size hardly increases further, because the pearling instability breaks up all drops above a certain volume. Only at very late stages of the simulation, a tendency to form large elongated drops can be noted.

To quantify the coarsening process, we use the total number of drops $N_{\mathrm{D}}(t)$ in the domain 222As the height profile $h(x,y,t)$ is a single, continuous field, we define an individual drop as a connected area $\Omega_{\mathrm{drop}}$ where the height exceeds a threshold slightly above the height of the adsorption layer (here $h_{\mathrm{thresh}}=1.05$)., as well as the drop size distribution $f(V,t)$ obtained by a Gaussian kernel density estimation (KDE) 333To obtain $f(V,t)$, at each step of the DNS, the volume $V$ of each drop is calculated by integrating the height profile $h(x,y,t)$ over the corresponding footprint $\Omega_{\mathrm{drop}}$. From the resulting list at time $t$ we calculate $f(V,t)$ using a KDE Scott (2015). In this description, $N_{\left[V,V+\mathrm{d}V\right]}=f(V)\mathrm{d}V$ is the number of drops with a volume in the interval $\left[V,V+\mathrm{d}V\right]$. Figure 2 shows the time-evolution of the normalized drop size distribution $\tilde{f}(V,t)=f(V,t)/N_{\mathrm{D}}(t)$ in the inclined case of Fig. 1 ($\alpha=1.0$), while the change in the total number of drops $N_{\mathrm{D}}(t)$ is presented for various inclinations in Fig. 5 (bottom panel, solid lines). We find, that the inclination-induced acceleration of coarsening results in a fast drop number decrease. In particular, this coarsening is always faster than the classical rigorous scaling law $N_{\mathrm{D}}(t)\sim t^{-\frac{3}{4}}$ in the horizontal case Otto et al. (2006), and further accelerates with increasing inclination. In the drop size distribution (cf. Fig. 2), the fast coarsening is visible as a strong broadening and shifting of the initially tightly peaked distribution towards larger volumes. Around $t\approx 1\cdot 10^{5}$, it then develops a second local maximum at $V\approx 6\cdot 10^{4}$, which grows in time at the cost of smaller values. Finally, after an inclination-dependent time $t_{\mathrm{c}}$, the coarsening almost stops, as indicated by a significant kink in the $N_{\mathrm{D}}(t)$ curve (Fig. 5 bottom) and a subsequent very slow decrease. In the DNS (see Fig. 1), this phase occurs for $t>t_{c}\approx 6.5\cdot 10^{4}$ and is characterized by drop ensembles consisting of similar-sized drops, in accordance with the quite uniform and almost stationary drop size distribution (Fig. 2). Therefore, the drops slide with small relative velocities, leading to only a few coalescence events. These mergings often result in large drops that are unstable w.r.t. pearling and break-up again. In this way, statistically an almost stationary state is reached and kept in which merging and break up of drops balance.

Stability Properties of Single Drops

The time evolution of the drop size distribution results from the interplay of drop interactions (dominated by their relative velocity) and stability properties of individual drops. Both information is presented in Fig. 3 in the form of a bifurcation diagram obtained by pseudo-arclength continuation within the PDE2Path framework Uecker et al. (2014). It shows for a single drop at fixed inclination the dependence of sliding velocity on drop volume $V$ (cf. Ref. Engelnkemper et al. (2016) for other cases and implementation details).

Figure 3 reveals the existence of a variety of different drop shapes, velocities and stability properties. For small drop volumes $V$, only simple, almost ellipsoidal cap-shaped drops exist (sub-branch (a)). Increasing $V$, this sub-branch terminates at a critical volume $V_{\mathrm{cr}}\approx 7.2\cdot 10^{4}$ in a saddle-node bifurcation, which also connects it to sub-branch (b), whose drops exhibit an elongated tail and are linearly unstable. Sub-branch (b) connects to the stable sub-branches (c,d) via another saddle-node bifurcation and a subsequent Hopf bifurcation (cf. Engelnkemper et al. (2016)). Although at small volumes only the drops of sub-branch (a) exist, from $V\approx 5.3\cdot 10^{6}$ onwards, we find a multistability of sub-branch (a) and the elongated drops of sub-branch (c). In the DNS, one mainly observes drops from sub-branch (a) because drops of larger volume than $V_{\mathrm{cr}}$ which are formed by merging are normally unstable and decay by pearling. However, sometimes the merged drop is elongated and linearly stable, i.e., on sub-branches (c,d).

Next, we connect the information gained from the bifurcation study of the individual drop to the ensemble dynamics. As at relatively low inclinations stable elongated drops are rarely formed, we focus on sub-branch (a): The bifurcation point at $V_{\mathrm{cr}}$ provides the stability limit for simple drops and, therefore, sets an upper limiting volume for the ensemble DNS. Figure 4 shows bifurcation curves together with the late-stage quasi-stationary drop size distributions obtained from DNS. Comparison shows, that the location of the main peak of the distribution is directly connected to the position of the saddle-node bifurcation at $V_{\mathrm{cr}}$: The number of drops with $V>V_{\mathrm{cr}}$ decreases significantly. Indeed, $V_{\mathrm{cr}}$ almost coincides with the r.h.s. inflection point of the size distribution. These observations equally hold for different inclinations $\alpha$ (see Fig. 4).

Further, we extract from bifurcation diagrams as the one in Fig. 3 power laws that relate (i) drop velocity and volume at fixed inclination

[TABLE]

and (ii) the critical $V_{\mathrm{cr}}$ and inclination,

[TABLE]

with $a_{0}=2.2\cdot 10^{-4}$, $\beta_{0}=0.569$, $a_{1}=7.18\cdot 10^{4}$ and $\beta_{1}=-1.40$ Engelnkemper et al. (2016).

Statistical Model

In the final step, the obtained ’single-drop’ information is employed to develop a minimal statistical model for the ensemble dynamics as characterized by the unnormalized drop size distribution $f(V,t)$. To capture the coarsening dynamics that is dominated by the interplay of collision-caused merging and instability-caused splitting of drops, we extend Smoluchowski’s continuous rate equation for coagulation Smoluchowski (1916), following the approach of Meakin et al. for breath figures (see Meakin (1992) and references therein). Thereby, loss and gain of drops of each volume are accounted for through continuous transition rate kernels for coalescence and fragmentation. The model conserves the total volume and reads

[TABLE]

[TABLE]

[TABLE]

The properties of the kernels $K$ and $J$ in this non-local evolution equation are crucial features of the coarse-grained model. We deduce them from the single-drop results (2) and (3) above and employ a minimum of free parameters and assumptions. In particular, the kernel $K(V_{a},V_{b})$ (cf. Eq. (5)) accounts for the coalescence of two drops with volumes $V_{a}$ and $V_{b}$. It sets the frequency of collisions as the ratio of the relative drop velocity and the mean distance $L/2$ between two drops on the domain. The drop velocities $U(V)$ are given by the obtained scaling law (2) with the only a priori unknown parameter being $k_{1}$. It is a measure for the reduction of the number of collisions because all drops slide in the same direction and therefore only interact with a subset of the other drops. The other kernel $J(V_{a},V_{b})$ (cf. Eq. (6)) with the sigmoid function $\sigma(V_{a},V_{\mathrm{cr}})$ 444We use $\sigma(V_{a},V_{\mathrm{cr}})=\frac{1}{2}\left(1+\tanh\left(\frac{V_{a}-(V_{\mathrm{cr}}+2b_{V})}{b_{V}}\right)\right)$ accounts for drop splitting and corresponds to the simplest implementation of the instability threshold obtained above. In particular, drops with $V_{a}>V_{\mathrm{cr}}$ fragment into two drops of volume $V_{b}$ and $V_{a}-V_{b}$, respectively, with equal probability for all $V_{b}<V_{\mathrm{cr}}$ (expressed by the Heaviside function $\Theta$) 555This is a simple rational choice as information about the specific instability timescales for specific drop sizes and particular fragmentation ratios is difficult to obtain and too complex for a simple model Engelnkemper et al. (2016).. Here, the free parameters are the smoothness $b_{V}$ of the transition to the unstable regime and the timescale ratio $j$ between fragmentation and coalescence processes.

The developed statistical model (Statistical Model) is solved numerically 666Discretized on the volume domain and with a 4th order Runge-Kutta time-stepping scheme. employing initial conditions corresponding to early stages of the DNS of Eq. (1) (e.g., in Fig. 1). As a result, the top row of Fig. 5 compares the two dynamics of the drop size distribution as measured in the DNS of the thin film equation (1) and in the simulation of the statistical model (Statistical Model). It shows a very good agreement of all main features, as e.g., the appearence of a second peak and the convergence to a quasi-stationary distribution.

Furthermore, the bottom panel of Fig. 5 compares the evolution of the drop number $N_{D}$ for different inclination angles $\alpha$ fixing in all cases parameters $k_{1}$ and $j$ such that a best fit results for simulations with $\alpha=1.0$. Nevertheless, the predictions of the statistical model for all inclinations agree very well with the DNS results. This gives clear evidence that the dynamics of the ensemble properties resulting from many individual complex coalescence and fragmentation processes can be rather well captured by our simple statistical model.

Conclusions

We have investigated the coarsening behavior of ensembles of interacting sliding drops employing a thin-film equation. We have shown in direct simulations that a balance of coalescence and fragmentation processes emerges that can be related to stability properties of individual drops as captured in a single-drop bifurcation diagram. Dynamically, the statistical ensemble properties converge to an almost stationary state. Further, based on the gained single-drop information, we have developed a minimal statistical model that faithfully captures the main ensemble dynamics and very well compares to the full direct numerical simulations. We believe that the proposed methodology of employing ’microscopic’ information in the form of bifurcation properties of individual entities (here drops), to derive coarse-grained ’macroscopic’ statistical models for the ensemble dynamics, represents a multi-scale approach that will prove useful in other nonlinear nonequilibrium systems.

We acknowledge partial support by DFG within the Sino-German Collaborative Research Center TRR 61, and GIF under grant I-1361-401.10/2016.