1 Introduction
The motion of multiple droplets contained within a moving immersing fluid in a Hele-Shaw geometry is relevant to various microfluidic systems, such as chemical and biological lab-on-a-chip devices (Stone, Stroock & Ajdari Reference Stone, Stroock and Ajdari2004; Squires & Quake Reference Squires and Quake2005; Zhao Reference Zhao2013). Generation of such configurations is commonly done by controlled injection, via T-junction, of one liquid into a second immiscible liquid (Garstecki & Whitesides Reference Garstecki and Whitesides2006; Christopher et al. Reference Christopher, Noharuddin, Taylor and Anna2008; Desreumaux et al. Reference Desreumaux, Caussin, Jeanneret, Lauga and Bartolo2013, among many others).
Since inertial effects and gravity are commonly negligible due to the small size associated with microfluidic configurations, the flow field is governed by surface related effects such as capillarity and viscosity. Modelling the dynamics of such systems requires calculation of the forces of interaction between multiple droplets moving relatively to each other and the surrounding fluid. These forces are commonly approximated by superposition of dipole solutions (Beatus, Tlusty & Bar-Ziv Reference Beatus, Tlusty and Bar-Ziv2006; Uspal & Doyle Reference Uspal and Doyle2012; Shani et al. Reference Shani, Beatus, Bar-Ziv and Tlusty2013; Fleury et al. Reference Fleury, Schiller, Thutupalli, Gompper and Seemann2014), which requires an assumption of a sufficiently large distance between the droplets. Given the lack of exact models of interaction forces between closely spaced droplets, models based on the dipole approximation are frequently applied for cases which violate the fundamental assumption of long distances between the droplets (e.g. Liu, Goree & Feng Reference Liu, Goree and Feng2012; Shen et al. Reference Shen, Leman, Reyssat and Tabeling2014). In addition, some systems involve droplet within a droplet configurations (He et al. Reference He, Edgar, Jeffries, Lorenz, Shelby and Chiu2005; Hindmarsh et al. Reference Hindmarsh, Su, Flanagan and Singh2005; Hanson et al. Reference Hanson, Chang, Graves, Li, Mason and Deming2008), which, to the best of our knowledge, have not been previously modelled in a Hele-Shaw geometry.
Exact solution of the interaction forces between closely spaced droplets in a Hele-Shaw cell requires applying no-penetration boundary conditions for circular droplets moving relative to one another, which due to the interaction experience a non-uniform external flow. For two droplets of identical radii Beatus, Bar-Ziv & Tlusty (Reference Beatus, Bar-Ziv and Tlusty2012) obtained such an exact solution by application of the method of images. The aim of this work is to apply transformation to a bipolar coordinate system to obtain an exact solution for the general case of two relatively moving droplets of different radii and viscosities, as well as a droplet within a droplet configuration, located within a moving immersing fluid in a Hele-Shaw cell. This paper is organized as follows. In § 2 we describe the microfluidic droplet configuration and relevant assumptions made in the model. In § 3 we obtain the pressure field created by two closely spaced, droplets, in the bipolar coordinate system, moving relative to one another. In § 4 we calculate the interaction forces between two droplets in the Cartesian coordinate system. We compare our model to experimental results, the existing dipole model as well as numerical computations. We then apply our results to examine the dynamics of two closely spaced droplets, and a droplet within a droplet, subject to an external oscillating velocity. In § 5 we give concluding remarks.
2 Problem formulation
We focus on the interaction forces between two closely spaced droplets (see figure 1 a), or a droplet within a droplet (see figure 1 b), positioned between two parallel flat plates separated by a small gap. The droplets are immersed within a different immiscible liquid flowing uniformly far from the droplets.
The Cartesian coordinate system is denoted by
$(x,y,z)$
and time is denoted by
$t$
. The coordinate
$x$
is defined as parallel to a line connecting the centres of the droplets. The upper and lower plates are parallel to the
$x{-}y$
plane. The gap between the parallel plates is denoted by
$g$
and is assumed to be small compared to all other length scales of the problem. The droplets are assumed circular in the
$x{-}y$
plane due to dominant surface tension effects. The centres of droplets
$1$
and
$2$
are denoted by
$(x_{1},y_{1})$
and
$(x_{2},y_{2})$
, respectively. The distance between the centres of the droplets is denoted by
$h$
. The radii of droplets
$1$
and
$2$
are denoted by
$r_{1}$
and
$r_{2}$
, respectively. The
$z$
-averaged liquid velocity in the
$x{-}y$
plane is denoted by
$(u,v)$
and the uniform
$z$
-averaged velocity of the surrounding liquid far from the location of the droplets is denoted by
$(u_{\infty },v_{\infty })$
. The velocities of droplets
$1$
and
$2$
are
$(u_{1},v_{1})$
and
$(u_{2},v_{2})$
, respectively. The viscosity of the surrounding liquid is
${\it\mu}_{s}$
and the viscosities of droplets
$1$
and
$2$
are denoted by
${\it\mu}_{1}$
and
${\it\mu}_{2}$
, respectively. We focus on configurations with negligible inertial effects.
Hereafter, normalized variables will be denoted by a capital letter and characteristic values will be denoted by an asterisk superscript. We define the normalized coordinates
$(X,Y,Z)$
and time
$T$
,
$$\begin{eqnarray}X=\frac{x}{l^{\ast }},\quad Y=\frac{y}{l^{\ast }},\quad Z=\frac{z}{g},\quad T=\frac{t}{l^{\ast }/u^{\ast }},\end{eqnarray}$$
where
$l^{\ast }$
and
$u^{\ast }$
are characteristic length scale and speed in the
$X{-}Y$
plane. The normalized
$Z$
-averaged velocity in the
$X{-}Y$
plane and normalized pressure
$P$
are,
$$\begin{eqnarray}U=\frac{u}{u^{\ast }},\quad V=\frac{v}{u^{\ast }},\quad P=\frac{p}{12{\it\mu}_{i}u^{\ast }l^{\ast }/g^{2}},\end{eqnarray}$$
where
${\it\mu}_{i}$
is the viscosity of the liquid in the region of calculation (see figure 2
b,c). The normalized droplet radii and distance between the droplets centres are
$(R_{1},R_{2})=(r_{1}/l^{\ast },r_{2}/l^{\ast })$
and
$H=h/l^{\ast }$
, respectively. The normalized velocities of the droplets and the velocity far from the droplets are
$(U_{1},V_{1})=(u_{1}/u^{\ast },v_{1}/u^{\ast })$
,
$(U_{2},V_{2})=(u_{2}/u^{\ast },v_{2}/u^{\ast })$
and
$(U_{\infty },V_{\infty })=(u_{\infty }/u^{\ast },v_{\infty }/u^{\ast })$
, respectively.
Figure 1. Schematic description of the droplet configuration and coordinate system. The coordinate
$x$
is defined as parallel to a line connecting the centres of the droplets. Panel (a) presents two closely spaced droplets and panel (b) presents a droplet within a larger droplet.
$r_{1},r_{2}$
are the radii of droplets
$1$
and
$2$
, respectively.
$h$
is the distance between the centres of the droplets,
$g$
is the gap between the upper and lower plates.
Using the assumptions of negligible inertial effects and a small gap between the plates, we apply the Hele-Shaw approximation, and thus the governing equation is the Laplace equation for the liquid pressure
${\rm\nabla}_{\Vert }^{2}P=0$
, where
$\boldsymbol{{\rm\nabla}}_{\Vert }$
is the two-dimensional nabla operator in the
$X{-}Y$
plane and
$\partial P/\partial Z=0$
. The relation between the
$Z$
-averaged normalized velocity,
$(U,V)$
, and the normalized pressure gradient is
$(U,V)=-\boldsymbol{{\rm\nabla}}_{\Vert }P$
. The Laplace equation is supplemented by the no-penetration condition at the boundary of the droplets,
$$\begin{eqnarray}(\boldsymbol{{\rm\nabla}}_{\Vert }P+(U_{i},V_{i}))\boldsymbol{\cdot }\hat{\boldsymbol{n}}=0,\end{eqnarray}$$
where
$\hat{\boldsymbol{n}}$
is a unit vector normal to the boundary of droplet
$i$
and pointing outward,
$i=1,2$
correspond to droplets
$1,2$
respectively, and the uniform velocity far from the droplets,
$(U,V)\rightarrow (U_{\infty },V_{\infty })$
.
3 Calculation of pressure distribution in a bipolar coordinate system
We define auxiliary coordinates, denoted by tildes, as
$(\tilde{X},{\tilde{Y}})$
$$\begin{eqnarray}\tilde{X}=X+U_{\infty }T,\quad {\tilde{Y}}=Y+V_{\infty }T,\end{eqnarray}$$
where the velocity of droplet
$i$
is
$(\tilde{U} _{i},\tilde{V}_{i})=(U_{i}-U_{\infty },V_{i}-V_{\infty })$
and
$(\tilde{U} _{\infty },\tilde{V}_{\infty })=(0,0)$
. We apply the bipolar coordinate transformation from the
$\tilde{X}{-}{\tilde{Y}}$
plane to the
${\it\sigma}{-}{\it\tau}$
plane as described in figure 2(a), where the curves of constant
${\it\sigma}$
(blue solid) and of constant
${\it\tau}$
(red dashed) are perpendicular and describe circles. The
$({\it\sigma},{\it\tau})$
coordinates have two focal points, denoted by
$F_{1}$
and
$F_{2}$
, located at
$(A,0)$
and
$(-A,0)$
, respectively, on the
$\tilde{X}$
-axis of the Cartesian coordinate system. The conformal transformation from the Cartesian coordinate system
$(\tilde{X},{\tilde{Y}})$
to the bipolar coordinates
$({\it\sigma},{\it\tau})$
is (Polyanin Reference Polyanin2001)
$$\begin{eqnarray}(\tilde{X},{\tilde{Y}})=A\left(\frac{\sinh {\it\tau}}{\cosh {\it\tau}-\cos {\it\sigma}},\frac{\sin {\it\sigma}}{\cosh {\it\tau}-\cos {\it\sigma}}\right),\end{eqnarray}$$
where
${\it\sigma}$
of the point
$L$
is the angle
$\angle F_{2}LF_{1}$
and
${\it\tau}$
is the natural logarithm of the ratio of distances
$|LF_{1}|$
and
$|LF_{2}|$
respectively,
${\it\tau}=\ln |LF_{1}|/|LF_{2}|$
. Equation (3.2) can be presented as
$$\begin{eqnarray}(\tilde{X}-A\cdot \coth {\it\tau})^{2}+{\tilde{Y}}^{2}=\frac{A^{2}}{\sinh ^{2}{\it\tau}},\end{eqnarray}$$
and thus the radius of droplet
$i$
is given by
$R_{i}=A/\sinh {\it\tau}$
and the droplet radius centre is located at
$A\coth {\it\tau}$
. From (3.2)
$R_{1}$
,
$R_{2}$
and
$H$
are related to
${\it\tau}_{1}$
,
${\it\tau}_{2}$
and
$A$
(see figure 2
a),
$$\begin{eqnarray}{\it\tau}_{1}=\sinh ^{-1}\left(\frac{A}{R_{1}}\right),\quad {\it\tau}_{2}=C\sinh ^{-1}\left(\frac{A}{R_{2}}\right)\end{eqnarray}$$
and
$$\begin{eqnarray}A=\frac{~\sqrt{(H^{2}-(R_{1}+R_{2})^{2})(H^{2}-(R_{1}-R_{2})^{2})}}{2H},\end{eqnarray}$$
where
$C=-1$
for two closely spaced droplets (figure 2
b) and
$C=1$
for a droplet within a droplet (figure 2
c).
Figure 2. Illustration of the bipolar coordinate system. Panel (a) presents isosurfaces of
${\it\tau}$
(red dashed), isosurfaces of
${\it\sigma}$
(blue smooth) and two focal points at
$(-A,0),(A,0)$
.
${\it\sigma}$
is the angle between the lines connecting point
$L$
and the two focal points and
${\it\tau}$
is the natural logarithm of the ratio of the lengths of the lines connecting
$L$
to the focal points. Panels (b) and (c) present two closely spaced droplets and a droplet within a droplet, respectively, in the
${\it\sigma}{-}{\it\tau}$
plane.
The governing equation in the
${\it\sigma}{-}{\it\tau}$
plane is
$$\begin{eqnarray}\frac{\partial ^{2}\tilde{P}}{\partial {\it\sigma}^{2}}+\frac{\partial ^{2}\tilde{P}}{\partial {\it\tau}^{2}}=0,\end{eqnarray}$$
in the rectangular domain
${\it\tau}_{1}\leqslant {\it\tau}\leqslant {\it\tau}_{2}$
and
$-{\rm\pi}\leqslant {\it\sigma}\leqslant {\rm\pi}$
. Since
${\it\sigma}\rightarrow {\rm\pi}$
and
${\it\sigma}\rightarrow -{\rm\pi}$
both represent the line connecting the centres of the circles, from opposing directions, equation (3.6) is supplemented by the periodic boundary conditions at
${\it\sigma}=-{\rm\pi},{\rm\pi}$
$$\begin{eqnarray}\tilde{P}({\it\sigma}=-{\rm\pi},{\it\tau})=\tilde{P}({\it\sigma}={\rm\pi},{\it\tau}),\quad \left.\frac{\partial \tilde{P}}{\partial {\it\sigma}}\right|_{{\it\sigma}=-{\rm\pi}}=\left.\frac{\partial \tilde{P}}{\partial {\it\sigma}}\right|_{{\it\sigma}={\rm\pi}},\end{eqnarray}$$
as well as no-penetration condition on the boundary of the droplets at
${\it\tau}_{i}$
(
$i=1,2$
corresponds to droplets
$1,2$
respectively).
$$\begin{eqnarray}\left.\frac{\partial P}{\partial {\it\tau}}\right|_{{\it\tau}={\it\tau}_{i}}+\tilde{U} _{i}A\frac{1-\cosh {\it\tau}_{i}\cos {\it\sigma}}{(\cosh {\it\tau}_{i}-\cos {\it\sigma})^{2}}-\tilde{V}_{i}A\frac{\sin {\it\sigma}\sinh {\it\tau}_{i}}{(\cosh {\it\tau}_{i}-\cos {\it\sigma})^{2}}=0.\end{eqnarray}$$
Applying the separation of variables and solving the eigenvalue problem yields the general form of the solution
$$\begin{eqnarray}\tilde{P}({\it\sigma},{\it\tau})=\mathop{\sum }_{n=1}^{\infty }({\it\alpha}_{n}\cos n{\it\sigma}+{\it\beta}_{n}\sin n{\it\sigma})({\it\gamma}_{n}\cosh n{\it\tau}+{\it\delta}_{n}\sinh n{\it\tau}).\end{eqnarray}$$
Utilizing (3.8) on
${\it\tau}={\it\tau}_{i}$
yields,
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\it\alpha}_{n}{\it\delta}_{n}=\frac{-I_{1}\sinh (n{\it\tau}_{2})+I_{2}\sinh (n{\it\tau}_{1})}{n\sinh (n({\it\tau}_{1}-{\it\tau}_{2}))},\quad {\it\alpha}_{n}{\it\gamma}_{n}=\frac{I_{1}\cosh (n{\it\tau}_{2})-I_{2}\cosh (n{\it\tau}_{1})}{n\sinh (n({\it\tau}_{1}-{\it\tau}_{2}))},\\ \displaystyle {\it\beta}_{n}{\it\delta}_{n}=\frac{-J_{1}\sinh (n{\it\tau}_{2})+J_{2}\sinh (n{\it\tau}_{1})}{n\sinh (n({\it\tau}_{1}-{\it\tau}_{2}))},\quad {\it\beta}_{n}{\it\gamma}_{n}=\frac{J_{1}\cosh (n{\it\tau}_{2})-J_{2}\cosh (n{\it\tau}_{1})}{n\sinh (n({\it\tau}_{1}-{\it\tau}_{2}))},\end{array}\right\}\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle & \displaystyle I_{i}=\frac{\tilde{U} _{i}A}{{\rm\pi}}\int _{-{\rm\pi}}^{{\rm\pi}}\cos n{\it\sigma}\frac{\cosh {\it\tau}_{i}\cos {\it\sigma}-1}{(\cosh {\it\tau}_{i}-\cos {\it\sigma})^{2}}\,\text{d}{\it\sigma}=2\tilde{U} _{i}An\text{e}^{-n|{\it\tau}_{i}|}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle J_{i}=\frac{\tilde{V}_{i}A}{{\rm\pi}}\int _{-{\rm\pi}}^{{\rm\pi}}\sin n{\it\sigma}\frac{\sinh {\it\tau}_{i}\sin {\it\sigma}}{(\cosh {\it\tau}_{i}-\cos {\it\sigma})^{2}}\,\text{d}{\it\sigma}=\left\{\begin{array}{@{}ll@{}}2\tilde{V}_{i}An\text{e}^{-n{\it\tau}_{i}}\quad & \text{if }{\it\tau}_{i}>0\\ -2\tilde{V}_{i}An\text{e}^{n{\it\tau}_{i}}\quad & \text{if }{\it\tau}_{i}<0.\end{array}\right. & \displaystyle\end{eqnarray}$$
$(X,Y)$
coordinates is thus obtained from (3.9) to (3.11) together with (3.2) and the transformation 
$$\begin{eqnarray}P(X,Y)=\tilde{P}(\tilde{X},{\tilde{Y}})+U_{\infty }X+V_{\infty }Y.\end{eqnarray}$$
4 Interaction forces and dynamics of two closely spaced droplets and a droplet within a droplet
In order to calculate the force acting on droplet
$i$
in the Cartesian coordinates
$(X,Y)$
we integrate
$$\begin{eqnarray}(F_{X,i},F_{Y,i})=-\oint _{droplet}P(X,Y)\left(\frac{\partial Y}{\partial S},\frac{\partial X}{\partial S}\right)\text{d}S,\end{eqnarray}$$
where
$S$
is a coordinate across the droplet boundary. In the
${\it\sigma}{-}{\it\tau}$
plane,
$S$
corresponds to the variation of
${\it\sigma}$
for the fixed value of
${\it\tau}$
. Thus the Jacobian simply reads
$\partial S/\partial {\it\sigma}$
. The dependence of the Jacobian in
$S$
eliminates its dependence in the reflection of the pressure in both
$X,Y$
directions leading to the simplified expression
$$\begin{eqnarray}(F_{X,i},F_{Y,i})=-\int _{-{\rm\pi}}^{{\rm\pi}}\tilde{P}({\it\sigma},{\it\tau}_{1})\left(\frac{\partial Y}{\partial {\it\sigma}},\frac{\partial X}{\partial {\it\sigma}}\right)\text{d}{\it\sigma}+{\rm\pi}R_{i}^{2}(U_{\infty },V_{\infty }),\end{eqnarray}$$
in the
${\it\sigma}{-}{\it\tau}$
plane. Substituting (3.9)–(3.12) and (3.2) into (4.2) yields the force that acts on droplet
$i$
due to the uniform flow and as well as the flow field interaction with the neighbouring droplet
$j$
$$\begin{eqnarray}\displaystyle (F_{X,i},F_{Y,i}) & = & \displaystyle -4{\rm\pi}A^{2}\mathop{\sum }_{n=1}^{\infty }n\text{e}^{-2n|{\it\tau}_{i}|}\coth (n(|{\it\tau}_{i}|+|{\it\tau}_{j}|))\cdot (U_{i}-U_{\infty },V_{i}-V_{\infty })\nonumber\\ \displaystyle & & \displaystyle +\,4{\rm\pi}A^{2}\mathop{\sum }_{n=1}^{\infty }n(\coth (n(|{\it\tau}_{i}|+|{\it\tau}_{j}|))-1)\cdot (U_{j}-U_{\infty },-V_{j}+V_{\infty })\nonumber\\ \displaystyle & & \displaystyle +\,{\rm\pi}R_{i}^{2}(U_{\infty },V_{\infty }).\end{eqnarray}$$
The first and the second terms of the expression (4.3) can be attributed to the interdroplet interaction forces exerted on the
$i\text{th}$
droplet which are induced by the instantaneous displacement (first term in (4.3)) and velocity (second term in (4.3)) of the second (
$j\text{th}$
) droplet, accordingly. The third term in the expression (4.3) is related to the force component exerted on a single droplet by the surrounding fluid (i.e. in the absence of interdroplet interactions).
Details of the derivation of (4.3) are given in appendix A. For the case of a droplet within a larger droplet we obtain
$$\begin{eqnarray}\displaystyle (F_{X,i},F_{Y,i}) & = & \displaystyle \displaystyle -4{\rm\pi}A^{2}\mathop{\sum }_{n=1}^{\infty }n\text{e}^{-2n{\it\tau}_{i}}\coth (n|{\it\tau}_{i}-{\it\tau}_{j}|)\cdot (U_{i},V_{i})\nonumber\\ \displaystyle & & \displaystyle \displaystyle -\,4{\rm\pi}A^{2}\mathop{\sum }_{n=1}^{\infty }\frac{n\text{e}^{-n({\it\tau}_{i}+{\it\tau}_{j})}}{\sinh (n|{\it\tau}_{i}-{\it\tau}_{j}|)}\cdot (U_{j},V_{j}),\end{eqnarray}$$
where
$(i,j)=(1,2)$
represents the external droplet and
$(i,j)=(2,1)$
represents the inner droplet. The forces
$(F_{X,i},F_{Y,i})$
are normalized by the characteristic force
$f^{\ast }=12{\it\mu}u^{\ast }l^{\ast 2}/g$
, where
${\it\mu}$
is defined by the liquid in the region of calculation (which is the immersing fluid for (4.3) and the viscosity of the external droplet for (4.4)). In the case of droplet within a droplet configuration the first and the second terms of the expression (4.4) can be again attributed to the interdroplet interaction forces exerted on droplet
$i$
, which are induced by the instantaneous displacement (first term in (4.4)) and velocity (second term in (4.4)) of droplet
$j$
.
Equations (4.3), (4.4) allow for exact (in the Hele-Shaw limit) solution of interaction forces for arbitrarily chosen droplets positions and velocities. However, obtaining the dynamics of the pair of droplets driven by external uniform flow (which may be closely spaced) requires computing their velocities from force balance between the forces applied by the surrounding fluid
$(f_{x,i},f_{y,i})$
and the droplets’ internal friction, denoted hereafter as
$(r_{x,i},r_{y,i})$
. Given the dominant effect of capillary forces in microconfigurations, surface tension significantly affects the friction force of the droplet. Experimental works commonly estimate the total friction of the droplet by measuring the ratio
${\it\beta}=u_{d}/u_{\infty }$
(e.g. Beatus et al.
Reference Beatus, Bar-Ziv and Tlusty2012; Liu et al.
Reference Liu, Goree and Feng2012; Shen et al.
Reference Shen, Leman, Reyssat and Tabeling2014), where
$u_{d}$
is the speed of an isolated non-wetting droplet immersed in a perfectly wetting fluid. Recently, Huerre et al. (Reference Huerre, Theodoly, Leshansky, Valignat, Cantat and Jullien2015) related
$u_{d}$
to the thickness of a thin lubrication film of the immersing fluid above and below the droplet, where for sufficiently small film thickness it was shown that
${\it\beta}\approx 0.25$
.
For simplicity, we hereafter use dimensional parameters and define
${\it\mu}_{1}$
,
${\it\mu}_{2}$
and
${\it\mu}_{s}$
are the viscosities of droplet
$1$
, droplet
$2$
and the surrounding liquid, respectively. The coefficient
${\it\beta}_{i}$
may be related to the internal friction force on droplet
$i$
, denoted as
$(r_{x,i},r_{y,i})=-12({\it\mu}_{i}{\rm\pi}r_{i}^{2}/g-c{\rm\pi}r_{i}^{2}/g)(u_{i},v_{i})$
, where
$c_{i}=2{\it\mu}_{s}/{\it\beta}_{i}-{\it\mu}_{i}-{\it\mu}_{s}$
. Applying the aforementioned force balance
$(r_{x,in},r_{y,in})=(f_{x,i},f_{y,i})$
yields the equations of motion for two closely spaced droplets
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{u_{i}}{u_{\infty }}=\frac{({\it\mu}_{s}g_{i}-{\it\mu}_{s}g_{12}+{\it\mu}_{s}r_{i}^{2})({\it\mu}_{j}r_{j}^{2}+c_{j}r_{j}^{2}+{\it\mu}_{s}g_{j})+{\it\mu}_{s}g_{12}({\it\mu}_{s}g_{j}-{\it\mu}_{s}g_{12}+{\it\mu}_{s}r_{j}^{2})}{{\it\mu}_{s}^{2}g_{12}^{2}+({\it\mu}_{1}r_{1}^{2}+c_{1}r_{1}^{2}+{\it\mu}_{s}g_{1})({\it\mu}_{2}r_{2}^{2}+c_{2}r_{2}^{2}+{\it\mu}_{s}g_{2})},\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{v_{i}}{v_{\infty }}=\frac{({\it\mu}_{s}g_{i}+{\it\mu}_{s}g_{12}+{\it\mu}_{s}r_{i}^{2})({\it\mu}_{j}r_{j}^{2}+c_{j}r_{j}^{2}+{\it\mu}_{s}g_{j})-{\it\mu}_{s}g_{12}({\it\mu}_{s}g_{j}+{\it\mu}_{s}g_{12}+{\it\mu}_{s}r_{j}^{2})}{{\it\mu}_{s}^{2}g_{12}^{2}+({\it\mu}_{1}r_{1}^{2}+c_{1}r_{1}^{2}+{\it\mu}_{s}g_{1})({\it\mu}_{2}r_{2}^{2}+c_{2}r_{2}^{2}+{\it\mu}_{s}g_{2})}.\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{u_{1}}{u_{\infty }}=\frac{v_{1}}{v_{\infty }}=\frac{2r_{1}^{2}{\it\mu}_{s}{\it\mu}_{1}(-c_{2}r_{2}^{2}-{\it\mu}_{2}r_{2}^{2}+{\it\mu}_{1}g_{2})}{{\it\mu}_{1}^{2}g_{12}^{2}+({\it\mu}_{s}r_{1}^{2}+c_{1}r_{1}^{2}-{\it\mu}_{1}g_{1})(-{\it\mu}_{2}r_{2}^{2}-c_{2}r_{2}^{2}+{\it\mu}_{1}g_{2})}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{u_{2}}{u_{\infty }}=\frac{v_{2}}{v_{\infty }}=\frac{-2r_{1}^{2}{\it\mu}_{s}{\it\mu}_{d}g_{12}}{{\it\mu}_{1}^{2}g_{12}^{2}+({\it\mu}_{s}r_{1}^{2}+c_{1}r_{1}^{2}-{\it\mu}_{1}g_{1})(-{\it\mu}_{2}r_{2}^{2}-c_{2}r_{2}^{2}+{\it\mu}_{1}g_{2})}, & \displaystyle\end{eqnarray}$$
$(u_{1},v_{1})$
represents the external droplet and
$(u_{2},v_{2})$
represents the inner droplet, 
$$\begin{eqnarray}\displaystyle g_{1}=4a^{2}\mathop{\sum }_{n=1}^{\infty }n\text{e}^{-2n{\it\tau}_{1}}\coth (n({\it\tau}_{1}-{\it\tau}_{2})),\quad g_{2}=4a^{2}\mathop{\sum }_{n=1}^{\infty }n\text{e}^{-2n|{\it\tau}_{2}|}\coth (n({\it\tau}_{1}-{\it\tau}_{2})),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle g_{12}=4a^{2}\mathop{\sum }_{n=1}^{\infty }\frac{-n\text{e}^{-2n({\it\tau}_{1}+|{\it\tau}_{2}|)}}{\sinh (n({\it\tau}_{1}-{\it\tau}_{2}))},\end{eqnarray}$$
and
$$\begin{eqnarray}a=\frac{\sqrt{(h^{2}-(r_{1}+r_{2})^{2})(h^{2}-(r_{1}-r_{2})^{2})}}{2h}.\end{eqnarray}$$
We note that (4.5)–(4.7) are the droplet velocities for a coordinate system where
$X$
is parallel to the line connecting the centres of the droplets. Thus, the coordinate system rotates with the motion of the droplets and hence
$(u_{\infty },v_{\infty })$
needs to be recalculated accordingly.
Figure 3. Comparison between (4.3) and (4.4) (smooth line), the dipole approximation (dashed line) and numerical calculations (dotted line). Panels (a,c,e) present force in the
$X$
direction versus the distance between droplets. Panels (b,d,f) present force in the
$Y$
direction versus the distance between droplets. Panels (a–d) present the interaction forces of two adjacent droplets (4.3). In panels (a,b)
$R_{1}=1$
,
$R_{2}=4$
,
$(U_{1},V_{1})=(1.92,-1.61)$
,
$(U_{2},V_{2})=(0,-1.5)$
,
$(U_{\infty },V_{\infty })=(2.5,0)$
. In panels (c,d)
$R_{1}=1$
,
$R_{2}=3$
,
$(U_{1},V_{1})=(0.52,1.93)$
,
$(U_{2},V_{2})=(1,0)$
and
$(U_{\infty },V_{\infty })=(1,0)$
. Panels (e,f) present the forces of a droplet within a droplet (4.4) where
$R_{1}=3$
,
$R_{2}=1$
,
$(U_{1},V_{1})=(1.414,1.414)$
and
$(U_{2},V_{2})=(-0.707,0.707)$
.
Figure 3 illustrates the forces acting on two closely spaced droplets given in (4.3) (smooth lines, panels a–d), as well as the force acting on a droplet within a droplet configuration given in (4.4) (smooth lines, panels e–f), where the positions and velocities of the droplets are prescribed. Both analytical solutions (4.3), (4.4) are compared to the forces computed from the numerical solution of the Laplace equation in the original (Cartesian) coordinates satisfying the boundary conditions (2.2). Here we note that the numerical solution of the Laplace equation has been performed in the commercially available software (COMSOL Multiphysics 4.5™) using a grid of
$10^{5}$
mesh points. Numerical computation of the force components (i.e.
$(F_{X,i},F_{Y,i})$
) applied on the droplets (for both configurations) is performed by the numerical integration of the pressure (calculated from the solution of the Laplace equation) along the boundary of each droplet projected into the horizontal and vertical directions. In the case of closely spaced droplets, the dipole approximation (dashed lines) is compared to the exact analytical solution (4.3) also. The forces in the
$x$
-direction (figure 3
a,c,e) and the
$y$
-direction (figure 3
b,d,f) are presented versus the distance between the droplets
$H{-}R_{1}{-}R_{2}$
in figure 3(a–d) and versus
$H$
in figure 3(e,f). In all cases, the analytical solutions given in (4.3), (4.4) are in the perfect agreement with the numerical solution of (3.6)–(3.8). In the case of
$H{-}R_{1}{-}R_{2}\ll 1$
, a significant difference exists between the dipole model and the exact solution (4.3). However, this difference vanishes for
$H{-}R_{1}{-}R_{2}\gg 1$
as
$H\rightarrow \infty$
and the force approaches the solution of a single isolated droplet,
${\rm\pi}R_{1}^{2}(2U_{\infty }-U_{1})$
. In figure 3(c,d) we illustrate the case of one droplet moving with the same velocity as the external flow. In this case
$H{-}R_{1}{-}R_{2}$
does not affect the interaction force computed by the dipole model, since the dipole is of zero strength. In comparison, the exact solution yields a significant increase in the interaction forces and thus the dipole model is this case fails both quantitatively and qualitatively. Figure 3(e,f) examines the forces acting on a droplet within a droplet, yielding approximately constant forces for
$H<1$
. Additionally, the interaction force increases with
$H$
as the small droplet approaches the boundary of the bigger droplet.
Figure 4. Accuracy of the dipole approximation compared with the exact solution (4.3). Panels (a,c) present error of the dipole approximation scaled by the total force,
$E_{T}=|F_{exact}-F_{dipole}|/|F_{d\infty }|$
. Panels (b,d) present the error for the dipole approximation scaled by the interaction force,
$E_{R}=|F_{exact}-F_{dipole}|/|F_{interaction}|$
. In panels (a,b) the droplet pair is parallel to the external flow. In panels (c,d) the droplet pair is perpendicular to the external flow.
$R_{1}/R_{2}=1$
(smooth line),
$R_{1}/R_{2}=2$
(dashed dotted line),
$R_{1}/R_{2}=4$
(dashed line), in all cases
$U_{\infty }=5$
.
Figure 5. Scaled speed of droplet pairs of the identical particles
$U/U_{d}$
versus the distance between the droplets
$H{-}R_{1}{-}R_{2}$
.
$U_{d}$
is the speed of a single isolated droplet. The line connecting the centres of the droplets is: (I) parallel to external flow (squares, dotted line and dotted–dashed line) or (II) perpendicular to the external flow (circles, smooth line and dashed line). The experimental data were taken from Shen et al. (Reference Shen, Leman, Reyssat and Tabeling2014) and are marked by squares and circles. Exact solution (4.3) is marked by smooth and dashed-dotted lines. Dipole approximation is marked by dashed and dotted lines.
Figure 6. Trajectories of the centres of two droplets due to an external rotating uniform velocity
$(U_{\infty },V_{\infty })=|U_{\infty }|(\cos ({\it\omega}T),\sin ({\it\omega}T))$
. Smooth lines are exact solution (4.6)–(4.7) and dashed lines mark dipole approximation. Panels (a–c),
$|U_{\infty }|=0.4$
,
${\it\omega}=0.025$
. Panel (a)
$R_{1}=R_{2}=1$
,
${\it\mu}_{1}/{\it\mu}_{2}=1$
. Panel (b)
$R_{1}=10$
,
$R_{2}=1$
and
${\it\mu}_{1}/{\it\mu}_{2}=1$
. Panel (c)
$R_{1}=R_{2}=1$
,
${\it\mu}_{1}/{\it\mu}_{2}=1000$
. Panel (d) presents a droplet within a droplet where
$|U_{\infty }|=0.5$
,
${\it\omega}=1.256$
,
$R_{1}=3$
,
$R_{2}=1$
and
${\it\mu}_{1}/{\it\mu}_{2}=0.5$
. Panels (e,f) present zoomed regions of the trajectories illustrated in panel (d).
Figure 4 illustrates the difference between (4.3) and the dipole approximation for the limit of
$H{-}R_{1}{-}R_{2}\rightarrow 0$
. In all cases
$U_{\infty }=5$
, and
$R_{1}/R_{2}=1,2,4$
, corresponding to smooth, dashed and dashed-dotted lines, respectively. The line connecting the centres of the droplets is parallel (panels a,b) or perpendicular (panels c,d) to the uniform flow. Since (4.3) is an exact solution (under the assumptions of the Hele-Shaw model), we define the error of the dipole approximation by
$E_{T}=|F_{exact}-F_{dipole}|/|F_{d\infty }|$
(panels a,c), where
$F_{exact}$
is the solution of (4.3),
$F_{dipole}$
is the force computed by the dipole approximation and
$F_{d\infty }$
is the force on an isolated droplet. For parallel configuration (panel a) clear minima are evident for
$U_{1}/U_{2}\approx 1$
where the error decreases to
${\approx}5\,\%$
. For the perpendicular configuration no minima are evident and the errors vary from
$25\,\%$
to
$40\,\%$
. In fact, various theoretical and experimental models concerned with the dynamics of microfluidic droplets (e.g. excitation of the phonon modes in a droplet lattice Beatus et al.
Reference Beatus, Tlusty and Bar-Ziv2006) require a knowledge of the interdroplet interaction forces. These interaction forces are represented by the difference between the total force and the force acting on an isolated droplet. In the present study we define the error of the interaction force as
$E_{R}=|F_{exact}-F_{dipole}|/|F_{interaction}|$
presented in (panels b,d), where
$F_{interaction}=F_{exact}-F_{d\infty }$
. The errors associated with the interaction forces are
${\approx}100\,\%$
for the parallel configuration (panel b) and
${\approx}60\,\%$
for the perpendicular configuration (panel d). As can be inferred from the results of figure 4, the dipole approximation of the total force is more accurate for the parallel configuration (panel a) in comparison to the perpendicular one (panel c). However, in the case of the interaction forces the opposite is true, namely the dipole approximation yields the more accurate approximation for the perpendicular configuration (panel d) compared to the parallel one (panel b). Thus, the minima in the error of the dipole model for the parallel configuration is due to a reduction of magnitude of the interaction forces for
$U_{1}/U_{2}\approx 1$
.
In figure 5 we compare the experimental results (squares, circles) presented in Shen et al. (Reference Shen, Leman, Reyssat and Tabeling2014) with the two analytical models, namely the classical dipole model (dotted line, dashed line) and the exact model (4.3) (smooth line, dash-dotted line) developed herein. In this figure the speed of the droplet pairs are plotted versus the distance between the droplets for the two distinct orientations: (I) the line connecting the centres of the droplets is parallel to external flow (squares, dotted line and dotted–dashed line) and (II) perpendicular to the external flow (circles, smooth line and dashed line). In case (I) the configuration is close to the region of minimal difference between the dipole model and the exact solution (4.3) (see figure 4). In this case the difference between the models is negligible compared to the experimental resolution. However, in case (II) (4.3) accurately captures the change in slope of the experimental data as
$H{-}R_{1}{-}R_{2}\rightarrow 0$
, and thus significantly outperforms the original dipole model.
Figure 6 presents the trajectories of the centres of two droplets subject to the uniform, rotating velocity field
$(U_{\infty },V_{\infty })=|U_{\infty }|(\cos ({\it\omega}T),\sin ({\it\omega}T))$
. Smooth lines denote calculations using the exact model (4.6), (4.7) and dashed lines correspond to the original dipole model. Figure 6(a) illustrates the motion of two identical droplets (
$R_{1}/R_{2}=1,{\it\mu}_{1}/{\it\mu}_{2}=1,|U_{\infty }|=0.4,{\it\omega}=0.025$
). By careful observation of the trajectories of both droplets shown in figure 6(a) one can clearly see that, due to the equal interdroplet interaction forces, the relative displacement between their centres remains permanently constant. However, the orbits of each droplet centre are perfectly periodic and circulate around the fixed points in the plane. It is worthwhile emphasising that the presence of the interaction forces brings about the modification of a single droplet trajectory. Figure 6(b) presents the trajectories of two droplets with different radii (
$R_{1}=1$
and
$R_{2}=0.1$
,
${\it\mu}_{1}/{\it\mu}_{2}=1,|U_{\infty }|=0.4,{\it\omega}=0.025$
). In contrast to the previously considered case (i.e. identical droplets), the mismatch in the radii of the two droplets introduces the asymmetry in the interaction forces applied to the droplets. In this case, the relative distance between the centres of the droplets is time varying and two distinct (slow and fast) rotational components can be clearly identified. The fast
${\it\omega}$
component is manifested by the local, high-frequency rotations drifting around the slowly evolving (slow component) outer orbit. In fact, the fast rotational component is induced by the external flow, whereas a slower-frequency component can be attributed to the relative displacement between the centres of the droplets. Importantly, the original dipole model yields a significant error in the approximation of the lower oscillatory mode. Trajectories shown in figure 6(c) correspond to the case where the asymmetry is induced by the mismatch in viscosities
${\it\mu}_{1}/{\it\mu}_{2}=1000$
(
$R_{1}/R_{2}=1$
, keeping all other parameters identical). Similarly to the case shown in figure 6(b), the two distinct (slow and fast) rotational components are evident in figure 6(c). Interestingly, the dipole approximation yields a slow drift component in the opposite direction to that of the exact solution. Finally, in figure 6(d–f) we considered a different droplet within a droplet configuration where a mismatch is introduced in both viscosities and radii
$R_{1}=3$
and
$R_{2}=1$
; again, two distinct frequencies are evident. To clearly illustrate the low-amplitude fast scale rotations of both the trajectories figure 6(d) we zoom in locally in 6(e) and figure 6(f). Zoomed regions are denoted with dashed rectangular frames in the inset of figure 6(d).
5 Concluding remarks
We applied a transformation to a bipolar coordinate system and obtained an exact solution for the interaction forces between two closely spaced droplets, as well as a droplet within a droplet, in a Hele-Shaw cell. The commonly used approximation by dipoles is shown to give significant errors (10 %–30 %) of the total force in the limit of small distance between the droplets. An exception is the case of two droplets moving at similar speeds, where the line connecting the centres of the droplets is parallel to the external uniform flow, yielding a small error of
${\approx}5\,\%$
or less. When we examine only the interaction forces, the errors are above
${\approx}50\,\%$
and the dipole-based model cannot be used to accurately capture the interaction forces.
The analysis is performed in the framework of the Hele-Shaw limit, assuming negligible inertial effects and shallow configurations
$g/l\rightarrow 0$
. For cases in which the viscosities of the droplets are negligible compared with the viscosity of the computed region we expect errors of
$(g/l)^{2}$
due to the neglected viscous terms. For cases in which the viscosities of the droplets are similar, or greater than the viscosity of the fluid in the computed region, the errors will increase to
$g/l$
due to a mismatch of the boundary conditions at the interface between the droplets and the surrounding fluid (see Gat, Frankel & Weihs Reference Gat, Frankel and Weihs2009).
Let us close the current discussion by identifying several interesting future extensions of the current work. Due to the dominant effect of surface tension in microscale flows, small changes in the geometry of the droplet may induce significant capillary forces, while keeping the droplet geometry in the
$X{-}Y$
plane approximately circular. Thus, the classic problem of multiple solutions for the motion of droplets and fingers in a Hele-Shaw cell (Kopf-Sill & Homsy Reference Kopf-Sill and Homsy1988) is still relevant to cases with dominant surface tension. Current experimental works resolve this by empirically estimating the speed of an isolated droplet. In addition, we believe that the results presented in this work provide a rather important platform for the future extensions of the analytical modelling of complex, long-range interactions emerging in dense one-dimensional and two-dimensional droplet media. Moreover, the newly proposed unit cell model of droplet within a droplet paves the way for the future theoretical and experimental investigations of the special acoustic properties of higher-dimensional droplet structures affected by the response of the inner droplets.
Appendix A. Calculation of interaction forces
In order to find the force acting on droplet
$i$
in the Cartesian coordinates
$(X,Y)$
we calculate (4.2).
$$\begin{eqnarray}(F_{X,i},F_{Y,i})=-\int _{-{\rm\pi}}^{{\rm\pi}}\tilde{P}({\it\sigma},{\it\tau}_{i})\left(\frac{\partial Y}{\partial {\it\sigma}},\frac{\partial X}{\partial {\it\sigma}}\right)\text{d}{\it\sigma}+{\rm\pi}R_{i}^{2}(U_{\infty },V_{\infty }).\end{eqnarray}$$
Substituting (3.9)–(3.11) into (A 1) yields,
$$\begin{eqnarray}\displaystyle (F_{X,i},F_{Y,i}) & = & \displaystyle {\rm\pi}R_{i}^{2}(U_{\infty },V_{\infty })-\mathop{\sum }_{n=1}^{\infty }\int _{-{\rm\pi}}^{{\rm\pi}}({\it\gamma}_{n}\cosh n{\it\tau}_{i}+{\it\delta}_{n}\sinh n{\it\tau}_{i})({\it\alpha}_{n}\cos n{\it\sigma}+{\it\beta}_{n}\sin n{\it\sigma})\nonumber\\ \displaystyle & & \displaystyle \times \,\left(\frac{\partial Y}{\partial {\it\sigma}},\frac{\partial X}{\partial {\it\sigma}}\right)\text{d}{\it\sigma},\end{eqnarray}$$
where
$$\begin{eqnarray}\frac{\partial Y}{\partial {\it\sigma}}=\frac{a(\cos ({\it\sigma})\cosh ({\it\tau}_{i})-1)}{(\cosh ({\it\tau}_{i})-\cos ({\it\sigma}))^{2}},\quad \frac{\partial X}{\partial {\it\sigma}}=\frac{a(\sin ({\it\sigma})\sinh ({\it\tau}_{i}))}{(\cosh ({\it\tau}_{i})-\cos ({\it\sigma}))^{2}}.\end{eqnarray}$$
The structure of (A 2) is simplified further by computing the following integrals,
$$\begin{eqnarray}\displaystyle & \displaystyle \int _{-{\rm\pi}}^{{\rm\pi}}\sin (n{\it\sigma})\frac{(\cos ({\it\sigma})\cosh ({\it\tau}_{i})-1)}{(\cosh ({\it\tau}_{i})-\cos ({\it\sigma}))^{2}}\,\text{d}{\it\sigma}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \int _{-{\rm\pi}}^{{\rm\pi}}\cos (n{\it\sigma})\frac{(\cos ({\it\sigma})\cosh ({\it\tau}_{i})-1)}{(\cosh ({\it\tau}_{i})-\cos ({\it\sigma}))^{2}}\,\text{d}{\it\sigma}=2{\rm\pi}n\text{e}^{(-n|{\it\tau}_{i}|)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \int _{-{\rm\pi}}^{{\rm\pi}}\cos (n{\it\sigma})\frac{(\sinh ({\it\tau}_{i})\sin ({\it\sigma}))}{(\cosh ({\it\tau}_{i})-\cos ({\it\sigma}))^{2}}\,\text{d}{\it\sigma}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \int _{-{\rm\pi}}^{{\rm\pi}}\sin (n{\it\sigma})\frac{(\sinh ({\it\tau}_{i})\sin ({\it\sigma}))}{(\cosh ({\it\tau}_{i})-\cos ({\it\sigma}))^{2}}\,\text{d}{\it\sigma}=\left\{\begin{array}{@{}ll@{}}n\text{e}^{-n{\it\tau}_{i}}\quad & \text{if }{\it\tau}_{i}>0\\ -n\text{e}^{n{\it\tau}_{i}}\quad & \text{if }{\it\tau}_{i}<0.\end{array}\right. & \displaystyle\end{eqnarray}$$
Substituting (A 4)–(A 7) into (A 3) and performing some rather simple algebraic manipulations yields the required expression (4.3).
