---
abstract: 'Droplet migration in a Hele–Shaw cell is a fundamental multiphase flow problem which is crucial for many microfluidics applications. We focus on the regime at low capillary number and three-dimensional direct numerical simulations are performed to investigate the problem. In order to reduce the computational cost, an adaptive mesh is employed and high mesh resolution is only used near the interface. Parametric studies are performed on the droplet horizontal radius and the capillary number. For droplets with an horizontal radius larger than half the channel height the droplet overfills the channel and exhibits a pancake shape. A lubrication film is formed between the droplet and the wall and particular attention is paid to the effect of the lubrication film on the droplet velocity. The computed velocity of the pancake droplet is shown to be lower than the average inflow velocity, which is in agreement with experimental measurements. The numerical results show that both the strong shear induced by the lubrication film and the three-dimensional flow structure contribute to the low mobility of the droplet. In this low-migration-velocity scenario the interfacial flow in the droplet reference frame moves toward the rear on the top and reverses direction moving to the front from the two side edges. The velocity of the pancake droplet and the thickness of the lubrication film are observed to decrease with capillary number. The droplet velocity and its dependence on capillary number cannot be captured by the classic Hele–Shaw equations, since the depth-averaged approximation neglects the effect of the lubrication film.'
author:
- Yue Ling
- 'Jose-Maria Fullana'
- Stéphane Popinet
- Christophe Josserand
title: 'Droplet migration in a Hele–Shaw cell: Effect of the lubrication film on the droplet dynamics'
---

INTRODUCTION {#sec:intro}
============

Droplet-based microfluidics is a promising tool for performing biochemical or chemical assays or to observe liquid-liquid extraction [@Stone_2004a; @Fair_2007a; @Teh_2008a]. Droplets are unit systems of controlled volume and content, within which mixing, reacting and/or transferring can be achieved. Therefore, a comprehensive understanding of droplet migration in microchannels is essential to many microfluidics applications. Rectangular channels with a width significantly larger than the height are commonly used in microfluidics devices. In such cases, the confinement of the droplet is only in the vertical direction and the droplet is exposed to a two-dimensional Poiseuille flow (as in a Hele–Shaw cell). The migration of the droplet is driven by the ambient fluid with an average inflow velocity $U_f$ and a constant terminal migration velocity $U_d$ will be reached.

The ratio between the horizontal radius of the droplet ($R$) and half of the channel height ($H$) has been shown to have a strong effect on the droplet migration velocity $U_d$ in a microchannel (hereafter we simply refer to $U_d$ as the droplet velocity since the transient process is not of interest in the present study). The overall variation of $U_d$ with $R/H$ can be seen in the experiments by Roberts [*et al.*]{}[@Roberts_2014a], (the experimental results will be shown later in section \[sec:3D\]) and can be divided into three different regimes according to the dominant mechanism:

-   Poiseuille-dominated regime (Regime P), $R/H\ll 1$. In this regime the droplet velocity is dominated by the Poiseuille flow in the Hele–Shaw cell. In the limit of $R/H\to 0$ (as the droplet is a point mass) the droplet at the centerline of the channel will move with the maximum fluid velocity ($U_{f,max}=1.5 U_f$). As $R$ increases the droplet is exposed to the parabolic fluid velocity profile with a decreasing velocity toward the top and bottom walls and the droplet slows down. Theoretical studies on the migration of droplets much smaller than the channel height have been conducted in many previous works, see for example Nadim and Stone [@Nadim_1991a] and Hudson [@Hudson_2010a].

-   Film-dominated regime (Regime F), $R/H\gg 1$: In this regime, the droplet overfills the channel and a thin lubrication film is formed between the droplet and the wall. The droplet loses its roughly spherical geometry and becomes “pancake-like”. Furthermore, the droplet velocity $U_d$ becomes independent of $R/H$ and is mainly dictated by the lubrication film. The film thickness $h$ is governed by the balance between surface tension and viscous stresses. As the former tends to minimize the interface curvature and pushes the droplet interface closer to the wall, the latter tends to open up the gap between the droplet and the wall. Therefore the capillary number ($\mathrm{Ca}=\mu_f U_d/\sigma$, where $\mu_f$ and $\sigma$ is the viscosity of the ambient fluid and the surface tension, respectively), which is the ratio between the viscous and surface tension forces, is critical in determining the film thickness and the droplet mobility in this regime [@Bretherton_1961a]. Beyond the capillary number, the viscosity ratio between the droplet and the ambient fluid has also been shown to have a significant influence on the lubrication film thickness and the flow pattern [@Hodges_2004a].

-   Transition regime (Regime T), $R/H \approx 1$: This regime connects the previous two regimes. As the droplet radius is comparable to half of the channel height, the effect of confinement on the droplet becomes significant. As a result, the droplet velocity decreases dramatically with $R/H$ in this regime.

The three droplets shown in Fig. \[fig:problem\_description\] schematically represent these three regimes. The axes corresponding to the channel length, height, and width are defined as $x$, $y$, and $z$, respectively. The experimental and numerical results of Roberts [*et al.*]{}[@Roberts_2014a] clearly show the variation of droplet velocity in the Poiseuille-dominated regime, and close agreement with theoretical predictions has been achieved. The experimental data in the film-dominated regime is less documented however. Furthermore, simulation results for the transition and film-dominated regimes are not provided, due to the numerical challenge in resolving the thin film [@Roberts_2014a]. The present study aims at improving the current understanding of droplet motion in microchannels, in particular the droplet dynamics in the transition and film-dominated regimes.

![Schematics of droplet migration in a Hele–Shaw cell for different droplet radii. (a) Side-view and (b) top-view. []{data-label="fig:problem_description"}](problem_description){width="0.6\columnwidth"}

In contrast with scenarios in which the droplet is much smaller than the channel height, the dynamics of droplets in the film-dominated regime ([*i.e.*]{}, $R/H\gg1$) are more complicated and less understood. Early efforts in predicting the velocity of a flat bubble in a Hele–Shaw cell go back to the pioneering theoretical work by Taylor and Saffman [@Taylor_1959a]. The Hele–Shaw equations are obtained by averaging the Stokes equations along the channel height (depth), and it can be shown that the bubble velocity $U_d$ is always larger than the average inflow velocity $U_f$. For low capillary numbers, the horizontal cross-section of the bubble is close to a circle (exhibiting a pancake shape similar to the present problem) [@Kopf-Sill_1988a], then $U_d/U_f=2$ according to the Taylor–Saffman theory [@Taylor_1959a; @Tanveer_1986a; @Tanveer_1987a]. In the later experiments by Kopf-Sill and Homsy [@Kopf-Sill_1988a], the bubble velocity was observed to be about an order of magnitude lower than the Taylor–Saffman prediction. Possible reasons for the discrepancy have been explored, such as an excessive drag due to a moving contact line within the film and Marangoni effects due to surfactant transport [@Saffman_1989a; @Park_1994a; @Maruvada_1996a], yet the underlying mechanisms responsible for the low mobility of the bubble is still not fully understood. In the limit when the 2D droplet is perfectly circular, the Hele–Shaw equations have also been solved analytically by Afkhami and Renardy [@Afkhami_2013a] and Gallaire [*et al.*]{}[@Gallaire_2014a] to obtain the droplet velocity for arbitrary viscosity ratios as $$\frac{U_d}{U_f}= \frac{2}{1+\mu_d/\mu_f}\, , \label{eq:HS_drop_vel}$$ where $\mu_d$ is the viscosity of the droplet fluid. Equation [(\[eq:HS\_drop\_vel\])]{} recovers the Taylor–Saffman solution $U_d/U_f=2$ for bubbles ($\mu_d/\mu_f=0$). Nevertheless, since the Hele–Shaw equations ignore the lubrication film, Eq. [(\[eq:HS\_drop\_vel\])]{} does not capture the effect of the film dynamics on the droplet velocity.

Due to the small scale of microchannels, experimental measurements of the detailed flow structure are very difficult. Recently Huerre [*et al.*]{}[@Huerre_2015a] presented quantitative optical measurements of lubrication film thicknesses as small as 20 nm for a channel height $2H=25 \mu$m (For such a thin film, the intermolecular forces become important for the film dynamics and droplet velocity. In the present study, we focus on the cases where the minimum film thickness is always larger than one micrometer. As a result, the effects of intermolecular forces can be ignored.) However it is very challenging to probe local values of velocity and pressure, and consequently the flow structure inside and outside of the droplet is hard to obtain. Numerical simulations which provide a complete description of the flow field can shed light on the present problem. It has been shown that numerical approximations of the Navier-Stokes equations can represent flow physics in many microfluidics applications, such as droplets/bubbles migration and the breakup of droplets at T-junctions [@Afkhami_2011a; @Bedram_2013a; @Hoang_2013a; @Hoang_2013b].

Nevertheless 3D direct numerical simulations of droplet migration in a microchannel are challenging as well. This is in particular true for droplets with large ratio $R/H$ and small Ca. Since capillary effects are dominant, the time step used in the temporal integration is controlled by the surface tension term in the momentum equation. For small capillary numbers the time step becomes very small and thus the number of steps for the terminal velocity of droplet to be reached is huge. In addition, at low capillary numbers the thickness of the lubrication film becomes very small (several orders of magnitude smaller than the channel height) [@Huerre_2015a]. Since the mesh resolution and the time step restriction are dictated by the smallest length scale in the problem, the total number of cells and integration steps in 3D simulations are both very large. Three-dimensional simulation results of droplet/bubble in microchannels are rarely seen until the very recent works by Hoang [*et al.*]{}[@Hoang_2013b], Roberts [*et al.*]{}[@Roberts_2014a] and Zhu and Gallaire [@Zhu_2016a]. However, as far as the authors are aware, a 3D direct numerical solution to the Navier-Stokes equations for droplets with a large $R/H$ and with a lubrication film formed has never been reported in the literature.

The goal of the present study is to investigate the droplet migration in a Hele–Shaw cell through 3D direct numerical simulations. The capillary number of the droplet is considered to be small enough so that the horizontal shape of the droplet is near-circular all the time. Particular attention will be paid to the effect of the lubrication film on the droplet dynamics. A two-phase flow solver, Gerris[@Popinet_2003a; @Popinet_2009a], is employed to solve the variable-density Navier-Stokes equation with surface tension. In the Gerris solver, the Volume-of-Fluid (VOF) method is used to resolve the interface between the two phases and an adaptive mesh is used to refine the mesh only in the lubrication film and near the interface, so that the computational cost can be reduced. The governing equations and the numerical methods are presented in section \[sec:equations\]. The simulation results will be shown in sections \[sec:2D\] and \[sec:3D\]. We will first perform simulations of 2D droplets migrating in a microchannel with parameters similar to the 3D case for the purpose of validation, see \[sec:2D\]. The simulation results and discussions for the 3D droplet are then presented in section \[sec:3D\]. Finally, we draw conclusions on the present study.

Governing equations and numerical methods {#sec:equations}
=========================================

Governing equations {#sec:gov_eqns}
-------------------

The one-fluid approach is employed to resolve the two-phase flow, where the phases corresponding to the droplet and the ambient fluid are treated as one fluid with material properties that change abruptly across the interface. The incompressible, variable-density, Navier-Stokes equations with surface tension can be written as $$\begin{aligned}
    \rho (\partial_t{\boldsymbol{u}} + {\boldsymbol{u}} \cdot \nabla {\boldsymbol{u}}) = -\nabla p 
    + \nabla \cdot (2\mu {\boldsymbol{D}}) + \sigma \kappa \delta_s {\boldsymbol{n}}\, , 
    \label{eq:mom} \\
    \nabla \cdot {\boldsymbol{u}} = 0 \, ,
    \label{eq:divfree}\end{aligned}$$ where $\rho$, $\mu$, ${\boldsymbol{u}}$, and $p$ represent density and viscosity, velocity and pressure, respectively. The deformation tensor is denoted by ${\boldsymbol{D}}$ with components $D_{ij}=(\partial_i u_j + \partial_j u_i)/2$. The third term on the right hand side of Eq. [(\[eq:mom\])]{} is a singular term, with a Dirac distribution function $\delta_s$ localized on the interface, and it represents the surface tension force. The surface tension coefficient is $\sigma$, and $\kappa$ and ${\boldsymbol{n}}$ are the local curvature and unit normal of the interface.

The volume fraction $C$ is introduced to distinguish the different phases, in particular $C=1$ in the computational cells with only the ambient fluid (respectively $C=0$ in the droplet), and its time evolution satisfies the advection equation $$\begin{aligned}
    \partial_t C + {\boldsymbol{u}} \cdot \nabla C =0 \, .
    \label{eq:color_func}\end{aligned}$$ The fluid density and viscosity are then defined by $$\begin{aligned}
    \rho  & = C \rho_f + (1-C) \rho_d \, , 
    \label{eq:density} \\
    \mu  & = C \mu_f + (1-C) \mu_d \, .
    \label{eq:viscosity}\end{aligned}$$ where the subscripts $f$ and $d$ represent the ambient fluid and the droplet respectively.

Numerical methods {#sec:numerics}
-----------------

The Navier-Stokes equations (Eqs. [(\[eq:mom\])]{} and [(\[eq:divfree\])]{}) are solved by the open-source package Gerris[@Popinet_2003a; @Popinet_2009a]. In Gerris, a finite volume approach based on a projection method is employed. A staggered-in-time discretisation of the volume-fraction/density and pressure leads to a formally second-order accurate time discretisation. The interface between the different fluids are tracked and followed using a Volume-of-Fluid (VOF) method [@Scardovelli_1999a; @Tryggvason_2011a]. A quad/octree spatial discretisation is used, which gives a very important flexibility allowing dynamic grid refinement into user-defined regions [@Popinet_2003a]. Finally the height-function (HF) method is utilized to calculate the local interface curvature, and a balanced-force surface tension discretization is used [@Francois_2006a; @Popinet_2009a].

Time step restriction {#sec:numerics}
---------------------

Both temporal and spatial schemes of integration are explicit. Therefore, numerical stability requires the time step for integration to satisfy the stability conditions corresponding to the advection, viscous, and surface tension terms, respectively. Due to the small capillary number in the present problem, the time step restriction is dominated by the contribution from the surface tension force. The conventional capillary time step restriction is given by Brackbill [*et al.*]{}[@Brackbill_1992a] as, $$\begin{aligned}
    \Delta t \le \Delta t_{BKZ} = \sqrt{ { \rho \Delta^3 \over \pi \sigma } }\, \label{eq:dt_surf_tension}\, ,\end{aligned}$$ where $\Delta$ represents the cell size. For droplets with large $R/H$, the thickness of the lubrication film is the smallest scale in the problem and thus dictates the smallest cell size and the time step to be used. Taking the conventional scaling of film thickness [@Landau_1942a; @Bretherton_1961a], the grid size can be estimated as $\Delta \sim h\sim H Ca^{2/3}$. (The scaling of the film thickness on capillary number will be discussed in detail in the sections \[sec:2D\] and \[sec:3D\].) Then Eq. [(\[eq:dt\_surf\_tension\])]{} becomes $$\begin{aligned}
    \Delta t_{BKZ} = \sqrt{ { \rho H^3 \mathrm{Ca}^2 \over \pi \sigma } }\, \label{eq:dt_surf_tension_BKZ}\, . \end{aligned}$$ If we assume that the terminal migration velocity is reached after a time of order $H/U_d$, then the total number of time steps can be estimated as $$\begin{aligned}
    \mathcal{N}_t \ge \sqrt{ {\pi \sigma  \over  \rho H \mathrm{Ca}^2 U_d^2} } \sim \mathrm{Ca}^{-3/2} \mathrm{Re}^{-1/2}\, \label{eq:Nt_surf_tension_BKZ}\, , \end{aligned}$$ where $\mathrm{Re}=\rho_f U_d H/ \mu_f$. As the total number of grid points (on a regular Cartesian grid) will scale like $\mathcal{N}_{\Delta} \sim (H/\Delta)^D \sim \mathrm{Ca}^{-2D/3}$ with $D$ the dimension of the problem. Then the total work (number of grid points times number of time steps) for a 3D problem will thus scale as $$\begin{aligned}
    \mathcal{N}_{\Delta} \mathcal{N}_t \sim \mathrm{Ca}^{-7/2} Re^{-1/2}\, \label{eq:totalwork_surf_tension_BKZ}\, . \end{aligned}$$

Given that capillary numbers for microfluidics applications can be as low as $10^{-5}$, we see that direct numerical simulations will be extremely expensive. The adaptive mesh refinement in Gerris allows to refine the mesh only in the lubrification film, which greatly decreases this cost, however this is still not sufficient to reach such low capillary numbers.

A different capillary time step restriction has been proposed by Galusinski and Vigneaux [@Galusinski_2008a]. They suggested that when both $\mathrm{Re}$ and $\mathrm{Ca}$ are very small, the capillary time step should satisfy $$\begin{aligned}
    \Delta t \le  \max(\Delta t_{BKZ}, \Delta t_{STK})\, \label{eq:dt_GV}\, ,\end{aligned}$$ where $\Delta t_{STK} = \mu_f \Delta / \sigma$ is the capillary time step related to the Stokes equation (when inertial effects are negligible). It can then easily be shown that the ratio $\Delta t_{BKZ}/\Delta t_{STK}\sim \mathrm{Re}^{-1/2}\mathrm{Ca}^{-1/6}$. Therefore, if $\mathrm{Re}< \mathrm{Ca}^{1/3}$ the larger time step $\Delta t_{STK}$ can be used. In the present study, however, we are interested in the cases where $\mathrm{Re}\sim O(1)$ and $\mathrm{Ca}$ varies from $10^{-3}$ to $10^{-1}$ (see the detailed values in section \[sec:3D\]). As a consequence, the criterion $\mathrm{Re}< \mathrm{Ca}^{1/3}$ is never satisfied and Eqs. [(\[eq:dt\_surf\_tension\])]{}-[(\[eq:totalwork\_surf\_tension\_BKZ\])]{} remain valid to estimate the time step and computational cost.

2D simulations of droplet migration in a microchannel {#sec:2D}
=====================================================

Before we present the 3D results, 2D simulations for a droplet moving in a microchannel are performed for validation purpose. The simulation setup is schematically shown in Fig. \[fig:problem\_description\_2d\]. The flow is symmetric in the $y$ direction with regard to the axis $y=0$, thus only the upper half of the channel is considered. The inflow from the left boundary is a parabolic velocity profile with mean value $U_f$ and maximum value $U_{f,max}=1.5 U_f$. An outflow boundary condition is applied at the right boundary of the domain; while the top boundary is taken as a no-slip wall condition. The length of the bubble is denoted by $2R$, and $R/H=2$. As the droplet is constrained by the channel, a lubrication film is formed between the droplet and the wall. The channel length is $L_x=16H$, which is long enough for the droplet to reach its terminal velocity $U_d$ before being affected by the outflow boundary.

The physical properties of the ambient fluid and the droplet are listed in Table \[tab:phys\_parameters\], which are similar to the recent experiments by Roberts [*et al.*]{}[@Roberts_2014a] The Reynolds number based on the channel height and average inflow velocity, $\mathrm{Re}_f=\rho_f U_f H/\mu_f$, is taken to be 1.35. The surface tension is varied from 0.0002 to 0.0159 N/m.

![Two-dimensional simulation setup for a droplet in a planar microchannel.[]{data-label="fig:problem_description_2d"}](problem_description_2d){width="0.8\columnwidth"}

   $\hat{\rho}_f$ (kg/m$^3$)   $\hat{\rho}_d$ (kg/m$^3$)   $\hat{\mu}_f$ (Pa s)   $\hat{\mu}_d$ (Pa s)   $\hat{H}$ ($\mu$m)   $\hat{U}_f$ (m/s)
  --------------------------- --------------------------- ---------------------- ---------------------- -------------------- -------------------
              750                        1000              $7.4\times 10^{-4}$         $ 10^{-3}$               $50$               0.0267

  : \[tab:phys\_parameters\] Physical properties of the ambient fluid and the droplet in the present study. [^1]

Bretherton [@Bretherton_1961a] studied theoretically and experimentally the dynamics of a long bubble in a cylindrical channel. In the limit of small capillary numbers, [*i.e.*]{}, $Ca=\mu_f U_d/\sigma \ll 1$, it is shown that the ratio between the film thickness $h$ and the half-height of the channel $H$ scales as $Ca^{2/3}$, as does the ratio between the excess of droplet velocity compared to the mean inflow velocity of the ambient fluid: $$\begin{aligned}
{ h \over H } = { (U_d-U_f) \over U_f} \sim \mathrm{Ca}^{2/3}\end{aligned}$$ the proportionality constant depends on the geometrical configuration (planar or cylindrical channels) and the viscosity ratio, although temperature and surfactants may also modify the surface tension. It can be shown that $h/H$ is identical to $(U_d-U_f)/U_f$ due to mass fluxes balance. In the original work of Bretherton [@Bretherton_1961a], the viscosity of the fluid within the bubble is considered to be zero. An extension to droplets of arbitrary viscosity was given later by Hodges [*et al.*]{}[@Hodges_2004a] The theory for a slender droplet in a cylindrical channel is readily modified for a 2D droplet confined by two parallel plates. The film thickness and the excess of droplet velocity for a 2D droplet in the limits of $\mu_d\to 0$ (bubbles) and $\mu_d\to \infty$ (very viscous droplets) are given by $$\frac{h}{H} = \frac{U_d - U_f }{U_f} \sim 0.643 \ (3 \mathrm{Ca})^{2/3}\, .
\label{eq:scaling_2d_inviscid}$$ and $$\frac{h}{H} = \frac{U_d - U_f }{U_f} \sim 0.51 \ (3 \mathrm{Ca})^{2/3}\, , 
\label{eq:scaling_2d_visc}$$ respectively [@Afkhami_2011a]. Similar scaling relations have also been derived for moving foams as well as bubble trains in micro-devices [@Baroud_2010a].

The multi-scale nature of the problem can be seen through Eqs. [(\[eq:scaling\_2d\_inviscid\])]{} and [(\[eq:scaling\_2d\_visc\])]{}. For a capillary number of $10^{-4}$, $h/H=0.00288$ and therefore the length scales involved in the system span over three orders of magnitudes. This fact highlights the importance of the adaptive mesh refinement technique for the present study, which reduces the overall computational costs by refining the mesh only in the lubrication film and near the interface. The shapes of the droplets and the corresponding streamlines in the droplet reference frame are shown in Fig. \[fig:2d\_results\](a) for different capillary numbers. The color of the streamlines denotes the velocity magnitude in the droplet reference frame. The droplet shape and flow pattern are in full agreement with the theory [@Bretherton_1961a; @Hodges_2004a]. Under the effect of the viscous stress, the radius of curvature at the front meniscus is smaller than that at the rear. A small bump is visible on the top-rear of the droplet. The four stagnation points on the droplet interface (A: top-rear; B: top-front; C: rear; D: front) and the other two inside the droplet (E: near the rear; F: near the front) are clearly seen. For small Ca, the major recirculation (A-E-F-B-A) is located at the center of the droplet, with two minor ones sitting on both sides (A-E-C-A and B-F-D-B). When Ca increases the major recirculation shifts toward the rear of the droplet and the two minor recirculation zones become smaller. When the capillary number reaches $\mathrm{Ca}=0.116$ the stagnation points C and E get very close and the minor recirculation near the rear disappears.

![(a) Droplet shape and the streamlines in the droplet reference frame. The color of the streamlines denotes the velocity magnitude in the droplet reference frame. (b) Relative droplet velocity and film thickness as functions of the capillary number. []{data-label="fig:2d_results"}](2d_results){width="0.99\columnwidth"}

Figure \[fig:2d\_results\](b) shows a quantitative comparison between the numerical and the theoretical results for the droplet velocity and the film thickness. The excess of droplet velocity $(U_d - U_f)/U_f$ (which is the same as $h/H$) are plotted as functions of the droplet capillary number $\mathrm{Ca}$ in log-log scale, for $\mathrm{Ca}$ varying from $1.24 \times 10^{-3}$ to $9.87 \times 10^{-2}$. The droplet velocity in the simulation is evaluated by a weighted-integral over all the cells occupied by the droplet fluid, [*i.e.*]{}, $ U_d =  \int_S u\ (1-C) \ dS / \int_S (1-C)\  dS$. It is seen that the simulation results match the theoretical results quite well. For large $\mathrm{Ca}$ the simulation results are slightly lower than the theoretical prediction. This observation is consistent with the previous studies by Afkhami [*et al.*]{}[@Afkhami_2011a]

3D simulations of droplet migration in a microchannel (Hele–Shaw cell) {#sec:3D}
======================================================================

In this section we present the 3D simulation results of droplet migration in a Hele–Shaw cell. The physical properties of the ambient fluid and the droplet are identical to those of the previous 2D simulations, but the droplet is 3D and the cross-section in the $x$-$z$ plane is close to a circle.

Simulation setup
----------------

The computational domain for a 3D droplet migrating in a Hele–Shaw cell is shown in Fig. \[fig:3d\_setup\]. It is assumed that the droplet motion and the corresponding flow is symmetrical in the $y$ and $z$ directions, with respect to the planes $y=0$ and $z=0$, respectively. Therefore, only a quarter of the overall domain is actually computed. Taking the half-height of the channel $H$ as the typical length scale, the length, height, and width of the computational domain are taken as $L_x=16H$, $L_y=H$, and $L_z=8H$, respectively.

![Three-dimensional simulation setup for a droplet in a Hele–Shaw cell.[]{data-label="fig:3d_setup"}](3d_setup){width="0.9\columnwidth"}

An inflow boundary condition is imposed on the left of the domain with a parabolic velocity profile. Correspondingly an outflow boundary condition is used on the right. The top of the domain is taken to be a no-slip wall moving with a velocity $U_w$ in the opposite direction with respect to the flow. Therefore the simulation is indeed conducted in a moving reference frame with a mean inflow velocity of $U_f-U_w$. The purpose of using a moving reference frame is to keep the droplet near the center of the domain and to avoid the influence of the boundaries. Ideally the terminal droplet velocity in the lab reference frame $U_d$ should be used for $U_w$, but since $U_d$ is unknown *a priori*, an initial guess of $U_d$ is used for $U_w$. Nevertheless as long as the droplet does not get too close to the inflow and outflow boundaries, the specific value of $U_w$ is immaterial to the simulation results. Finally the other boundaries are set to symmetric boundary conditions.

The droplet is taken to be initially stationary in the simulation frame. Driven by the ambient fluid flow the droplet approaches the terminal droplet velocity as time evolves. The initial shape of droplets with $R/H < 1$ is taken to be a sphere with radius $R$. For $R/H \gtrsim 1$ the cross-section at $y=0$ is a circle with radius $R_0$. The top of the droplet is initially flat with a constant distance of $0.1H$ to the top wall. The “edge" of the droplet is formed by a quarter of a circle of radius of $0.9H$ rotated with respect to the $y$ axis, see Fig. \[fig:3d\_setup\]. Similar to the droplet velocity, the droplet will evolve to the final shape for which the surface tension and viscous forces at the interface are in balance. The final value of $R$, based on the horizontal cross-section area $S_{y=0}$ (at the plane $y=0$) as $R=\sqrt{2S_{y=0}/\pi}$, is slightly different from its initial value $R_0$. Due to the initial configuration we impose, the difference between $R$ and $R_0$ is generally small. (Therefore from here on we do not distinguish $R_0$ and $R$ unless otherwise indicated.) Similar to the velocity of the top wall, the specific details of the initial shape of the droplet are irrelevant to the final results as long as the terminal droplet velocity and shape are reached.

An adaptive mesh is used to discretize the spatial domain. The finest cell size used is $\Delta_{\min}=H/128$. The mesh is refined at the interface based on the gradient of the volume fraction in each cell. Mesh refinement is also conducted using an *a posteriori* error estimate of all three velocity components as a cost function for adaptation. A mesh refinement study was performed in the 2D simulations, showing that $\Delta_{\min}$ is sufficient for the range of capillary numbers considered in the present study.

   $r=\rho_f$/$\rho_d$   $m=\mu_f/\mu_d$   $\mathrm{Re}_f=\rho_f U_f H/\mu_f$    $R/H$          $\mathrm{Ca}_f=\mu_f U_f/\sigma$
  --------------------- ----------------- ------------------------------------ ---------- --------------------------------------------
          0.75                0.74                        1.35                  0.5 to 2   $6.17\times10^{-3}$ to $9.87\times10^{-2}$

  : \[tab:dmls\_parameters\] Key dimensionless parameters in 3D simulations.

The physical parameters are identical to the 2D simulation presented in section \[sec:2D\] (see Table \[tab:phys\_parameters\]). Parametric studies are performed by varying the surface tension coefficient $\sigma$ and the droplet radius $R$. Taking the ambient fluid density $\rho_f$, the half of the channel height $H$, and the mean inflow velocity $U_f$ as characteristic scales, the unknown variables and physical parameters can be non-dimensionalized and the governing equations (Eqs. [(\[eq:mom\])]{}–[(\[eq:color\_func\])]{}) are solved in dimensionless form. The key dimensionless parameters of the present problem are given in Table \[tab:dmls\_parameters\]. As shown in the table, the ratio $R/H$ and the capillary number based on the average inflow velocity $\mathrm{Ca}_f$ are chosen as the tuning parameters. The capillary numbers in the 3D simulation are significantly larger than that in experiments by Roberts [*et al.*]{}[@Roberts_2014a], in which $\mathrm{Ca}_{f,\mathrm{Exp}}=9.87\times 10^{-4}$. The smallest $\mathrm{Ca}_f$ considered here is about an order of magnitude larger than this experimental value. As illustrated in the previous 2D results, the lubrication film formed is very small for small capillary numbers. It is prohibitively difficult to simulate droplets with a larger $R/H\gtrsim 1$ and a small $\mathrm{Ca}_f$ comparable to the experimental value with current computational resources. Nevertheless it is expected that as long as the thickness of the lubrication film does not reach the regime where van der Waals forces are important, the effect of the lubrication film on the droplet dynamics can still be captured with simulation results at larger $\mathrm{Ca}_f$. At the end of this section an extensive discussion on the effect of surface tension will be given.

Effect of droplet horizontal radius on droplet dynamics
-------------------------------------------------------

![(a) Temporal evolution of mean droplet velocity. (b) Terminal mean droplet velocity $U_d$, (c) thickness (height) of the gap (film) between the droplet interface and the top channel wall, denoted by $h$, and (d) relative $u$-velocities at the top and the side of the droplet ($U_y$ and $U_z$) as a function of the ratio between the droplet radius $R$ and half of the channel height $H$. The experimental and theoretical results are shown for comparison [@Roberts_2014a; @Nadim_1991a; @Hudson_2010a] . []{data-label="fig:udvsR"}](udvsR){width="0.96\columnwidth"}

The temporal evolution of the droplet velocity is shown in Fig. \[fig:udvsR\](a) for $\mathrm{Ca}_f=9.87\times 10^{-2}$ and $R/H$ varying from 1 to 2. When $R/H$ increases, $U_d$ decreases and in contrast, the time it takes for the droplet to reach its terminal velocity increases (while identical initial conditions are imposed). From $R/H=1$ to 1.25, $U_d$ decreases significantly (about 9%). The decrease of $U_d$ with $R/H$ becomes more gradual for larger $R/H$. It is seen that the difference between $U_d$ for $R/H=1.75$ and $R/H=2$ becomes invisible as $t>2$.

The simulation results of the droplet velocity as a function of $R/H$ are shown in Fig. \[fig:udvsR\](b). The experimental and simulation results of Roberts [*et al.*]{}[@Roberts_2014a] and the theoretical results of Nadim and Stone [@Nadim_1991a] for small $R/H$ are also shown for comparison. The present predictions of $U_d$ are in good agreement with both the theoretical and experimental results for $R/H<1$.

For $R/H>1$ the simulation results of $U_d$ are larger than those measured in the experiment. The discrepancy is expected since smaller surface tension coefficients (or larger capillary numbers) are considered in the simulation. Nevertheless the simulation results capture well the trend of variation of $U_d$ with $R/H$, and the three different regimes (Poiseuille-dominated, transition, and film-dominated regimes) can be clearly seen. In particular, the variation of $U_d$ in the transition and the film-dominated regimes is much better illustrated. In the transition regime ($0.75\lesssim R/H \lesssim 1.25$), the droplet just overfills the channel forming a lubrication film. Correspondingly, the droplet velocity decreases rapidly. In the film-dominated regime, for $R/H>1.5$, $U_d$ is observed to reach a plateau gradually. The plateau value of $U_d$ decreases with $\mathrm{Ca}_f$, approaching the experimental value which corresponds to a much smaller $\mathrm{Ca}_f$.

The three different regimes can also be identified from the variation of film thickness with $R/H$ as shown in Fig. \[fig:udvsR\](c). In the Poiseuille-dominated regime, as the gap between the droplet and the top wall is relatively large, its thickness is not affected by the wall. The droplet remains close to a sphere, so $h/H$ closely follows $1-R/H$. As $R/H$ approaches 1, (in the transition regime,) the film thickness starts to deviate from $1-R/H$ and becomes governed by the viscous stress and the surface tension. When $R/H\gg 1$, (in the film-dominated regime,) the film thickness becomes completely independent from the horizontal radius and becomes dependent on $\mathrm{Ca}_f$.

At last the droplet circulation velocities in the three different regimes are shown in Fig. \[fig:udvsR\](d). The $u$-velocities at the top (the point with maximum $y$ on the droplet) and the side of the droplet (the point with minimum $z$ on the droplet) are denoted as $U_y$ and $U_z$, respectively. The simulation results of Roberts [*et al.*]{}[@Roberts_2014a] and the theoretical results of Hudson [@Hudson_2010a] are shown for comparison. For $R/H \ll 1$, namely in the Poiseuille-dominated regime, $U_y$ and $U_z$ vary as parabolic functions with respect to $R/H$. The present predictions of $U_y$ are in good agreement with both the previous numerical results and the theoretical prediction. In the film-dominated regime ($R/H \gg 1$), it is observed that both $U_y$ and $U_z$ reach constant values, similar to $U_d$. It is also interesting to notice that the relative velocities $U_y-U_d$ and $U_z-U_d$ for different $\mathrm{Ca}_f$ seem to overlap.

A simple model for droplet dynamics
-----------------------------------

An interesting observation from the simulation results is that the velocity of a 3D droplet can vary in a much wider range than the velocity of a 2D droplet. In particular, the droplet velocity in the film-dominated regime is observed to be generally smaller than the average inflow velocity. This low droplet velocity has also been observed in experiments [@Roberts_2014a; @Huerre_2015a].

From Fig. \[fig:udvsR\](b) it can be seen that when $R/H$ increases from 0.5 to 2, the droplet velocity $U_d$ for $\mathrm{Ca}_f=2.47\times 10^{-2}$ decreases from 1.4 to 0.8. In this section, we propose a simple model of mass fluxes balance to qualitatively illustrate the droplet mobility in the film-dominated regime and the resulting 3D flow patterns.

\[sec:model\]

![Schematics of mass fluxes budget for a droplet in (a) 2D and (b) 3D channels. []{data-label="fig:mass_cons"}](mass_cons){width="0.85\columnwidth"}

In the 2D planar (or axisymmetric) geometry the ambient fluid which enters the channel can only move along the flow direction, and thus must push the droplet to move faster than the average inflow velocity to ensure mass conservation. If $h$ is given the corresponding value of $U_d$ can be obtained by the balance of mass fluxes [@Bretherton_1961a]. Taking a control volume (of zero thickness) near the rear of the droplet as shown in Fig. \[fig:mass\_cons\](a), the mass fluxes are identical to zero on both the $y$ and $z$ boundaries and thus the fluxes from the left and right must cancel. Considering a bubble with large $R/H$ and small $h/H$, with a stress-free interface (and thus a constant velocity within the film), the relation between $h$ and $U_d$ can be expressed as $$\frac{U_d-U_f}{U_d}\frac{H}{h} = 1\, . 
    \label{eq:2D_mass_cons}$$ It is seen that $U_d > U_f$ for non-zero $h$. For small $h/H$ the excess of droplet velocity is small and the droplet velocity $U_d$ is close to the average inflow velocity $U_f$. In the 3D analogue, since the ambient fluid can flow around the droplet from the sides (or edge) of the droplet, the mass fluxes balance and the flow pattern becomes more complex. Similarly to the 2D case we look at a thin layer of small thickness $\epsilon \ll h$ along the symmetric plane ($z\in (0,\epsilon)$) near the rear of the droplet. As there is zero flux on the symmetric plane $z=0$, the transverse mass flux is only non-zero at the plane $z=\epsilon$, normalized by the flux from the right, [*i.e.*]{}, $U_d h \epsilon$, which can be expressed as $$F_{z,rear} = 1 - \frac{U_d-U_f}{U_d}\frac{H}{h} \, , 
    \label{eq:3D_mass_cons}$$ assuming the variation of $U_f$ in $z\in (0,\epsilon)$ is negligible. The transverse flux $F_{z,rear}$ is not uniformly distributed on the plane $z=\epsilon$, instead it is concentrated near the stagnation point at the rear of the droplet. A similar analysis can be performed for the front of the droplet and it can be shown that $F_{z,rear}$ and $F_{z,front}$ have the same magnitude but different signs, which creates a circulation around the edge of the droplet.

For a given $h$, $F_{z,rear}$ is related to $U_d$ through mass conservation, however the fluxes balance is not sufficient to determine $U_d$. This extra spatial degree of freedom of mass flux thus allows a much wider variation of droplet velocity. In contrast to the 2D case, where $U_d$ is always slightly larger than $U_f$, in the 3D case $U_d/U_f$ has been observed to vary from as low as $0.1$ to about 2 in experiments [@Park_1994a; @Huerre_2015a].

Furthermore the two scenarios of $U_d>U_f$ and $U_d<U_f$ correspond to two qualitatively different three-dimensional flow patterns over the droplet. Assuming $h\ll H$, it can be shown that $F_{z,rear}<0$ when $U_d>U_f$. The “pure" bubble experiment by Park [*et al.*]{}[@Park_1994a] falls in this scenario of high-migration-velocity (HMV). In the droplet reference frame the sweeping motion of the ambient fluid on the edge of the droplet is from the front to the rear. Combined with the fact that the flow on the top of the droplet is always from the front to the rear, this means that the overall flow in the droplet reference frame is in the rear-to-front direction.

On the other hand the scenario of low-migration-velocity (LMV), [*i.e.*]{}, $U_d<U_f$ and thus $F_{z,rear} > 0$, is very different. The droplets in the present simulation and in the experiments of Roberts [*et al.*]{}[@Roberts_2014a] and Huerre [*et al.*]{}[@Huerre_2015a] fall in this category. The positive transverse flux near the droplet rear ($F_{z,rear}>0$) indicates that the fluid actually goes around the droplet from the rear to the front in the droplet reference frame. As a consequence, the relative flows on the top and the side of the droplets are in opposite directions, (which are also illustrated in Fig. \[fig:udvsR\](d)). A complex interfacial flow sweeping over the droplet is then formed, in which the flow (in the droplet reference frame) first goes from the front to the rear on the top of the droplet then returns from the rear back to the front on the edge of the droplet. (A visualization of this interfacial flow will be presented in Fig. \[fig:interfacial\_vel-vector\_dropframe\].)

Effect of three-dimensional flow structure on droplet dynamics
--------------------------------------------------------------

### Streamlines and flow circulation

![Streamlines inside and outside of the droplet in the droplet reference frame for $\mathrm{Ca}_f=9.87\times 10^{-2}$ and $R/H=0.5$, 1 and 2. The color of the streamlines denotes the corresponding relative velocity magnitude. The pink line indicates the droplet interface.[]{data-label="fig:streamlines"}](streamlines){width="1\columnwidth"}

![Streamlines for $\mathrm{Ca}_f=9.87\times 10^{-2}$ and $R/H=2$ in the droplet reference frame: (a) Internal circulation and external circulation near the front of the droplet. (b) Flow within the film between the droplet and the wall and the internal circulation inside the droplet. The color of the streamlines denotes the magnitude of the relative velocity. The angular coordinate $\theta$ is defined in the droplet reference frame with the origin located at the center of the horizontal cross-section.[]{data-label="fig:streamlines_3D"}](streamlines_3D){width="1.0\columnwidth"}

The three-dimensional flow structure inside and outside of the droplet is shown in Fig. \[fig:streamlines\] for $\mathrm{Ca}_f=9.87\times 10^{-2}$ and $R/H=0.5$, 1 and 2. The streamlines on the two symmetric planes $z=0$ and $y=0$ are plotted in the droplet reference frame and colored by the relative velocity magnitude. As the droplet moves faster for smaller droplets the relative velocity magnitude is larger for smaller $R/H$. The flow pattern for $R/H<1$ here is also consistent with the simulation results of Roberts [*et al.*]{}[@Roberts_2014a].

It is interesting to compare the 3D flow pattern for the droplet with $R/H=2$ to its 2D analogue shown in Fig. \[fig:2d\_results\]. In the vertical plane $z=0$ the ambient fluid is moving from right to left at the top ($y=H$) and in a reverse direction at the center plane ($y=0$). As a result two external circulations develop outside of the droplet, similar to the 2D analogue shown in Fig. \[fig:2d\_results\](a). The external flow will also drive internal circulation inside the droplet. Similarly to the 2D case, six stagnation points (A, B, C, D on the droplet interface and E, F inside the droplet) can be seen, however the stagnation flow around these points becomes three dimensional. For example in the 2D case, the flow reaches the rear stagnation point C and can only turn toward point A and then turn again to the left and leave the domain, forming a counter-clockwise external circulation, see Fig. \[fig:2d\_results\](a). In 3D the fluid has an extra degree of freedom and can go in the $z$ direction around the droplet toward the stagnation point D, see the streamlines on the $x$-$z$ plane in Fig. \[fig:streamlines\](c). This transverse fluid motion near the stagnation point C is characterized by $ F_{z,rear} $ in the model presented in section \[sec:model\]. Due to the 3D flow structure the minor internal circulations (A-E-C-A and B-D-F-B) in the 2D case are invisible here.

The transverse fluid motion significantly affects the external and internal circulating flow. As shown in Fig. \[fig:streamlines\_3D\](a), the 3D streamlines inside and outside of the droplet become very complicated and are very dependent on the angular coordinate $\theta$. (Here $\theta$ is defined as the angle between the normal to the interface in the $x$-$z$ plane and the $x$ axis.) For $\theta=0$ the circulation is perfectly aligned with the symmetric $x$-$y$ plane, and the flow pattern is similar to the 2D case. As $\theta$ increases the external circulation is twisted in a complex manner. The incoming flows for nonzero $\theta$ near the top of the channel are observed to turn toward the center plane $z=0$, (indicated by black arrows in Fig. \[fig:streamlines\_3D\](a)). As the flow approach the droplet interface, it turns to the side to adjust to the droplet shape, (see the magenta arrow). At about $\theta \approx 60$ degrees the streamlines outside of the droplet are almost tangential to the droplet interface. Therefore, in contrast with the typical quasi-2D flows in a Hele–Shaw cell, droplet migration in the microchannel generates a fully 3D flow.

### Velocity field on the droplet interface

![(a) Top- (b) front- and (c) rear-views of relative velocity in the droplet reference frame $\{u-U_d,v,w\}$ on the droplet interface for $\mathrm{Ca}_f=9.87\times 10^{-2}$ and $R/H=0.5$, 1 and 2. The color denotes the magnitude of the relative velocity. []{data-label="fig:interfacial_vel-vector_dropframe"}](interfacial_vel-vector_dropframe){width="0.99\columnwidth"}

The vectors of the interfacial velocity in the droplet reference frame are shown in Fig. \[fig:interfacial\_vel-vector\_dropframe\]. The color of the interface indicates the relative velocity magnitude. Regions with white color indicate zero velocity, [*i.e.*]{}, stagnation regions on the droplet interface. The four stagnation points (A, B, C, D) on the droplet interface and 3D stagnation flows in the vicinity can be clearly seen.

Driven by the flow within the film, the interfacial flow on the top of the droplet moves to the left; while the flow on the edge near the center plane $y=0$ moves in the opposite direction. An interesting and complex interfacial flow is then formed on the droplet interface, which is highlighted with yellow dashed lines in Fig. \[fig:interfacial\_vel-vector\_dropframe\]. At first the flow starts from the vicinity of stagnation point B (indicated by the blue dashed line $\mathrm{B}_1$$\mathrm{B}_2$ ), and goes toward the rear. It is seen that the velocity vectors on the top interface are not parallel to the $x$ axis, instead they turn toward the side of the droplet. The angle between the top interfacial flow and the $-x$ axis is denoted by $\alpha$ here, and it can be seen that $\alpha$ varies significantly with the angle $\theta$ of the relative flow near the wall (in the droplet reference frame). For $\theta=0$, $\alpha=0$, but then $\alpha$ increases significantly with $\theta$ up to about 45 degrees. Then top interfacial flow is observed to curve back toward the center plane ($z=0$) before reaching stagnation point A. After reaching point A the flow reverses its direction and moves toward the front stagnation point D along the edge of the droplet. For the top interfacial flow leaving from $\mathrm{B}_1$$\mathrm{B}_2$ with a large $\theta$, the flow may not get a chance to turn back to stagnation point A before returning to the front. As a result a ring of zero interfacial velocity (the region in white color) is formed on the edge of the top where the interfacial flow reverses direction. Near the rear of the droplet the expected up-going interfacial flow from stagnation points C to A is not seen due to the strong interfacial flow along the edge. The stagnation flows near C and A seem to merge, forming a stagnation area at the droplet rear.

Furthermore, the streamlines in Fig. \[fig:streamlines\_3D\](b) show that the flow within the lubrication film is influenced by the complex interfacial flow on the droplet interface. Similar to the interfacial flow on the top of droplet (see Fig. \[fig:interfacial\_vel-vector\_dropframe\]), when the flow passes through the lubrication film, the streamlines off the center plane $z=0$ also first turn toward the side and then return back to the streamwise direction.

Configuration and dynamics of the lubrication film
--------------------------------------------------

![(a) Thickness of the lubrication film between the interface and the wall $h/H$ and (b) shear $\partial u/\partial y$ on the interface for $\mathrm{Ca}_f=9.87\times 10^{-2}$ and $R/H=0.5$, 1 and 2. The film for $R/H=2$ is highlighted by pink dotted lines. []{data-label="fig:interface_height_shear"}](interface_height_shear){width="0.85\columnwidth"}

As the film between the droplet and the wall serves to “lubricate" the droplet, when the film thickness decreases (for example due to increasing surface tension) it is expected the droplet will slow down correspondingly. The increase of droplet velocity with the film thickness has been confirmed in the experiments by Huerre [*et al.*]{}[@Huerre_2015a]. To better understand the effect of the lubrication film on the droplet dynamics, we plot the thickness of the lubrication film and the shear on the droplet top interface in Fig. \[fig:interface\_height\_shear\] for $\mathrm{Ca}_f=9.87\times 10^{-2}$ and $R/H=0.5$, 1 and 2.

As $R/H$ is small the gap between the droplet and the wall is large and thus strictly speaking there exists no film, the term is retained only for convenience. When $R/H$ is larger or comparable to unity a thin lubrication film with small thickness can be clearly seen. In particular for $R/H=2$ the film area (highlighted by the pink dotted line) becomes quite large compared to the cross-section of the droplet at the $y=0$ plane. Instead of being centered, the lubrication film is shifted toward the rear (see Fig. \[fig:interfacial\_vel-vector\_dropframe\]). The overall film exhibits a “peach" shape.

Values of the shear $\partial u/\partial y$ on the droplet interface are shown in Fig. \[fig:interface\_height\_shear\](b). When $R/H$ is small the droplet is away from the top wall and the shear on the top of the droplet is indeed quite small. The shear increases dramatically as $R/H$ becomes comparable or larger than unity. For $R/H=2$ it is observed that the shear on the droplet is not uniformly distributed and the large negative shear is located near the rear. The region of high shear has a croissant shape, and minimum shear rate reaches about $(\partial u/\partial y)_{\min}=-8$ (in dimensionless unit), the magnitude of which is more than twice larger than the minimum shear in the ambient 2D Poiseuille flow, [*i.e.*]{}$(\partial u/\partial y)_{\min, \mathrm{Poiseuille}}=-3$. It is obvious that the strong shear stress induced by the film is closely related to the low mobility of droplets with larger $R/H$.

Effect of capillary number on film and droplet dynamics
-------------------------------------------------------

As mentioned at the beginning of this section, the capillary numbers considered in the present 3D simulations are larger than that in experiments. To justify that the simulation results presented above at least qualitatively reflect the correct flow physics of droplet migration in a Hele–Shaw cell, the capillary number is varied to investigate its effect on the results. From the droplet velocity previously shown in Fig. \[fig:udvsR\](b), the change in $U_d$ is invisibly small for $R/H=0.5$ as $\mathrm{Ca}_f$ varies from $9.87\times 10^{-2}$ to $6.17\times 10^{-3}$. Though smaller capillary numbers are used, the surface tension is still sufficient to keep the droplet spherical and the droplet velocity is then mainly determined by the viscous stress. The droplet velocity obtained in simulations for $\mathrm{Ca}_f=9.87\times 10^{-2}$, which is two orders of magnitude larger than that in experiments, still agrees quite well with the experimental and theoretical results. On the other hand, as $R/H\gtrsim 1$, the effect of the capillary number on the droplet velocity becomes much stronger. As shown in Fig. \[fig:udvsR\](b) and (c), the droplet velocity and film thickness decrease by about $12.2\%$ and $53.3\%$ for $R/H=2$, when $\mathrm{Ca}$ decreases from $9.87\times 10^{-2}$ to $2.47\times 10^{-2}$.

The film thickness and shear on the droplet interface for $\mathrm{Ca}_f=2.47\times 10^{-2}$ and $R/H=0.5$, 1 and 2 are shown in Fig. \[fig:interface\_height\_shear\_SmallST\]. The surface tension here is four times that in Fig. \[fig:interface\_height\_shear\] while the other parameters remain the same. (Correspondingly $\mathrm{Ca}_f$ in Fig. \[fig:interface\_height\_shear\_SmallST\] is a quarter of that in Fig. \[fig:interface\_height\_shear\].) When $R/H=0.5$ there is no difference in the gap thickness. For $R/H=2$ the film becomes significantly thinner as surface tension increases. (Therefore different legends for film thickness are used in in Figs. \[fig:interface\_height\_shear\](a) and \[fig:interface\_height\_shear\_SmallST\](a) for better visualization.) The minimum film thickness for $\mathrm{Ca}_f=2.47\times 10^{-2}$ is about $h_{\min}=0.015H$ while for $\mathrm{Ca}_f=9.87\times 10^{-2}$, $h_{\min}=0.048H$. It is also observed that when $\mathrm{Ca}_f$ decreases, the lubrication film still shifts to the rear, however the shift becomes smaller and the droplet shape becomes more symmetric with regard to the $y$-$z$ plane. Furthermore, the overall “peach" shape of the film remains, but the smallest film thickness for $R/H=2$ seems to move to the two sides of the film.

When $\mathrm{Ca}_f$ decreases the general distribution of the shear rate within the film is similar. The high-shear region still holds a croissant shape. The minimum shear near the rear decreases to about $-10$ which is about 25% larger than that for $\mathrm{Ca}_f=9.87\times 10^{-2}$ in magnitude. This indicates that the decreases of $\mathrm{Ca}_f$ results in a decrease of $h$ which in turn introduces a stronger shear drag on the droplet and reduces the droplet migration velocity.

![(a) Thickness of the lubrication film between the interface and the wall $h/H$ and (b) shear $\partial u/\partial y$ on the interface for $\mathrm{Ca}_f=2.47\times 10^{-2}$ and $R/H=0.5$, 1 and 2. []{data-label="fig:interface_height_shear_SmallST"}](interface_height_shear_SmallST){width="0.85\columnwidth"}

The minimal film thickness on the symmetric $x$-$y$ plane is plotted as a function of $\mathrm{Ca}$ in Fig. \[fig:ud-h\_vs\_Ca\](a), where $\mathrm{Ca}=\mu_f U_d/\sigma$ is the capillary number based on the droplet velocity $U_d$. The 2D theory, [*i.e.*]{}, Eq. [(\[eq:scaling\_2d\_inviscid\])]{}, is also plotted for comparison. It is seen that $h_{z=0,\min}/H$ in the 3D case is much smaller than in the 2D case. For $\mathrm{Ca}=9.87\times 10^{-2}$, $h_{z=0,\min}/H$ in 3D is an order of magnitude smaller than in 2D, as already visualized in Figs. \[fig:2d\_results\] and \[fig:streamlines\]. The difference between 3D and 2D film thickness seems to decrease with $\mathrm{Ca}$. The film thickness variation with $\mathrm{Ca}$ seems to follow a power law but the power coefficient is about 0.4 instead of 2/3 as in 2D.

The droplet velocity $U_d$ decreases with $\mathrm{Ca}_f$ as shown in Fig. \[fig:ud-h\_vs\_Ca\](b). Here the capillary number based on the average inflow velocity is used following previous works [@Kopf-Sill_1988a; @Park_1994a]. The droplet velocity obtained by the Hele–Shaw equations (Eq. [(\[eq:HS\_drop\_vel\])]{}), which is about 0.85, is also plotted for comparison. The depth-averaged approximation misses the effect of the lubrication film on the droplet dynamics, and thus Eq. [(\[eq:HS\_drop\_vel\])]{} fails to capture the variation of the droplet velocity with capillary number. It is seen that the Hele–Shaw solution under-predicts $U_d/U_f$ for larger $\mathrm{Ca}_f$ and over-predicts $U_d/U_f$ for small $\mathrm{Ca}_f$. The droplet velocity only matches the Hele–Shaw solution at $\mathrm{Ca}_f\approx 0.05$. Similar observations have also been made by Zhu and Gallaire [@Zhu_2016a]. The range of $\mathrm{Ca}_f$ considered in the present simulation is an order of magnitude larger than the experimental value. It is hard to tell if the simulation results would match the experimental measurements [@Roberts_2014a] if the same $\mathrm{Ca}_f$ is used. Nevertheless it seems like the computed $U_d$ would approach the experimental value as $\mathrm{Ca}_f$ decreases.

As discussed previously the mass fluxes balance in 3D is substantially different from 2D. In the 2D case it is known that $(U_d-U_f)/U_f=h/H$ (for $h/H\ll1$). In 3D, $(U_d-U_f)/U_f$ becomes negative. It is shown that $U_d$ also decreases with $h/H$, but the decrease is non-linear and follows a power law about $(h/H)^{1/4}$ as indicated in Fig. \[fig:ud-h\_vs\_Ca\](c).

The simple model of transverse flux near the rear of the droplet (Eq. [(\[eq:3D\_mass\_cons\])]{}) can serve as a measure of the three-dimensional “strength” of the flow. The variation of $F_{z,rear}$ is plotted as a function of the film thickness in Fig. \[fig:ud-h\_vs\_Ca\](d). When $h$ decreases it is observed that $F_{z,rear}$ increases rapidly. As the film becomes thinner, the drag induced by the lubrication film becomes stronger and it is more difficult for the ambient fluid to push the droplet forward. As a result the ambient fluid entering the channel tends to go around the droplet from the transverse edge instead of pushing the droplet. Therefore it is quite clear that both the 3D flow effect and the lubrication film dynamics control the mobility of a droplet in a Hele–Shaw cell.

![(a) The minimum film thickness at the $z=0$ plane and (b) droplet velocity as functions of the capillary number. The droplet velocity obtained by the Hele–Shaw equations (Eq. [(\[eq:HS\_drop\_vel\])]{}) is also plotted for comparison. (c) Excess of droplet velocity compared to the average inflow velocity and (d) Normalized transverse flux near the rear of the droplet as a function of the film thickness. []{data-label="fig:ud-h_vs_Ca"}](ud-h_vs_Ca){width="0.99\columnwidth"}

Discussions of computational costs
----------------------------------

The present simulations are performed on 8 to 24 processors for 3 to 25 days depending on cases. The adaptive meshes used for $\mathrm{Ca}_f=9.87\times 10^{-2}$ and $2.47\times 10^{-2}$ and $R/H=2$ at the plane $z=0$ are shown in Fig. \[fig:2D-interface\_with-mesh\]. The main constraint on the computational efficiency comes from the small time step which is in turn restricted by the surface tension (or capillary number), see Eq. [(\[eq:dt\_surf\_tension\])]{}. When surface tension increases four times, the time step is divided by two and namely the computational time doubles assuming the mesh is unchanged. However, as indicated in Fig. \[fig:2D-interface\_with-mesh\](b) that if $\mathrm{Ca}_f$ is further decreased, the minimum mesh size in the adaptive mesh, [*i.e.*]{}, $\Delta_{min}=H/128$, will not be sufficient and thus must be refined. Therefore, even with an adaptive mesh, it is too expensive to simulate droplets at capillary number comparable to experiments for large $R/H$. To perform 3D direct simulation of droplets at low capillary number, novel numerical schemes to compute the surface tension are required. If the surface tension term in Eq. [(\[eq:mom\])]{} can be treated in an implicit manner similar to what is done for the viscous term in Stokes flow (at low Reynolds number), then the time-step restriction due to the explicit integration of the surface tension term can be lifted (the authors are not aware of this type of scheme in the literature). Nevertheless the development of an implicit treatment of surface tension is out of the scope of the present paper and will be relegated to future work.

![The adaptive mesh on the plane $z=0$ for $R/H=2$ and (a) $\mathrm{Ca}_f=9.87\times 10^{-2}$ and (b) $\mathrm{Ca}_f=2.47\times 10^{-2}$. []{data-label="fig:2D-interface_with-mesh"}](2D-interface_with-mesh){width="1.0\columnwidth"}

CONCLUSIONS
===========

The migration of a droplet in a Hele–Shaw cell has been investigated by 3D direct numerical simulations. We focused on the regime at low capillary number where the droplet remains circular in the horizontal plane. Parametric studies were performed by varying the droplet horizontal radius and the capillary number. For droplets with an horizontal radius $R$ larger than the half-height of the channel $H$, the droplet overfills the channel and a thin lubrication film is formed between the droplet and the wall. An adaptive two-phase flow solver is utilized for the simulations and mesh is locally refined to accurately resolve the thin film. The ratio $R/H$ has a significant impact on the droplet velocity and controls three regimes: the Poiseuille-dominated regime, the film-dominated regime, and the transition regime. The droplet velocity decreases dramatically with $R/H$ in the transition regime. In the film-dominated regime, the droplet horizontal radius is much larger than the channel height and the droplet holds a pancake-like shape. Then the migration velocity is found to become independent of the radius and is mainly dictated by the capillary number. The simulation also shows that the droplet velocity in the film-dominated regime decreases with the capillary number and is generally lower than the average inflow velocity. The lubrication film dynamics and the three-dimensional flow structure seem to both contribute to the low mobility of the droplet. In this low-migration-velocity scenario the interfacial flow on the top of the droplet moves from front toward the rear in the droplet reference frame, but reverses its direction moving toward the front from the two sides of the the droplet. As the depth-averaged Hele–Shaw equations ignore the effect of the lubrication film on the droplet dynamics, their solution fails to capture the dependence of droplet velocity on capillary number. As the shear induced by the lubrication film has a strong influence on the droplet dynamics, the viscosity ratio can be an important parameter for the present problem, and we intend to investigate its influence in future work.

Acknowledgements {#acknowledgements .unnumbered}
================

This project has been supported by the the ANR TRAM project (ANR-13-BS09-0011). The simulations of this paper are conducted on our laboratory cluster and on the CINES Occigen machine for which we gratefully acknowledge grant x20152b7325 from GENCI.

[32]{}ifxundefined \[1\][ ifx[\#1]{} ]{}ifnum \[1\][ \#1firstoftwo secondoftwo ]{}ifx \[1\][ \#1firstoftwo secondoftwo ]{}““\#1””@noop \[0\][secondoftwo]{}sanitize@url \[0\][‘\
12‘\$12 ‘&12‘\#12‘12‘\_12‘%12]{}@startlink\[1\]@endlink\[0\]@bib@innerbibempty @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****, ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [“,” ]{} (),  @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [**]{} (, ) @noop [****,  ()]{} @noop [****, ()]{} @noop [****, ()]{} @noop [****,  ()]{} @noop [****,  ()]{}

[^1]: Symbol with $\hat{}$ denotes dimensional parameters.
