2.1 Fluid dynamics

The physical system considered is a droplet with a viscous interface, centred in and surrounded by an abruptly started infinite Newtonian flow. In experiments on droplets and emulsions (Mun & McClements Reference Mun and McClements2006; Feigl et al. Reference Feigl, Megias-Alguacil, Fischer and Windhab2007), the characteristic length is droplet radius $a\sim 10{-}100~{\rm\mu}\text{m}$ and the characteristic velocity is $u_{0}\sim 100~{\rm\mu}\text{m}~\text{s}^{-1}$ . In combination with density ${\it\rho}\sim 10^{3}~\text{kg}~\text{m}^{-3}$ and dynamic viscosity ${\rm\mu}\sim 100~\text{mPa}~\text{s}$ , an approximate Reynolds number is not larger than $Re\sim 10^{-4}$ for previous experiments and reach $Re\sim 10^{-2}$ considering water. As a result, the dynamics of the fluid inside and outside of the drop are appropriately described by the Stokes equations

(2.1a,b ) $$\begin{eqnarray}-\boldsymbol{{\rm\nabla}}p+{\it\mu}{\rm\Delta}\boldsymbol{v}+\boldsymbol{f}=0,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}=0\end{eqnarray}$$

for pressure $p$ and applied force $\boldsymbol{f}$ . In the following the inner and outer viscosities will be noted as ${\it\mu}^{int}$ and ${\it\mu}^{ext}$ , respectively. Due to the linearity of the Stokes equations, the velocity $\boldsymbol{v}$ at the interface is the sum of the ambient flow $\boldsymbol{v}^{\infty }$ and a ‘perturbation’ resulting from the force density on the interface ${\it\Gamma}$ . For shear flow, $\boldsymbol{v}^{\infty }=(\dot{{\it\epsilon}}y,0,0)$ , for shear rate $\dot{{\it\epsilon}}$ (figure 1). For Poiseuille flow, $\boldsymbol{v}^{\infty }$ is described later in (3.6). Far from the drop interface, as $|\boldsymbol{x}|\rightarrow \infty$ , the influence of the interfacial forces vanishes and the boundary conditions are

(2.2a,b ) $$\begin{eqnarray}\boldsymbol{v}\rightarrow \boldsymbol{v}^{\infty }\quad \text{and}\quad p\rightarrow 0.\end{eqnarray}$$

The droplet and ambient fluid are assumed to be immiscible, but their dynamics are coupled. Velocity is continuous across the droplet surface ${\it\Gamma}$ . Thus, at point $\boldsymbol{x}$ on ${\it\Gamma}$ ,

(2.3) $$\begin{eqnarray}\boldsymbol{v}^{ext}(\boldsymbol{x})=\boldsymbol{v}^{int}(\boldsymbol{x})=\boldsymbol{v}^{s}(\boldsymbol{x})\end{eqnarray}$$

for external and internal fluid velocities $\boldsymbol{v}^{ext}$ and $\boldsymbol{v}^{int}$ . As a result of the immiscibility, fluid does not flow across the droplet surface and the droplet volume does not change. As a result, the fluid velocities are identical to the velocity of the membrane,

(2.4) $$\begin{eqnarray}\frac{\text{D}\boldsymbol{x}}{\text{D}t}=\boldsymbol{v}^{s}(\boldsymbol{x})\end{eqnarray}$$

for $\boldsymbol{x}$ on ${\it\Gamma}$ and material derivative $\text{D}/\text{D}t$ . Additionally, the droplet surface is assumed to be at mechanical equilibrium with hydrodynamic forces provided by the hydrodynamic stress tensor $\bar{{\bf\sigma}}$ . This quasistatic equilibrium requires a surface force density $\boldsymbol{f}$ , as

(2.5) $$\begin{eqnarray}[[\bar{{\bf\sigma}}]]\boldsymbol{\cdot }\boldsymbol{n}+\boldsymbol{f}=0\end{eqnarray}$$

in which $[[\bar{{\bf\sigma}}]]=\bar{{\bf\sigma}}^{ext}-\bar{{\bf\sigma}}^{int}$ is the traction jump at the droplet surface ${\it\Gamma}$ . $\boldsymbol{n}$ is the unit outward normal vector to the surface.

2.2 Mechanical properties of the surface

The mechanical properties of the droplet surface are represented by the Boussinesq–Scriven stress tensor, originally described by Boussinesq (Reference Boussinesq1913) and generalized by Scriven (Reference Scriven1960). A Boussinesq–Scriven surface is characterized by a linear constitutive law with surface tension ${\it\gamma}$ , shear surface viscosity ${\it\mu}_{s}$ and dilational surface viscosity ${\it\mu}_{d}$ . In these terms, the surface stress tensor $\bar{{\bf\sigma}}_{s}$ may be formulated as

(2.6) $$\begin{eqnarray}\bar{{\bf\sigma}}_{s}={\it\gamma}\bar{\unicode[STIX]{x1D64B}}+({\it\mu}_{d}-{\it\mu}_{s}){\it\Theta}\bar{\unicode[STIX]{x1D64B}}+2{\it\mu}_{s}\bar{\unicode[STIX]{x1D65A}}\end{eqnarray}$$

for surface projection tensor $\bar{\unicode[STIX]{x1D64B}}$ , surface rate of dilation ${\it\Theta}$ and surface rate of deformation tensor $\bar{\unicode[STIX]{x1D65A}}$ (Secomb & Skalak Reference Secomb and Skalak1982). $\bar{\unicode[STIX]{x1D64B}}$ is defined using the unit outward normal vector $\boldsymbol{n}$ to the surface and the unit tensor $\bar{\unicode[STIX]{x1D644}}$ , as

(2.7) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D64B}}=\bar{\unicode[STIX]{x1D644}}-\boldsymbol{n}\boldsymbol{n}^{\text{T}}.\end{eqnarray}$$

In terms of the surface gradient $\boldsymbol{{\rm\nabla}}_{s}$ and surface velocity $\boldsymbol{v}^{s}$ , the surface rate of deformation tensor $\bar{\unicode[STIX]{x1D65A}}$ is given by

(2.8) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D65A}}={\textstyle \frac{1}{2}}\bar{\unicode[STIX]{x1D64B}}(\boldsymbol{{\rm\nabla}}_{s}\boldsymbol{v}^{s}+(\boldsymbol{{\rm\nabla}}_{s}\boldsymbol{v}^{s})^{\text{T}})\bar{\unicode[STIX]{x1D64B}}\end{eqnarray}$$

and the surface rate of dilation ${\it\Theta}$ by

(2.9) $$\begin{eqnarray}{\it\Theta}=\bar{\unicode[STIX]{x1D64B}}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}_{s}\boldsymbol{v}^{s}.\end{eqnarray}$$

In this paper, both shear and dilational surface viscosities ${\it\mu}_{s}$ and ${\it\mu}_{d}$ are constant. Likewise, surface tension ${\it\gamma}$ is also constant and, in the absence of an external flow, the droplet’s equilibrium shape is spherical. As a result, droplets in subsequent simulations are initially spherical.

2.3 Numerical model of the surface

The numerical model represents the droplet interface with Loop elements as a subdivision surface (Loop Reference Loop1987; Cirak, Ortiz & Schroder Reference Cirak, Ortiz and Schroder2000), a powerful method to computer aided design. Using a triangular mesh based on successive refinements of an icosahedron, the subdivision method guarantees $C^{2}$ continuity almost everywhere on the surface (except on the 12 vertices of the initial icosahedron, where the surface is only $C^{1}$ ). However, these vertices do not necessitate special treatment, since membrane forces are computed by a finite element method which relaxes the continuity requirement compared to a direct computation (see below for more details). The characteristic feature of a subdivision surface is that the displacement field of an element depends on the nodal displacement of nearest-neighbour elements, in addition to its own nodal displacement. These nearest-neighbour elements are said to comprise the $1$ -ring about the element; regular elements have $12$ elements in their $1$ -ring. For a given element $e$ , position $\boldsymbol{x}$ is computed as

(2.10) $$\begin{eqnarray}\boldsymbol{x}({\it\xi},{\it\eta})=\mathop{\sum }_{n\in E_{n}}N_{n}^{e}({\it\xi},{\it\eta})\boldsymbol{x}_{n}\end{eqnarray}$$

for the curvilinear coordinates $({\it\xi},{\it\eta})$ on element $e$ , node $n$ in the $1$ -ring $E_{n}$ about element $e$ , box-spline shape functions $N_{n}^{e}$ and the nodal value $\boldsymbol{x}_{n}$ . The velocity $\boldsymbol{v}({\it\xi},{\it\eta})$ on the surface is computed in the same way as $\boldsymbol{x}$ . The box-spline shape functions used for both regular and irregular elements are discussed in the appendix of Cirak et al. (Reference Cirak, Ortiz and Schroder2000). This surface representation provides accurate derivatives with respect to ${\it\xi}$ and ${\it\eta}$ , which are necessary to compute the surface stress tensor. This method provides an improvement compared to our previous code and notably for the surface representation (Boedec, Leonetti & Jaeger Reference Boedec, Leonetti and Jaeger2011; Boedec, Jaeger & Leonetti Reference Boedec, Jaeger and Leonetti2012). Two examples of meshing of droplets with dilational viscosity in shear flow for two numbers of elements are shown in figure 2. Loop subdivision elements have been used in previous studies on capsules (Le Reference Le2010; Huang et al. Reference Huang, Chang and Sung2012) and vesicles (Spann, Zhao & Shaqfeh Reference Spann, Zhao and Shaqfeh2014). However, the numerical methods to compute membrane stress and the fluid solver are different.

Figure 2. Two cases of meshing characterized by the number $N$ of elements of a capsule with dilational viscosity in shear flow: $Ca=0.3$ , $Bq_{s}=0$ and $Bq_{d}=1$ . The colour code corresponds to the magnitude of the interface velocity. (a) $N=320$ , (b) $N=1280$ .

With this framework established, the surface velocity gradient is computed on each element as

(2.11) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}_{s}\boldsymbol{v}=\unicode[STIX]{x1D622}^{{\it\alpha}{\it\beta}}\frac{\partial v_{i}}{\partial s^{{\it\alpha}}}\frac{\partial x_{j}}{\partial s^{{\it\beta}}}\unicode[STIX]{x1D626}_{i}\otimes \unicode[STIX]{x1D626}_{j},\end{eqnarray}$$

with Cartesian components of velocity $v_{i}$ and position $x_{j}$ , local curvilinear coordinates $s^{{\it\alpha}}=({\it\xi},{\it\eta})$ and the local contravariant metric tensor $\unicode[STIX]{x1D622}^{{\it\alpha}{\it\beta}}$ . With (2.11) substituted into the components of (2.6), the Boussinesq–Scriven stress tensor may be computed. Transformed to curvilinear coordinates, the stress tensor takes the form

(2.12) $$\begin{eqnarray}{\it\sigma}^{{\it\alpha}{\it\beta}}=\unicode[STIX]{x1D622}^{{\it\alpha}{\it\gamma}}\unicode[STIX]{x1D622}^{{\it\beta}{\it\delta}}\frac{\partial x_{i}}{\partial s^{{\it\gamma}}}\frac{\partial x_{j}}{\partial s^{{\it\delta}}}{\it\sigma}_{ij}.\end{eqnarray}$$

We follow the finite element method described by Cirak et al. (Reference Cirak, Ortiz and Schroder2000) to solve the weak form of the equation

(2.13) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}_{s}\boldsymbol{\cdot }{\bf\sigma}^{{\it\alpha}{\it\beta}}-\boldsymbol{f}=0\end{eqnarray}$$

for surface force density $\boldsymbol{f}$ . Introducing a virtual displacement field ${\it\delta}\boldsymbol{x}$ , (2.13) reads:

(2.14) $$\begin{eqnarray}\int _{S}\left(\frac{1}{2}{\it\sigma}^{{\it\alpha}{\it\beta}}{\it\delta}\unicode[STIX]{x1D622}_{{\it\alpha}{\it\beta}}+\boldsymbol{f}\boldsymbol{\cdot }{\it\delta}\boldsymbol{x}\right)\hspace{-3.0pt}\,\text{d}S=0,\quad \forall {\it\delta}\boldsymbol{x},\end{eqnarray}$$

where ${\it\delta}\unicode[STIX]{x1D622}_{{\it\alpha}{\it\beta}}={\it\delta}\boldsymbol{x}_{,{\it\alpha}}\boldsymbol{\cdot }\boldsymbol{x}_{,{\it\beta}}+\boldsymbol{x}_{,{\it\alpha}}\boldsymbol{\cdot }{\it\delta}\boldsymbol{x}_{,{\it\beta}}$ is the virtual variation of the metric. Thus, with the weak form of membrane equilibrium, only first derivatives of shape and virtual displacement are needed to compute membrane forces. Using the same interpolation basis (Loop elements) for the virtual displacement ${\it\delta}\boldsymbol{x}$ and for the force $\boldsymbol{f}$ , this equation writes:

(2.15) $$\begin{eqnarray}\mathop{\sum }_{e=1}^{N_{el}}\int _{S^{e}}\mathop{\sum }_{n=1}^{N^{e}}\mathop{\sum }_{m=1}^{N^{e}}\left(\frac{1}{2}{\it\sigma}^{{\it\alpha}{\it\beta}}[(N_{n,{\it\alpha}}N_{m,{\it\beta}}+N_{n,{\it\beta}}N_{m,{\it\alpha}}){\it\delta}x_{i}^{n}x_{i}^{m}]+N_{n}N_{m}f_{i}^{n}{\it\delta}x_{i}^{m}\right)\hspace{-3.0pt}\,\text{d}S^{e}=0,\end{eqnarray}$$

where the integral over the surface $S$ has been split as a sum of integration over element $e$ surface $S^{e}$ , with $N_{e}$ elements in the mesh and the interpolation (2.10) by Loop element has been used, with $N_{n}$ as a notation for the $n$ th shape function of element $e$ , evaluated at coordinates $({\it\xi},{\it\eta}):N_{n}=N_{n}^{e}({\it\xi},{\it\eta})$ and $x_{i}^{n}$ stands for the nodal value $n$ of the $i$ th Cartesian component of position. The integration over surfaces $S^{e}$ is computed using Gauss quadrature points. Thus, after rearranging the terms, and solving the linear system with the nodal values of Cartesian component of force $f_{i}^{n}$ , the computation of $\boldsymbol{f}$ may be then formulated symbolically in terms of linear operators $F^{{\it\gamma}}$ and $F^{{\it\nu}}$ acting on ${\it\gamma}$ and $\boldsymbol{v}$ , respectively, as

(2.16) $$\begin{eqnarray}\boldsymbol{f}=F^{{\it\gamma}}({\it\gamma})+F^{{\it\nu}}(\boldsymbol{v}).\end{eqnarray}$$

2.4 Boundary element model and update scheme

The Stokes equations are solved using the standard boundary element method (Pozrikidis Reference Pozrikidis1992). The boundary element method is evaluated over the same mesh used to describe the droplet surface and velocity is assumed to be continuous across this interface. For Cartesian component $i$ of the velocity $\boldsymbol{v}$ at $\boldsymbol{x}$ , the boundary integral is

(2.17) $$\begin{eqnarray}\displaystyle v_{i}(\boldsymbol{x}) & = & \displaystyle v_{i}^{\infty }(\boldsymbol{x})+\frac{1}{8{\rm\pi}{\it\mu}^{ext}}\int _{S}G_{ij}(\boldsymbol{x}_{0},\boldsymbol{x})f_{j}(\boldsymbol{x}_{0})\,\text{d}S(\boldsymbol{x}_{0})\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1-{\it\lambda}}{8{\rm\pi}}\int _{S}\left[v_{i}(\boldsymbol{x}_{\mathbf{0}})-v_{i}(\boldsymbol{x})\right]T_{ijk}(\boldsymbol{x},\boldsymbol{x}_{\mathbf{0}})n_{k}(\boldsymbol{x}_{\mathbf{0}})\,\text{d}S(\boldsymbol{x}_{\mathbf{0}}),\end{eqnarray}$$

with free space Stokeslet $G_{ij}=({\it\delta}_{ij}/r)+(X_{i}X_{j}/r^{3})$ and Stresslet $T_{ijk}=-6(X_{i}X_{j}X_{k}/r^{5})$ , point $\boldsymbol{x}_{0}$ on ${\it\Gamma}$ , vector $\boldsymbol{X}=\boldsymbol{x}_{0}-\boldsymbol{x}$ and $r=\Vert \boldsymbol{X}\Vert$ . The viscosity contrast has been introduced as follows:

(2.18) $$\begin{eqnarray}{\it\lambda}={\it\mu}^{int}/{\it\mu}^{ext}.\end{eqnarray}$$

As a result, at point $\boldsymbol{x}$ on ${\it\Gamma}$ , the discretized version of (2.17) may be written as

(2.19) $$\begin{eqnarray}\boldsymbol{v}=\boldsymbol{v}^{\infty }+\unicode[STIX]{x1D642}\boldsymbol{f}+(1-{\it\lambda})\unicode[STIX]{x1D64F}\boldsymbol{v}.\end{eqnarray}$$

Taking into account the representation of the interfacial load in (2.16), the velocity of the droplet interface $\boldsymbol{v}(t)$ is computed by solving the equation

(2.20) $$\begin{eqnarray}\boldsymbol{v}(\boldsymbol{x})=\boldsymbol{v}^{\infty }(\boldsymbol{x})+\unicode[STIX]{x1D642}(F^{{\it\gamma}}({\it\gamma})+F^{{\it\nu}}(\boldsymbol{v}))+(1-{\it\lambda})\unicode[STIX]{x1D64F}\boldsymbol{v}.\end{eqnarray}$$

This equation could be integrated explicitly in time, using the previous velocity field $\boldsymbol{v}^{\boldsymbol{t}}$ to compute $\boldsymbol{x}^{\boldsymbol{t}+\text{d}\boldsymbol{t}}$ . Doing so would result in a constraint on the time step depending both on surface tension and surface viscosities. To avoid the constraint due to surface viscosity, we choose instead to solve (2.20) at each substep of the time stepping algorithm. Thus, as the equation is implicit for velocity $\boldsymbol{v}(t)$ , it is recast as

(2.21) $$\begin{eqnarray}(I-\unicode[STIX]{x1D642}F^{{\it\nu}}-(1-{\it\lambda})\unicode[STIX]{x1D64F})\boldsymbol{v}(\boldsymbol{x})=\boldsymbol{v}^{\infty }(\boldsymbol{x})+\unicode[STIX]{x1D642}F^{{\it\gamma}}({\it\gamma}).\end{eqnarray}$$

A generalized minimum residual (GMRES) solver is used to solve (2.21) at each time step, providing a fully implicit solution for $\boldsymbol{v}(t)$ . As a result, this approach differs from similar methods employing Boussinesq–Scriven stress, such as Rodrigues et al. (Reference Rodrigues, Ausas, Mut and Buscaglia2015), which provide only semi-implicit solutions to the fully coupled problem. The cost of an implicit solution and iterative method, of course, is that multiple linear systems of equations are solved at each time step. Nonetheless, the iterative approach ensures that the primary restriction on time step $\text{d}t$ is surface tension, not surface viscosity. According to the no-slip boundary condition in (2.4), updated position $\boldsymbol{x}(t+\text{d}t)$ on the droplet is computed using the Runge–Kutta–Fehlberg method RK45, with an adaptive step size $\text{d}t$ . The same method has been recently used to calculate the dynamics of an elastic capsule without bending resistance in a planar elongation flow (de Loubens et al. Reference de Loubens, Deschamps, Boedec and Leonetti2015a ).

As the droplet elongates in flow, the quality of the elements composing the mesh might degrade, especially for near-breakup simulations. Note that the inclusion of surface viscosities tends generally to smooth the elongation of elements on the whole mesh, thus leading to a better overall mesh quality. When needed, we use a remeshing algorithm which updates the position of the nodes while preserving the deformed shape. To do so, we slide the nodes tangentially on the surface, in order to have more nodes in regions of high curvature, while also keeping area of elements locally homogeneous.

2.5 Dimensionless parameters and variables

Several dimensionless parameters describe the surface tension and viscosities. The ratio of hydrodynamic stress to the resistance of surface tension ${\it\gamma}$ is reflected in the capillary number (or Weber number)

(2.22) $$\begin{eqnarray}Ca=\frac{{\it\mu}^{ext}\dot{{\it\epsilon}}a}{{\it\gamma}}\end{eqnarray}$$

for ambient fluid viscosity ${\it\mu}^{ext}$ , shear rate $\dot{{\it\epsilon}}$ and initial droplet radius $a$ . Borrowing the terminology of Erni (Reference Erni2011), surface viscosity strengths are measured by the dimensionless Boussinesq numbers,

(2.23) $$\begin{eqnarray}\displaystyle Bq_{s} & = & \displaystyle \frac{{\it\mu}_{s}}{{\it\mu}^{ext}a}\end{eqnarray}$$

(2.24) $$\begin{eqnarray}\displaystyle Bq_{d} & = & \displaystyle \frac{{\it\mu}_{d}}{{\it\mu}^{ext}a}\end{eqnarray}$$

$Bq_{s}=Bq_{d}$

$Bq$

$Bq=0$

${\it\beta}$

(2.25) $$\begin{eqnarray}{\it\beta}=(Bq)(Ca).\end{eqnarray}$$

Introduced by Barthès-Biesel & Sgaier (Reference Barthès-Biesel and Sgaier1985) for capsules with surface viscosity, ${\it\beta}$ is comparable to the Weissenberg number in rheology. In the part on the results of this paper, we discuss clean drops with a contrast of viscosity ${\it\lambda}$ and drops with a viscous interface.

Additionally, physical quantities are non-dimensionalized: the lengths by the droplet radius a, the time by the inverse of the shear rate $1/\dot{{\it\epsilon}}$ , the pressure and membrane stress by ${\it\mu}^{ext}\dot{{\it\epsilon}}$ .

3.1 Clean droplet in shear flow

The shape of a clean droplet ( $Bq=0$ ) in shear flow has been extensively studied numerically in the past using the boundary element method (Rallison Reference Rallison1981; Kennedy et al. Reference Kennedy, Pozrikidis and Skalak1994), advancing front method (Kwak & Pozrikidis Reference Kwak and Pozrikidis1998) and volume-of-fluid method (Li et al. Reference Li, Renardy and Renardy2000), among many others. Such computational work has been shown to be consistent with experimental results (Rumscheidt & Mason Reference Rumscheidt and Mason1961; Kwak & Pozrikidis Reference Kwak and Pozrikidis1998). Many authors have also determined well known analytical results based on the main assumption of small deformation theory: (Taylor Reference Taylor1934; Chaffey & Brenner Reference Chaffey and Brenner1967; Cox Reference Cox1969; Barthès-Biesel & Acrivos Reference Barthès-Biesel and Acrivos1973; Rallison Reference Rallison1980; Vlahovska, Blawzdziewicz & Loewenberg Reference Vlahovska, Blawzdziewicz and Loewenberg2005, Reference Vlahovska, Blawzdziewicz and Loewenberg2009a ). These approaches have their own domain of validity (Acrivos Reference Acrivos1983; Rallison Reference Rallison1984). We recall below the analytical results well discussed in the literature which are sufficient for the validation of the code.

The shape in shear flow is mainly characterized by the so-called Taylor deformation parameter $D=(L-B)/(L+B)$ , in which $L$ and $B$ are the major and minor axes of the ellipsoid having the same moment of inertia as the droplet. These axes remain in the $z=0$ plane for the duration of the simulations. Another salient parameter is the inclination ${\it\theta}$ which represents the angle between the direction of flow and the major axis of the droplet.

In his seminal work, Taylor (Reference Taylor1934) has determined $D$ at the first order in $Ca$ and inclination ${\it\theta}$ at the zeroth order in $Ca$ in both limits of low and high viscosity contrast  ${\it\lambda}$ :

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\lambda}\gg 1;\quad Ca=O(1)\rightarrow D_{T}=\frac{5}{4{\it\lambda}};\quad {\it\theta}=0 & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\lambda}=O(1);\quad Ca\ll 1\rightarrow D_{T}=Ca\frac{19{\it\lambda}+16}{16{\it\lambda}+16};\quad {\it\theta}=\frac{{\rm\pi}}{4}. & \displaystyle\end{eqnarray}$$

Cox (Reference Cox1969) has proposed compact analytical equations which cover a large range of viscosity contrast assuming slight deformations from the spherical shape which can be quantified for example by $D\ll 1$ :

(3.3a,b ) $$\begin{eqnarray}D_{Cox}=\frac{Ca{\displaystyle \frac{19{\it\lambda}+16}{16{\it\lambda}+16}}}{\sqrt{1+\left({\displaystyle \frac{19{\it\lambda}Ca}{20}}\right)^{2}}};\quad {\it\theta}_{Cox}=\frac{{\rm\pi}}{4}-\frac{1}{2}\text{atan}\left(\frac{19{\it\lambda}Ca}{20}\right).\end{eqnarray}$$

If the viscous stress is large ${\it\lambda}Ca\gg 1$ (or very small), (3.1) (or (3.2)) is obtained. Further analysis has shown that the equations of Cox (Reference Cox1969) are in fact valid in the limit of large viscosity contrast (Rallison Reference Rallison1980). However, these equations are useful to validate our code, to highlight the global correct agreement between theory and numerics, the consistence of the high viscosity contrast limit and the need of high numerical accuracy to distinguish numerically models even if theory has already established their respective validities.

Figure 3. Clean droplet. Steady-state deformation: the Taylor parameter $D=(L-B)/(L+B)$ of a clean droplet in the limit of small capillary number $Ca=10^{-3}$ . Comparisons are made with analytical results of Taylor (Reference Taylor1934) and Cox (Reference Cox1969).

Figure 4. Clean droplet. Steady-state inclination angle ${\it\theta}$ of a clean droplet in the limit of small capillary number $Ca=10^{-3}$ . Comparisons are made with analytical results of Taylor (Reference Taylor1934) and Cox (Reference Cox1969).

We performed numerical simulations of a clean droplet in shear flow covering 12 decades of viscosity contrast from $10^{-6}$ to $10^{6}$ . The capillary number was chosen as $Ca=10^{-3}$ to be consistent with the assumption of small deformation of models. The number of elements were 5120 elements without remeshing and 1280 elements with remeshing without notable differences. Both limits of Taylor (Reference Taylor1934) are recovered as shown in figures 3 for the deformation and 4 for the inclination angle at weak and high viscosity contrasts. In the case of the Taylor parameter, the error $|D-D_{Cox}|/D_{Cox}$ on the Taylor parameter $D$ is approximately $3.0\times 10^{-5}$ for ${\it\lambda}<10$ , less than $3.0\times 10^{-6}$ for ${\it\lambda}\geqslant 10^{5}$ and reaches the maximum $2.0\times 10^{-4}$ for ${\it\lambda}=100$ . This very good agreement could be due to the choice of a very small capillary number which is true in part. Then, the capillary number has also been varied from $0.001$ to $0.4$ for a range of viscosity contrast from $2$ to $100$ : figure 5. We have checked that the slopes of numerical simulations of $D$ versus the capillary number are in excellent agreement with the theory of Taylor (Reference Taylor1934). It is quite surprising that the results are in good agreement with the equations of Cox (Reference Cox1969) up to $Ca=0.1$ for a large range of viscosity contrast. But they become not satisfactory for $Ca\geqslant 0.3$ . To conclude on the Taylor parameter of deformation $D$ , the comparisons between numerics, analytical results of Taylor (Reference Taylor1934) and Cox (Reference Cox1969) are very satisfactory at high viscosity contrast as expected but also at small capillary numbers. However, the Cox (Reference Cox1969) model should fail to describe accurately the shapes of droplets at moderate and small capillary numbers (Rallison Reference Rallison1980). Our numerical code should be able to highlight this discrepancy on the linear variations of deformation parameters and provides additional information on the domains of validities of models.

Figure 5. Clean droplet. Steady-state deformation for a large range of capillary numbers $Ca$ from $0.001$ to $0.4$ and for various viscosity contrasts ${\it\lambda}$ equal to 2, 5, 10 and 100.

At $Ca=10^{-3}$ , the error on the inclination angle is larger, of the order of $10^{-3}$ . Thus, we can expect a larger error at higher capillary numbers such as $0.1$ for example which is not the case for $D$ . In the limit of small capillary number and small deformation, Chaffey & Brenner (Reference Chaffey and Brenner1967) derived another equation for the inclination angle as a function of viscosity contrast ${\it\lambda}$ and capillary number $Ca$ :

(3.4) $$\begin{eqnarray}{\it\theta}_{CB}=\frac{{\rm\pi}}{4}-Ca\frac{(19{\it\lambda}+16)(2{\it\lambda}+3)}{80(1+{\it\lambda})}.\end{eqnarray}$$

This equation has been recovered by several authors (Barthès-Biesel & Acrivos Reference Barthès-Biesel and Acrivos1973; Rallison Reference Rallison1980; Vlahovska et al. Reference Vlahovska, Blawzdziewicz and Loewenberg2005, Reference Vlahovska, Blawzdziewicz and Loewenberg2009a ) who calculate the power expansion in a small parameter, mainly the capillary number $Ca$ , of all the physical quantities: shape deformation, curvature, pressure, velocities and stress tensor. It means in particular that the viscous stress must be small, ${\it\lambda}Ca\ll 1$ contrary to the model of Cox (Reference Cox1969). Note that this result has also been validated experimentally (Guido & Villone Reference Guido and Villone1998; Guido, Greco & Villone Reference Guido, Greco and Villone1999).

Figure 6. Clean droplet. Inclination angle as a function of the capillary number $Ca$ and the viscosity contrast: dashed lines (full Cox equation), line (linear Chaffey–Brenner equation), symbols (numerical simulations). Colour code corresponds to ${\it\lambda}=1$ (black), ${\it\lambda}=5$ (blue), ${\it\lambda}=10$ (green) and ${\it\lambda}=100$ (red, $x$ ). Symbol code corresponds to ${\it\lambda}=1$ (disc), ${\it\lambda}=5$ (square), ${\it\lambda}=10$ (asterisk) and ${\it\lambda}=100$ (cross).

This relation is completely different from (3.3). Indeed, for a zero viscosity contrast, the result of Cox (Reference Cox1969) is ${\rm\pi}/4$ whatever the capillary number in agreement with Taylor (Reference Taylor1934) but in contradiction with the expression of Chaffey & Brenner (Reference Chaffey and Brenner1967): ${\it\theta}_{CB}={\rm\pi}/4-(3/5)Ca$ . Moreover, it is possible to develop (3.3) in the limit of small viscous stress ${\it\lambda}Ca\ll 1$ providing the linear variation of the inclination with the capillary number, highlighting the difference between the theoretical results. Thus, ${\it\theta}_{Cox}={\rm\pi}/4-Ca(19{\it\lambda}/40)$ . This is especially important to validate our numerical code. With a fixed viscosity contrast, the slopes are different meaning that the code must be able to differentiate clearly the models by a careful measurement of the inclination angle. Moreover, a linear variation as (3.4) valid for small viscous stress and $Ca$ is always an opportunity to validate a numerical code. It is sufficient to decrease the capillary number up to a clear linear variation and to measure the slope. To perform this, we have varied the capillary number from $0.001$ to $0.4$ and the viscosity contrast from 1 to 100. The number of elements was 1280 or 5120. The residual of GMRES was fixed to $10^{-9}$ for the resolution of (2.21) and the precision to $10^{-8}$ for the selection of the maximal step size in the RK45 time integrator. As expected, for ${\it\lambda}Ca>1$ , agreement is not reached with the (3.4) but there is good agreement with the Cox results as shown in figure 6(a). On the contrary, for ${\it\lambda}Ca<1$ , the agreement is excellent with (3.4): figure 6(b) and the curve ${\it\lambda}=1$ for figure 6(a). If this comparison allows us to validate our code for a clean droplet in Stokes flow, our calculations also permit us to extract the limit of validity of (3.4) which is ${\it\lambda}Ca\leqslant 0.2$ . Indeed, with this criterion and regardless of the viscosity, the agreement between the result of Chaffey & Brenner (Reference Chaffey and Brenner1967) and our numerical results are excellent. Numerical studies on drops and capsules focus often on the Taylor parameter of deformation $D$ , and more scarcely on the inclination angle while in this study, the latter is however essential to validate our code by differentiating unambiguously the linear models. Indeed, without deformation, the droplet radius is known while the inclination angle is not definite at this zeroth order. When the droplet deforms slightly from the spherical shape, the first order corresponds to the equations of Taylor (Reference Taylor1934) which provides an inclination angle of ${\rm\pi}/4$ and a linear variation of $D$ with $Ca$ . Thus, equation of Chaffey & Brenner (Reference Chaffey and Brenner1967) provides an upper order of the inclination angle which is more dependent on models while remaining linear.

Figure 7. Clean droplet with ${\it\lambda}=1$ . Comparisons with previous numerical studies. (a) The Taylor parameter $D$ versus the capillary number $Ca$ . (b) The inclination angle versus $Ca$ .

In the literature, to our knowledge, the case ${\it\lambda}=1$ has mainly been studied for comparisons with analytical expressions and to validate numerical codes: figure 7. Our results depart from linear theory earlier than Kwak & Pozrikidis (Reference Kwak and Pozrikidis1998), at $Ca=0.25$ , but match very well Li et al. (Reference Li, Renardy and Renardy2000), in doing so considering the Taylor parameter $D$ . The inclination angle ${\it\theta}$ differs more among the numerical studies in the literature. The results of Kennedy et al. (Reference Kennedy, Pozrikidis and Skalak1994) and Li et al. (Reference Li, Renardy and Renardy2000) do not really match the linear theory of Chaffey & Brenner (Reference Chaffey and Brenner1967) but provide a good approximation and in any case, are largely better than the linear result of Cox (Reference Cox1969). Indeed, the slopes at ${\it\lambda}=1$ are $35/32\approx 1.09$ (Chaffey & Brenner Reference Chaffey and Brenner1967) compared with $0.5\ast (19/20)=0.475$ (Cox Reference Cox1969). While the previous numerical results are performed for a capillary number larger than 0.1, the present work extends the study to a minimum capillary number of $10^{-3}$ , allowing a full comparison with linear theory, notably in the more difficult case of an inclination angle. An additional concern for this line of inquiry is the critical capillary number $Ca_{c}$ , beyond which a steady-state deformation does not exist for a clean droplet with ${\it\lambda}=1$ . We find the critical value $Ca_{c}\approx 0.43$ , which is consistent with the range $0.37-0.43$ found in previous computational and experimental studies (Rallison Reference Rallison1981; Kennedy et al. Reference Kennedy, Pozrikidis and Skalak1994; Li et al. Reference Li, Renardy and Renardy2000; Cristini et al. Reference Cristini, Guido, Alfani, Bławzdziewicz and Loewenberg2003; Fischer & Erni Reference Fischer and Erni2007).

3.2 Droplet with viscous interface in shear flow

First, to evaluate the accuracy and convergence of the surface viscosity model, we consider a spherical droplet with $Ca=\infty$ and $Bq=Bq_{s}=Bq_{d}=1$ placed in shear flow. The initial viscous force distribution on the surface is computed and compared with the analytical values for each mesh. If the velocity far from the droplet is $\boldsymbol{V}=\dot{{\it\epsilon}}y\unicode[STIX]{x1D626}_{\boldsymbol{x}}$ , the analytical viscous force is $\boldsymbol{f}=(-3y+8yx^{2})\unicode[STIX]{x1D626}_{\boldsymbol{x}}+(-3x+8xy^{2})\unicode[STIX]{x1D626}_{\boldsymbol{y}}+8xyz\unicode[STIX]{x1D626}_{\boldsymbol{z}}$ . Error is computed by calculating the $\ell ^{2}$ norm of the relative error at each node and using the Voronoi region about each node to integrate over the surface of the sphere. As seen in figure 8(a), first-order convergence of the error with respect to the number of elements $N$ is observed.

Figure 8. Droplet with viscous interface and ${\it\lambda}=1$ . (a) Error in viscous force for a spherical droplet with $Bq=Bq_{s}=Bq_{d}=1$ and infinite $Ca$ . The number of elements $N$ varies from 80 to 20 480. The order of convergence is one (solid line). (b) Variation with dimensionless time $\dot{{\it\epsilon}}t$ of the Taylor parameter $D$ with two sets of parameters $(Ca;{\it\lambda};Bq)=(0.5;1;10)$ and $(Ca;{\it\lambda};Bq)=(0.35;0.1;5)$ for three numbers $N$ of elements.

Second, to evaluate the numerical convergence over the duration of a simulation, we consider a droplet with a viscous surface in shear flow. The dependence of the deformation $D$ on grid refinement (or number of elements) is shown in figure 8(b), for two droplets with different viscosity contrast and Boussinesq number in shear flow. While there is noticeable error for the coarsest mesh (using only 80 elements), finer meshes agree very well for both sets of parameters. Subsequent simulations are conducted with 320 elements except for the breakup which can necessitate 1280 elements with remeshing or more.

3.3 Droplet with viscous interface in Poiseuille flow

A third validation setting of a viscous interface is provided by Poiseuille flow. Recent studies of droplets with viscous surfaces, Schwalbe et al. (Reference Schwalbe, Phelan, Vlahovska and Hudson2011) and Reusken & Zhang (Reference Reusken and Zhang2013), consider the migration velocity of a droplet in Poiseuille flow and propose it as a benchmark problem. Indeed, Schwalbe’s analysis establishes that the droplet’s migration velocity for ${\it\lambda}=1$ is

(3.5) $$\begin{eqnarray}\frac{\boldsymbol{U}_{analytic}}{{\it\alpha}a^{2}}=-\frac{2Bq_{d}+3}{3(2Bq_{d}+5)}\unicode[STIX]{x1D626}_{y}\end{eqnarray}$$

for a droplet with radius $a$ . Thus, as $Bq_{d}\rightarrow \infty$ , $\boldsymbol{U}_{m}/{\it\alpha}a^{2}$ tends toward $-(1/3)\unicode[STIX]{x1D626}_{y}$ .

Thus, a spherical droplet is centred in a planar Poiseuille flow $v^{\infty }$ ,

(3.6) $$\begin{eqnarray}\boldsymbol{v}^{\infty }=(V-{\it\alpha}y^{2})\unicode[STIX]{x1D626}_{x},\end{eqnarray}$$

for speed $V$ at the centre line and ${\it\alpha}$ being proportional to the flow profile’s curvature (Schwalbe et al. Reference Schwalbe, Phelan, Vlahovska and Hudson2011). The migration velocity $\boldsymbol{U}_{m}$ may then be defined as

(3.7) $$\begin{eqnarray}\boldsymbol{U}_{m}(t)=\frac{1}{{\it\Omega}(t)}\int _{{\it\Omega}(t)}(\boldsymbol{v}(\boldsymbol{x},t)-\boldsymbol{v}^{\infty }(0))\,\text{d}\boldsymbol{x}\end{eqnarray}$$

for velocity $\boldsymbol{v}(\boldsymbol{x},t)$ on the droplet, undisturbed velocity $\boldsymbol{v}^{\infty }(0)$ at the centre line and droplet volume ${\it\Omega}$ (Reusken & Zhang Reference Reusken and Zhang2013). We consider a capillary number $Ca=0.01$ to maintain the deformations below $O(10^{-4})$ . The results were checked by another approach taking $Ca=1$ but a very small time step of $10^{-5}$ . The velocity is measured at short times preventing the appearance of any deformation (no inertia). There are no differences between the two calculations.

Figure 9. Droplet with viscous interface. (a) Error for $\boldsymbol{U}_{m}/{\it\alpha}a^{2}$ for $Bq_{d}=0$ , $1$ and $10$ with a second-order convergence. (b) Computed $\boldsymbol{U}_{m}/{\it\alpha}a^{2}$ compared with Schwalbe et al. (Reference Schwalbe, Phelan, Vlahovska and Hudson2011).

The instantaneous $\boldsymbol{U}_{m}$ was computed and the error is defined as $\Vert \boldsymbol{U}_{m}-\boldsymbol{U}_{analytic}\Vert _{2}$ . A grid refinement study is shown in figure 9(a), for $Bq_{d}=0$ , $1$ and $10$ . The convergence is of second order with respect to the number of triangular elements $N$ . A comparison with analytic results for a broad range of $Bq_{d}$ is given in figure 9(b), using a mesh with $N=320$ triangles. Again, excellent agreement occurs over the entire set of $Bq_{d}$ considered.