Film thickness distribution in gravity-driven pancake-shaped droplets rising in a Hele-Shaw cell

TL;DR

This study investigates the shape, velocity, and film thickness distribution of a buoyant droplet rising in a Hele-Shaw cell through experiments, simulations, and lubrication theory, revealing key insights into the film behavior and droplet shape.

Contribution

It provides a comprehensive analysis of droplet dynamics and film thickness in a Hele-Shaw cell, combining experimental, numerical, and theoretical approaches, and introduces a model for the film distribution.

Findings

Droplet velocity is unaffected by film thickness in the studied regime.

Experimental and numerical film thickness distributions agree well for iso-viscous cases.

The film thickness follows the Aussillous & Quéré model with fitted parameters.

Abstract

We study here experimentally, numerically and using a lubrication approach; the shape, velocity and lubrication film thickness distribution of a droplet rising in a vertical Hele-Shaw cell. The droplet is surrounded by a stationary immiscible fluid and moves purely due to buoyancy. A low density difference between the two mediums helps to operate in a regime with capillary number $Ca$ lying between $0.03-0.35$, where $Ca=\mu_o U_d /\gamma$ is built with the surrounding oil viscosity $\mu_o$, the droplet velocity $U_d$ and surface tension $\gamma$. The experimental data shows that in this regime the droplet velocity is not influenced by the thickness of the thin lubricating film and the dynamic meniscus. For iso-viscous cases, experimental and three-dimensional numerical results of the film thickness distribution agree well with each other. The mean film thickness is well captured by the Aussillous & Qu\'er\'e (2000) model with fitting parameters. The droplet also exhibits the ''catamaran'' shape that has been identified experimentally for a pressure-driven counterpart (Huerre $\textit{et al}$. 2015). This pattern has been rationalized using a two-dimensional lubrication equation. In particular, we show that this peculiar film thickness distribution is intrinsically related to the anisotropy of the fluxes induced by the droplet's motion.

Figures

Tables

Table 1. Table 1 : Dynamic viscosity μ i subscript 𝜇 𝑖 \mu_{i} and density ρ i subscript 𝜌 𝑖 \rho_{i} of the inner oils. Those of the outer oil are μ o = 560 subscript 𝜇 𝑜 560 \mu_{o}=560  mPa s and ρ o = 972 subscript 𝜌 𝑜 972 \rho_{o}=972  kg m -3 except for the case of Ricin (type 2) + 10 % percent 10 10\% ethanol for which the outer oil has μ o = 319 subscript 𝜇 𝑜 319 \mu_{o}=319 mPa s, ρ o = 970.5 subscript 𝜌 𝑜 970.5 \rho_{o}=970.5  kg m -3 . The interfacial surface tension between the inner-outer oils is γ 𝛾 \gamma . These properties were obtained at a room temperature 20 ∘ C.

Inner oil	${\mu }_i$ (mPa s)	${\rho }_i$ (kg m-3)	$\gamma$ (mN m-1)
Linseed oil	49	929	2.88
Sunflower oil	69.5	921.5	2.73
Sesame oil	71	919.1	2.64
Olive oil	79.3	913.3	2.55
Peanut oil	83.8	913.3	3.11
Ricin oil (type 1) + $10\%$ ethanol	302	943.1	4.97
Ricin oil (type 2) + $10\%$ ethanol	322	943.3	4.46

Table 2. Table 2 : List of notation, definition and working range.

Symbol	Definition	Expression	Working range
$W$	cell gap	-	4.59-4.8 mm
$U_d$	drop velocity	-	0.4-1.6 mm s-1
$Q$	injected drop volume	-	0.44-1.5 ml
$a$	pancake equivalent radius	$\sqrt{Q/\pi W}$	5.8-10.3 mm
$\alpha$	aspect ratio	$a/W$	1.25-2.27
$\Delta \rho$	density difference	$\left|{\rho }_i-{\rho }_o\right|$	27.3-58.8 kg m-3
${\mu }_i$	dynamic viscosity - droplet	-	49-322 mPa s
${\mu }_o$	dynamic viscosity - outer medium	-	319-560 mPa s
$\gamma$	interfacial surface tension	-	2.5-5.0 mN m-1
$\lambda$	dynamic viscosity ratio	${\mu }_i/{\mu }_o$	0.09-1.01
$Ca$	capillary number	${\mu }_o{U}_d/\gamma$	0.03-0.35
$Bo$	Bond number	$\Delta \rho g{a}^2/\gamma$	1.8-23.3

Equations

$$\frac{H_{\infty }}{W}=\frac{1}{2}\frac{P(3Ca{)}^{2/3}}{1+PQ(3Ca{)}^{2/3}},$$

$$U_d\sim \frac{\Delta \rho g{W}^2}{({\mu }_i+{\mu }_o)},$$

$$U_d=\frac{\Delta \rho g{W}^2}{12{\mu }_o(\lambda +1)}.$$

$$12{\alpha }^2(\lambda +1)Ca=Bo.$$

$$\frac{H_{\infty }}{W}=\frac{1}{2}\frac{P(3Ca{)}^{2/3}}{1+PQ(3Ca{)}^{2/3}},$$

$$\frac{x^2}{(L/2{)}^2}+\frac{y^2}{(T/2{)}^2}{e}^{cx}=0,$$

$$\frac{\partial }{\partial \bar{x}}\left[{\bar{H}}^3\left(\frac{1}{3Ca}{\bar{\kappa }}_x-\frac{Bo}{3Ca}\right)-\bar{H}\right]+\frac{\partial }{\partial \bar{y}}\left({\bar{H}}^3\frac{1}{3Ca}{\bar{\kappa }}_y\right)=0,$$

$$\mathbf{n}=\frac{{\left(-{\bar{H}}_x,-{\bar{H}}_y,1\right)}^T}{\sqrt{1+{\bar{H}}_x^2+{\bar{H}}_y^2}}.$$

$$\bar{\kappa}_{f,r}(\bar{r},\theta)=\underbrace{\frac{2a}{W}\Bigg{(}\frac{1+T_{f,r}(3Ca|\cos\theta|)^{2/3}}{1+Z_{f,r}(3Ca|\cos\theta|)^{2/3}}}_{\text{$(\bar{r},\bar{z})$-plane}}\Bigg{)}+\underbrace{\frac{\pi}{4}\frac{1}{\bar{r}}}_{\text{$(\bar{r},\theta)$-plane}},$$

$$\frac{{\bar{H}}_{\infty }^3}{3Ca}(\underset{{\Delta }^{2}h}{\underbrace{h_{\bar{x}\bar{x}\bar{x}\bar{x}}+2h_{\bar{x}\bar{x}\bar{y}\bar{y}}+h_{\bar{y}\bar{y}\bar{y}\bar{y}}}})-\left(1+\frac{{\bar{H}}_{\infty }^{2}Bo}{Ca}\right)h_{\bar{x}}=0.$$

$$\frac{{\bar{H}}_{\infty }^3}{3Ca}{h}_{\bar{r}\bar{r}\bar{r}\bar{r}}-\left(1+\frac{{\bar{H}}_{\infty }^{2}Bo}{Ca}\right)\left(\cos \theta h_{\bar{r}}-\frac{\sin \theta }{\bar{r}}{h}_{\theta }\right)=0.$$

$$h_{\bar{r}\bar{r}\bar{r}}=K_{p}h,$$

$$K_p=\left(\frac{3Ca}{{\bar{H}}_{\infty }^3}+\frac{3Bo}{{\bar{H}}_{\infty }}\right)\cos \theta ,$$

$$0=$$

$$u=$$

$$v=$$

$$Q^u=$$

$$Q^v=$$

$$\frac{\partial H}{\partial t}+\frac{\partial Q^u}{\partial x}+\frac{\partial Q^v}{\partial y}=0.$$

$$\frac{\partial H}{\partial t}+\frac{\partial}{\partial x}\Bigg{[}\frac{H^{3}}{3\mu_{o}}(\underbrace{\gamma\kappa_{x}}_{\mathbf{I}}-\underbrace{\Delta\rho g}_{\mathbf{II}})-\underbrace{U_{d}H}_{\mathbf{III}}\Bigg{]}+\frac{\partial}{\partial y}\Bigg{[}\frac{H^{3}}{3\mu_{o}}\underbrace{\gamma\kappa_{y}}_{\mathbf{I}}\Bigg{]}=0.$$

Full text

1 Introduction

Transport of droplets and bubbles in confined environments is a common process in engineering applications, such as microscale heat transfer and cooling using a slug flow (Kandlikar, 2012; Magnini et al., 2013), enhanced oil recovery based on foam injections where bubbles move in porous media (Farajzadeh et al., 2009) and microfluidic engineering using, droplets as micro-reactors (Song et al., 2006) and cell-encapsulating micro-compartments (He et al., 2005), to name a few. The study of transported droplets in confined dimensions also extends to biological science where red blood cells traversing passages with non-axisymmetric geometries were analysed (Halpern & Secomb, 1992).

Pioneering work has been initiated for a long bubble translating inside a straight cylindrical tube by Taylor (1961) conducting experiments and Bretherton (1961) combining experiments and asymptotic analysis. The analysis of Bretherton showed that the lubrication equations, at a very small capillary number $Ca$, were similar to their one-dimensional version assuming spanwise invariance. He established the famous asymptotic relation between the uniform film thickness $H_{\infty}$ and the capillary number in the $Ca<10^{-3}$ regime, namely, $H_{\infty}/W=P(3Ca)^{2/3}/2$, where $W$ is the tube diameter and $P$ a coefficient. The capillary number $Ca=\mu_{o}U_{d}/\gamma$ is built with the carrier phase dynamic viscosity $\mu_{o}$, the droplet velocity $U_{d}$ and the surface tension $\gamma$ between the two fluids. Aussillous & Quéré (2000) proposed

[TABLE]

as the Taylor’s law including an empirical coefficient $Q=2.5$, with the coefficient $P$ inherited from Bretherton (1961); this law was validated against the experimental data of Taylor (1961) for $Ca<2$. The empirical relation was rationalised by incorporating into the analysis of Bretherton the so-called “tube-fitting”condition, namely, that the bubble-film combination should fit inside the tube (Klaseboer et al., 2014). Besides those work considering the steady translation, Yu et al. (2018) has recently investigated how the lubrication film evolves between two steady states of a Bretherton bubble by combining theory, experiments and simulations.

Contrary to the translating bubble in a capillary tube, a bubble moving in a Hele-Shaw cell (two closely gapped parallel plates) resembles a flattened pancake. This configuration is relevant to microfluidic applications (Baroud et al., 2010) where the thickness of the microfluidic chips is much smaller than their horizontal dimension. Owing to the mathematical similarity between the governing equations of the depth-averaged Hele-Shaw flow and those of the two-dimensional (2D) irrotational flow as proved by Stokes (1898) and commented by Lamb (1993), potential flow theory was adopted to study the motion of a Hele-Shaw bubble theoretically (Taylor & Saffman, 1959) and numerically (Tanveer, 1986). Based on the stress jump derived by Bretherton (1961) and Park & Homsy (1984), 2D depth-averaged simulations including the leading-order effects of the dynamic meniscus were also carried out (Meiburg, 1989).

Motivated by the applications of droplet-based microfluidics, several works have been recently conducted to investigate the dynamics of a pressure-driven Hele-Shaw droplet. Huerre et al. (2015) and Reichert et al. (2018) performed high-precision experiments using reflection interference contrast microscopy technique to study the pressure-driven droplets, observing the so-called “catamaran” droplet shape. Simulations based on a finite volume method (Ling et al., 2016) and a boundary integral method (BIM) (Zhu & Gallaire, 2016) were carried out, confirming such a peculiar interfacial feature. It has to be mentioned that the much earlier work of Burgess & Foster (1990) performing a multi-region asymptotic analysis subtly revealed this feature for a Hele-Shaw bubble, which was rather unnoticed.

Limited work has been conducted for the gravity-driven droplets in a Hele-Shaw cell. Eri & Okumura (2011) and Yahashi et al. (2016) studied experimentally such configurations, trying to build up the scaling laws for the viscous drag friction of the Hele-Shaw droplets. Recently, Keiser et al. (2018) conducted experiments to study a sedimenting Hele-Shaw droplet, focusing on its velocity as a function of confinement, viscosity contrast and the lubrication capacity of the carrier phase.

In this work, we combine experiments, simulations and a lubrication model solved numerically to study the buoyancy-driven translation of a droplet inside a vertical Hele-Shaw cell. We examine the droplet velocity, film thickness and how they vary with the density and viscosity difference between the droplet and carrier phase. We introduce the experimental setup in §2, followed by the experimental results of the droplet mean velocity and film thickness in §3.1 and §3.2 respectively. The comparison between the three-dimensional (3D) BIM simulations and the experiments is shown in §3.3. The lubrication equation employed to model the problem is presented in §4 where the numerical solution of the lubrication equation is compared to the 3D simulation results in §4.1. The film thickness pattern is rationalised by solving the linearised 2D lubrication equation, which is presented in §4.2. We finally summarise our results in §5 with some discussions.

2 Experimental setup

A vertical Hele-Shaw cell made of two parallel glass plates, separated by a gap $W$, is filled with silicone oil of dynamic viscosity 560 mPa s and density 972 kg m*-3*, measured at $20^{\circ}$C. An oil drop is injected into the silicone oil medium from the bottom using a syringe as shown in figure 1(a). The drop moves as a result of buoyancy. The higher the density difference between the inner and outer medium, the higher the drop velocity $U_{d}$. The spanwise and streamwise cell dimensions are sufficiently large compared to the drop size to avoid any finite size effects from the lateral walls. On the other hand, the droplet is highly confined in the wall normal direction. The droplet radius $a$ is always larger than the cell gap $W$. Given the compliance of the glass walls, the thickness of which is bounded by our optical measurement tools, the cell gap $W$ lies in the range of $[4.59-4.8]$ mm and is recorded every time before the drop injection (see table 2).

The injected oils are tested beforehand to ensure non-wetting conditions for the oil droplet on the cell plates. The outer silicone oil totally wets the glass plate and forms a thin film of thickness $H$, between the drop and the glass’s interface (see figure 1(b)), which is measured using a CCI optical pen (see details in Appendix A). The pen is either placed fixed such that it measures the film thickness only along the centreline $L$ (centreline film thickness) or is mounted on a linear translational stage to perform lateral scans while the drop moves longitudinally. With an acquisition frequency of 200-500 Hz, scanning amplitude of 20-30 mm and frequency of 2-3 Hz, the obtained experimental data are interpolated in Matlab to obtain the film thickness maps for the entire drop. Droplet size and velocity determine the optimal acquisition frequency for the thickness sensor, and the scanning amplitude and frequency for the linear translational stage.

We observe that for the chosen inner oil volumetric range, the droplet in-plane shape is no longer a circle but closer to an oval, hence we refer to the drop longitudinal length (along the direction of gravity) as $L$ and to the transverse length as $T$, as shown in figure 1(b). The drop motion is captured using a Phantom Miro M310 camera with a Nikon 105 mm macro lens. The spatio-temporal analysis of the movie ensures uniform drop velocity as the drop moves along a longitudinal distance of $5L$ or more.

The drop volume $Q$ is expressed as a pancake of radius $a$ and height $W$-$2H_{\infty}$, where $H_{\infty}$ is the mean film thickness. We can simplify $Q$ as $\pi a^{2}W$ when $H_{\infty}\ll W$. For the volumetric range used for the inner oils, we found that the longitudinal and transverse lengths, $L$ and $T$, scale as the pancake radius $a$. The aspect ratio $\alpha$ is expressed as the ratio $a/W$. Keeping the cell gap $W$ fixed, data for different aspect ratios are obtained using three different volumes (0.5ml, 1ml, 1.5ml) for each oil.

Six oils with physical properties as mentioned in table 1 are used. The surface tension $\gamma$ between the inner and outer medium is measured using Teclis tensiometer and the oil viscosity and density are measured using Anton Paar SVM${}^{\text{TM}}$ 3000 viscometer. The experiment is performed at $20^{\circ}$C-$22^{\circ}$C.

The ratio $\lambda$ between the dynamic viscosity of the inner and outer phase lies between $[0.09-0.54]$. In addition to this range, another set of experiments is performed with $\lambda=1.01$, where the outer medium is silicone oil ($\mu_{o}=319$ mPa s, $\rho_{o}=970.5$ kg m*-3*) and the inner medium is a mixture of ricin oil and $10\%$ ethanol ($\mu_{i}=322$ mPa s, $\rho_{i}=943.3$ kg m*-3*) for three different drop volumes. The interfacial surface tension between these oils is 4.46 mN m*-1*. Notations for physical parameters and their definitions are detailed in table 2.

3 Experimental acquisition of the drop characteristics and their comparison with 3D BIM simulations

3.1 Experimental results for drop velocity

Considering the bulk dissipation only, the resulting viscous drag force acting on the drop scales as $F_{d}\sim(\mu_{i}+\mu_{o})U_{d}\pi a^{2}{W}^{-1}$. Unlike Okumura (2018), we consider both the inner and outer viscosities since they are of the same order. Balancing the total drag force with the buoyancy force, $F_{g}\sim\Delta\rho g\pi a^{2}W$, we obtain a scaling for the droplet mean velocity as

[TABLE]

where $\Delta\rho$ is the density difference and $g=9.81$ m s*-2*.

Under the assumption of cylindrical penny-shaped wetting drops, a theoretical expression for the drop velocity can be obtained from Maxworthy (1986), Bush (1997) and Gallaire et al. (2014). Gallaire et al. (2014) deduced the drop velocity in a Hele-Shaw cell, subjected to both buoyancy and Marangoni flow, using depth-averaged Stokes equations, called the Brinkman equations. In the absence of Marangoni effect and at leading order, Bush (1997) and Gallaire et al. (2014) predicted the mean drop velocity as

[TABLE]

Introducing the Bond number $Bo$ (refer table 2) we can rewrite equation (3) using the aspect ratio $\alpha$ as,

[TABLE]

The experimental data are plotted against the theoretical equation (4) in figure 2. Following the trend predicted by equation (4), figure 2 signifies the dominant forces in play are buoyancy and viscous drag due to the volume of fluid displaced by the drop. The dissipation induced in the thin film as well as the one in the dynamic meniscus region are found not to play a role in the selected parameter range. However it has been observed that for low $Ca$ ranges, the dissipation in the thin film (Keiser et al., 2018) and in the dynamic meniscus (Reyssat, 2014) have to be taken into account.

3.2 Experimental results for film thickness

Film thickness maps were measured for different droplet velocities. Since the thickness sensor fails to capture the data in presence of high thickness gradient, no data is acquired along the drop edges, as shown in figure 3(d), where the black curve represents the drop in-plane boundary. For different aspect ratios, qualitatively similar thickness maps were obtained, with a high film thickness on the front edge, a constant film thickness in the centre and very low film thickness along the lateral edges of the drop, overall resembling a catamaran-like shape. The spanwise and streamwise cut made along the film thickness are shown in figure 3(b-c). The centreline film thickness indicated by the streamwise cut at $y=0$ (figure 3(b)) clearly shows a monotonically decreasing film thickness pattern, followed by a region of constant film thickness $H_{\infty}$ which then reaches a minimum value of $H_{min}$. At the rear of the droplet, the strong thickness gradient reverses the direction to have an increasing thickness profile close to the drop receding edge, thus posing technical issues to capture the film thickness.

A similar centreline film thickness profile was obtained for all the droplets with a distinct value of $H_{\infty}$ and $H_{min}$. These profiles are very similar to the ones of Bretherton (1961) for pressure-driven droplets, and as already noted in other works of pancakes (Huerre et al., 2015, Zhu & Gallaire, 2016 and Reichert et al., 2018). Nondimensionalising the mean and minimum values along the centreline using the cell gap $W$ and plotting them as a function of $Ca$ shows a saturating trend for higher $Ca$ (figure 4). The experimental data are fitted based on the Taylor’s law model (Taylor, 1961; Aussillous & Quéré, 2000), according to which apart from the static and dynamic meniscus regions, the lubrication film has a constant thickness of $H_{\infty}$ given as:

[TABLE]

where the coefficients $P=0.544$ and $Q=2.061$ are obtained from the best fit curve for the experimental data. The nonlinear equation (5) is fitted using the Matlab function sseval such that the objective function, defined as the sum of squared errors between the real data of $H_{\infty}/W$ and the one predicted by equation (5), using any pair of parameters $P$ and $Q$, is the minimum. The L2 error norm for $H_{\infty}/W$ between the fitted and actual data is $0.02$.

The fitting coefficients compare well with Aussillous & Quéré (2000) and Klaseboer et al. (2014), where the mean film thickness model for bubbles ($\lambda=0$) is based on the Taylor’s law with coefficient $P=0.643$, and $Q=2.5$ and $2.79$ respectively. Fitting coefficients obtained from a 3D BIM simulation of Zhu & Gallaire (2016) for pressure-driven flows and $\lambda=1$ show the same order of magnitude as the experimental ones, with $P=0.6$ and $Q=1.5$.

In figure 4(a), we see that our experimental data for mean film thickness are bounded by the predicted values for the two extreme viscosity ratios of $\lambda=0$ and $\lambda=1$. Comparing the thickness predictions by Klaseboer et al. (2014) and Zhu & Gallaire (2016) for $Ca=0.1$ we see that the thickness variation is $20\%$ as $\lambda$ increases from 0 to 1. This is consistent with the $18\%$ (approximately) increase as reported in Martinez & Udell (1990) for pressure-driven drops in an axisymmetric tube. Further, this variation in thickness reduces to a merely $11\%$ for $Ca=0.05$, as $\lambda$ changes from 0 to 1.

The same model when used for fitting the minimum film thickness profile $H_{min}/W$ gives fitting coefficients $P=0.372,Q=1.247$ with an L2 error norm between the fitted and actual data as $0.025$.

Motivated by the qualitative agreement for the mean film thickness value between the experimental data and the 3D BIM simulations for pressure-driven droplets (figure 4(a)), we perform a 3D BIM simulation using the solver developed in Zhu & Gallaire (2016), suitably adapted for gravity-driven droplets. Details of the numerical scheme are referred to that paper.

3.3 Comparison with 3D simulations

The current numerical simulations only address the cases where the inner and outer viscosities are the same, namely $\lambda=1$. To realize it experimentally, three different drop volumes, 0.44ml, 1ml and 1.5ml, were injected in the Hele-Shaw cell resulting in $Ca$=0.032, 0.043 and $0.046$ with the corresponding $Bo$=1.81, 4.04 and 6.2. The error in volume injected decreased from $10\%$ to $3\%$ as we moved from the smallest to the largest drop volume.

The experimental film thickness maps for the chosen $Ca$ range show that the precise shape of the pancake in-plane shape is close to an oval. Hence, the experimental in-plane drop shape is obtained by fitting the following equation

[TABLE]

on an instantaneous image of the drop where $c$ is a fitting parameter (figure 3(a)).

Experimentally, due to large thickness gradient along the drop edges, the CCI sensor fails to capture the thickness in these regions. Thus the map is obtained for an area smaller than the in-plane shape of the drop (black curve in figure 5). On the contrary, the numerical simulations are capable of retrieving the complete film thickness map, but for making a visually effective comparison between the experiments and numerics, only the part of the numerical result, with the same area as the experimental data, is shown in figure 5. It’s top/bottom half corresponds to the numerical/experimental data. The red dashed curve refers to the numerical in-plane shape of the drop.

Both the experiments and simulations capture the formation of catamarans at the lateral transition regions, a uniform film thickness in the centre and a very high film thickness at the front edge of the drop. The agreement is almost quantitative. The relative error in the uniform film thickness $H_{\infty}$ for drop volumes 0.44 ml, 1 ml and 1.5 ml, is $5\%,3\%$ and $2\%$ with absolute values as $12\ \mu$m, $7\ \mu$m and $6\ \mu$m, respectively.

The numerical solution is further validated by making several streamwise (figure 6) and spanwise (figure 7) cuts along the largest drop of volume $1.5$ ml. Along the centreline, the 3D BIM simulation captures precisely the lubrication film variation: large film thickness at front edge, followed by a constant thickness profile, ending in a small oscillation before posing an increasing trend at the rear edge. There is a good quantitative comparison between the experiments and numerics, with a slight variation in the film thickness along the advancing meniscus.

4 Analysis of the film thickness pattern

In order to rationalise the film thickness pattern observed in §3.3, we model hereafter the problem using a lubrication approach. For simplicity, we formulate the 2D lubrication equation assuming the drop dynamic viscosity $\mu_{i}=0$.

4.1 Formulating the nonlinear 2D lubrication equation

Applying the long-wavelength assumption (Oron et al., 1997) and by neglecting inertia, the 2D nonlinear lubrication equation (see details in Appendix B) for the film thickness ${H}$ separating the interface from the wall, in the reference frame moving at the drop velocity $U_{d}$, can be derived. Using the pancake radius $a$ as the characteristic length and $a/U_{d}$ as the characteristic time, the dimensionless lubrication equation for the steady profile in the dimensionless coordinate system $\bar{x},\bar{y}$ is written as

[TABLE]

where ${\bar{\kappa}}$ is the mean curvature of the interface, given by $\bar{\kappa}=\nabla\cdot{\mathbf{n}}$, where the unit normal vector ${\mathbf{n}}$ on the interface is given by

[TABLE]

Note the anisotropy of the fluxes in equation (7): both the buoyancy and the motion in the $\bar{x}$ direction do not affect the flux in the $\bar{y}$ direction, breaking the isotropy induced by the capillary pressure gradient.

The nonlinear equation (7) together with the equation for the interface curvature $\bar{\kappa}$ are solved numerically by the commercial finite element solver COMSOL Multiphysics. The two variables for this coupled system of partial differential equations are $\bar{H}$ and $\bar{\kappa}$. As boundary conditions we impose the film thickness $\bar{H}=W/2a$ and the mean curvature $\bar{\kappa}=\bar{\kappa}_{f,r}$ at the droplet mid height. The mean curvature boundary condition in the static meniscus is composed by a component in the $(\bar{r},\theta)$-plane and a component in the $(\bar{r},\bar{z})$-plane. In the spirit of Meiburg (1989) and Nagel (2014), we consider the local capillary number defined with the normal velocity to the static cap for the mean curvature boundary condition model:

[TABLE]

where the coefficients with subscript $f$ have to be used for $\theta\in[-\pi/2,\pi/2]$ and the ones with subscript $r$ for $\theta\in[\pi/2,3\pi/2]$, where $\theta$ and $\bar{r}$ are defined as $\theta=\arctan(\bar{y}/\bar{x})$ and $\bar{r}=(\bar{x}^{2}+\bar{y}^{2})^{1/2}$, respectively. The values of the coefficients are $T_{f}=2.285$, $T_{r}=-0.5067$, $Z_{f}=0.4075$ and $Z_{r}=-0.1062$. The curvature boundary condition model in the $(\bar{r},\bar{z})$-plane is inspired by the equivalent model of Balestra et al. (2018), which has been developed by an extensive study for the 2D planar Stokes problem. The validity of this model has recently been corroborated by Atasi et al. (2018) for pancake bubbles. The correction $\pi/4$ for the in-plane curvature $1/\bar{r}$ in the $(\bar{r},\theta)$-plane, where $\bar{r}=1$ for a circular geometry, has been derived asymptotically by Park & Homsy (1984). Note that a more involved model could be used to describe the out-of-plane curvature (($\bar{r},\bar{z}$)-plane) in the lateral transition regions (Burgess & Foster, 1990).

In the present work we extract the pancake shape from the results of the 3D BIM simulations for $\lambda=1$. As explained in §3.3, the in-plane boundaries of the deformed pancake in the $(\bar{r},\theta)$-plane can be well described by equation (6).

It has to be stressed that the used lubrication equation should not, a priori, be used in the static meniscus region close to the boundary, where the interface slope is large. However, we have found that such an approach gives surprisingly good results if one uses the model for the static rim curvature (9) for the curvature boundary condition (see Balestra (2018) for a discussion of the axisymmetric case), which directly sets the film thickness profile in that region. Hence, such an approach can be used to numerically obtain the film thickness profile over the entire domain, also behind its validity range.

The comparison between the film thickness profile obtained by the solution of the nonlinear lubrication equation using the model equation (9) for the static cap mean curvature $\bar{\kappa}_{f,r}(r,\theta)$, with the one obtained by the 3D BIM simulations, is shown in figure 8. One can observe that both methods predict the formation of catamarans at the lateral transition regions, a uniform film thickness in the centre and oscillations at the back. In spite of the strong assumptions made for this model, the agreement is surprisingly good, even with an iso-viscous drop ($\mu_{i}=\mu_{o}$). The relative error in the uniform film thickness is of $10\%$ and its absolute value is $30\,\mu$m. The thin-film pattern shown by both approaches, as well as by the experiments, is therefore indeed a robust feature. Supported by this agreement, we investigate the thin-film pattern using the linearized version of this simple 2D lubrication model, which is computationally much cheaper than the 3D Stokes simulations.

4.2 Qualitative analysis of thickness pattern using the linearized 2D lubrication equation

The qualitative nature of the film thickness pattern can be inferred by performing a linear analysis of the 2D lubrication equation (7). With the use of the film thickness decomposition $\bar{H}=\bar{H}_{\infty}+\varepsilon h$, where $\bar{H}_{\infty}=H_{\infty}/a$, the linear equation for the first-order film thickness correction reads:

[TABLE]

The film thickness $\bar{H}_{\infty}$ is expressed using the empirical model (5) (Taylor (1961), Aussillous & Quéré (2000) and Balestra et al. (2018)), with $P=0.643$ and $Q=2.2$.

Equation (10) for the film thickness correction around the uniform film thickness $\bar{H}_{\infty}$ can be solved as a boundary-value problem, as recently conducted by Atasi et al. (2018). In contrast to the nonlinear solution of Sec. 4.1, here we only solve the lubrication equation from the thin-film region up to the beginning of the dynamic meniscus region. This is equivalent to looking at the first-order correction of the uniform thin film region due to the matching of the film thickness in the dynamic meniscus region to a larger value. In the present context, we impose a film thickness correction $h=A$ and a mean curvature of the order $\Delta h=1/A+1/\bar{r}$ on the perimeter, with $A$ as a constant value of $10^{-3}\times\bar{H}_{\infty}$. This boundary condition does not have to be understood as a rigorous matching approach, but rather as a way to find the structure of the film thickness profile in the region where it is close to be uniform. A rigorous matching for the limit $Ca\ll 1$, can be found in Park & Homsy (1984). The maps of the film thickness correction $h$, together with some profiles along the streamwise and spanwise directions, are shown in figure 9.

First, it can be clearly observed that the linear lubrication equation with a perturbed film thickness and curvature along the domain boundary is able to reproduce the catamaran-like pattern observed in pancake droplets as seen in §3.3 and §4.1. The film thickness is the smallest in the lateral part of the pancake (see figure 9(a)), so that its 3D shape resembles the hull of a catamaran. Therefore, we can conclude that this pattern is the generalization of the one-dimensional oscillations found by Bretherton (1961) at the rear of an axisymmetric bubble for a 2D concave structure, like a pancake droplet, and is intrinsically related to the anisotropy of the equation.

Second, the film thickness correction along the streamwise direction $\bar{x}$ (see figure 9(c)) deviates from a uniform profile as expected from Bretherton’s theory (Bretherton, 1961). The film thickness oscillates at the rear meniscus and increases monotonically at the front meniscus. Note that the film thickness correction in the uniform film region of a pancake is not vanishing as the base film thickness $\bar{H}_{\infty}$ is given by equation (5), which is an asymptotic estimate for $Bo=0$ but not an exact solution of the lubrication equation with $Bo\neq 0$. Furthermore, it can be observed that the more one moves away from $\bar{y}=0$, the more the thickness of the film is reduced. Therefore, the thickness of the film left by the front meniscus is not uniform.

To better highlight this crucial point, we show in figure 9(d) the normalized difference between the film thickness correction and its value at $\bar{y}=0$. The film thickness decreases as $|\bar{y}|$ increases, before increasing again close to the edge to match the boundary condition.

These qualitative observations can be rationalised by simplifying the linear lubrication equation (10) for the different regions of the domain (see figure 9(b)). The lubrication equation (10) in polar coordinates can be simplified to

[TABLE]

For small polar angles $\theta$, the contribution $(\sin\theta/\bar{r})h_{\theta}$, which corresponds to the flux in the tangential direction, can be neglected so that the linear lubrication equation becomes, after integration along $\bar{r}$:

[TABLE]

with

[TABLE]

which is the linearised one dimensional equation for the Landau-Levich-Derjaguin-Bretherton problem (Landau & Levich (1942), Derjaguin (1943) and Bretherton (1961)) in the radial direction $\bar{r}$ projected on the streamwise direction. Therefore, we know from the solution of Bretherton (1961) that the film thickness is oscillating at the rear meniscus and monotonically increasing at the front one. Focusing for now on $Bo=0$, we know that the thickness deposited by a front meniscus depends on the velocity normal to the interface. In this case, one has therefore $\bar{H}_{\infty}\sim Ca_{p}^{2/3}$ with $Ca_{p}=Ca\cos\theta$ as the local capillary number at a given polar angle $\theta$. Hence, the film thickness in the central region of the pancake varies like $(Ca\cos\theta)^{2/3}$. Similar results have been reported for a pressure-driven red blood cell traversing in a non-axisymmetric passage (Halpern & Secomb, 1992) and in pancake droplets by Reichert et al. (2018). Once a given film thickness is set by the front meniscus, the same thickness will be present over the entire thin film region at the corresponding spanwise location $\bar{y}$. The good agreement between the dependence in $({\cos\theta})^{2/3}$ of the film thickness and the profile along the spanwise direction obtained by resolving the 2D lubrication equation is shown in figure 9(d).

Similarly, the oscillations at the rear meniscus depend on the polar angle. For a pancake droplet, due to the film thickness nonuniformity resulting from the nonuniform deposition at the front, the wavelength of the oscillations at the back scales as $\lambda_{osc}\sim(Ca\cos\theta)^{1/3}$. Given the $1/3$ power-law dependence, the wavelength is almost unchanged, before rapidly reducing to [math] when $\theta\rightarrow\pm\pi/2$ (see figure 9(a)).

It is important to note that a plane cut of the film thickness at a given angle $\theta$ does not present a region of constant film thickness. A pancake droplet cannot be seen just as the collection of different one-dimensional profiles obtained by the solution to equation (12) for different polar angles $\theta$. In fact, the film thickness at any spanwise location $\bar{y}$ is set by the front meniscus at the corresponding polar angle $\theta$ and equation (12) only indicates the scaling of this film thickness as well as the oscillations at the back.

For $\theta\rightarrow\pm\pi/2$, which corresponds to the lateral meniscus region (see figure 9(b)), the tangential flux term $(\sin\theta/\bar{r})h_{\theta}$ in equation (11) can no longer be neglected. Burgess & Foster (1990) performed an involved analysis of the lubrication equation in this region for a pancake droplet at low capillary numbers and found that the local film thickness in the so-called lateral transition regions scales as $Ca^{4/5}$ rather than as $Ca^{2/3}$. Therefore, for $Ca\ll 1$, the film thickness in these lateral regions is much smaller than the one in the other regions. This explains why one observes catamaran-like structures in the lateral regions of pancake droplets. Note that for the $Ca$-range considered in the present study, the film thickness in the lateral transition regions is still sufficiently large so that the viscous dissipation can be neglected also in these regions, as confirmed by the results of Sec. 3. Furthermore, Burgess & Foster (1990) have also shown that the polar extent of these lateral regions scales as $Ca^{1/5}$, whereas their radial extent scales as $Ca^{2/5}$ instead of as $Ca^{1/3}$ that one has for the length of rear and front dynamic menisci of axisymmetric droplets (see also Hodges et al. (2004)).

5 Conclusions and perspectives

We report the velocity, mean film thickness and thickness map for a droplet moving due to buoyancy in a vertical Hele-Shaw cell. The mean drop velocity compares well with the leading order velocity expression of Gallaire et al. (2014). This signifies that buoyancy and viscous drag force are the dominant forces in our experimental parameter range with the viscous dissipation in the film thickness and in the dynamical meniscus having a negligible effect on the droplet velocity. On the contrary, the dimensionless mean film thickness data was dependent on the dimensionless droplet velocity, expressed as $Ca$, and was fitted with the Taylor’s law model (Aussillous & Quéré, 2000).

We also obtained the complete film thickness maps using a CCI optical pen mounted on a linear translation stage. Based on a boundary integral method, 3D Stokes equations were solved. These numerical results for $\lambda=1$ were in very good agreement with our experimental results. The thickness pattern had a distinct catamaran-like shape as experimentally observed for pressure-driven flows in Huerre et al. (2015) and Reichert et al. (2018).

To understand the nature of the thickness patterns observed experimentally and numerically, the problem was approached using a lubrication equation, which was solved as a boundary value problem, rather than as an initial value problem, as recently conducted by Atasi et al. (2018). In spite of all the crude assumptions done for developing the model, its nonlinear solution for the film thickness profile of a pancake bubble compared surprisingly well with the results of 3D BIM simulations, evidencing the robustness of the thin-film pattern.

In order to unravel the structure of the film thickness profile, we linearized the lubrication model and solved for the linear thickness corrections around a uniformly-thick film. We have been able to show that not only the oscillations at the rear meniscus, but also the catamaran-like pattern can be directly retrieved by solving the linear 2D lubrication equation when perturbing the film thickness at the boundaries, which mimics the presence of a meniscus of greater film thickness. In particular, the catamaran-like structure results from the anisotropic flux induced by the motion of the walls with respect to the pancake and the need to match the film thickness to larger values in the dynamic meniscus region surrounding the region where the thin film is rather uniform. This pattern is therefore independent of the force driving the motion. In fact, in totally different contexts, the same pattern is also found in drops levitating on a moving substrate (Hodges et al., 2004; Lhuissier et al., 2013) as well as in oleoplaning drops (Daniel et al., 2017). In the central part of the pancake droplet, the thickness left by the front meniscus scales as $(Ca\cos\theta)^{2/3}$, and depends therefore on the velocity normal to the interface. This scaling no longer holds in the lateral transition region, where the component of the flux tangential to the interface becomes important and the thickness of the film is much smaller, resulting in the formation of the catamaran-like structure.

Finally, we would like to highlight a contrasting feature seen between drops moving in a cylindrical tube and that in a Hele-Shaw cell. In cylindrical tubes the main difference between the pressure-driven and buoyancy-driven bubble motion, as found by Bretherton (1961), is that buoyant bubbles may remain stuck if the capillary radius is less than $0.918l_{c}$, where the capillary length $l_{c}=\sqrt{\gamma/\Delta\rho g}$. This failure results from the impossibility to match the static gravity-corrected meniscus shape with the flat thin film region. A similar result was obtained recently in the planar geometry by Lamstaes & Eggers (2017) with a prefactor of $0.847$. Interestingly enough, we suspect that there is no such bubble arrest in Hele-Shaw cells, as a consequence of the additional direction which adds a degree of freedom in the curvature. As shown in figure 10, our experimental results show marked drop motion when the half cell gap is below $0.847l_{c}$.

I. S. thanks the Swiss National Science Foundation (grant no. 200021-159957). The computer time is provided by the Swiss National Supercomputing Centre (CSCS) under project ID s603. An ERC starting grant “SimCoMiCs 280117” is gratefully acknowledged. L. Z. thanks a VR International Postdoc Grant “2015- 06334” from the Swedish Research Council. We thank Ludovic Keiser for fruitful discussions.

Appendix A CCI working principle

We hereby describe the principle of the Confocal Chromatic Imaging technique. An achromatic lens decomposes the incident white light into a continuum of monochromatic images which constitutes the measurement range. The light reflected by a sample surface put inside this range is collected by a beam splitter. A pinhole allows then to block the defocused light that does not come from the sample surface. Eventually, the spectral repartition of the collected light is analyzed by a spectrometer. The wavelength of maximum intensity is detected and the distance value is deduced from a calibration curve. Several reflecting interfaces can be detected at the same time, allowing thickness measurement of thin transparent layers.

Appendix B Derivation of two-dimensional nonlinear lubrication equation for pancake droplets

The derivation of the model equation presented in §4.2 is briefly outlined hereafter. Considering the same physical properties for the droplet and the outer medium as outlined in §1, and under the assumption of negligible inertia (Oron et al., 1997) with $\rho_{i}\ll\rho_{o}$, and $\mu_{i}\ll\mu_{o}$ the three dimensional momentum equations reads

[TABLE]

Using $L$ as the characteristic length scale in $x$ and $y$ direction and the film thickness $H$ as the characteristic length scale in $z$ direction, the film aspect ratio $\epsilon$ is defined as $\epsilon=H/L$. The long wavelength approximation is employed since $\epsilon\ll 1$. Mass conservation indicates that the characteristic velocity in $z$ direction ($W$) is much smaller than the other two components ($U$ in $x$ and $V$ in $y$ direction), $W\sim\epsilon U\ll U$ and $W\sim\epsilon V\ll V$. The Stokes equation simplifies as

[TABLE]

Integrating equation (19) in $z$ and applying dynamic boundary conditions yields $p=p_{0}-\gamma\kappa$ where ${\kappa}$ is the mean curvature of the interface. Integrating equation (17) and (18) twice in $z$ and considering $u(z=0)=-U_{d}$ and the zero-slip boundary condition $v(z=0)=0$ as well as the zero-shear-stress interface $\partial u(z=H)/\partial z=0$ and $\partial v(z=H)/\partial z=0$ yields the velocity components:

[TABLE]

Since the inner medium density $\rho_{i}\ll\rho_{o}$, we replace $\rho_{o}$ by $\Delta\rho$ in equation (20) where $\Delta\rho=\rho_{o}-\rho_{i}$ represents the density difference between the inner and outer fluid. Integrating equation (20) and (21) in $z$ from [math] to $H$ yields the flux in $x$ and $y$ direction as

[TABLE]

Finally integrating the continuity equation and applying the Leibniz Integral rule and the kinematic boundary condition at the interface yield the mass conservation equation expressed as

[TABLE]

Introducing equation (22) and (23) in equation (24) finally yields the lubrication equation:

[TABLE]

The terms $\mathbf{I}$ in the spatial variation of the flux corresponds to the surface tension effects, term $\mathbf{II}$ to the variation due to the buoyancy force and term $\mathbf{III}$ accounts for the reference frame moving with the drop. Note the anisotropy of the fluxes: both the buoyancy and the motion in the $x$ direction do not affect the flux in the $y$ direction, breaking the isotropy induced by the capillary pressure gradient.

Using the pancake radius $a$ as the characteristic length and $a/U_{d}$ as the characteristic time, the dimensionless lubrication equation for the steady profile in the dimensionless coordinate system $\bar{x},\bar{y}$ is written as equation (7).