2.1. Averaged equations for dispersed two-phase flows

To clarify the respective roles of forces and higher-order moments in dispersed two-phase flow modelling, we introduce the averaged mass and momentum equations for the dispersed and continuous phases, obtained through an ensemble-averaging procedure. For a monodisperse suspension, the governing mass and momentum equations for each phase can be expressed as (Fintzi & Pierson Reference Fintzi and Pierson2026b )

(2.1) \begin{align} \phi _{\!f} + \phi &\approx 1, \\[-12pt]\nonumber \end{align}

(2.2) \begin{align} \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{U} &= 0 , \\[-12pt]\nonumber \end{align}

(2.3) \begin{align} (\partial _t + \boldsymbol{U}_{\!p} \boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }})\phi &= - \phi \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{U}_{\!p}, \\[-12pt]\nonumber \end{align}

(2.4) \begin{align} \rho _p \phi (\partial _t + \boldsymbol{U}_{\!p} \boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }})\boldsymbol{U}_{\!p} &= \phi (\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol \varSigma + \rho _p \boldsymbol{g}) + \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol \varSigma ^p + \boldsymbol{F} , \\[-12pt]\nonumber \end{align}

(2.5) \begin{align} \rho _{\!f} \phi _{\!f} (\partial _t + \boldsymbol{U}_{\!f} \boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }}) \boldsymbol{U}_{\!f} &= \phi _{\!f} \left (\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{\varSigma } + \rho _{\!f} \boldsymbol{g}\right ) + \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol \varSigma ^f - \boldsymbol{F}, \end{align}

respectively. Note that in contrast to Fintzi & Pierson (Reference Fintzi and Pierson2026b ), all ensemble averaged quantities are written using capital letters. The subscripts $_{\!f}$ and $_p$ refer to the continuous phase and the dispersed phase, respectively. The vector $\boldsymbol{g}$ is the acceleration of gravity, $\rho _k$ is the density, $\mu _k$ is the viscosity and $\boldsymbol{U}_k$ is the averaged velocity of phase $k$ . We use the notation $\phi = n_{\!p} v_p$ , where $n_{\!p}$ is the droplet number density and $v_p$ is the volume of a single droplet. Here $\phi _{\!f}$ is the volume fraction of the continuous phase. Here $\boldsymbol \varSigma = - P_{\!f} \boldsymbol \delta + 2 \mu _{\!f} \boldsymbol{E}$ , with $ \boldsymbol{E} = ( {1}/{2})[\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{U} + (\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{U})^\dagger ]$ , is the mean Newtonian stress of the mixture, with $P_{\!f}$ being the mean hydrodynamic pressure, and $\boldsymbol{U} = \phi _{\!f} \boldsymbol{U}_{\!f} + \phi \boldsymbol{U}_{\!p}$ the averaged velocity of the mixture (note that the definition of U is only exact in the homogeneous regime (Fintzi & Pierson Reference Fintzi and Pierson2026b )) $, \boldsymbol{\varSigma }^p$ and $\boldsymbol{\varSigma }^f$ are the effective stresses of the dispersed and continuous phase, respectively. Finally, $\boldsymbol{F}$ represents the interphase momentum exchange.

The number density, continuous phase volume fraction, interphase force and effective stresses of the dispersed and continuous phase can be expressed formally as (Fintzi & Pierson Reference Fintzi and Pierson2026b ),

(2.6) \begin{align} n_{\!p} &= \big \langle \delta _p\big \rangle , \\[-12pt]\nonumber \end{align}

(2.7) \begin{align} \phi _{\!f} &= \big \langle \chi _{\!f}\big \rangle , \\[-12pt]\nonumber \end{align}

(2.8) \begin{align} \boldsymbol{F} &= \left \langle \delta _p \oint _{\varGamma _{\alpha }} \boldsymbol \sigma ^*_{\!f}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma \right \rangle \!, \\[-12pt]\nonumber \end{align}

(2.9) \begin{align} \boldsymbol{\varSigma }^p &= -\rho _pv_p\big \langle \delta _p \boldsymbol{u}_\alpha '\boldsymbol{u}_\alpha '\big \rangle , \\[-12pt]\nonumber \end{align}

(2.10) \begin{align} \boldsymbol{\varSigma }^f &= - \big \langle \chi _{\!f}\rho _{\!f} \boldsymbol{u}_{\!f}'\boldsymbol{u}_{\!f}'\big \rangle + \left \langle \delta _p{\oint _{\varGamma _{\alpha }} \boldsymbol{r}\boldsymbol \sigma ^{*}_{\!f}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \delta _p\int _{\varOmega _{\alpha }} 2\mu _{\!f}\boldsymbol{e}_d^* {\rm d}\varOmega }\right \rangle \nonumber \\ &- \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\left \langle \delta _p{ \frac {1}{2}\oint _{\varGamma _{\alpha }} \boldsymbol{rr}\boldsymbol \sigma ^{*}_{\!f}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \delta _p\int _{\varOmega _{\alpha }} 2\mu _{\!f} \boldsymbol{r} \boldsymbol{e}_d^* {\rm d}\varOmega }\right \rangle + \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\boldsymbol{\nabla }} (\ldots ), \end{align}

respectively. Here the operator $\langle \ldots \rangle$ denotes ensemble averaging, $\chi _{\!f}$ is the phase indicator function of the continuous phase, $\delta _p = \sum _\alpha \delta (\boldsymbol{x}-\boldsymbol{x}_\alpha )$ is the distribution pointing on the particles centre of mass (denoted $\boldsymbol{x}_\alpha$ ) and n the unit vector normal (outward) to the droplets surface. Here $\boldsymbol{u}_\alpha$ is the centre of mass velocity of a particle labelled $\alpha$ . Here $\varOmega _\alpha$ and $\varGamma _\alpha$ represent the volume and surface, respectively, of the droplets $\alpha$ . The superscript $'$ indicates the relative values of a quantity with respect to its phase-or particle-averaged value. Specifically, $p_{\!f}' = p_{\!f}^0 - P_{\!f}$ , $\boldsymbol{u}_\alpha ' = \boldsymbol{u}_\alpha - \boldsymbol{U}_{\!p}$ and $\boldsymbol{u}_{\!f}' = \boldsymbol{u}_{\!f}^0 -\boldsymbol{U}_{\!f}$ , with $\boldsymbol{u}_{\!f}^0$ and $p_{\!f}^0$ , the local (non-averaged) velocity and pressure fields of the continuous phase. The averaged velocities and pressure are defined as $\boldsymbol{U}_{\!p} = \langle \delta _p{\boldsymbol{u}_\alpha }\rangle / n_{\!p}$ , $\boldsymbol{U}_{\!f} = \langle {\chi _{\!f} \boldsymbol{u}_{\!f}^0 }\rangle / \phi _{\!f}$ and $P_{\!f} = \langle {\chi _{\!f} p_{\!f}^0}\rangle / \phi _{\!f}$ , respectively. The superscript $^*$ represents the relative values of a quantity with respect to the mixture ensemble averaged value, such that $\boldsymbol{\sigma }_{\!f}^* = \boldsymbol{\sigma }_{\!f}^0 - \boldsymbol{\varSigma }$ and $\boldsymbol{e}_d^* = \boldsymbol{e}_d^0 - \boldsymbol{E}$ , where $\boldsymbol{\sigma }_{\!f} = -p_{\!f}^0 + \mu _{\!f} [\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\!f}^0 + (\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\!f}^0)^\dagger ]$ is the local stress of the continuous phase, and $\boldsymbol{e}_d^0 =( {1}/{2})(\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_d^0+^\dagger \boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_d^0)$ , is the internal shear rate within the droplet phase. Note that $\boldsymbol{u}_d^0$ is the internal velocity field of the fluid within the droplets, in opposition to $\boldsymbol{u}_\alpha$ which is the centre of mass velocity.

In this work, our focus is specifically on $\boldsymbol{F}$ and on the second and third contributions appearing in the effective stress tensor $\boldsymbol \varSigma ^f$ .

2.2. The closure problem

In the present context, the closure problem consists in expressing all ensemble-averaged terms of the form $\langle {\ldots }\rangle$ in (2.8) to (2.10) solely in terms of the macroscopic variables namely $n_{\!p}$ , $\phi _{\!f}$ , $P_{\!f}$ , $\boldsymbol{U}_{\!p}$ and $\boldsymbol{U}_{\!f}$ . We now assume that the droplets are spherical with radius $a$ . The surface exchange terms may be rewritten in terms of integrals of conditional averaged quantities (Lhuillier Reference Lhuillier1992; Zhang & Prosperetti Reference Zhang and Prosperetti1997; Fintzi Reference Fintzi2025). If we take the example of the drag force term, we obtain

(2.11) \begin{equation} \left \langle \delta _p \oint _{\varGamma _{\alpha }} \boldsymbol \sigma ^*_{\!f}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma \right \rangle = \int _{\mathbb{R}^3} P[\boldsymbol{w}|\boldsymbol{x},t] n_{\!p}[\boldsymbol{x},t]\oint _{\varGamma _{p}} \Big \{\! -p_{\!f}^1 \boldsymbol \delta + \mu _{\!f} [\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}^1 + ^\dagger \boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}^1] \Big \}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma (\boldsymbol{r}) d^3 \boldsymbol{w}, \end{equation}

where $p_{\!f}^1[\boldsymbol{r}|\boldsymbol{x},\boldsymbol{w}]$ and $\boldsymbol{u}^1[\boldsymbol{r}|\boldsymbol{x},\boldsymbol{w}]$ , are the pressure and velocity fields ensemble averaged on every configuration where a droplet is located at $\boldsymbol{x}$ with centre of mass velocity w, minus the ensemble averaged pressure and velocity fields, i.e. $P_{\!f}[\boldsymbol{r},t]$ , and $\boldsymbol{U}[\boldsymbol{r},t]$ , r being the points on the droplets surface. Here $\varGamma _p$ denotes the domain of the surface of the test droplet, whose surface is described by the equation $|\boldsymbol{x}- \boldsymbol{y}| - a=r - a = 0$ . Here, $ P[\boldsymbol{w}|\boldsymbol{x},t]$ is the probability of finding a droplet with centre of mass velocity w knowing that a droplet is present at x at time $t$ . The equations governing $p_{\!f}^1$ and $\boldsymbol{u}^1$ can be obtained by conditionally averaging the local mass and momentum equations (Fintzi Reference Fintzi2025). The kinematics boundary conditions at the interface of the test droplet are well-defined only if the conditional average is performed on a sample of configurations where the droplets possess the same centre of mass velocity and shape. Hence, conditionally averaging on the centre of mass position x, and on the centre of mass velocity w, is a necessary step to reduce the degree of freedom of the conditionally averaged problem.

In the following, we consider arbitrary values of density and viscosity ratios, defined as

(2.12) \begin{align} \lambda = \mu _p/\mu _{\!f}, && \zeta =\rho _p /\rho _{\!f}, \end{align}

respectively. Additionally, we limit this theoretical investigation to dilute suspensions ( $\phi \ll 1$ ). We defined the droplet Reynolds number as

(2.13) \begin{equation} \textit{Re} = \frac {U\rho _{\!f} a}{\mu _{\!f}}, \end{equation}

where $a$ is the radius of the droplets, and $U = |\boldsymbol{w}_r[\boldsymbol{y}]|$ , with $\boldsymbol{w}_r[\boldsymbol{x}] = \boldsymbol{U}[\boldsymbol{x}] - \boldsymbol{w}$ , the magnitude of the relative velocity between the centre of mass velocity of the test droplet and the ensemble averaged velocity of the mixture evaluated at the centre of mass position y. Note that the value of the Reynolds number is considered arbitrary up to § 4.2 where we consider a small but finite Reynolds number. Finally, the capillary number,

(2.14) \begin{equation} \textit{Ca} = \mu _{\!f} U/\gamma , \end{equation}

is assumed to be small enough so that the droplets remain spherical.

2.3. Conditionally averaged equations

At the leading order in droplet volume fraction, the governing equations for $\boldsymbol{u}^1$ and $p_{\!f}^1$ are equivalent to those of an isolated droplet immersed in an unbounded domain (Hinch Reference Hinch1977; Koch Reference Koch1993). Hence, we consider the problem of an isolated test droplet immersed in an arbitrary flow where only the uniform relative motions induce inertia. The disturbance pressure and velocity fields ( $\boldsymbol{u}^1$ and $p_{\!f}^1$ ), relative to the position of a test droplet are denoted $\boldsymbol{u}_{o}$ , $\boldsymbol{u}_{i}$ , $p_{o}$ and $p_{i}$ , for the velocity outside, and inside, the pressure outside, and inside the volume of the test droplet, respectively. Figure 1 displays a schematic representation of the problem.

Figure 1.

Sketch of the problem. Here $\boldsymbol{u}_{o/i}[\boldsymbol{x}|\boldsymbol{y},\boldsymbol{w}]$ is the disturbance velocity field evaluated at x , generated by a test droplet positioned at y with centre of mass velocity w . Here n is the surface normal (outward) vector of the droplet. The exterior domain of the test droplet is denoted $\varOmega _o$ , and the interior domain $\varOmega _i$ .

Outside and inside the volume of the test droplet we can write the mass and momentum equations for the disturbance fields ( $p_{o/i}$ and $\boldsymbol{u}_{o/i}$ ) by conditionally averaging the Navier–Stokes equations (Koch Reference Koch1993; Fintzi Reference Fintzi2025). This leads to the equations for the disturbance fields originally derived by Maxey & Riley (Reference Maxey and Riley1983) and Gatignol (Reference Gatignol1983). However, it is important to emphasize that – even in dilute suspensions – conditional averaging introduces additional contributions in (2.16) and (2.18) (see Koch (Reference Koch1993) equation (9), or Fintzi (Reference Fintzi2025)), which are not present in Maxey & Riley (Reference Maxey and Riley1983) and Gatignol (Reference Gatignol1983)’s equations. In particular, one extra term corresponds to the velocity-variance of the conditionally averaged velocity field (see Koch (Reference Koch1993) or chapter 4 in Fintzi (Reference Fintzi2025)). In dilute flow, this term turns out to be negligible (Koch Reference Koch1993, Appendix A). Consequently, in the present context, we consider that the conditionally averaged equations are equivalent to the equations for an isolated droplet. The only difference is that conditionally averaged equations are explicitly related to ensemble-averaged fields (U and $P_{\!f}$ ) through their boundary conditions (Fintzi Reference Fintzi2025). In dimensionless form, they read

(2.15) \begin{align} \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{u}_{o} &= 0, \\[-12pt]\nonumber \end{align}

(2.16) \begin{align} \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol \sigma _{\!o} &= \textit{Re} [ \partial _t \boldsymbol{u}_{o} + \boldsymbol{u}_{o}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{o} + \boldsymbol{u}_{o}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{w}_r + \boldsymbol{w}_r \boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{o}] \nonumber \\ &= \textit{Re}\! \boldsymbol{f}_{\!o}, \\[-12pt]\nonumber \end{align}

for all $|\boldsymbol{x}- \boldsymbol{y}| = |\boldsymbol{r}| = r \gt 1$ , and

(2.17) \begin{align} \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{u}_{i} &= 0, \\[-12pt]\nonumber \end{align}

(2.18) \begin{align} \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol \sigma _{i} &= (\zeta /\lambda -1) \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol \varSigma + \frac {\zeta \textit{Re}}{\lambda } [ \partial _t \boldsymbol{u}_{i} + \boldsymbol{u}_{i}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{i} + \boldsymbol{u}_{i}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{w}_r + \boldsymbol{w}_r \boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{i}]\nonumber \\ &= (\zeta /\lambda -1) \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol \varSigma +\frac {\zeta }{\lambda }{\textit{Re}}\boldsymbol{f}_{i} = \boldsymbol{f}_{i}^{\textit{tot}} , \\[-12pt]\nonumber \end{align}

for $r\lt 1$ . The terms $\boldsymbol{f}_{o/i}$ denote dimensionless effective body force incorporating the inertial terms from the Navier–Stokes equations. Additionally, $\boldsymbol \sigma _{i/o} = -p_{i/o}\boldsymbol \delta + 2\boldsymbol{e}_{i/o}$ with $\boldsymbol{e}_{i/o} = ( {1}/{2})[\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{i/o} + ^{\dagger }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{i/o}]$ is the dimensionless Newtonian stresses in the interior and exterior of the test droplet. Note that the distances have been made dimensionless using the radius $a$ , and the velocities using the velocity scale $U$ . Likewise, the stresses in the exterior of the test droplet $\boldsymbol \sigma _o$ , and the ensemble averaged stress $\boldsymbol \varSigma$ , have been made dimensionless using the viscous scale $U \! \mu _{\!f} /a$ , and the stress in the interior of the test droplet $\boldsymbol \sigma _i$ , using the scale $U \mu _d /a$ . Recall that $\boldsymbol \sigma _{i}$ is a disturbance field defined with respect to the mean velocity $\boldsymbol{U}$ , and the mean stress $\boldsymbol \varSigma$ . Hence, because of the difference in viscosity and density between the dispersed and continuous phases (i.e. between $\boldsymbol \varSigma$ and $\boldsymbol \sigma _i$ ), $\boldsymbol{f}_{i}^{\textit{tot}}$ includes in addition the divergence of the ensemble-averaged stress times the ratio of the material properties: $(\zeta /\lambda - 1)$ . Note that in (2.16) and (2.18), and in the following, $\boldsymbol{w}_r$ , $\boldsymbol{U}$ and $\boldsymbol{E}$ are dimensionless. In particular, they correspond to $\boldsymbol{w}_r/U$ , $\boldsymbol{U}/U$ and $\boldsymbol{E}a/U$ , respectively.

At the surface of the droplet ( $r = 1$ ) the continuity of velocity and the non-deformability of the test droplet ( ${\textit{Ca}} \ll 1$ ), and the imposed velocity of the test droplet centre of mass, imposes

(2.19) \begin{align} \boldsymbol{u}_{i} - \boldsymbol{u}_{o}= 0, && (\boldsymbol{u}_{i}+\boldsymbol{w}_r) \boldsymbol{\cdot }\boldsymbol{n} = 0. \end{align}

Here $\boldsymbol{u}_{i}+\boldsymbol{w}_r$ can be interpreted as the velocity observed in the reference frame moving with the test droplet. The correct boundary condition for the disturbance shear rate is (at $r=1$ )

(2.20) \begin{equation} \boldsymbol{n}\boldsymbol{\cdot }2 [\boldsymbol{e}_{o} - \lambda \boldsymbol{e}_{i} + (1-\lambda )\boldsymbol{E} ]\boldsymbol{\cdot }(\boldsymbol \delta - \boldsymbol{nn}) = \boldsymbol{b}\boldsymbol{\cdot }(\boldsymbol \delta - \boldsymbol{nn}), \end{equation}

with, $\boldsymbol{b}[\boldsymbol{r},t]$ a tangential stress jump at the interface normalized by $\mu _{\!f} U / a$ . Here we have considered a stress jump in order to stay general when deriving the reciprocal relation (see § 3), nevertheless, $\boldsymbol{b}$ will be taken to be zero in the final calculation of the moments of forces § 4. In particular, the contribution of Marangoni stresses to the rheology of the droplet suspension is not included in the present analysis (see Fintzi & Pierson (Reference Fintzi and Pierson2026b ) for further discussion). Note that in general the ensemble averaged quantities $\boldsymbol{U}$ and $\boldsymbol{E}$ as well as $\boldsymbol{b}$ depend on the position along the surface of the test droplet and time, so that $\boldsymbol{w}_r[\boldsymbol{r},t]$ , $\boldsymbol{E}[\boldsymbol{r},t]$ and $\boldsymbol{b}[\boldsymbol{r},t]$ in (2.20) and (2.19) are not spatially constant. Hence, it is important to note that although $\boldsymbol{w}_r[\boldsymbol{y}]$ is a unit vector, $\boldsymbol{w}_r[\boldsymbol{x}]$ is not, because in the general case $\boldsymbol{U}[\boldsymbol{x}] \neq \boldsymbol{U}[\boldsymbol{y}]$ .

For sufficiently large $|\boldsymbol{y}-\boldsymbol{x}|$ , it is assumed that the conditionally averaged velocity and pressure fields are not affected by the presence of the test droplet, hence they are equal to their respective ensemble averaged values. Hence, far from the test droplet the disturbance fields, which are defined as the difference between conditionally and ensemble averaged fields, vanish. It reads

(2.21) \begin{align} \lim _{r \to \infty }\boldsymbol{u}_{o}[\boldsymbol{r}|\boldsymbol{y},\boldsymbol{w}] = 0, && \lim _{r \to \infty }p_{o}[\boldsymbol{r}|\boldsymbol{y},\boldsymbol{w}]= 0. \end{align}

3.1. Auxiliary solutions

The reciprocal theorem requires the use of a known solution. We call that solution the auxiliary solution, since its only purpose is to compute the real solution of the equations introduced above. The auxiliary velocity and stress fields are denoted by $\hat {\boldsymbol{u}}_{i/o}$ and $\hat {\boldsymbol \sigma }_{i/o} = - \hat {p}_{i/o}\boldsymbol \delta + \boldsymbol{\boldsymbol{\nabla }} \hat {\boldsymbol{u}}_{i/o} + ^\dagger \boldsymbol{\boldsymbol{\nabla }} \hat {\boldsymbol{u}}_{i/o}$ , respectively. In general, all auxiliary quantities are written with a hat.

In the auxiliary problem, we neglect all inertial effects and consider an arbitrary quadratic flow; thus $\hat {p}_{o/i}$ and $\hat {\boldsymbol{u}}_{o/i}$ satisfy

(3.1) \begin{align} \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\hat {\boldsymbol{u}}_{o} = 0 , &&\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\hat {\boldsymbol \sigma }_{o} = 0 , \\[-12pt]\nonumber \end{align}

(3.2) \begin{align} \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\hat {\boldsymbol{u}}_{i} = 0 , && \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\hat {\boldsymbol \sigma }_{i} = 0 , \end{align}

for $r\gt 1$ and $r\lt 1$ , respectively. The boundary conditions at $r\to \infty$ read

(3.3) \begin{align} \lim _{r \to \infty }\hat {\boldsymbol{u}}_{o}[\boldsymbol{r}|\boldsymbol{y},\boldsymbol{w}] = 0, && \lim _{r \to \infty }\hat {p}_{o}[\boldsymbol{r}|\boldsymbol{y},\boldsymbol{w}]= 0, \end{align}

and at $r=1$ ,

(3.4) \begin{align} \hat {\boldsymbol{u}}_{i} &= \hat {\boldsymbol{u}}_{o}, \\[-12pt]\nonumber \end{align}

(3.5) \begin{align} \hat {\boldsymbol{u}}_{i} \boldsymbol{\cdot }\boldsymbol{n} &= - \hat {\boldsymbol{w}}_r \boldsymbol{\cdot }\boldsymbol{n}, \\[-12pt]\nonumber \end{align}

(3.6) \begin{align} \boldsymbol{n}\boldsymbol{\cdot }2[\hat {\boldsymbol{e}}_{o} - \lambda \hat {\boldsymbol{e}}_{i} + (1-\lambda ) \hat {\boldsymbol{E}} ]\boldsymbol{\cdot }(\boldsymbol \delta - \boldsymbol{nn}) &= \hat {\boldsymbol{b}}\boldsymbol{\cdot }(\boldsymbol \delta - \boldsymbol{nn}), \end{align}

where we introduced $\hat {\boldsymbol{w}}_r = \hat {\boldsymbol{U}} - \boldsymbol{w}$ as the undisturbed relative velocity field of the auxiliary problem, similarly $\hat {\boldsymbol{E}}$ is the undisturbed shear rate of the auxiliary problem. Note that for the test problem it is important that $\hat {\boldsymbol{b}}\neq 0$ , this will be useful to write a reciprocal relation in § 3. Indeed, this stress discontinuity is a mathematical artefact in order to evaluate the surface integral of the shear stress at the interface of the test droplet (Raja et al. Reference Raja, Subramanian and Koch2010) but has no physical meaning here. For our purposes, we only need $\hat {\boldsymbol{b}}$ and $\hat {\boldsymbol{U}}$ to be quadratic fields. Hence, we consider the relations

(3.7) \begin{align} \hat {\boldsymbol{w}}_r(\boldsymbol{y} + \boldsymbol{r}) &= \hat {\boldsymbol{U}} - \boldsymbol{w} + \boldsymbol{r} \boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }}\hat {\boldsymbol{U}}|_{\boldsymbol{r}=0} + \frac {1}{2}\boldsymbol{rr} : \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\boldsymbol{\nabla }}\hat {\boldsymbol{U}}|_{\boldsymbol{r}=0} + \ldots , \\[-12pt]\nonumber \end{align}

(3.8) \begin{align} \hat {\boldsymbol{E}}(\boldsymbol{y} + \boldsymbol{r}) &= \hat {\boldsymbol{E}}|_{\boldsymbol{r}=0} + \boldsymbol{r} \boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \hat {\boldsymbol{E}}|_{\boldsymbol{r}=0}+ \ldots , \\[-12pt]\nonumber \end{align}

(3.9) \begin{align} \hat {\boldsymbol{b}}(\boldsymbol{y} + \boldsymbol{r}) &= \hat {\boldsymbol{b}}|_{\boldsymbol{r}=0} + \boldsymbol{r} \boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \hat {\boldsymbol{b}}|_{\boldsymbol{r}=0} + \frac {1}{2}\boldsymbol{rr} : \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\boldsymbol{\nabla }} \hat {\boldsymbol{b}}|_{\boldsymbol{r}=0}+ \ldots . \end{align}

According to the linearity of the Stokes equations, (3.1) and (3.2), and of the boundary conditions, (3.4) to (3.6), we deduce that $\hat {\boldsymbol{u}}_{i/o}$ and $\hat {p}_{i/o}$ must be linear combination of spherical harmonics proportional to $\hat {\boldsymbol{w}}_r|_{\boldsymbol{r}=0}$ , $\hat {\boldsymbol{b}}|_{\boldsymbol{r}=0}$ and their derivatives (Brenner Reference Brenner1964). Hence, we can write the general solution of $\boldsymbol{u}_{i/o}$ and $p_{i/o}$ as

(3.10) \begin{align} \begin{pmatrix} \hat {\boldsymbol{u}}_{o}\\[7pt] \hat {p}_{o}\\[7pt] \hat {\boldsymbol{u}}_{i}\\[7pt] \hat {p}_{i} \end{pmatrix} = \begin{pmatrix} \mathbb{U}_{o}^{(1)} + \mathbb{U}_{o}^{(2)}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} + \mathbb{U}_{o}^{(3)} :\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\boldsymbol{\nabla }} & \mathbb{U}_{o}^{\text{(1-b)}} + \mathbb{U}_{o}^{\text{(2-b)}}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} + \mathbb{U}_{o}^{\text{(3-b)}} :\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\boldsymbol{\nabla }} \\[7pt] \mathbb{P}_{o}^{(1)} + \mathbb{P}_{o}^{(2)}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} + \mathbb{P}_{o}^{(3)} :\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\boldsymbol{\nabla }} & \mathbb{P}_{o}^{\text{(1-b)}} + \mathbb{P}_{o}^{\text{(2-b)}}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} + \mathbb{P}_{o}^{\text{(3-b)}} :\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\boldsymbol{\nabla }} \\[7pt] \mathbb{U}_{i}^{(1)} + \mathbb{U}_{i}^{(2)}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} + \mathbb{U}_{i}^{(3)} :\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\boldsymbol{\nabla }} & \mathbb{U}_{i}^{\text{(1-b)}} + \mathbb{U}_{i}^{\text{(2-b)}}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} + \mathbb{U}_{i}^{\text{(3-b)}} :\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\boldsymbol{\nabla }} \\[7pt] \mathbb{P}_{i}^{(1)} + \mathbb{P}_{i}^{(2)}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} + \mathbb{P}_{i}^{(3)} :\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\boldsymbol{\nabla }} & \mathbb{P}_{i}^{\text{(1-b)}} + \mathbb{P}_{i}^{\text{(2-b)}}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} + \mathbb{P}_{i}^{\text{(3-b)}} :\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\boldsymbol{\nabla }} \\ \end{pmatrix} \boldsymbol{\cdot }\begin{pmatrix} \hat {\boldsymbol{w}}_r\\[7pt] \hat {\boldsymbol{b}} \end{pmatrix}\!. \end{align}

The tensor $\mathbb{U}_{i/o}^{(n)}$ and $\mathbb{P}_{i/o}^{(n)}$ are $(n+1)$ - and $n$ -order tensors given in Leal (Reference Leal2007) and in Appendix F of Fintzi & Pierson (Reference Fintzi and Pierson2026b ), that depend only on the relative coordinate $\boldsymbol{r} = \boldsymbol{x}-\boldsymbol{y}$ and the viscosity ratio $\lambda$ . We also define the $(n+2)$ -order tensors,

(3.11) \begin{equation} \mathbb{S}_{i/o}^{(n)} = - \boldsymbol \delta \mathbb{P}_{i/o}^{(n)} + \boldsymbol{\boldsymbol{\nabla }} \mathbb{U}_{i/o}^{(n)} + ^\dagger \boldsymbol{\boldsymbol{\nabla }} \mathbb{U}_{i/o}^{(n)}, \end{equation}

where $n =1,2,3$ or $1\text{-}b,2\text{-}b,3\text{-}b$ . (The $^\dagger$ used in (3.11) represents the transpose operator which acts on the two closest indices. For example, the transpose of the triadic $\boldsymbol{abc}$ may be written $^\dagger \boldsymbol{abc} = \boldsymbol{bac}$ or $\boldsymbol{abc}^\dagger = \boldsymbol{acb}$ .) Equation (3.11) provides the stresses fields of the auxiliary problem upon contraction with $\hat {\boldsymbol{w}}_r$ , and $\hat {\boldsymbol{b}}$ , and the higher derivatives.

3.2. Singularity solutions

Following Stone et al. (Reference Stone, Brady and Lovalenti2001) and Raja et al. (Reference Raja, Subramanian and Koch2010), we also consider the auxiliary problem of a point source located at the origin. The reason why this solution is necessary is that, in the previous auxiliary problem $\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\hat {\boldsymbol{U}}= 0$ . Hence, as demonstrated below the previous auxiliary solutions are unable to provide a formula for the trace of the first moment using the reciprocal theorem. Because we also aim to compute the traces of the second moment of forces, one also needs the solutions of a point force and a potential dipole. These two classes of solutions, i.e. point source and point force (and their derivatives), obey the non-homogeneous Stokes equations

(3.12) \begin{align} {{\nabla} }^2 \hat {\boldsymbol{u}}_o = \boldsymbol{\boldsymbol{\nabla }} \hat {p}_o + (\boldsymbol{Q}^{(n)}\odot {{\nabla} }^{(n)})\boldsymbol{\boldsymbol{\nabla }} \delta (\boldsymbol{r}), && {{\nabla} }^2 \hat {\boldsymbol{u}}_o = \boldsymbol{\boldsymbol{\nabla }} \hat {p}_o + (\boldsymbol{R}^{(n+1)}\odot{{\nabla} }^{(n)}) \delta (\boldsymbol{r}), \end{align}

respectively. Where $\boldsymbol{Q}^{(n)}$ and $\boldsymbol{R}^{(n+1)}$ are two arbitrary nth and $(n+1)$ th order tensors, respectively. Similarly, the ${{\nabla} }^{(n)}$ denote nth order outer product of gradients vectors. The operator $\odot$ represents the contraction product on the nth common indices. Both equations are completed by the continuity equation $\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{u}_{o}$ . The solution of the point source and point force read

(3.13) \begin{align} \hat {\boldsymbol{u}}_o = -\frac {1}{4\pi }{{\nabla} }^{(n+1)}(1/r)\odot \boldsymbol{Q}^{(n)}, && \hat {\boldsymbol{u}}_o = \frac {1}{8\pi }{{\nabla} }^{(n)}[{{\nabla} }^{(2)} - \boldsymbol \delta {{\nabla} }^2]r \odot \boldsymbol{R}^{(n+1)}, \end{align}

respectively. Note that $( {1}/{8\pi })[{{\nabla} }^{(2)} - \boldsymbol \delta {{\nabla} }^2]r = - (\boldsymbol{nn} + \boldsymbol \delta )/(8\pi r)$ is the free-space Green function of Stokes equations (Pozrikidis et al. Reference Pozrikidis1992). The stress tensors read

(3.14) \begin{align} \hat {\boldsymbol \sigma }_o &= -\frac {1}{2\pi }{{\nabla} }^{(n+2)}(1/r)\odot \boldsymbol{Q}^{(n)}, \\[-12pt]\nonumber \end{align}

(3.15) \begin{align} \hat {\boldsymbol \sigma }_o &= \frac {1}{8\pi }{{\nabla} }^{(n)}[ 2{{\nabla} }^{(3)} -({\boldsymbol{\nabla }} \boldsymbol \delta + ^\dagger \boldsymbol{\boldsymbol{\nabla }}\boldsymbol \delta +\boldsymbol \delta \boldsymbol{\boldsymbol{\nabla }}){{\nabla} }^2 ]r\odot \boldsymbol{R}^{(n+1)}, \end{align}

for the point source and point forces solutions, respectively. More details on how to obtain those solutions are given in Chapter 6 of Pozrikidis (Reference Pozrikidis2011) or in Appendix A.

As an example, the solution of a point source (i.e. (3.13) (first equation) with $n=0$ ) used in the calculation of the first moment of force reads

(3.16) \begin{align} \hat {\boldsymbol{u}}_{o} = \frac {Q^{(0)}}{4\pi } \boldsymbol{n}r^{-2}, && \hat {\boldsymbol \sigma }_{o} = \frac {Q^{(0)}}{2\pi }\left ( \boldsymbol \delta - 3 \boldsymbol{nn} \right )r^{-3}. \end{align}

Note that this expression is valid throughout the domain excluding the point $\boldsymbol{r} = \boldsymbol{x} - \boldsymbol{y} = 0$ , hence we may use either the subscript $o$ or $i$ . The solution for a point source dipole $(\boldsymbol{Q}^{(1)}$ ) and point force $(\boldsymbol{R}^{(1)})$ will also be used to derive relations for the traces of the second moment of force. In summary, as the values of $\hat {\boldsymbol{w}}_r$ , $\hat {\boldsymbol{b}}$ , $\boldsymbol{R}^{(n)}$ and $\boldsymbol{Q}^{(n)}$ , are entirely arbitrary, the solutions given by (3.10), (3.13) and (3.16) can be used as a tool to derive expressions for the first and second moments of the forces by means of the reciprocal theorem derived in the next section.

3.3. Reciprocal theorem

The derivation of the general formula proceeds in three main steps: (i) we write a reciprocal relation in the exterior of the test droplet; (ii) we then write the corresponding relation in the interior of the test droplet and (iii) finally, using the boundary conditions at the droplet interface, we combine the two expressions to obtain the final reciprocal relation.

We first take the dot product of (2.16) with $\hat {\boldsymbol{u}}_{o}$ , and the dot product of (3.1) with $\boldsymbol{u}_{o}$ , subtracting both expressions gives

(3.17) \begin{equation} \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }(\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\hat {\boldsymbol{u}}_o) = \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }(\hat {\boldsymbol \sigma }_{o}\boldsymbol{\cdot }\boldsymbol{u}_o) + \textit{Re} (\hat {\boldsymbol{u}}_{o}\boldsymbol{\cdot }\boldsymbol{f}_{\!o}). \end{equation}

To derive this relation we used the identities $\boldsymbol \sigma _o :\boldsymbol{\boldsymbol{\nabla }} \hat {\boldsymbol{u}}_o = 2\boldsymbol{e}_o :\boldsymbol{\boldsymbol{\nabla }} \hat {\boldsymbol{u}}_o = 2{\boldsymbol e}_o : \hat {\boldsymbol e}_o$ , and $\hat {\boldsymbol \sigma }_o :\boldsymbol{\boldsymbol{\nabla }} {\boldsymbol{u}}_o = 2\hat {\boldsymbol{e}}_o : \boldsymbol{e}_o$ , which imply $\boldsymbol \sigma _o :\boldsymbol{\boldsymbol{\nabla }} \hat {\boldsymbol{u}}_o = \hat {\boldsymbol \sigma }_o :\boldsymbol{\boldsymbol{\nabla }} {\boldsymbol{u}}_o$ . This simplification explains why complete knowledge of the velocity field in the exterior of the droplet becomes unnecessary when $\boldsymbol{f}_{\!o} = 0$ . Consequently, the power of the reciprocal theorem is rooted in the commutativity of the operation $\textbf e_o:\hat {\textbf e}_o$ . Integrating this expression on the whole domain $\varOmega _{o}$ , and using the divergence theorem we obtain

(3.18) \begin{equation} - \oint _{\varGamma _{p}} \hat {\boldsymbol{u}}_{o} \boldsymbol{\cdot }\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma = -\oint _{\varGamma _{p}} \boldsymbol{u}_{o}\boldsymbol{\cdot }\hat {\boldsymbol \sigma }_{o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma + \textit{Re}\int _{\varOmega _{o}} \hat {\boldsymbol{u}}_{o}\boldsymbol{\cdot }\boldsymbol{f}_{\!o} {\rm d}\varOmega . \end{equation}

The surface integrals appear only on the boundary of the test droplet because far from the test droplet ( $r\to \infty$ ), the fields $\hat {\boldsymbol{u}}_o$ , $\hat {\boldsymbol \sigma }_o$ , $\boldsymbol{u}_o$ and $\boldsymbol \sigma _o$ vanish (see (2.21) and (3.3)). Equation (3.18) already provides a sufficient reciprocal relation when considering solid particles (Masoud & Stone Reference Masoud and Stone2019). Indeed, in that case $\boldsymbol{u}_{o} = -\boldsymbol{w}_r$ , and $\hat {\boldsymbol{u}}_{o} = -\hat {\boldsymbol{w}}_r$ at $r=1$ , because of the no-slip boundary condition at the particle interface, hence leading directly to a formula for the force traction on the surface of the solid particle (first term on the left-hand side of (3.19)). Of course, one still needs to provide an approximation for $\boldsymbol{f}_{\!o}$ if inertia is considered. For fluid particles, a few more steps are necessary. Inserting the relation $\boldsymbol{u}_o = \boldsymbol{u}_o+\boldsymbol{w}_r-\boldsymbol{w}_r$ and $\hat {\boldsymbol{u}}_o = \hat {\boldsymbol{u}}_o+\hat {\boldsymbol{w}}_r-\hat {\boldsymbol{w}}_r$ gives

(3.19) \begin{align} &\oint _{\varGamma _{p}} \hat {\boldsymbol{w}}_{r} \boldsymbol{\cdot }\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - 2\oint _{\varGamma _{p}} (\hat {\boldsymbol{u}}_{o} + \hat {\boldsymbol{w}}_r) \boldsymbol{\cdot }\boldsymbol e_{o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma \nonumber \\ &= \oint _{\varGamma _{p}} \boldsymbol{w}_{r}\boldsymbol{\cdot }\hat {\boldsymbol \sigma }_{o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - 2\oint _{\varGamma _{p}} (\boldsymbol{u}_{o} + \boldsymbol{w}_r)\boldsymbol{\cdot }\hat {\boldsymbol e}_{o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma + \textit{Re}\int _{\varOmega _{o}} \hat {\boldsymbol{u}}_{o}\boldsymbol{\cdot }\boldsymbol{f}_{\!o} {\rm d}\varOmega . \end{align}

Where we used the relation, $(\hat {\boldsymbol{u}}_{o} + \hat {\boldsymbol{w}}_r) \boldsymbol{\cdot }\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} = 2 (\hat {\boldsymbol{u}}_{o} + \hat {\boldsymbol{w}}_r) \boldsymbol{\cdot }\boldsymbol e_{o}\boldsymbol{\cdot }\boldsymbol{n}$ , allowed because $(\hat {\boldsymbol{u}}_{o} + \hat {\boldsymbol{w}}_r) = (\boldsymbol \delta - \boldsymbol{nn})\boldsymbol{\cdot }(\hat {\boldsymbol{u}}_{o} + \hat {\boldsymbol{w}}_r)$ on the surface of the test droplet, see (3.5). Similar considerations lead to $({\boldsymbol{u}}_{o} + {\boldsymbol{w}_r}) \boldsymbol{\cdot }\hat {\boldsymbol \sigma }_{o}\boldsymbol{\cdot }\boldsymbol{n} =2 ({\boldsymbol{u}}_{o} + {\boldsymbol{w}_r}) \boldsymbol{\cdot }\hat {\boldsymbol e}_{o}\boldsymbol{\cdot }\boldsymbol{n}$ .

The equations governing the flow inside the droplets have not yet been considered. We now take the dot product of (2.18) with $(\hat {\boldsymbol{u}}_i + \hat {\boldsymbol{w}}_r)$ and of (3.2) with $(\boldsymbol{u}_i+ \boldsymbol{w}_r)$ , subtracting both expression leads to

(3.20) \begin{equation} \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }[\boldsymbol \sigma _{i}\boldsymbol{\cdot }(\hat {\boldsymbol{u}}_i+\hat {\boldsymbol{w}}_r)] - 2\boldsymbol{e}_i : \hat {\boldsymbol{E}} = \boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }[\hat {\boldsymbol \sigma }_{i}\boldsymbol{\cdot }(\boldsymbol{u}_i+\boldsymbol{w}_r)] - 2\hat {\boldsymbol{e}}_i :\boldsymbol{E} + (\hat {\boldsymbol{u}}_i+\hat {\boldsymbol{w}}_r)\boldsymbol{\cdot }\boldsymbol{f}_{i}^{\textit{tot}}, \end{equation}

where we used the relation, $- \boldsymbol \sigma _i :\boldsymbol{\boldsymbol{\nabla }} (\hat {\boldsymbol{u}}_i+\hat {\boldsymbol{w}}_r)+ \hat {\boldsymbol \sigma }_i :\boldsymbol{\boldsymbol{\nabla }} ({\boldsymbol{u}}_i+{\boldsymbol{w}}_r) = - 2\boldsymbol{e}_i :\hat {\boldsymbol{E}} + 2\hat {\boldsymbol{e}}_i :\boldsymbol{E}$ . Integrating this relation over the domain $\varOmega _i$ , and using the divergence theorem, as well as similar simplifications used in (3.19), gives

(3.21) \begin{align} &\oint _{\varGamma _{p}} (\hat {\boldsymbol{u}}_i+\hat {\boldsymbol{w}}_r)\boldsymbol{\cdot }\boldsymbol{e}_{i}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \int _{\varOmega _{i}} \boldsymbol{e}_i : \hat {\boldsymbol{E}} {\rm d}\varOmega \nonumber \\ &= \oint _{\varGamma _{p}} (\boldsymbol{u}_i+\boldsymbol{w}_r)\boldsymbol{\cdot }\hat {\boldsymbol{e}}_{i}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \int _{\varOmega _{i}} \hat {\boldsymbol{e}}_i :\boldsymbol{E} {\rm d}\varOmega + \frac {1}{2} \int _{\varOmega _{i}} (\hat {\boldsymbol{u}}_i+\hat {\boldsymbol{w}}_r)\boldsymbol{\cdot }\boldsymbol{f}_{i}^{\textit{tot}} {\rm d}\varOmega. \end{align}

The first terms on both the left- and right-hand sides of (3.21) may be obtained using (2.20) and (3.21). Indeed, since $(\hat {\boldsymbol{u}}_i+\hat {\boldsymbol{w}}_r)\boldsymbol{\cdot }(\boldsymbol \delta -\boldsymbol{nn}) = (\hat {\boldsymbol{u}}_i+\hat {\boldsymbol{w}}_r)$ , multiplying (2.20) by $(\hat {\boldsymbol{u}}_i+\hat {\boldsymbol{w}}_r)$ yields

(3.22) \begin{equation} -2 \lambda (\hat {\boldsymbol{u}}_{i/o}+\hat {\boldsymbol{w}}_r) \boldsymbol{\cdot }\boldsymbol{e}_i \boldsymbol{\cdot }\boldsymbol{n} = \boldsymbol{b}\boldsymbol{\cdot }(\hat {\boldsymbol{u}}_{i/o}+\hat {\boldsymbol{w}}_r) -2 (\hat {\boldsymbol{u}}_{i/o}+\hat {\boldsymbol{w}}_r)\boldsymbol{\cdot }\boldsymbol{e}_{o}\boldsymbol{\cdot }\boldsymbol{n} -2 (1 -\lambda ) (\hat {\boldsymbol{u}}_{i/o}+\hat {\boldsymbol{w}}_r)\boldsymbol{\cdot }\boldsymbol{E} \boldsymbol{\cdot }\boldsymbol{n}, \end{equation}

with a similar relation for the test problem boundary condition. Adding (3.19) and (3.21) times $2\lambda$ , while using the boundary condition given by (3.22) leads to the expression

(3.23) \begin{align} &\oint _{\varGamma _{p}} \hat {\boldsymbol{w}}_{r} \boldsymbol{\cdot }\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma -\lambda \int _{\varOmega _{i}} 2\boldsymbol{e}_i : \hat {\boldsymbol{E}} {\rm d}\varOmega - (1-\lambda ) \oint _{\varGamma _{p}} 2(\boldsymbol{u}_{i} + \boldsymbol{w}_r)\boldsymbol{\cdot }\hat {\boldsymbol{E}}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma + \oint _{\varGamma _{p}} (\boldsymbol{u}_{i} + \boldsymbol{w}_r)\boldsymbol{\cdot }\hat {\boldsymbol{b}} {\rm d}\varGamma \nonumber \\ &\qquad \quad = \oint _{\varGamma _{p}} \boldsymbol{w}_{r}\boldsymbol{\cdot }\hat {\boldsymbol \sigma }_{o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \lambda \int _{\varOmega _{i}} 2\hat {\boldsymbol{e}}_i :\boldsymbol{E} {\rm d}\varOmega - (1-\lambda ) \oint _{\varGamma _{p}} 2(\hat {\boldsymbol{u}}_{i} + \hat {\boldsymbol{w}}_r) \boldsymbol{\cdot }\boldsymbol{E} \boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma \nonumber \\ &\qquad \qquad \quad + \oint _{\varGamma _{p}} (\hat {\boldsymbol{u}}_{i} + \hat {\boldsymbol{w}}_r) \boldsymbol{\cdot }\boldsymbol{b} {\rm d}\varGamma + \lambda \int _{\varOmega _{i}} (\hat {\boldsymbol{u}}_i+\hat {\boldsymbol{w}}_r)\boldsymbol{\cdot }\boldsymbol{f}_{i}^{\textit{tot}} {\rm d}\varOmega + \textit{Re}\int _{\varOmega _{o}} \hat {\boldsymbol{u}}_{o}\boldsymbol{\cdot }\boldsymbol{f}_{\!o} {\rm d}\varOmega . \end{align}

To obtain the final form of the reciprocal theorem, we proceed by noting that the second and third terms on the left-hand side of (3.23) may be combined upon using the divergence theorem on the third term, and using the relation

(3.24) \begin{equation} - \lambda \boldsymbol{e}_i : \hat {\boldsymbol{E}} - (1-\lambda )\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }[\hat {\boldsymbol{E}} \boldsymbol{\cdot }(\boldsymbol{u}_i + \boldsymbol{w}_r)] = - \boldsymbol{e}_i : \hat {\boldsymbol{E}} - (1-\lambda ) \hat {\boldsymbol{E}}:\boldsymbol{E} - (1-\lambda )\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\hat {\boldsymbol{E}} \boldsymbol{\cdot }(\boldsymbol{u}_i + \boldsymbol{w}_r). \end{equation}

One can also derive the same manipulation on the second and third term on the right-hand side of (3.23), namely

(3.25) \begin{equation} - \lambda \hat {\boldsymbol{e}}_i : {\boldsymbol{E}} - (1-\lambda )\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }[{\boldsymbol{E}} \boldsymbol{\cdot }(\hat {\boldsymbol{u}}_i + \hat {\boldsymbol{w}}_r)] = - \hat {\boldsymbol{e}_i} : {\boldsymbol{E}} - (1-\lambda ) \hat {\boldsymbol{E}}:\boldsymbol{E} - (1-\lambda )\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }{\boldsymbol{E}} \boldsymbol{\cdot }(\hat {\boldsymbol{u}}_i + \hat {\boldsymbol{w}}_r). \end{equation}

Because the product $\hat {\boldsymbol{E}} : {\boldsymbol{E}}$ is commutative, the second term on the right-hand side of (3.24) and (3.25) cancel each other in the final expression. Injecting both (3.25) and (3.24) in (3.23) leads to the final form of the reciprocal theorem for droplets, namely

(3.26) \begin{align} \oint _{\varGamma _{p}} &\hat {\boldsymbol{w}}_{r} \boldsymbol{\cdot }\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \int _{\varOmega _{i}} 2\boldsymbol{e}_i : \boldsymbol{\boldsymbol{\nabla }}\hat {\boldsymbol{U}} {\rm d}\varOmega - (1-\lambda ) \int _{\varOmega _{i}} (\boldsymbol{u}_{i} + \boldsymbol{w}_r)\boldsymbol{\cdot }{{\nabla} }^2 \hat {\boldsymbol{U}} {\rm d}\varOmega + \oint _{\varGamma _{p}} (\boldsymbol{u}_{i} + \boldsymbol{w}_r)\boldsymbol{\cdot }\hat {\boldsymbol{b}} {\rm d}\varGamma \nonumber \\ &\quad = \oint _{\varGamma _{p}} \boldsymbol{w}_{r}\boldsymbol{\cdot }\hat {\boldsymbol \sigma }_{o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \int _{\varOmega _{i}} 2\hat {\boldsymbol{e}}_i :\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{U} {\rm d}\varOmega - (1-\lambda ) \int _{\varOmega _{i}} (\hat {\boldsymbol{u}}_{i} + \hat {\boldsymbol{w}}_r) \boldsymbol{\cdot }{{\nabla} }^2 \boldsymbol{U} {\rm d}\varOmega \nonumber \\ &\qquad + \oint _{\varGamma _{p}} (\hat {\boldsymbol{u}}_{i} + \hat {\boldsymbol{w}}_r) \boldsymbol{\cdot }\boldsymbol{b} {\rm d}\varGamma + \lambda \int _{\varOmega _{i}} (\hat {\boldsymbol{u}}_i+\hat {\boldsymbol{w}}_r)\boldsymbol{\cdot }\boldsymbol{f}_{\!i}^{\textit{tot}} {\rm d}\varOmega + \textit{Re}\int _{\varOmega _{o}} \hat {\boldsymbol{u}}_{o}\boldsymbol{\cdot }\boldsymbol{f}_{\!o} {\rm d}\varOmega . \end{align}

Note that we used the relations, $2{\boldsymbol{e}_i} : \hat {\boldsymbol{E}} = 2{\boldsymbol{e}_i} : \hat {\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{U}}$ and $\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{E} = ( {1}/{2}){{\nabla} }^2 \boldsymbol{U}$ because the background flow is divergence free equation (2.2). Upon making the good choice for $\hat {\boldsymbol{w}}_r$ and $\hat {\boldsymbol{b}}$ (i.e. whether it is constant, linear, or quadratic, see (3.7)), one recovers on the left-hand side of (3.26) the expression of the force, first and second moments applied on the droplet in a yet arbitrary flow described by $\boldsymbol{w}_r$ and $\boldsymbol{b}$ , and their gradients.

In (3.26) recall that $\boldsymbol{n}$ is the outward normal to the test-droplet surface ( $\varGamma _p$ ), while $\varOmega _o$ and $\varOmega _i$ represent the domains outside and inside the test droplet, respectively. On the right-hand side of (3.26) all the terms are either known quantities from the auxiliary problem ( $\hat {\boldsymbol \sigma }_{o},\hat {\boldsymbol{e}}_i, \hat {\boldsymbol{u}}_{i/o/r}$ ), or boundary conditions of the ‘real’ problem ( $\boldsymbol{U},\boldsymbol{w}_r$ and $\boldsymbol{b}$ ). However, one exception remains, that is the inertial term $\boldsymbol{f}_{\!o}$ (respectively, $\boldsymbol{f}_{\!i}^{\textit{tot}}$ ), which is a function of the unknown velocity field $\boldsymbol{u}_{o}$ (respectively, $\boldsymbol{u}_i$ ). Hence, to compute the right-hand side of (3.26), we must use some approximations allowing us to compute the last two terms of (3.26) without the knowledge of the complete solution. This is the subject of the next section.

4.1. Formulae for the force, first and second moments

Let us consider the configurations of the auxiliary problem listed in table 1.

Table 1.

Configurations of the auxiliary problem. Configuration 1 is used to compute the force, configurations 2 and 3 to compute the first moment of the force, configurations 4–6 to compute the second moment of the force and configuration 7 to compute the droplet internal shear rate.

Table 1. Long description

A table with three columns and seven rows. The columns are labeled Configurations, w_r_hat(x), b(x), and Point source / force solutions. Each row lists different configurations and their corresponding values. Row 1: Configuration 1, w_r_hat[r] = w_r_hat, b[r] = 0, Q^(n) = 0, R^(n) = 0. Row 2: Configuration 2, w_r_hat[r] = r · ∇U, b[r] = 0, Q^(n) = 0, R^(n) = 0. Row 3: Configuration 3, w_r_hat[r] = 0, b[r] = 0, Q^(0). Row 4: Configuration 4, w_r_hat[r] = 1/2 rr · ∇∇U, b[r] = 0, Q^(n) = 0, R^(n) = 0. Row 5: Configuration 5, w_r_hat[r] = 0, b[r] = 0, Q^(1). Row 6: Configuration 6, w_r_hat[r] = 0, b[r] = 0, R^(1). Row 7: Configuration 7, w_r_hat[r] = 0, b[r] = b, Q^(n) = 0, R^(n) = 0.

For any of configurations 1 to 7, the expressions for $\hat {\boldsymbol{\sigma }}_{i/o}$ and $\hat {\boldsymbol{u}}_{i/o}$ required in (3.26) can be obtained from (3.10), (3.11) and (3.16).

Substituting configuration 1 into (3.26) and factoring out $\hat {\boldsymbol{w}}_r$ yields

(4.1) \begin{align} \oint _{\varGamma _{p}} \boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma &= \oint _{\varGamma _{p}} \boldsymbol{w}_{r}\boldsymbol{\cdot }\mathbb{S}_{o}^{(1)}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \int _{\varOmega _{i}} \mathbb{S}_{i}^{(1)} :\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{U} {\rm d}\varOmega - (1-\lambda ) \int _{\varOmega _{i}} (\mathbb{U}_{i}^{(1)} + \boldsymbol \delta ) \boldsymbol{\cdot }{{\nabla} }^2 \boldsymbol{U} {\rm d}\varOmega \nonumber \\ &\quad + \oint _{\varGamma _{p}} (\mathbb{U}_{i/o}^{(1)} + \boldsymbol \delta ) \boldsymbol{\cdot }\boldsymbol{b} {\rm d}\varGamma + \lambda \int _{\varOmega _{i}} (\mathbb{U}_{i}^{(1)} + \boldsymbol \delta )\boldsymbol{\cdot }\boldsymbol{f}^{\textit{tot}}_i {\rm d}\varOmega + \textit{Re}\int _{\varOmega _{o}} \mathbb{U}_{o}^{(1)}\boldsymbol{\cdot }\boldsymbol{f}_{\!o} {\rm d}\varOmega , \end{align}

which is a general formula for the hydrodynamic drag forces. Based on this formula, one can rederive Faxen law, for example (Stone et al. Reference Stone, Brady and Lovalenti2001; Pozrikidis Reference Pozrikidis2011).

The general formula for the first moment of force may be obtained by setting configuration 2 into (3.26), it yields

(4.2) \begin{align} &\oint _{\varGamma _{p}} r_j(\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n})_i {\rm d}\varGamma - \int _{\varOmega _{i}} 2(\boldsymbol{e}_i)_{ij} {\rm d}\varOmega \overset {i\neq j}{=} \oint _{\varGamma _{p}} (\boldsymbol{w}_{r})_l (\mathbb{S}^{(2)}_{o})_{ijkl} n_k {\rm d}\varGamma - \int _{\varOmega _{i}} (\mathbb{S}^{(2)}_{i})_{ijkl}(\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{U})_{kl}{\rm d}\varOmega \nonumber \\ &\qquad \qquad + (\lambda -1) \int _{\varOmega _{i}} (\mathbb{U}_{i}^{(2)} + \boldsymbol{r}\boldsymbol \delta )_{\textit{ijk}} ({{\nabla} }^2 \boldsymbol{U})_k {\rm d}\varOmega + \oint _{\varGamma _{p}} (\mathbb{U}_{i}^{(2)} + \boldsymbol{r}\boldsymbol \delta )_{\textit{ijk}} \text{b}_k {\rm d}\varGamma \nonumber \\ &\qquad \qquad + \lambda \int _{\varOmega _{i}} (\mathbb{U}_{i}^{(2)} + \boldsymbol{r}\boldsymbol \delta )_{\textit{ijk}} (\boldsymbol{f}^{\textit{tot}}_i)_k {\rm d}\varOmega + \textit{Re}\int _{\varOmega _{o}} (\mathbb{U}_{o}^{(2)})_{\textit{ijk}}(\boldsymbol{f}_{\!o})_k {\rm d}\varOmega . \end{align}

It is important to note that to derive (4.2) we have factorized the equations by the tensor $(\boldsymbol{\boldsymbol{\nabla }} \hat {\boldsymbol{U}})_{ij}$ . Because $(\boldsymbol{\nabla }\hat {\boldsymbol{U}})_{kk} = \boldsymbol{\nabla }\boldsymbol{\cdot }\hat {\boldsymbol{U}} = 0$ , taking the trace of (4.2) leads to an incorrect expression, since it would mean factoring by zero. Therefore, it is important to note that this equality holds only for $i \neq j$ . Hence, if we note $M_{ij}$ the whole first moment, then (4.2) only provides its traceless part, namely $M_{ij} - 1/3\delta _{ij}M_{kk}$ . Nevertheless, using configuration 3 in (3.18), we obtain a formula for the trace of the first moment ( $M_{kk}$ ), which reads

(4.3) \begin{equation} \oint _{\varGamma _{p}} \boldsymbol{r} \boldsymbol{\cdot }\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n}{\rm d}\varGamma = - \textit{Re}\int _{\varOmega _{o}} \frac {\boldsymbol{r}}{r^3}\boldsymbol{\cdot }\boldsymbol{f}_{\!o} {\rm d}\varOmega . \end{equation}

Indeed, the integral of $\boldsymbol{u}_o \boldsymbol{\cdot }\hat {\boldsymbol{\sigma }}_o$ over the droplet surface on the right-hand side of (3.17) vanishes because $\boldsymbol{u}_o \boldsymbol{\cdot }\hat {\boldsymbol{\sigma }}_o \propto \boldsymbol{u}_o \boldsymbol{\cdot }\boldsymbol{n}$ at $r=1$ , and $\boldsymbol{u}_o$ is divergence-free. By adding $\delta _{ij}$ times (4.3) to (4.2), one obtains an expression valid for every pair of indices $i$ and $j$ , that is, the complete first-moment tensor.

Using configuration 4 (see table 1), we obtain the general relation for the second moment of the force (in index notation),

(4.4) \begin{align} &\frac {1}{2}\oint _{\varGamma _{p}} (\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n})_i r_jr_k {\rm d}\varGamma - \int _{\varOmega _{i}} 2(\boldsymbol{e}_i )_{ij} r_k {\rm d}\varOmega \overset {i\neq j,k}{=} \oint _{\varGamma _{p}} (\boldsymbol{w}_{r})_l (\mathbb{S}^{(3)}_{o})_{ijklm} n_m {\rm d}\varGamma - \int _{\varOmega _{i}} (\mathbb{S}^{(3)}_{i})_{ijklm} (\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{U})_{lm} {\rm d}\varOmega \nonumber \\ &\qquad \qquad - (1-\lambda ) \int _{\varOmega _{i}} \left((\mathbb{U}^{(3)}_{i})_{ijkl} + \frac {1}{2}r_kr_j\delta _{il}\right)({{\nabla} }^2 \boldsymbol{U})_l {\rm d}\varOmega + \oint _{\varGamma _{p}} \left((\mathbb{U}^{(3)}_{i})_{ijkl} + \frac {1}{2}r_kr_j\delta _{il}\right) \text{b}_l {\rm d}\varGamma \nonumber \\ &\qquad \qquad + \lambda \int _{\varOmega _{i}} \left((\mathbb{U}^{(3)}_{i})_{ijkl} + \frac {1}{2}r_kr_j\delta _{il}\right) (\boldsymbol{f}_{i}^{\textit{tot}})_l {\rm d}\varOmega + \textit{Re}\int _{\varOmega _{o}} (\mathbb{U}^{(3)}_{o})_{ijkl}\boldsymbol{\cdot }\boldsymbol{f}_{\!o} {\rm d}\varOmega . \end{align}

Here we used the relation $\int _{\varOmega _{i}} (\boldsymbol{u}_i+\boldsymbol{w}_r) {\rm d}\varOmega =v_p(\boldsymbol{w}- \boldsymbol{w}) = 0$ . This formula is valid for a droplet immersed in an otherwise arbitrary background flow with velocity $\boldsymbol{U}$ and an interfacial jump $\boldsymbol{b}$ . However, it is only valid for $i \neq j,k$ . Indeed, we have factored out $(\boldsymbol{\nabla }\boldsymbol{\nabla }\hat {\boldsymbol{U}})_{kji}$ , and $(\boldsymbol{\nabla }\boldsymbol{\nabla }\hat {\boldsymbol{U}})_{kii} = (\boldsymbol{\nabla }\boldsymbol{\nabla }\hat {\boldsymbol{U}})_{iji} = 0$ . Consequently, if we denote by $K_{\textit{ijk}}$ the complete second moment, (4.4) provides only the traceless part of $K_{\textit{ijk}}$ under any contraction over $i,j$ or $i,k$ . Therefore, a separate expression is still required for the trace of (4.4), together with a procedure to combine this trace with the previously undefined traceless part (with respect to $i,j$ and $i,k$ ) of $K_{\textit{ijk}}$ . This procedure is detailed in (Appendix B.1) and yields

(4.5) \begin{equation} G_{\textit{ijk}} = K_{\textit{ijk}} + \frac {1}{8} (K_{\textit{lkl}} - 3K_{\textit{llk}})\delta _{ij} + \frac {1}{8} (K_{\textit{llj}} -3 K_{\textit{ljl}})\delta _{ik}. \end{equation}

One can verify that taking the trace of this expression over $(i,k)$ or $(i,j)$ yields zero. To obtain the whole second moment from $G_{\textit{ijk}}$ (i.e. from (4.4)) we use the relation,

(4.6) \begin{equation} K_{\textit{ijk}} = G_{\textit{ijk}} + \frac {1}{8}(3K_{\textit{llk}}-K_{\textit{lkl}}) \delta _{ij} + \frac {1}{8}(3K_{\textit{ljl}} - K_{\textit{llj}}) \delta _{ik}, \end{equation}

where $\boldsymbol{G}$ is defined as the traceless part of $\boldsymbol{K}$ under contractions over the index pairs $ij$ and $ik$ . To obtain the traces of the second moment, i.e. $K_{\textit{lkl}}$ and $K_{\textit{llk}}$ , we can use a procedure similar to the one used for (4.3). For the moment, we delay this procedure to the next section, where the real problem will be defined in greater detail.

Recall that $\boldsymbol \sigma _o$ and $\boldsymbol{e}_i$ are the local stress and shear rate, relative to the bulk stress ( $\boldsymbol \varSigma$ ), and bulk shear rate ( $\boldsymbol{E}$ ), respectively. Consequently, the left-hand sides of (4.1) to (4.4) are exactly the terms we seek in (2.8) to (2.10).

Although these formulae are sufficient to close the averaged momentum equations, one may be interested in the expression for the integral of $\boldsymbol{e}_i$ over the droplet volume, independently of the integral of $\boldsymbol{r}\boldsymbol \sigma _o\boldsymbol{\cdot }\boldsymbol{n}$ , in contrast to (4.2) which provides the sum of these two terms. Indeed, this term can be used for the calculation of the droplet deformation, see for example Fintzi & Pierson (Reference Fintzi and Pierson2026b ). Using configuration 7 of table 1 we obtain

(4.7) \begin{align} \oint _{\varGamma _{p}} (\boldsymbol{u}_{i} + \boldsymbol{w}_r)\boldsymbol{r} {\rm d}\varGamma & = \oint _{\varGamma _{p}} \boldsymbol{w}_{r}\boldsymbol{\cdot }\mathbb{S}_{o}^{{(2-b)}}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \int _{\varOmega _{i}} \mathbb{S}_{i}^{{(2-b)}} :\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{U} {\rm d}\varOmega + (\lambda -1) \int _{\varOmega _{i}} \mathbb{U}_{i}^{\text{(2-b)}} \boldsymbol{\cdot }{{\nabla} }^2 \boldsymbol{U} {\rm d}\varOmega \nonumber \\ &\quad + \oint _{\varGamma _{p}} \mathbb{U}_{i}^{{(2-b)}} \boldsymbol{\cdot }\boldsymbol{b} {\rm d}\varGamma + \lambda \int _{\varOmega _{i}} \mathbb{U}_{i}^{{(2-b)}}\boldsymbol{\cdot }\boldsymbol{f}_{i}^{\textit{tot}} {\rm d}\varOmega + \textit{Re}\int _{\varOmega _{o}} \mathbb{U}_{o}^{{(2-b)}}\boldsymbol{\cdot }\boldsymbol{f}_{\!o} {\rm d}\varOmega . \end{align}

Upon using the divergence theorem on the left-hand side term, one can obtain a formula for the integral of the velocity gradient inside the droplet, whose symmetric part corresponds to the integral of $\boldsymbol{e}_i$ over the droplet volume. Hence, the left-hand side of (4.7) corresponds rigorously to the integral of $\boldsymbol{e}_i$ over the droplet volume.

4.2. Force moments on a droplet immersed in a steady uniform flow at low Reynolds number

This subsection considers only a uniform steady state relative motion in the real problem, such that $\boldsymbol{w}_r= \boldsymbol{U}-\boldsymbol{w}$ is not a function of space and time. Furthermore, we assume clean interfaces, such that $\boldsymbol{b}=\boldsymbol{0}$ .

4.2.1. Oseen drag force

In this section, we show how to obtain the $O(\textit{Re})$ correction to the drag force on a translating spherical droplet using the reciprocal theorem. Although this result has been known for a long time (Taylor & Acrivos Reference Taylor and Acrivos1964), the present derivation allows us to introduce useful notation and to clarify the methodology employed in the subsequent sections. Using (4.1), we obtain

(4.8) \begin{equation} \oint _{\varGamma _{p}} \boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma =\boldsymbol{w}_{r}\boldsymbol{\cdot }\oint _{\varGamma _{p}} \mathbb{S}_{o}^{(1)}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma + \zeta \textit{Re} \int _{\varOmega _{i}} (\mathbb{U}_{i}^{(1)} + \boldsymbol \delta )\boldsymbol{\cdot }\boldsymbol{f}_{i} {\rm d}\varOmega + \textit{Re}\int _{\varOmega _{o}} \mathbb{U}_{o}^{(1)}\boldsymbol{\cdot }\boldsymbol{f}_{\!o} {\rm d}\varOmega . \end{equation}

The term involving the divergence of mean stress cancelled out because the latter is a constant vector in this configuration and $\int _{\varOmega _{i}} \boldsymbol{w}_r +\boldsymbol{u}_i {\rm d}\varOmega =0$ . The first term on the right-hand side is by definition the Stokes flow contribution to the drag, since $\boldsymbol{w}_{r}\boldsymbol{\cdot }\mathbb{S}_{o}^{(1)}$ represents the stress field of a translating spherical droplet in Stokes flow conditions (3.10). Thus,

(4.9) \begin{equation} - \boldsymbol{w}_{r}\boldsymbol{\cdot }\oint _{\varGamma _{p}} \mathbb{S}_{o}^{(1)}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma = - 2 \pi \frac {2+3\lambda }{\lambda +1} \boldsymbol{w}_r = \boldsymbol{F}_{s}, \end{equation}

where we introduced the vector $\boldsymbol{F}_{\!s}$ which represents the Stokes force from the droplet on the continuous phase.

The real challenge is to compute the last two integrals on the right-hand side of (4.8), which require expressions for the forcing terms $\boldsymbol{f}_{\!o}$ and $\boldsymbol{f}_{\!i}$ . In a situation where only steady and uniform relative motions are present, we deduce from (2.16) to (2.18) that

(4.10) \begin{equation} \boldsymbol{f}_{i/o} = (\boldsymbol{u}_{i/o} + \boldsymbol{w}_r)\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{i/o}, \end{equation}

where we recall that $\boldsymbol{u}_{i/o}$ are the yet unknown velocity fields inside or outside the translating droplet. To obtain the $O(\textit{Re})$ correction to the drag force, one only needs to obtain a uniformly valid approximation of $\boldsymbol{f}_{i/o}$ accurate at $O(1)$ in $\textit{Re}$ . It is well known that the fields $\boldsymbol{u}_{i/o}$ , can be approximated by $\boldsymbol{u}_{i/o} \approx \mathbb{U}^{(1)}_{i/o}\boldsymbol{\cdot }\boldsymbol{w}_r$ only at distance $r \lt O(\textit{Re}^{-1})$ (Kaplun & Lagerstrom Reference Kaplun and Lagerstrom1957; Proudman & Pearson Reference Proudman and Pearson1957). Otherwise, for $r \gt O(\textit{Re}^{-1})$ , one may use the velocity field generated by a point force in the Oseen equation, i.e. the solution of (Pozrikidis Reference Pozrikidis2011)

(4.11) \begin{align} -\boldsymbol{\boldsymbol{\nabla }}\! p_{\textit{out}} +{{\nabla} }^2 \boldsymbol{u}_{\textit{out}} + \boldsymbol{F}_{\!s}\delta (\boldsymbol{r}) &= \textit{Re} \boldsymbol{w}_r \boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\textit{out}}, \end{align}

where $p_{\textit{out}}$ and $\boldsymbol{u}_{\textit{out}}$ are the outer pressure and velocity fields. This equation may be solved in the sense of generalized functions (Pozrikidis (Reference Pozrikidis2011), chapter 6), or with the use of Fourier transforms (Candelier & Mehlig Reference Candelier and Mehlig2016), it gives $\boldsymbol{u}_{\textit{out}} = \mathbb{U}^{(\textit{out})}\boldsymbol{\cdot }\boldsymbol{F}_{\!s}$ , with

(4.12) \begin{equation} \mathbb{U}^{(\textit{out})} = \frac {e^X}{4\pi r}\boldsymbol \delta + \boldsymbol{\boldsymbol{\nabla }} \left ( \frac {e^X-1}{8\pi X} (\boldsymbol{w}_r - \boldsymbol{n}) \right )\!, \end{equation}

where $X = ({r \textit{Re}}/{2}) (\boldsymbol{w}_r \boldsymbol{\cdot }\boldsymbol{n} - 1)$ . Recall that $\boldsymbol{w}_r$ is the dimensionless relative velocity vector, hence, it is a unit vector. A uniformly valid solution for the disturbance velocity field outside the test droplet and for arbitrary $r$ may be obtained by combining (4.12) and (3.10). It yields

(4.13) \begin{equation} \boldsymbol{u}_o = \boldsymbol{u}_{\textit{in}} + \textit{Re} \boldsymbol{u}_{\textit{out}}^*, \end{equation}

where $\boldsymbol{u}_{\textit{in}} = \mathbb{U}^{(1)}\boldsymbol{\cdot }\boldsymbol{w}_r$ , $\textit{Re} \boldsymbol{u}_{\textit{out}}^* = [\mathbb{U}^{(\textit{out})} - \mathbb{G}]\boldsymbol{\cdot }\boldsymbol{F}_{\!s}$ and $\mathbb{G} = (\boldsymbol \delta + \boldsymbol{nn})/(8\pi r)$ denotes the free-space Green’s function of Stokes flows. Note that expanding (4.12) and using the expression of $\mathbb{G}$ one directly gets

(4.14) \begin{align} \boldsymbol{u}_{\textit{out}}^* & = \left \{ \left [\frac {e^X -1}{8\pi X} \boldsymbol \delta + \frac {e^X-X-1}{16 \pi X^2} (\boldsymbol{nn} -\boldsymbol \delta )\right ](\boldsymbol{w}_r\boldsymbol{\cdot }\boldsymbol{n} - 1)\right .\nonumber \\[3pt] &\quad \quad \left . + \frac {(X-1) e^X + 1}{X^2 16 \pi } (\boldsymbol{w}_r - \boldsymbol{n}) (\boldsymbol{w}_r - \boldsymbol{n}) \right \}\boldsymbol{\cdot }\boldsymbol{F}_{\!s}. \end{align}

From (4.14), one observes that in the limit $X\to 0$ , or $r\to 0$ at $\textit{Re}$ fixed, or equivalently $\textit{Re} \to 0$ at fixed $r$ , the quantity $\boldsymbol{u}_{\textit{out}}^*$ remains bounded at $O(1)$ . For $r \to \infty$ or $X\to -\infty$ (since $\boldsymbol{w}_r\boldsymbol{\cdot }\boldsymbol{n}-1 \leq 0$ ), this function behaves at most as $\sim 1/({\textit{Re}}r)$ depending on the values of $\boldsymbol{w}_r\boldsymbol{\cdot }\boldsymbol{n}$ . Using (4.14) and (4.10), we write $\boldsymbol{f}_{\!o}$ as

(4.15) \begin{equation} \boldsymbol{f}_{\!o} = (\boldsymbol{u}_{\textit{in}}+\boldsymbol{w}_r)\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\textit{in}} + \textit{Re} \boldsymbol{w}_r\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\textit{out}}^\ast + \ldots \end{equation}

where the ellipsis denotes the remaining contributions such as $\textit{Re} \boldsymbol{u}_{\textit{in}}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\textit{out}}$ and $\textit{Re}^2 \boldsymbol{u}_{\textit{out}}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\textit{out}}$ , which will be shown to be of order $o(\textit{Re})$ in the final result.

The last term on the right-hand side of (4.1) can now be evaluated as

(4.16) \begin{align} \int _{\varOmega _{o}} \mathbb{U}_{o}^{(1)}\boldsymbol{\cdot }\boldsymbol{f}_{\!o} {\rm d}\varOmega &= \int _{1\lt |\boldsymbol{r}|\lt \infty }{ \mathbb{U}_{o}^{(1)}\boldsymbol{\cdot }(\boldsymbol{u}_{\textit{in}} + \boldsymbol{w}_r)\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\textit{in}} }d^3\boldsymbol{r} \nonumber \\ &+ \textit{Re} \int _{1\lt |\boldsymbol{r}|\lt \infty }{ \mathbb{U}_{o}^{(1)}\boldsymbol{\cdot }\boldsymbol{w}_r \boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\textit{out}}^\ast }d^3\boldsymbol{r}. \end{align}

Although all these integrals could be computed directly, it is useful to proceed to some simplifications. Because $\mathbb{U}_o^{(1)}$ and therefore $\boldsymbol{u}_{\textit{in}}$ are even functions of $\boldsymbol{n}$ , and consequently $\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\textit{in}}$ an odd function of n, the product, $\mathbb{U}_{o}^{(1)}\boldsymbol{\cdot }(\boldsymbol{u}_{\textit{in}} + \boldsymbol \delta )\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\textit{in}}$ is odd in n. Thus, this term vanishes upon integration over the spherical unit surface centred at the origin. Hence, the inner solution does not contribute to the $O(\textit{Re})$ correction to the drag force. The integral of $ (\mathbb{U}_i^{(1)}+\boldsymbol \delta )\boldsymbol{\cdot }\boldsymbol{f}_{\!i}$ over $\varOmega _i$ , in (4.1), vanishes for similar reasons.

The second integral of (4.16) represents the contribution of the outer solution. To evaluate the scalings with respect to $\textit{Re}$ of the integrals involving $\boldsymbol{u}_{\textit{out}}^\ast$ , one must consider separately the inner region where $\boldsymbol{u}_{\textit{out}}^\ast \sim O(1)$ and the outer region where $\boldsymbol{u}_{\textit{out}}^\ast \sim O(1/{\textit{Re}}r)$ . For example, the term in the second integral of (4.16) scales as $\textit{Re} O(r^{-2}) r^2 {\rm d}r$ in the inner region, yielding a contribution of $O(1)$ , after integration from $r=1$ to $r=\textit{Re}^{-1}$ , which is of course not negligible. However, if we only consider the part of $\mathbb{U}_o^{(1)}$ which decays as $\sim r^{-3}$ , then this contribution scales as $\textit{Re} O(r^{-4}) r^2 {\rm d}r= \textit{Re} r^{-2}$ , which yields a contribution of $O(\textit{Re})$ after integration. In the outer region now, one sees that the term which scales in $O(r^{-3})$ in $\mathbb{U}_o^{(1)}$ contributes to $\textit{Re} O(\textit{Re}^{-1} r^{-5}) r^2 {\rm d}r= r^{-3}$ , which yields a contribution of $O(\textit{Re}^2)$ after integration between $\textit{Re}^{-1}$ and $r\to \infty$ . Consequently, in this integral, one may use the approximation

(4.17) \begin{equation} \mathbb{U}_o^{(1)} =-|\boldsymbol{F}_{\!s}| \boldsymbol{\cdot }\mathbb{G} + O(r^{-3}), \end{equation}

which will therefore yield an error of $O(\textit{Re}^2)$ in the final result. Using the same reasoning we see that the term $\textit{Re} \mathbb{U}_o^{(1)} \boldsymbol{\cdot }\boldsymbol{u}_{\textit{in}}\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\textit{out}}^\ast r^2 {\rm d}r$ scales as $\sim \textit{Re} r^{-1}$ in the inner region, which integrates to $O({\textit{Re}}\ln \textit{Re})$ , and as $\sim r^{-2}$ in the outer region which integrates to $O(\textit{Re})$ , so this term and the others remain negligible. Using these approximations, we only need to evaluate the second integral of (4.16) to obtain the drag force, namely

(4.18) \begin{align} \textit{Re} \boldsymbol{F}_{\!s} \boldsymbol{\cdot }\int _{1\lt r\lt \infty }{ \boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}^*_{(out)}\boldsymbol{\cdot }\mathbb{G} } d^3 \boldsymbol{r} = \boldsymbol{w}_r\frac {(\boldsymbol{F}_{\!s}\boldsymbol{\cdot }\boldsymbol{F}_{\!s})}{16\pi } + O(\textit{Re}). \end{align}

The integral in (4.18) can be evaluated using the Fourier convolution theorem (see Masoud & Stone Reference Masoud and Stone2019), or by direct integration in spherical coordinates, which is the method adopted here, see Appendix D. Finally, by using (4.18) and (4.9) in (4.8) we obtain the drag force on the droplet as

(4.19) \begin{equation} \oint _{\varGamma _{p}} \boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma = - \boldsymbol{F}_{\!s}+ \frac {{\textit{Re}}}{16\pi }(\boldsymbol{F}_{\!s}\boldsymbol{\cdot }\boldsymbol{F}_{\!s})\boldsymbol{w}_r = 2\pi \frac {2+3\lambda }{\lambda +1}\boldsymbol{w}_r+ \pi \frac {(2+3\lambda )^2}{(\lambda +1)^2} \frac {{\textit{Re}}}{4} \boldsymbol{w}_r. \end{equation}

We recovered the classic results derived by Taylor & Acrivos (Reference Taylor and Acrivos1964).

4.2.2. First moment of force

The first moment may be computed directly from (4.1) and (4.3). The former equation reads

(4.20) \begin{align} \oint _{\varGamma _{p}} r_j(\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n})_i {\rm d}\varGamma - \int _{\varOmega _{i}} 2(\boldsymbol{e}_i)_{ij} {\rm d}\varOmega &\overset {i\neq j}{=} \zeta \textit{Re} \int _{\varOmega _{i}} (\mathbb{U}_{i}^{(2)} + \boldsymbol{r}\boldsymbol \delta )_{\textit{ijk}} (\boldsymbol{f}_{i})_k {\rm d}\varOmega + \textit{Re}\int _{\varOmega _{o}} (\mathbb{U}_{o}^{(2)})_{\textit{ijk}}(\boldsymbol{f}_{\!o})_k {\rm d}\varOmega , \end{align}

where we noticed that $\oint _{\varGamma _{p}} \mathbb{S}^{(2)}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma = 0$ , whence the right-hand side of (4.3) and (4.20) only contains inertial contributions. This is because the stresslet on a translating droplet in pure Stokes flow is zero (Kim & Karrila Reference Kim and Karrila2013).

We now focus on computing the last inertial term on the right-hand side of (4.20), which involves the forcing term $\boldsymbol{f}_{\!o}$ , still given by (4.10). It is found that the outer contribution to this integral is negligible at leading order in $\textit{Re}$ , and in this case the $O(\textit{Re})$ inertial effects can be obtained through a regular perturbation procedure. Indeed, for $r\lt O(\textit{Re}^{-1})$ , the term $\textit{Re} \boldsymbol{w}_r\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\textit{out}}^\ast$ scales as $\sim \textit{Re} r^{-1}$ , while $\mathbb{U}_o^{(2)} \sim r^{-2}$ . Consequently, the integrand associated with the stresslet scales as $O({\textit{Re}}r^{-3})\times r^2$ , which when integrated from $1$ to $\textit{Re}^{-1}$ yields a term of $O({\textit{Re}}\ln \textit{Re})$ . Likewise, in the outer region, $\textit{Re} \mathbb{U}_o^{(2)} \boldsymbol{\cdot }\boldsymbol{w}_r\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \boldsymbol{u}_{\textit{out}}^\ast r^2 {\rm d}r$ scales as $O(r^{-2})$ therefore after integration from $\textit{Re}^{-1}$ to $\infty$ , the outer contribution to the integral is of $O(\textit{Re})$ . The convergence of the volume integrals indicates the regular nature of the inertial contribution (Dabade et al. Reference Dabade, Marath and Subramanian2015). Similar arguments are presented in Stone et al. (Reference Stone, Brady and Lovalenti2001) and Raja et al. (Reference Raja, Subramanian and Koch2010) when computing the inertial correction to the stresslet in shear-dominated flows. Since the droplet lies entirely within the Stokes region, the Stokes flow field may be used to compute $\boldsymbol{f}_{\!i}$ , and consequently the inertial terms involving $\boldsymbol{f}_{\!i}$ in (4.20).

Based on this argument and (3.10) and (4.10), one may rewrite the forcing terms as

(4.21) \begin{equation} \boldsymbol{f}_{i/o} = \boldsymbol{w}_r \boldsymbol{\cdot }(\mathbb{U}_{i/o}^{(1)} +\boldsymbol \delta )\boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }} \mathbb{U}_{i/o}^{(1)}\boldsymbol{\cdot }\boldsymbol{w}_r, \end{equation}

and compute the last integrals on the right-hand side of (4.20) and (4.3), yielding

(4.22) \begin{align} \oint _{\varGamma _{p}} \boldsymbol{r}\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \int _{\varOmega _{i}} 2\boldsymbol{e}_i {\rm d}\varOmega &\overset {i\neq j}{=} -\pi \textit{Re} \frac {63 \lambda ^{3} + 150 \lambda ^{2} + 112 \lambda + 28}{60 (\lambda +1)^{3}} \boldsymbol{w}_r \boldsymbol{w}_r, \\[-12pt]\nonumber \end{align}

(4.23) \begin{align} \oint _{\varGamma _{p}} \boldsymbol{r}\boldsymbol{\cdot }\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma &= \pi \textit{Re} \frac {3 \lambda ^2+6\lambda + 4}{12 (\lambda +1)^2} \boldsymbol{w}_r\boldsymbol{\cdot }\boldsymbol{w}_r. \end{align}

Note that the integrals only involve tensor products of the position vector $\boldsymbol{r}$ with the distance function $r = |\boldsymbol{r}|$ . We have used the open-source library sympy of Python to compute these integrals and the subsequent ones in this work. Specifically, the integration over the unit sphere is done using the methodology described in Barthès-Biesel & Acrivos (Reference Barthès-Biesel and Acrivos1973) and Pozrikidis et al. (Reference Pozrikidis1992).

To obtain a formula for the whole first moment, we sum the deviatoric part of (4.22) with (4.23) times $\delta _{ij}/3$ . It gives

(4.24) \begin{align} \oint _{\varGamma _{p}} \boldsymbol{r}\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \int _{\varOmega _{i}} 2\boldsymbol{e}_i {\rm d}\varOmega = -\pi \textit{Re} \frac {63 \lambda ^{3} + 150 \lambda ^{2} + 112 \lambda + 28}{60 (\lambda +1)^{3}} \boldsymbol{w}_r \boldsymbol{w}_r\nonumber \\[3pt] + \pi \textit{Re} \frac {26 \lambda ^{3} + 65 \lambda ^2+54\lambda + 16}{60 (\lambda +1)^3} (\boldsymbol{w}_r\boldsymbol{\cdot }\boldsymbol{w}_r) \boldsymbol \delta . \end{align}

Additionally, using (4.7) and (3.10) and similar arguments for the estimation of $\boldsymbol{f}_{i/o}$ one finds

(4.25) \begin{equation} \oint _{\varGamma _{p}} \boldsymbol{u}_{i} \boldsymbol{r} {\rm d}\varGamma = \int _{\varOmega _{i}} \boldsymbol{e}_{i} {\rm d}\varOmega = -\pi \textit{Re} \frac {12\lambda ^2+23\lambda +10}{100(\lambda +1)^3}\left[\boldsymbol{w}_r\boldsymbol{w}_r-\frac {1}{3} (\boldsymbol{w}_r\boldsymbol{\cdot }\boldsymbol{w}_r )\boldsymbol \delta \right] \end{equation}

which may be used to compute the droplet deformation in that particular situation following the procedure outlined in Fintzi & Pierson (Reference Fintzi and Pierson2026b ).

4.2.3. Second moment of force

Now, let us turn our attention to the second moment of forces. Since $\boldsymbol{w}_r$ is uniform and $\boldsymbol{b}=0$ , (4.4) reads

(4.26) \begin{align} & \frac {1}{2}\oint _{\varGamma _{p}} r_kr_j ( \boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n})_i {\rm d}\varGamma - 2\int _{\varOmega _{i}} r_k (\boldsymbol{e}_i)_{ji} d\varOmega \overset {i\neq j,k}{=} (\boldsymbol{w}_{r})_l\oint _{\varGamma _{p}} (\mathbb{S}^{(3)}_{o})_{ijklm}n_m {\rm d}\varGamma \nonumber \\ &\qquad \qquad + \zeta \textit{Re} \int _{\varOmega _{i}} \left( (\mathbb{U}_{i}^{(3)})_{ijkl} + \frac {1}{2} r_kr_j\delta _{il} \right) (\boldsymbol{f}_{i})_l {\rm d}\varOmega + \textit{Re}\int _{\varOmega _{o}} (\mathbb{U}_{o}^{(3)})_{ijkl} (\boldsymbol{f}_{\!o})_l {\rm d}\varOmega . \end{align}

This formula gives the second moment of forces (left-hand side of (4.26)) for a droplet translating in a uniform flow. As noted, previously we have factored out by $(\boldsymbol{\boldsymbol{\nabla }}\boldsymbol{\boldsymbol{\nabla }}\hat {\boldsymbol{U}})_{kji}$ , to derive this formula, hence (4.26) is only valid when $i\neq j,k$ . To obtain the traces of the second moment, i.e. $K_{\textit{lkl}}$ and $K_{\textit{llk}}$ , we use a procedure similar to the one used for (4.3) (see Appendix B.2). Indeed, based on the singularity solutions given by (3.13) we may show that

(4.27) \begin{align} \frac {1}{2}\oint _{\varGamma _{p}} \boldsymbol{nn} \boldsymbol{\cdot }\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma = - \frac {1}{2}\oint _{\varGamma _{p}} \boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma + 3\boldsymbol{w}_r \int _{\varOmega _{i}} {\rm d}\varOmega - \frac {\textit{Re}}{2}\int _{\varOmega _{o}} (\boldsymbol \delta + \boldsymbol{nn})r^{-1}\boldsymbol{\cdot }\boldsymbol{f}_{\!o} {\rm d}\varOmega , \\[-12pt]\nonumber \end{align}

(4.28) \begin{align} \frac {1}{2}\oint _{\varGamma _{p}} \boldsymbol{nn} \boldsymbol{\cdot }\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \int _{\varOmega_{i}} 2 \boldsymbol{r}\boldsymbol{\cdot }\boldsymbol{e}_i {\rm d}\varOmega = \frac {1}{6}\oint _{\varGamma _{p}} \boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma + \frac {\textit{Re}}{6}\int _{\varOmega _{o}} (\boldsymbol \delta - 3\boldsymbol{nn})r^{-3}\boldsymbol{\cdot }\boldsymbol{f}_{\!o} d\varOmega . \end{align}

It is interesting to note that the drag force term appears explicitly on the right-hand side of these expressions. Here, it should be understood that (4.27) corresponds to $\mathbb{K}_{\textit{llk}}$ , while (4.28) corresponds to $\mathbb{K}_{\textit{lkl}}$ .

Whether, it is (4.26) and (4.27), or (4.28), the vectors $\boldsymbol{f}_{i/o}$ are still needed, and given by (4.10), with the approximation given by $\boldsymbol{u}_{i/o}= \boldsymbol{w}_r\boldsymbol{\cdot }\mathbb{U}_{i/o}^{(1)}$ close from the test droplet, and $\boldsymbol{u}_{\textit{out}}= \boldsymbol{F}_{\!s}\boldsymbol{\cdot }\mathbb{U}^{(\textit{out})}$ in the Oseen region. Additionally, similarly than with (4.17) we introduce the approximation (Nadim & Stone Reference Nadim and Stone1991)

(4.29) \begin{equation} (\mathbb{U}^{(3)}_o)_{ijkl} = \pi \frac {\lambda }{(\lambda +1)}(\mathbb{G})_{l i}\delta _{jk} + O(r^{-3}). \end{equation}

Hence, the last integrals on the right-hand side of (4.26) to (4.28) will be evaluated using the same approach as for (4.8) and (4.20). Using the same symmetry arguments as in the two previous sections, we arrive at the conclusion that in (4.26) to (4.28) all the internal contributions integrated over the volume of the test droplet vanish except the second term of (4.27). By making use of the same symmetry argument as in the drag force calculation, and of (4.29), we deduce that in (4.26) and (4.27), only the far field contribution is relevant. The result of the integral presented in (4.18) is used to compute these two contributions. In the last integral on the right-hand side of (4.28), only the inner velocity field could in principle contribute, due to the rapid $r^{-3}$ decay of the integrand. However, the contribution from the inner velocity field vanishes upon integration, since the final expression is odd in $\boldsymbol{n}$ . Overall, the methodology for computing the integrals requires nothing beyond what has already been presented in the previous two subsections. This yields

(4.30) \begin{align} \frac {1}{2}\oint _{\varGamma _{p}} r_kr_j ( \boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n})_i {\rm d}\varGamma - 2\int _{\varOmega _{i}} r_k (\boldsymbol{e}_i)_{ji} {\rm d}\varOmega &\overset {i\neq j,k}{=} \pi \frac {\lambda }{\lambda +1}(\boldsymbol{w}_r)_i \delta _{jk} + \textit{Re} \pi \frac {\lambda (2+3\lambda )}{8(\lambda +1)^2} \delta _{jk} (\boldsymbol{w}_r)_i, \end{align}

(4.31) \begin{align} &\hspace{6.8pc} \frac {1}{2}\oint _{\varGamma _{p}} \boldsymbol{nn} \boldsymbol{\cdot }\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma = \pi \frac {\lambda +2}{\lambda +1}\boldsymbol{w}_r + \textit{Re} \pi \frac {(2+\lambda )(2+3\lambda )}{8(\lambda +1)^2}\boldsymbol{w}_r, \end{align}

(4.32) \begin{align} \frac {1}{2}\oint _{\varGamma _{p}} \boldsymbol{nn} \boldsymbol{\cdot }\boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n} {\rm d}\varGamma - \int _{\varOmega _{i}} 2 \boldsymbol{r}\boldsymbol{\cdot }\boldsymbol{e}_i {\rm d}\varOmega &= \pi \frac {3\lambda +2}{3( \lambda +1) } \boldsymbol{w}_{r} + Re\pi \frac {(3\lambda +2)^2}{24( \lambda +1)^2 } \boldsymbol{w}_{r}. \end{align}

Combining all contributions and using (4.5) and (4.6) yields the expression for the complete second moment of the forces, namely

(4.33) \begin{align} \frac {1}{2}\oint _{\varGamma _{p}} r_kr_j ( \boldsymbol \sigma _{\!o}\boldsymbol{\cdot }\boldsymbol{n})_i {\rm d}\varGamma - 2\int _{\varOmega _{i}} r_k (\boldsymbol{e}_i)_{ji} {\rm d}\varOmega & = \frac {\pi }{\lambda +1} \left[ \frac {2}{3}\delta _{ij} (\boldsymbol{w}_r)_k + \lambda \delta _{jk} (\boldsymbol{w}_r)_i \right]\nonumber \\ &\quad + \frac {\pi \textit{Re}(3\lambda +2)}{8(\lambda +1)^2} \left[ \frac {2}{3}\delta _{ij} (\boldsymbol{w}_r)_k + \lambda \delta _{jk} (\boldsymbol{w}_r)_i \right]\!. \end{align}

One may note that the Stokes-flow contribution derived in Zhang & Prosperetti (Reference Zhang and Prosperetti1997) and Fintzi & Pierson (Reference Fintzi and Pierson2026b ) is recovered in the first line of (4.33).

In Appendix C we derived the $O(\textit{Re})$ inner pressure and velocity fields around a solid sphere in the same configuration, and derived the force moments from these fields, thereby validating the present methodology since we recover the same expressions for the force moments, at least in the limit $\lambda \to \infty$ .

The formulae given by (4.24) and (4.33) constitute the main results of the paper.

Figure 2.

Plots of the coefficients of the first moment of force divided by $\textit{Re}$ as a function of the viscosity ratio $\lambda$ . Here solid line is the coefficient in front of the isotropic part $C_{\textit{Iso}} = (3\lambda ^2 +6\lambda +4)/(\lambda +1)^2 /12$ (see (4.23)); the dashed line is the coefficient of the deviatoric part of the first moment $C_{{Dev}} = - (63\lambda ^3 +150\lambda ^2 +112\lambda +28)/ (\lambda +1)^3 /60$ (see (4.22)).

We have plotted in figure 2 the coefficients of the first moment as functions of the viscosity ratio. In this figure, $C_{\textit{Iso}}$ denotes the coefficient associated with the isotropic contribution to the first moment equation (4.23), while $C_{{Dev}}$ corresponds to the coefficient of the deviatoric contribution equation (4.22). One observes that $C_{\textit{Iso}}$ exhibits only a weak dependence on $\lambda$ , whereas $C_{{Dev}}$ is significantly larger for bubbles ( $\lambda =0$ ) than for solid particles ( $\lambda \to \infty$ ). Notably, $C_{{Dev}}$ is closely related to droplet deformation (Fintzi Reference Fintzi2025), and it is shown in Fintzi & Pierson (Reference Fintzi and Pierson2026a ) how the classical results of Taylor & Acrivos (Reference Taylor and Acrivos1964) for the inertial deformation of a translating droplet can be recovered. One may note the absence of $\zeta$ in these expressions, indicating that the droplet density has no effect on the first and second moments, at least for $\textit{Re} \ll 1$ , ${\textit{Ca}} \ll 1$ and in the quasi-steady-state configuration. This arises because, upon direct calculation, the first integral on the right-hand side of (4.20) vanishes. Moreover, as shown previously, the inertial contributions integrated over the volume of the test droplet also vanish in the expressions for the second moment; consequently, the droplet density has no effect on the second moment. A similar observation was made by Magnaudet et al. (Reference Magnaudet, Takagi and Legendre2003), who showed that the inertial lift force acting on a droplet translating parallel to a wall does not depend on the droplet density. However, Taylor & Acrivos (Reference Taylor and Acrivos1964) demonstrated that the density ratio enters the force balance through droplet deformation. Therefore, one may expect the density ratio to influence the force moments for a slightly deformed droplet.