2.1. Experimental set-up

Figure 1(a) describes the experimental set-up. Here, a glycerol/water pendant drop of radius $R=1$ mm comes in contact with a quiescent silicone oil film of thickness $H$. The radial expansion of the neck of the liquid bridge formed as the oil spreads on the aqueous drop during engulfment was imaged using a Photron Fastcam SA5 with a Nikon 200 mm f/4 AF-D Macro lens and two teleconverters (Nikon $2\times$ and $1.7\times$) at 3000 fps. This neck radius $r$ was tracked in Matlab using partial area subpixel edge detection techniques affording resolutions below $6\ \mathrm {\mu }{\rm m}\ {\rm pixel}^{-1}$ (Trujillo-Pino et al. Reference Trujillo-Pino, Krissian, Alemán-Flores and Santana-Cedrés2013; Trujillo-Pino Reference Trujillo-Pino2019). A complete list of the materials used is included in A.2.

Figure 1.

(a) Experimental set-up. Water/glycerol pendant drops of radius $R\approx 1$ mm and viscosity $\eta _w=[0.035, 0.068, 0.154, 0.705]\ \textrm {Pa}\ \textrm {s}$ are brought into contact with a silicone oil film of thickness $H$ and viscosity $\eta _o = [0.33, 0.73, 1.5, 7.6]\ \textrm {Pa}\ \textrm {s}$. Viscosity ratios are held at $\eta _w/\eta _o = 1/10$. The radius of the post-contact liquid bridge $r$ is quantified using high-speed photography. (b) Initialization of the three-phase flow simulations. The red phase represents the aqueous drop with density $\rho _w$, viscosity $\eta _w$ and radius $R$. The blue phase represents the oil film with density $\rho _o$, viscosity $\eta _o$ and film thickness $H$. The colourless background represents air with density $\rho _a$ and viscosity $\eta _a$.

The drop's initial radius $R$ is held constant, and $Oh$ is varied by changing the oil film viscosity $\eta _o=[0.33, 0.73, 1.5, 7.6]\ \textrm {Pa}\ {\rm s}$. This range of oil viscosities was created by mixing 20 cSt and 1000 St silicone oils. Drops consisted of deionized (DI) water and glycerol mixtures with viscosities $\eta _w=[0.035, 0.068, 0.154, 0.705]\ {\rm Pa}\ {\rm s}$ chosen to produce drop–film viscosity ratios $\eta _w/\eta _o\approx 0.1$. All viscosities were characterized on an Anton-Paar 702 Rheometer with parallel plate geometry from shear rates $0.1$–$100$ s$^{-1}$, and all mixtures were confirmed to behave as Newtonian fluids.

Film thickness ratios were varied as $H/R = [0.2\unicode{x2013}4.4]$ by varying the film's thickness $H$. For $H < 1$ mm, a glass slide cleaned with DI water and isopropyl alcohol was placed on a laboratory balance, and a prescribed mass of silicone oil was deposited on the slide. The slide was then placed on a level surface, and the oil was allowed to spread fully until it became pinned by the slide's corners, a phenomenon that occurs even for fluids that wet a surface perfectly (Oliver, Huh & Mason Reference Oliver, Huh and Mason1977), forming a level oil–air interface of known thickness (Brochard-Wyart et al. Reference Brochard-Wyart, Hervet, Redon and Rondelez1991; Noblin, Buguin & Brochard-Wyart Reference Noblin, Buguin and Brochard-Wyart2004). For $H \geqslant 1$ mm, glass slides were glued together with a hot-glue gun to form a rectangular well. The well depth was controlled using machinery shims. A prescribed mass of oil was again added using a laboratory balance. To avoid complications in side view imaging from the concave meniscus formed by the oil on the glass wall, care was taken to ensure that the desired film thickness was slightly larger than the well depth, so imaging would occur on a level convex meniscus not obstructed by the oil's contact line.

The surface tension between the oil phase and the glycerol/water phase was measured to be $\sigma _{ow} \approx 0.02\ \textrm {N}\ \textrm {m}^{-1}$ by comparing the change in contact angle $\alpha$ of a sessile glycerol/water drop on a clean glass slide that is surrounded by either air or silicone oil. The Young–Dupré equation was then applied, assuming that the surface tension between the air and the glass is $0$ and that $\alpha = 0^\circ$ for oil on glass.

2.2. Simulation set-up and parameters

Numerical simulations of the drop engulfment process were performed in a rectangular domain; see figure 1. A water/glycerol drop with radius $R$, density $\rho _w$ and viscosity $\eta _w$ is placed in contact with an oil film. The centre of the drop is placed at $H+R+\delta$, where $H$ is the film thickness, $\delta$ is the diffuse interface thickness, $\rho _o$ is the density, and $\eta _o$ is the oil's viscosity. Air, with density $\rho _ a$ and viscosity $\eta _ a$, makes up the background phase. Table 1 provides the relative magnitudes of the parameters and the range of parameters explored in the simulations and experiments. Densities and viscosities are fixed in the simulations: $\rho _w/\rho _a=\rho _o/\rho _a=100$ and $\eta _o/\eta _a=1000$, $\eta _w/\eta _a=100$. The principal interest is to characterize the viscous and inertial engulfment dynamics for different ratios $H/R$, which is achieved by exploring the parameter space $Oh=[0.01\unicode{x2013}40]$ and $H/R=[0.03\unicode{x2013}4]$; see table 1.

Table 1.

Physical properties of the liquids used in the ternary flow system.

We consider a ternary flow system that is composed of a water/glycerol drop, an oil film and air, where the spreading coefficients determine the equilibrium contact angles (Pannacci et al. Reference Pannacci, Bruus, Bartolo, Etchart, Lockhart, Hennequin, Willaime and Tabeling2008; Guzowski et al. Reference Guzowski, Korczyk, Jakiela and Garstecki2012). The spreading factor can be calculated as $S_o=\sigma _{wa}-\sigma _{wo}-\sigma _{ao}$. Since the surface tension of oil is typically lower than that of water/glycerol, the silicone oil will perfectly wet the water/glycerol surface. Thus we consider only spreading coefficients $S_o\geqslant 0$, when the oil fully engulfs the water/glycerol drop at equilibrium.

Inherent to the phase field model are diffuse interfaces with finite thicknesses. Here, we ensure that the interface thickness is small by using a Cahn number $Cn=\delta /R\leqslant 0.05$, which achieves the sharp interface limit by performed convergence tests shown in § A.3. The convergence tests also indicate that the oil film thickness does not affect the results, even for $H/R=0.1$. Since we are interested in drops with radii smaller than the capillary length, we are justified in neglecting the gravitational force in our simulations. From the simulations, we extract the fluid velocity fields, the interfacial shapes, the position of the apparent contact line, and the drop's centre of mass calculated as $C_m=\sum \phi _w y/({\rm \pi} R^2 )$, where $\phi _w$ denotes the order parameter of the water droplet, and $y$ is the vertical displacement scaled by the drop diameter $2R$.

A detailed comparison between two-dimensional (2-D) and three-dimensional (3-D) simulations is conducted in § A.1, which shows that there are only minor differences between the two. For the 3-D simulations, the initial movement of the contact line is also restricted by another curvature, which makes the evolution slower at short times than in the 2-D case. Nevertheless, from the comparison, it does not affect the power law for the contact liquid bridge growth. Thus we focus the article on 2-D simulations using a $D2Q9$ lattice, as these give a good qualitative description of the flow.

Although numerous validation tests of the conservative phase field method have been conducted, we admit that the method still has drawbacks. For example, simulations of large interface deformations in the inertial regime may produce non-physical ghost droplets of oil along the water droplet's interface near the water/oil/air contact line. The reason for these ghost droplets might be pressure oscillations due to the inertial regime's low viscosities and large capillary forces. In addition, the initial conditions are such that there is a small overlap between the droplet and film interfaces, which causes inaccuracy in the surface tension calculation. As the interfaces evolve, this overlap will quickly disappear. Other possibilities include the enforcement of a constant mobility, which is a non-physical parameter in the conservative phase field method (Sun & Beckermann Reference Sun and Beckermann2007). In order to mitigate non-physical droplets that are not present in the data presented here, we consider an oil spreading factor $S_o=0$ in the inertial regime in our simulations, and use the experimental value $S_o=0.5\sigma _{wa}$ in the viscous regime; see table 1.

2.3. Conservative phase field lattice Boltzmann method

The velocity–pressure LBM (Baroudi & Lee Reference Baroudi and Lee2021) is utilized to solve the mass and momentum transport equations, and the conservative phase field LBM is employed to capture the interfaces in order to simulate the three-phase engulfment process in two and three dimensions. In the following equations, the order parameter $0\leqslant \phi _i\leqslant 1$ is used to separate different phases, where $\phi _i=1$ represents the bulk region of component $i$, and $0<\phi _i<1$ represents the interface region. In our simulation, we have three order parameters: $\phi _o$ for the oil film, $\phi _w$ for the water droplet, and $\phi _a$ for the air. Since the summation of three order parameters is $\sum _i\phi _i=1$, we here solve for only $\phi _w$ and $\phi _o$, and obtain $\phi _a$ from $\phi _a=1-\phi _o-\phi _w$. For order parameter $\phi _ i$, the discrete Boltzmann equation can be expressed as

(2.1)\begin{equation} \left(\frac{\partial}{\partial t}+\boldsymbol{e}_\alpha\boldsymbol{\cdot}\boldsymbol{\nabla}\right)h_\alpha^i={-}\frac{1}{\lambda_\phi}\,(h_\alpha^i-h_\alpha^{i,eq})+{S_{\alpha}^i}, \end{equation}

where $h_\alpha ^i$ represents the particle distribution function for order parameter $\phi _i$ along the $\boldsymbol {e}_\alpha$ direction. Here, $\boldsymbol {e}_\alpha$ denotes the discrete velocity vector in a $D2Q9$ or $D3Q27$ lattice (Lee & Lin Reference Lee and Lin2005), and $\lambda _\phi$ is the relaxation time. The equilibrium distribution function $h_\alpha ^{i,eq}$ can be calculated as

(2.2)\begin{equation} h_\alpha^{i,eq}=\phi_i\varGamma_\alpha, \end{equation}

where $\varGamma _\alpha$ is the Maxwell–Boltzmann equilibrium distribution function presented by the method of Hermite polynomials and discretized by Gauss quadrature (He & Luo Reference He and Luo1997), and $\varGamma _\alpha$ takes the form

(2.3)\begin{equation} \varGamma_\alpha=t_\alpha\left[1+\left(\frac{\boldsymbol{e}_\alpha \boldsymbol{\cdot}\boldsymbol{u}}{c_s^2}+\frac{\left(\boldsymbol{e}_\alpha \boldsymbol{\cdot}\boldsymbol{u}\right)^2}{2c_s^4}-\frac{\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{u}}{2c_s^2}\right)\right], \end{equation}

where $t_\alpha$ is the weight, $c_s$ is the speed of sound, and $\boldsymbol {u}$ is the macroscopic velocity vector of the flow field (He & Luo Reference He and Luo1997). The last term, $S_\alpha ^i$, is the intermolecular forcing term

(2.4)\begin{equation} {S}_\alpha^i=\varGamma_\alpha(\boldsymbol{e}_\alpha-\boldsymbol{u})\boldsymbol{\cdot}\left(\frac{4 }{\delta}\,\frac{\boldsymbol{\nabla}\phi_i}{|\boldsymbol{\nabla}\phi_i|}\,\phi_i(1-\phi_i)- \frac{\phi_i^2}{\sum_{j=1}^3\phi_j^2}\sum_{j=1}^3\frac{4 }{\delta}\,\frac{\boldsymbol{\nabla}\phi_j}{|\boldsymbol{\nabla}\phi_j|}\,\phi_j(1-\phi_j)\right), \end{equation}

where $\delta$ denotes the finite interface thickness.

The discrete Boltzmann equation for mass and momentum transport can be expressed as

(2.5)\begin{equation} \left(\frac{\partial }{\partial t}+\boldsymbol{e}_\alpha \boldsymbol{\cdot}\boldsymbol{\nabla}\right)g_\alpha ={-}\frac{1}{\lambda}\,(g_\alpha-g_\alpha^{eq})+F_\alpha, \end{equation}

where $g_\alpha$ represents the particle distribution function for the pressure, with the function

(2.6)\begin{equation} g_\alpha^{eq}=t_\alpha\bar{p}+\varGamma_\alpha c_s^2 -t_\alpha c_s^2, \end{equation}

where $\bar {p}=P/\rho$. Here, $P$ represents the dynamic pressure, and $\rho$ is the mixture density. Also, $F_\alpha$ represents the forcing term

(2.7)\begin{align} F_\alpha&={-}\varGamma_\alpha (\boldsymbol{e}_\alpha-\boldsymbol{u})\boldsymbol{\cdot}\left(\frac{1}{\rho}\,\boldsymbol{\nabla} P\right)+\varGamma_\alpha(0)\,(\boldsymbol{e}_\alpha-\boldsymbol{u})\boldsymbol{\cdot}\left(\boldsymbol{\nabla} \bar{p}\right)\nonumber\\ &\quad +\varGamma_\alpha (\boldsymbol{e}_\alpha-\boldsymbol{u})\boldsymbol{\cdot}\left[\frac{\nu}{\rho}\,( \boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{\rm T})\, \boldsymbol{\nabla}\rho+\frac{1}{\rho}\,\boldsymbol{F}_s\right], \end{align}

where $\nu$ is the kinematic viscosity, and $\boldsymbol {F}_s$ is the surface tension force. In our simulations, the continuum surface tension force model is employed (Kim Reference Kim2005):

(2.8)\begin{equation} \boldsymbol{F}_{s}={-}\frac{3\delta}{2}\sum_i\sigma_i\,\boldsymbol{\nabla}\boldsymbol{\cdot} \left(\frac{\boldsymbol{\nabla}\phi_i}{|\boldsymbol{\nabla}\phi_i|}\right) |\boldsymbol{\nabla}\phi_i|\,\boldsymbol{\nabla}\phi_i, \end{equation}

where $\sigma _i=(\sigma _{ij}+\sigma _{ik}-\sigma _{jk})/2$. The gradient terms of the macroscopic values that appear above can be calculated by the second-order isotropic finite difference method (Lee & Lin Reference Lee and Lin2005).

The governing equations recovered through the Chapman–Enskog expansion are

(2.9)$$\begin{gather} \frac{\partial \bar{p}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \bar{p}+ c_s^2\, \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$

(2.10)$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{uu})={-}\frac{1}{\rho}\,\boldsymbol{\nabla} P+\frac{1}{\rho}\,\boldsymbol{\nabla}\boldsymbol{\cdot} \eta\left(\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^ {\rm T}\right)\nonumber\\ \quad -\frac{3\delta}{2\rho}\sum_{i=1}^{3}\sigma_i\, \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\boldsymbol{\nabla}\phi_i}{|\boldsymbol{\nabla}\phi_i|}\right) |\boldsymbol{\nabla}\phi_i|\,\boldsymbol{\nabla}\phi_i, \end{gather}$$

(2.11)$$\begin{gather}\frac{\partial\phi_i}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\phi_i\boldsymbol{u})= \boldsymbol{\nabla}\boldsymbol{\cdot}M\left(\vphantom{\frac{\phi_i^2}{\sum_{j=1}^3\phi_j^2}\sum_{j=1}^3}\boldsymbol{\nabla}\phi_i-\frac{4}{\delta}\, \frac{\boldsymbol{\nabla}\phi_i}{|\boldsymbol{\nabla}\phi_i|}\,\phi_i(1-\phi_i)\right.\nonumber\\ \quad \left.+\frac{\phi_i^2}{\sum_{j=1}^3\phi_j^2}\sum_{j=1}^3\frac{4}{\delta}\,\frac{ \boldsymbol{\nabla}\phi_j}{|\boldsymbol{\nabla}\phi_j|}\,\phi_j(1-\phi_j)\right), \end{gather}$$

and the density and viscosity defined for the three phases are

(2.12)$$\begin{gather} \rho=\sum_{i=1}^3\rho_i\phi_i, \end{gather}$$

(2.13)$$\begin{gather}\eta=\sum_{i=1}^3\eta_i\phi_i. \end{gather}$$

At the boundaries of the domain, we apply symmetric boundary conditions to the left and right boundaries, while the no-slip boundary condition is applied to the top and bottom boundaries. The details of the LBM derivation and Chapman–Enskog expansion can be found in Zhao & Lee (Reference Zhao and Lee2023).