Drag on a spheroidal particle at clean and surfactant-laden interfaces: effects of particle aspect ratio, contact angle and surface viscosities

Abstract

Translation of a non-spherical particle trapped at a membrane or at a complex interface between fluids is a relevant situation occurring in biological, technological and everyday life systems. Here, we consider prolate spheroidal particles, translating at a clean or (insoluble) surfactant-laden planar interface located between a viscous fluid and air, protruding into the surrounding subphase. Both the subphase and monolayer contribute to the total resistance experienced by the particle, which in turn is a function of interface and bulk viscosities, particle size and aspect ratio as well as the immersion depth of the particle. We explore how the drag on a spheroidal particle at a viscous interface can both rise or decrease with particle size depending on the dimensionless Boussinesq and Marangoni numbers. For a surfactant-laden interface, the surfactant distribution in the vicinity of the moving spheroid is significantly affected by the particle's immersion depth. When a particle sinks more in the viscous fluid, as determined by the involved surface tensions, the difference in surfactant concentration between front and rear of the particle decreases. For the drag coefficient of a spherical particle at an incompressible interface at low shear Boussinesq numbers, we propose a correction to previously reported analytic expressions. We probe both parallel and perpendicular friction coefficients as they are significantly different depending on particle shape, qualitatively different depending on surface shear viscosity, and we resolve the full three-dimensional distortion of the flow field around the moving particle.

1. Introduction

At vanishing Reynolds numbers, also known as creeping flow or the Stokes flow regime, the flow field around a translating sphere fully immersed in an unbounded fluid is one of the best characterized in fluid mechanics (Leal 2007). The analytical solutions of the flow field around asymmetric prolate and oblate particles in an unbounded fluid with no-slip boundary conditions on the particles are also known (Brenner 1963). In the case of a sphere moving along the interface between two fluids, however, the problem becomes more complex and cannot be solved analytically (Chisholm & Stebe 2021). This implies that the drag force on an adsorbed particle moving tangentially to an interface cannot be calculated analytically for arbitrary contact angles. Adding surfactant to the interface brings additional complications (Manikantan & Squires 2020). The excess of surfactant molecules leads to an interface viscosity, which increases the resistance of, and consequently the drag on, the particle. A translating particle disturbs the surface concentration of the surfactant at the interface. The variation in the surfactant concentration causes a gradient in interface tension and gives rise to an extra force, known as the Marangoni force, in the opposite direction of the particle motion. A variation in the surfactant concentration also generates a corresponding Marangoni flow, from the high concentration (low interface tension) region to the low concentration (high interface tension) region that counterbalances the Marangoni force (Pourali et al. 2021).

The aforementioned phenomena constituting a complex interaction between the translating particle, bulk fluid, interface rheology and surfactant transport are such that numerical or experimental methods are required to reveal the nature of this phenomena and determine the frictional forces on particles moving along an interface (Jaensson, Anderson & Vermant 2021). The drag coefficient is an important quantity in predicting and analysing interfacial rheological properties or trajectories obtained from particle tracking experiments (Bonales et al. 2007; Maestro et al. 2011). Knowledge about the motion and diffusion of an anisotropic particle at the interface is crucial for the understanding of biological systems (Ding, Warriner & Zasadzinski 2001), micro-organism motion (Lauga et al. 2006; Masoud & Stone 2014; Shaik & Ardekani 2017), micrometre sized ‘Marangoni surfers’ (Dietrich et al. 2020), and inclusions such as proteins, or other membrane-bound particles, in biological membranes, and the design of artificial membranes (Ally & Amirfazli 2010). Small particle probes also give microrheological methods increased sensitivity compared with macroscopic methods (Samaniuk & Vermant 2014).

A starting point for the analysis and quantitative understanding of these systems is the hydrodynamic model for the translation of a cylindrical particle in a slab of a viscous, incompressible membrane with thickness $d$, which mimics the protein motion in bilayers presented by Saffman & Delbrück (1975). According to Saffman's model, the drag force should be a logarithmic function of particle size. The bulk-surface coupling is described by a hydrodynamic length scale, the Saffman–Delbrück length $\mathcal {L}$ (Saffman & Delbrück 1975; Saffman 1976),

(1.1)\begin{equation} \mathcal{L}=\frac{\eta^s}{\eta_a + \eta}, \end{equation}

where $\eta ^s$ is the membrane surface viscosity, usually considered proportional to membrane thickness $d$, and $\eta _a$ and $\eta$ are bulk viscosities of the adjacent fluids, here air (vanishing $\eta _a$) and a viscous fluid. Consider a situation where the particle axis of symmetry is in the direction normal to the membrane and translating with a constant velocity normal to its axis of rotational symmetry. The model predicts the following relation for the drag coefficient $f$ for a non-protruding cylindrical particle with radius $R$ whose length $l$ is identical with the membrane thickness (Saffman 1976; Dimova et al. 2000):

(1.2)\begin{equation} f = \frac{4{\rm \pi} \eta^s}{\ln ( \mathcal{L}/R) - C}. \end{equation}

Here $C = 0.5772257$ is the Euler–Masceroni constant and the ratio $\mathcal {L}/R$ is called the Boussinesq number. Equation (1.2) assumes $\mathcal {L}\gg R$. In general, for a particle with characteristic linear size $a$, the Boussinesq number $\mathcal {L}/a$ quantifies the relative importance of the interfacial shear stress to bulk stress. In the rest of this work we will denote this number by ${\textit {Bq}_{1}}$, as we are going to consider a dilatational surface viscosity as well, giving rise to an additional dimensionless number ${\textit {Bq}_{2}}$.

Numerous biological studies refer to Saffman's continuum approach and some of them show disagreement with the original study. For example, Gambin et al. (2006) measured translational diffusion coefficients $D=k_{B}T/f$ of cylindrical peptides in a surfactant bilayer and showed that the drag coefficient $f$ is proportional to the radius $R$ of the diffusing object, $f\propto \eta ^s R$, and, therefore, there is a sharp qualitative disagreement with the Saffman–Delbrück model concerning the dependence on $R$. The effects of inclusion size in the membrane was studied by Levine & MacKintosh (2002) and Levine, Liverpool & MacKintosh (2004). They calculated the drag coefficient of a non-protruding rigid cylindrical rod ($l\ll R$) by solving the coupled equations for in-plane and out-of-plane fluid motions, assuming incompressibility of both the bulk and the membrane. In their work, the rod's axis of symmetry is parallel to the interface. They showed that (i) for small objects (specifically, $l \ll \mathcal {L}$ ), the drag coefficients become independent of both the rod orientation and aspect ratio; and (ii) for larger rods ($l>\mathcal {L}$), with high aspect ratio, the drag coefficient in perpendicular motion $f^{\perp }$ becomes purely linear in the rod length $l$. These results are qualitatively different from the motion of a rod in a three-dimensional bulk fluid with viscosity $\eta$. In the limit of low Reynolds number, the viscous drag on a rod is anisotropic and exhibits a logarithmic length dependency,

(1.3)\begin{equation} f = \frac{2{\rm \pi} \eta l} {\ln (Al/R)}, \quad f^{{\perp}} = 2f, \end{equation}

where $f$ is the drag coefficient in parallel motion (the particle's centre moves in the direction of the particle's symmetry axis), and $A$ is a numerical factor of order of unity (Kirkwood & Auer 1951; Klopp, Stannarius & Eremin 2017). Comparing the results in two and three dimensions shows that the relation in (1.3) breaks down in two dimensions.

Exploiting a similar approach, Fischer (2004b) derived the drag force on an ideal needle of vanishing thickness moving in a surface film overlying a fluid of depth $H$. He showed that for a parallel motion at high Boussinesq numbers, at similar viscosities and water depths, the drag on a needle equals that on a disk if its length is 3.3 (for $H \gg R$) or 10.9 (for $H \ll R$) times longer than the diameter of the disk. For perpendicular motion at high Boussinesq number, a needle experiences the same drag as for parallel motion, if it is shorter by the factor $1/e$ than the edge-on moving needle.

There have been attempts to solve the Stokes equation for a protruding particle, at finite contact angles. Danov et al. (1995) were the first to solve the Stokes equation numerically for a three-dimensional particle that protrudes into the subphase. They modelled a compressible interface characterized by interface shear and dilatational viscosity. They reported the results for particles with contact angles between $20^{\circ }$ and $90^{\circ }$. In this model, the interface surface tension was assumed to be a constant, therefore, the effect of the Marangoni force was neglected. Dimova et al. (2000) later considered the Marangoni effect and surfactant diffusion by using the Gibbs elasticity of the interface in their modified model, but they still neglected a coupling between interface and subflow. These calculations are only valid for small deviations in the surfactant excess concentration, i.e. for small Péclet numbers. Fischer, Dhar & Heinig (2006) used an approach different from Levine et al. (2004), solving for the stresses due to the subphase and at the contact line separately for interfaces with a shear viscosity, under the assumption that the interface is incompressible. They presented solutions for contact angles between $0^{\circ }$ and $180^{\circ }$, as well as for immersed particles in the liquid near to the interface. Stone & Masoud (2015) studied the protrusion of oblate particles at the interface. They used a perturbation expansion for the velocity and drag force on the particle as a function of particle protrusion. Using the Saffman drag force as a zero order term in the expansion and applying the reciprocal theorem they calculated the first-order term in the expansion which is due to the protrusion. We are not aware of previous works that considered a spheroid at a viscous, compressible interface that generates Marangoni flows.

To fill this gap of insight, we here investigate the effect of contact angle on the Marangoni flow and the drag coefficient of spherical and prolate spheroidal particles (ellipsoids of revolution) translating with a constant velocity, both parallel and tangential to their principle axis, at the flat interface between fluid and air, while the interface is possibly carrying insoluble surfactant. The interface is characterized by both shear and dilatational viscosities. This work extends our previous study (Pourali et al. 2021) in two directions: to particles with different contact angles and to non-spherical particles.

For those readers mainly interested in real-world applications, it important to stress that the present work does not take into account any specific coupling between surface viscosities and surfactant concentration; all three are treated as independent parameters. This artificial setting allows us to explore the effects of the individual contributions to the flow fields and drag coefficients. A majority of the past literature suggests that insoluble surfactant monolayers are almost inherently incompressible due to Marangoni flow (Klingler & McConnell 1993; Steffen et al. 2001; Wurlitzer, Schmiedel & Fischer 2002; Fischer 2004a; Fischer et al. 2006; Manikantan & Squires 2020), while it has been recently clarified that the interface can be incompressible by dilatational viscosity, Marangoni effects or a combination of both (Pourali et al. 2021). The same work had shown that the drag coefficient of a spherical particle, symmetrically immersed at an incompressible interface, within the limit of vanishing interface shear viscosity, exhibited the same value regardless of the origin of incompressibility. The aforementioned empirical findings thus do not allow us to draw conclusions about the relevance of Marangoni effects. In fact, for the special case of inviscid interfaces carrying surfactants, the product between Marangoni (${\textit {Ma}}$) and Péclet ($\textit {Pe}_s$) numbers, but not just ${\textit {Ma}}$ alone, plays a crucial role in determining the dilatation of the interface (Pourali et al. 2021). At high values of ${\textit {Ma}}\,\textit {Pe}_s$ the interface is incompressible. A high value of interface dilatational viscosity makes the interface incompressible regardless of ${\textit {Ma}}\,\textit {Pe}_s$. When we are going to study drag coefficients at incompressible interfaces we are thus free to choose one out of the two routes. Treating the surface viscosities as independent variables is furthermore supported by recent experiments where the rheological (extra) stresses are large with respect to the thermodynamic ones (Peppicelli et al. 2019). Since there are different routes to interface (in)compressibility, inferring the exact nature of the interface using particle probes was shown to possibly lead to inconsistent results, e.g. when changing the particle size or aspect ratio (Samaniuk & Vermant 2014), giving further motivation to the current work.

After presenting the governing equations, definitions of dimensionless numbers like ${\textit {Ma}}$ and $\textit {Pe}_s$, and the numerical scheme in § 2, within the results § 3 we are going to investigate various extremal and less extremal situations, focus on the drag coefficient, discuss the effect of the flow and surfactant concentration fields, compare with theoretical results, for both spherical and spheroidal particles. Assumptions to be made, concerning the interfacial equation of state or the neglect of the interface deformations are discussed in §§ 2.3 and 3.1, respectively. We here choose the length of the axis of rotational symmetry of the spheroid to define dimensionless quantities. All results to be presented then equally apply, after suitable scaling with $\mathcal {D}$, to spheroids with constant volume or constant surface area, as explained in § 3.4. Conclusions are provided in § 4.

2. Model and methods

A sketch of the system under study is shown in figure 1. A particle is translating with a constant velocity $U$ in the $x$-direction at the interface between a viscous fluid ($y<0$) and air ($y>0$). The particle shape is defined by the lengths $a$, $b$ and $c$ of the three principal semi-axes. Here we study prolate spheroids with half-axes $a\ge b=c$, which includes the sphere as a special case ($a=b=c \equiv R$). Spheroids translate either parallel or perpendicular to their axis of uniaxial symmetry, so that the system is torque-free, and can reach a steady state. The submergence of the particle is defined by the coordinate of its centre at $y=-h$, while the interface is located at $y=0$ and separates a viscous fluid with viscosity $\eta$ at $y < 0$ and the inviscid fluid at $y > 0$ (figure 1).

Figure 1. Prolate spheroid at the interface between air (white) and a viscous fluid (blue). The interface is eventually laden with surfactant (not shown here, but if so, we add a thick black line representing surfactant). The particle geometry is specified by the three principal semi-axes $a, b$ and $c$ with $a\ge b=c$ for a prolate spheroid (and $a\le b=c$ for an oblate spheroid, so that $a$ is the length of the main, not necessarily longest, axis). Their ratio is $\mathcal {D}=a/b\in [0,\infty ]$ with $\mathcal {D}>1$ for prolate spheroids, $\mathcal {D}=1$ for a sphere and $\mathcal {D}\in (0,1)$ for oblate spheroids. The particle translates with constant velocity $U$ in a positive $x$-direction, while its symmetry axis resides within the interfacial $x$–$z$-plane. For parallel and perpendicular motions, the axis of uniaxial symmetry is either aligned in the $x$- or $z$-direction, respectively. The submergence of the particle is defined by the $y$-coordinate of its centre, denoted by $h$, giving rise to dimensionless negative immersion depth $\mathcal {H}=h/b$, i.e. $\mathcal {H}=-1$ if the particle is fully immersed in water, in grazing contact with the interface. We are going to introduce dimensionless quantities using the principal length $a$ as the length unit (§ 2.4). All results then equally apply, after suitable scaling with $\mathcal {D}$, to spheroids with constant volume or constant surface area (§ 3.4).

2.1. Hydrodynamic equations

The fluid is considered incompressible. Its dynamics can thus be modelled with the Stokes equations

(2.1)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\rm \pi}=\boldsymbol{0}, \end{gather}

(2.2)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

where $\boldsymbol {u}$ is the fluid velocity field, $p$ is the pressure field and $\boldsymbol {{\rm \pi} } = -p \boldsymbol{\mathsf{I}} + \boldsymbol {\tau }$ the stress tensor for a Newtonian fluid modelled by $\boldsymbol {\tau } = 2 \eta \boldsymbol{\mathsf{D}}$. Here $\boldsymbol{\mathsf{D}} = [\boldsymbol {\nabla }\boldsymbol {u} + (\boldsymbol {\nabla }\boldsymbol {u})^\textrm {T}]/2$ is the rate of deformation tensor. All fields appearing in our equations are spatio-temporal fields that depend on $x$, $y$, $z$ and time $t$. A similar decomposition $\boldsymbol {{\rm \pi} }^s = \gamma \boldsymbol{\mathsf{I}}_s + \boldsymbol {\tau }^s$ can also be written for the interface stress tensor, where $\gamma$ is the surface tension of the interface, $\boldsymbol{\mathsf{I}}_s=\boldsymbol{\mathsf{I}}-\boldsymbol {nn}$ the surface or tangential projection tensor with surface normal $\boldsymbol {n}$ and $\boldsymbol {\tau }^s$ the extra surface stress tensor (Brenner 1991; Jaensson & Vermant 2018; Venerus & Öttinger 2018). It is assumed to be given by the Boussinesq–Scriven constitutive law (Boussinesq 1913; Scriven 1960)

(2.3)\begin{equation} \boldsymbol{\tau}^s = (\kappa^s - \eta^s)(\boldsymbol{\mathsf{I}}_s:\boldsymbol{\mathsf{D}}_s)\boldsymbol{\mathsf{I}}_s + 2\eta^s \boldsymbol{\mathsf{D}}_s, \end{equation}

where $\eta ^s$ and $\kappa ^s$ are the shear and dilatational viscosities of the interface, respectively. The surface rate of deformation tensor $\boldsymbol{\mathsf{D}}_s$ appearing in (2.3) is defined as $2\boldsymbol{\mathsf{D}}_s = (\boldsymbol {\nabla }_s \boldsymbol {u}_s) \boldsymbol {\cdot } \boldsymbol{\mathsf{I}}_s + \boldsymbol{\mathsf{I}}_s \boldsymbol {\cdot } (\boldsymbol {\nabla }_s \boldsymbol {u}_s)^\textrm {T}$, where $\boldsymbol {u}_s$ is the velocity $\boldsymbol {u}$ evaluated at the interface, and $\boldsymbol {\nabla }_s=\boldsymbol{\mathsf{I}}_s\boldsymbol {\cdot }\boldsymbol {\nabla }$ is the surface gradient operator (Brenner 1991).

The velocity $\boldsymbol {u}$ is assumed to vanish on the simulation box surface. Boundary conditions for $\boldsymbol {u}$ at the fluid interface are defined by continuity of velocity in the tangential direction $\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {t} = \boldsymbol {u}_s \boldsymbol {\cdot } \boldsymbol {t}$ (Venerus & Öttinger 2018), where $\boldsymbol {t}$ is a unit tangent vector residing in the $x$–$z$-plane, and a vanishing normal velocity $\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {n} = 0$. Moreover, conservation of momentum yields a momentum jump balance in the tangential direction,

(2.4)\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\rm \pi} \boldsymbol{\cdot} \boldsymbol{t} ={-}(\boldsymbol{\nabla}_s \!\boldsymbol{\cdot} \boldsymbol{\tau}^s) \boldsymbol{\cdot} \boldsymbol{t}+ K_{\rm \pi} (\boldsymbol{\nabla}_s \ln \varGamma) \boldsymbol{\cdot} \boldsymbol{t}, \end{equation}

where the Gibbs–Marangoni modulus $K_{{\rm \pi} }=\varGamma \partial \varPi ^s /\partial \varGamma$ allows one to relate the gradient in surface pressure to the gradient in surfactant concentration $\varGamma$ as $\boldsymbol {\nabla }_s \varPi ^s = K_{{\rm \pi} } \boldsymbol {\nabla }_s \ln \varGamma$. Surface pressure $\varPi ^s$ is defined as a difference between surface tension of a clean interface $\gamma _0$ and surface tension in the presence of surfactant, $\varPi ^s(\varGamma ) = \gamma _0-\gamma (\varGamma )$. The velocity and pressure fields thus receive their time dependency through $\varGamma$. The evolution of the surfactant concentration $\varGamma$ is governed by the unsteady surface convection–diffusion (SCD) equation (Brenner & Leal 1978; Stone 1990; Brenner 1991)

(2.5)\begin{equation} \frac{\partial \varGamma}{\partial t} + \boldsymbol{\nabla}_s\boldsymbol{\cdot} (\varGamma \boldsymbol{u}_s)= D_s \nabla_s^2 \varGamma, \end{equation}

where $D_s$ is the surface diffusivity of the surfactant at the planar interface. The Stokes equations supplemented by (2.5) are the governing equations for $\boldsymbol {u}$ and $\varGamma$ as a function of position and time. These equations are solved with an initial condition of $\varGamma =\varGamma _0$ and subject to vanishing surfactant flux from the interface boundary.

2.2. Drag force

For a spheroidal particle translating with a constant velocity $\boldsymbol {U}$ at the interface, the relation between the drag force on the particle $\boldsymbol {F}$ and the velocity, in general, is a complex function of particle geometry, bulk and interface rheological properties, the interfacial equation of state and the transport properties of the surfactant. It moreover depends on the orientation of the anisotropic particle with respect to $\boldsymbol {U}$. The drag force on the spheroidal particle embedded at the interface is the sum of three contributions: (i) the bulk force

(2.6)\begin{equation} \boldsymbol{F}^b=\int_{S_p} \boldsymbol{n}_p \boldsymbol{\cdot} \boldsymbol{\rm \pi} \, \textrm{d}S, \end{equation}

where $\boldsymbol {n}_p$ is the unit normal vector to the surface of the particle; (ii) the interface viscous force

(2.7)\begin{equation} \boldsymbol{F}^s=\int_{\partial S_p}\boldsymbol{n}_p \boldsymbol{\cdot} \boldsymbol{\tau}^s\,\textrm{d}l \end{equation}

due to the extra surface stress tensor, where $\partial S_p$ is the elliptic perimeter of the particle at the interface; and (iii) the elastic or Marangoni contribution to the drag force

(2.8)\begin{equation} \boldsymbol{F}^M=\int_{\partial S_p}\gamma (\varGamma) \boldsymbol{n}_p \boldsymbol{\cdot} \boldsymbol{\mathsf{I}}_s\,\textrm{d}l \end{equation}

due to the non-uniform distribution of the surfactant in the contact line.

So far, we presented the general formulation. Because in our set-up the particle translates in the $x$-direction with its main axis aligned in either the $x$- or $z$-direction, the force has only an $x$-component, $F_x = -fU$, the remaining two components vanish for symmetry reasons.

2.3. Interfacial equation of state

For an isothermal system, the surface tension is solely a (typically nonlinear) function of the surfactant concentration $\varGamma$, i.e. $\gamma = \gamma (\varGamma )$. In this work we assume a linearized equation of state: $\gamma _L = \gamma _* - \varGamma k _{B}T_*$, where $\gamma _*=\gamma (\varGamma _0)+\varGamma _0k_{B}T_*$ is a constant, and $k_{B}T_*$ represents the negative of the slope of $\gamma$ with respect to $\varGamma$, taken at $\varGamma _0$, the equilibrium surfactant surface number density in the absence of the particle. The linearized Gibbs modulus is thus $K_{\rm \pi} = \varGamma k_{B}T_*$. For the special case of sufficiently small $\varGamma _0\rightarrow 0$, one has $T_*\rightarrow T$ and the equation of state reduces to the ‘ideal gas’ form $\gamma =\gamma _0-\varGamma k_{B}T$ due to remaining kinetic contributions, and where $\gamma _0=\gamma (0)$ is the surface tension of the clean interface.

Nonlinearities are likely to set in with increasing concentration, and nonlinear equations of state have been discussed in the literature (Lopez & Hirsa 2000; Manikantan & Squires 2020). For cases where the concentration does not vary significantly across the interface, an equation of state linearized about a certain $\varGamma$ might still serve as a good approximation, so that $\gamma _0$ and $k_{B}T$ just receive new interpretations. For this reason, the linearized form does not restrict this study to systems at very low surfactant concentrations. We are not aware of an established and essentially parameter-free nonlinear interfacial equation of state that could have been used instead of the linearized one, for the purpose of the present study. Similarly, for the current analysis, surface viscosities are assumed to be system parameters, and their dependence on $\varGamma$ is neglected (Scriven 1960; Ortega, Ritacco & Rubio 2010). A clean, surfactant-free interface has $\varGamma =0$. In most practical cases one needs some surface active species at the interface to get significant extra (viscous) stresses, but surface viscosities for clean interfaces have also been reported (Earnshaw 1981).

2.4. Dimensionless numbers and equations

For the problem at hand, we use the constant velocity $U$, the semi-axis $a$ of the particle, the viscosity $\eta$ of the fluid and the initial surfactant concentration $\varGamma _0$ to introduce reduced units, and to come up with a number of dimensionless parameters such as two Bousinessq numbers ${\textit {Bq}_{1}}$ and ${\textit {Bq}_{2}}$, ratio of surface shear viscosity to bulk viscosity, and ratio between surface dilatation viscosity to bulk viscosity, and the surface Péclet number $\textit {Pe}_s$, a ratio of diffusion time, $a^2/D_s$, to the characteristic time for particle motion, $a/U$, i.e.

(2.9a–c)\begin{equation} {\textit{Bq}_{1}}=\frac{\eta^s}{a\eta}, \quad {\textit{Bq}_{2}}=\frac{\kappa^s}{a\eta}, \quad \textit {Pe}_s =\frac{aU}{D_s}. \end{equation}

The particle aspect ratio $\mathcal {D}$ is defined by the ratio between half-axes $a$ and $b$, and its dimensionless immersion $\mathcal {H}$ is specified by the ratio between $h$ and $b$, i.e.

(2.10a,b)\begin{equation} \mathcal{D} = \frac{a}{b}, \quad \mathcal{H} = \frac{h}{b}, \end{equation}

so that $\mathcal {H}\le -1$ at complete immersion and $\mathcal {H}=0$ for a particle symmetrically located at the interface. Making the interfacial equation of state dimensionless, we obtain $\gamma (\varGamma )/\eta U = {\textit {Ca}}^{-1} - {\textit {Ma}} \, \varGamma /\varGamma _0$, giving rise to the dimensionless capillary number ${\textit {Ca}}=\eta U/\gamma _*$ (ratio between bulk stress and interface tension force) and Marangoni number

(2.11)\begin{equation} {\textit{Ma}}= \frac{\varGamma_0 k_{B}T_*}{\eta U}, \end{equation}

the ratio between the force tending to deform the interface and surface elasticity which tends to restore the original shape of the film.

Dimensionless variables, marked by a tilde (only in the current sentence), are therefore defined uniquely in terms of their dimensional counterparts, such as $\tilde {\gamma }=\gamma /\eta U$, $\tilde {t}=t U/a$, $\tilde {x}=x/a$, $\tilde {\boldsymbol {u}}=\boldsymbol {u}/U$, $\tilde {\varGamma }=\varGamma /\varGamma _0$, $\tilde {\boldsymbol {{\rm \pi} }}=\boldsymbol {{\rm \pi} } a/(\eta U)$, $\tilde {\boldsymbol {F}}=\boldsymbol {F}/(\eta U a)$ and $\tilde {f}=f/(\eta a)$. For the rest of this work, variables are meant to be dimensionless, and we omit asterisks for brevity. The dimensionless equation of state thus reads as $\gamma ={\textit {Ca}}^{-1}-{\textit {Ma}}\,\varGamma$. Similarly, $\varPi ^s=K_{{\rm \pi} }={\textit {Ma}}\,\varGamma$ in the dimensionless form, which shows that the Marangoni contribution to the force is proportional to ${\textit {Ma}}$. The Stokes equations are replaced by $-\boldsymbol {\nabla } p + \nabla ^2 \boldsymbol {u}=0$ and $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = 0$, the dimensionless extra surface stress tensor (2.3) now involves the Bousinessq numbers,

(2.12)\begin{equation} \boldsymbol{\tau}^s = ({\textit{Bq}_{2}} - {\textit{Bq}_{1}})(\boldsymbol{\mathsf{I}}_s:\boldsymbol{\mathsf{D}}_s)\boldsymbol{\mathsf{I}}_s + 2 {\textit{Bq}_{1}} \boldsymbol{\mathsf{D}}_s, \end{equation}

and the SCD equation (2.5) becomes

(2.13)\begin{equation} \frac{\partial \varGamma}{\partial t} + \boldsymbol{\nabla}_s\boldsymbol{\cdot} (\varGamma \boldsymbol{u}_s)= \frac{1}{\textit {Pe}_s}\nabla_s^2 \varGamma, \end{equation}

featuring the Péclet number $\textit {Pe}_s$. The tangential momentum jump balance, (2.4), remains unchanged.

We use the finite element method (FEM), implemented in an in-house code, to solve the full system of equations, where the dimensionless parameters ${\textit {Bq}_ {1}}$, ${\textit {Bq}_ {2}}$, $\textit {Pe}_s$, $\mathcal {D}$, $\mathcal {H}$ and ${\textit {Ma}}$ defined by (2.9a–c)–(2.11) completely describe the problem. The Marangoni number enters only in the presence of surfactant. In the limit of $\textit {Pe}_s \to 0$, the surfactant concentration remains uniform and regaining the uniform distribution is instantaneous, therefore, the surface pressure remains constant. In this surface viscosity-dominated regime both the ${\textit {Ma}}$ and $\textit {Pe}_s$ numbers do not enter as parameters.

2.5. Numerical implementation

The numerical implementation is similar to the one found in Pourali et al. (2021), with the notable difference that we explicitly track the motion of the particle in the current work. In order to diminish the boundary effects a large, cubic simulation box is used with box size $L=400$ (figure 2). An extensive study on the influence of the size of the bounding box for a spherical particle at an incompressible interface was conducted in our previous work (Pourali et al. 2021). There, we showed that a box size of $L=400$ yields an error of less than 1 % for the drag coefficient (compared with a much larger box). Moreover, due to symmetry of the problem in the $x$–$y$-plane, only half of the domain is meshed, and appropriate symmetry conditions are employed. As in this previous work, we denote the simulation domain containing fluid and air by $\varOmega$, the box boundary by $\partial \varOmega$, the air–liquid interface by $S_I$, the particle surface by $S_p$ and the intersection of the particle and the air–liquid interface by $\partial S_p$. The stationary regime is reached in each case at $t=30$, as we will demonstrate in § 3.6.

Figure 2. Finite element mesh used in the simulation for a prolate particle with $\mathcal {D} = 3$ at $\mathcal {H}=-0.6$. The blue surface shows the air–liquid interface $S_I$. The box is cubic with box sizes $L_x=L_y=L_z=L=400$, which is 400 times the longer half-axis of this spheroid.

2.5.1. Arbitrary Lagrange Euler formulation

To track the moving particle, the arbitrary Lagrange Euler (ALE) formulation is adopted (Hu, Patankar & Zhu 2001). As explained later, the mesh is moved exactly with the particle velocity on $S_p$, whereas it is stationary on the boundaries $\partial \varOmega$ of the enclosing box. To compensate for the motion of the mesh, the mesh velocity $\boldsymbol {u}_m$ is subtracted from the convective terms. For our problem, this only affects the SCD, which becomes

(2.14)\begin{equation} \left.\frac{\partial \varGamma}{\partial t}\right|_{\boldsymbol{x}_m} + (\boldsymbol{u}_s - \boldsymbol{u}_\text{m}) \boldsymbol{\cdot} \boldsymbol{\nabla}_s \varGamma + ( \boldsymbol{\nabla}_s \boldsymbol{\cdot} \boldsymbol{u}_s)\varGamma= D_s \nabla_s^2 \varGamma, \end{equation}

where we note that the partial derivative on the left-hand side is for a fixed nodal coordinate $\boldsymbol {x}_m$.

2.5.2. Weak forms

The weak form of the balance of momentum (2.1) and balance of mass (2.2) amounts to find $\boldsymbol {u}$ and $p$ such that (Carrozza, Hulsen & Anderson 2020; Balemans, Hulsen & Anderson 2016; Jaensson, Hulsen & Anderson 2017)

(2.15)\begin{gather} \int_{\varOmega} \boldsymbol{D}_v : 2\eta\boldsymbol{D} \, \textrm{d}V- \int_{\varOmega} (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}) p \, \textrm{d}V + \int_{S_I} (\boldsymbol{\nabla} \boldsymbol{v})^\textrm{T} : \boldsymbol{\rm \pi}^s \, \textrm{d}A = 0, \end{gather}

(2.16)\begin{gather}\int_{\varOmega} q (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}) \, \textrm{d}V= 0, \end{gather}

for all admissible test functions $\boldsymbol {v}$ and $q$. Here $\boldsymbol{D}_v = [\boldsymbol {\nabla } \boldsymbol {v} + (\boldsymbol {\nabla } \boldsymbol {v})^\textrm {T}]/2$ and $\boldsymbol {v}$ vanishes on the Dirichlet boundaries where the velocity $\boldsymbol {u}$ is prescribed, i.e. the box boundary $\partial \varOmega$ and particle surface $S_p$ (Donea & Huerta 2003).

The weak form of the SCD equation amounts to find $\varGamma$ such that

(2.17)\begin{equation} \int_{S_I} r \left[ \frac{\partial \varGamma}{\partial t} + (\boldsymbol{u}_s - \boldsymbol{u}_\text{m}) \boldsymbol{\cdot} \boldsymbol{\nabla}_s \varGamma + ( \boldsymbol{\nabla}_s \boldsymbol{\cdot} \boldsymbol{u}_s)\varGamma \right] \, \textrm{d}A={-}\int_{S_I} D_s (\boldsymbol{\nabla}_s r) \boldsymbol{\cdot} (\boldsymbol{\nabla}_s \varGamma) \,\textrm{d}A \end{equation}

for all admissible test functions $r$. Here we used the zero-flux boundary condition on $S_I\cap \partial \varOmega$.

2.6. Spatio-temporal discretization

The (2.15)–(2.17) are solved using the Galerkin FEM on meshes that are boundary fitted to the particle and the interface, and which are moved in time to track the motion of the particle. Tetrahedral $P_2P_1$ elements (Taylor & Hood 1973) are used for $\boldsymbol {u}$ and $p$ whereas triangular $P_2$ elements are used for $\varGamma$. Mesh generation is done using Gmsh (Geuzaine & Remacle 2009), which allows for great control over the local element size.

2.6.1. Time integration

Time integration commences by generating a mesh with nodal coordinates $\boldsymbol {x}_m = [x_m,y_m,z_m]^\textrm {T}$, based on the initial particle's centre position $\boldsymbol {X}_0 = [X_0,0,0]^\textrm {T}$. Then, at the beginning of a step at time $n \Delta t + \Delta t$ (for which we use the short-hand notation $n+1$, e.g. $X(n \Delta t + \Delta t)=X_{n+1}$), the exact new particle position is determined from $X_{n+1} = X_n + \Delta t \, U$, where we note that only the $X$ location has to be updated. The particle displacement is thus given by $\Delta X = X_{n+1} - X_n = \Delta t\, U$. To ensure a smooth deformation of the mesh, the new particle position is used to update the nodal coordinates $\boldsymbol {x}_m$ of the mesh by solving the following Laplace equation:

(2.18)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} (K_e \boldsymbol{\nabla} (\Delta x_m)) = 0 \quad (\boldsymbol{x}_m \in \varOmega), \end{gather}

(2.19)\begin{gather}\Delta x_m = 0 \quad (\boldsymbol{x}_m \in \partial \varOmega), \end{gather}

(2.20)\begin{gather}\Delta x_m = \Delta X \quad (\boldsymbol{x}_m \in S_\text{p}), \end{gather}

(2.21)\begin{gather}\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\nabla}(\Delta x_m) = 0 \quad (\boldsymbol{x}_m \in S_\text{I}). \end{gather}

Here $K_e$ is a diffusion coefficient which varies per element and which is chosen as the inverse of the element size. This approach ensures that most of the mesh deformation takes place in larger elements, minimizing mesh distortion (Hu et al. 2001). The mesh displacement field is then used to update the $x$-components of the node coordinates via $(x_m)_{n+1} =(x_m)_{n} + \Delta x_m$, while the $y$- and $z$-components remain unchanged.

With the mesh coordinates known at the new time step, the mesh velocity is found using a first-order backward differencing scheme for the first step,

(2.22)\begin{equation} (\boldsymbol{u}_m)_1 =\frac{(\boldsymbol{x}_m)_1-(\boldsymbol{x}_m)_0}{\Delta t}, \end{equation}

whereas a second-order backward differencing scheme is used for subsequent steps $n\ge 1$,

(2.23)\begin{equation} (\boldsymbol{u}_m)_{n+1} = \frac{\frac{3}{2}(\boldsymbol{x}_m)_{n+1}-2 (\boldsymbol{x}_m)_n+\frac{1}{2} (\boldsymbol{x}_m)_{n-1}}{\Delta t}. \end{equation}

On the updated mesh, the weak form of the governing equations, (2.15)–(2.17), are solved following an explicit scheme, as will be explained next. For comparison, we have also implemented an implicit scheme, as described in Appendix B. The comparison shows that results are basically unaffected by the choice of scheme.

2.6.2. Explicit scheme

In the explicit scheme the weak form of the balance equations, (2.15) and (2.16), is solved first to obtain the velocity and pressure at step $n+1$. For the surface concentration $\varGamma$, which is needed to evaluate the surface tension $\gamma =\gamma (\varGamma _{n+1})$ at step $n+1$, a first-order extrapolation is used for the first step $\gamma (\varGamma _{n+1})=\gamma (\varGamma _0)$ with $\varGamma _0=1$, whereas a second-order extrapolation is used for subsequent steps $\gamma (\varGamma _{n+1})=\gamma (2\varGamma _n - \varGamma _{n-1})$.

After solving this system for $\boldsymbol {u}_{n+1}$ and $p_{n+1}$, the surface velocity $(\boldsymbol {u}_s)_{n+1}$ is readily extracted, and used in the weak form of the SCD equation (2.17). Using a first-order semi-implicit Euler scheme for the first time step, we obtain

(2.24)\begin{equation} \frac{\varGamma_{1} - \varGamma_0}{\Delta t} + [(\boldsymbol{u}_s)_1 - (\boldsymbol{u}_m)_{1}] \boldsymbol{\cdot} \boldsymbol{\nabla}_s\varGamma_{1} + [\boldsymbol{\nabla}_s\boldsymbol{\cdot} (\boldsymbol{u}_s)_{n+1}]\varGamma_{1} = \boldsymbol{\nabla}_s \boldsymbol{\cdot} \left( D_s \boldsymbol{\nabla}_s\varGamma_{1} \right). \end{equation}

For subsequent time steps, a second-order, semi-implicit Gear scheme is used to evaluate $\varGamma _{n+1}$,

(2.25)$$\begin{gather} \frac{\frac{3}{2}\varGamma_{n+1} - 2\varGamma_n+\frac{1}{2}\varGamma_{n-1}}{\Delta t} + [(\boldsymbol{u}_s)_{n+1} - (\boldsymbol{u}_m)_{n+1} ]\boldsymbol{\cdot} \boldsymbol{\nabla}_s\varGamma_{n+1} + [\boldsymbol{\nabla}_s\boldsymbol{\cdot} (\boldsymbol{u}_s)_{n+1}]\varGamma_{n+1} \nonumber\\ = \boldsymbol{\nabla}_s \boldsymbol{\cdot} \left( D_s \boldsymbol{\nabla}_s\varGamma_{n+1} \right). \end{gather}$$

After solving the resulting system for $\varGamma _{n+1}$, all variables are now known at step $n+1$, and time integration can continue to the next time step.

3. Results and discussion

Within the following sections we are going to investigate (i) a spherical particle enclosing a certain contact angle with a clean, i.e. surfactant-free, fully compressible and inviscid interface; (ii) a spherical particle at a viscous interface dominated by surface viscosities; (iii) a prolate spheroidal particle symmetrically located at a surfactant-free, fully compressible and inviscid interface; (iv) a prolate spheroid at a surface viscosity-dominated regime; and (v) a prolate spheroid at a sufactant-laden, partially compressible and viscous interface. While in (i)–(iii) the focus is on the drag coefficient and comparison with existing theoretical results, in (iv) and (v) we are going to further investigate the flow and surfactant concentration fields at the interface, for both spherical and spheroidal particles. All results to be presented are obtained for particles moving at constant speed, deep in the stationary regime. The case of incompressible interfaces, realized at large ${\textit {Bq}_ {2}}$, or large ${\textit {Ma}}\,\textit {Pe}_s$, is captured by (ii) for the case of a sphere, and by (iv) for the case of a prolate spheroid.

3.1. Spherical particle at a clean, surfactant-free, fully compressible and inviscid interface

For a spherical particle with radius $R$ and, thus, $\mathcal {D}=1$ at a clean interface, we solve the above equations in the absence of surfactant, $\varGamma =0$, at vanishing Bousinessq numbers ${\textit {Bq}_ {1}}={\textit {Bq}_ {2}} = 0$ for a given dimensionless immersion level $\mathcal {H}=h/R$. This level $\mathcal {H}$ is determined by three involved interfacial tensions, and the drag coefficient is solely determined by the bulk force $\boldsymbol {F}^b$.

The surface activity of particles does not resemble surfactant molecules, due to their amphiphilic nature. A liquid at rest intersects a solid particle at a unique angle, which is called the contact angle, and is a key parameter when dealing with solid particles at the fluid interfaces. This concept has similarity with wetting phenomena. In wetting phenomena a droplet spreads on the surface of a solid surface, in contact with air. In the wetting process a new surface (between liquid and solid or liquid and air) is formed. Creating a new surface changes the energy of the system. Surface (or interface) energy is the work per unit area needed (or generated) to create a new surface. In the wetting process the balance of the interfacial energies between the solid, the liquid and the air, determines whether the liquid spreads or not. Paraphrasing from Zanini & Isa (2016), the extent by which a droplet spreads is determined by the point at which the energy gained in reducing the interface between the solid and air equals the energy penalty paid in creating new liquid–solid and liquid–air interfaces. This translates into the mechanical equilibrium of the three interfacial tensions at the three-phase contact line. Similarly, we can define the contact angle of a solid particle at a fluid interface. The presence of the particle at the interface is energetically favourable if this configuration has a lower energy than the situation where the particle is completely immersed in one of the fluids. Figure 3 shows a spherical particle with radius $R$ at the interface between air and a viscous fluid.

Figure 3. Schematic representation of a single spherical particle at the interface between air and a viscous fluid, at a macroscopic contact angle $\theta$.

The total energy of this system in the absence of flow, for the case of $|\mathcal {H}|\le 1$, can be written as (Davies et al. 2014)

(3.1)\begin{equation} E = {\rm \pi}R^2 \gamma\left[\mathcal{H}^2 + 2\mathcal{H}\left(\frac{\gamma_{pa}-\gamma_{pf}}{\gamma}\right) + \frac{2\gamma_{pf}+2\gamma_{pa}-\gamma}{\gamma}\right], \end{equation}

so that $E=4{\rm \pi} R^2 \gamma _{pf}$ if the particle is fully dissolved in the liquid ($\mathcal {H}<-1$), and where $\gamma _{pa}$, $\gamma _{pf}$ and $\gamma$ denote the interfacial tensions between particle–air, particle–fluid and fluid–air respectively. The equilibrium position of the particle is obtained from the condition $\partial E/\partial \mathcal {H} = 0$. This yields

(3.2)\begin{equation} \mathcal{H} ={-}\frac{\gamma_{pa} - \gamma_{pf}}{\gamma}, \end{equation}

where $\mathcal {H}$ is from now on the dimensionless equilibrium position. From trigonometry, as long as $|\bar {\mathcal {H}}|\le 1$, the particle prefers to reside at the interface and its equilibrium contact angle $\theta$ is given by Young's equation

(3.3)\begin{equation} \cos \theta ={-}\mathcal{H}. \end{equation}

Using this definition of the contact angle, the energy of the equilibrium configuration is given by

(3.4)\begin{equation} \bar{E} ={-}\gamma {\rm \pi}R^2 (1-\cos \theta)^2. \end{equation}

These energy considerations, which are based on the values of surface tensions of the different phases, therefore determine whether a floating spherical particle at the interface is energetically favourable or not. While the vertical particle motion is strongly suppressed or confined about the reduced equilibrium position $\bar {\mathcal {H}}$, the particle can freely move tangential to the interface. Knowing the dimensionless drag coefficient $f=6{\rm \pi}$ on a translating sphere in an unbounded fluid, one might expect that the parallel viscous drag experienced by a particle partially immersed in (or partially wetted by) two fluids is an average between the drags of the two fluids, weighted by the particle immersion depth in each phase.

It is straightforward to generalize the above expressions for the contact angles and immersion depth as a function of the surface tension to spheroids (Appendix A). For the equilibrium contact angle, measured in the $x$–$y$-plane (figure 1), we obtain

(3.5)\begin{equation} \cos\theta ={-}\frac{\mathcal{D}\mathcal{H}}{\sqrt{1-(1-\mathcal{D}^2)\mathcal{H}}}, \end{equation}

while the more important aspect is that (3.2) approximately holds also for spheroids (figure 20). We therefore proceed using the dimensionless immersion level $\mathcal {H}$ to present results, while keeping in mind that $\mathcal {H}$ can basically be replaced by a dimensionless combination of surface tensions. As for the sphere, $\mathcal {H}\in [-1,1]$ and $\mathcal {H}=-1$ represents the case of a spheroid that is completely immersed in the viscous liquid, in grazing contact with the interface.

For some particular cases, when one phase is highly viscous such as for the air–water system considered in this work, the contribution of the less viscous phase can be neglected. In these cases, the drag coefficient of a half-immersed particle at the interface, i.e. when $|\bar {\mathcal {H}}|\ll 1$ and $\theta \approx 90^{\circ }$, is $f\approx 3{\rm \pi}$ (Ranger 1978). With some additional assumptions, the drag coefficient of a particle at an interface has been quantitatively evaluated by hydrodynamic calculations. Zabarankin (2007) obtained the solution for hydrophilic particles ($\theta < 90^{\circ }$) by applying the same symmetry argument to a pair of fused spheres. The flow is computed for this new body, obtained by reflection of the immersed section of the sphere, and the numerical solutions were derived for a few contact angle values. Analytical expressions were later given by Dörr et al. (2016) and Villa et al. (2020) for hydrophilic particles, $\theta < 90^{\circ }$ ($\mathcal {H}\in [-1,0]$),

(3.6)\begin{equation} f_{{\rm \pi}/2} = \frac{1}{2}\left [1 + \frac{9}{16}\cos \theta - 0.139 \cos^2\theta + {{O}}(\cos^3 \theta)\right], \end{equation}

and for $90^{\circ } <\theta <180^{\circ }$ or equivalently, $\mathcal {H}\in [0,1]$,

(3.7)\begin{equation} f \approx 2\left (1 - \frac{\theta}{\rm \pi} \right) f_{{\rm \pi}/2} + 2 \left(\frac{\theta}{\rm \pi} - \frac{1}{2} \right) f_{\rm \pi}, \end{equation}

where $f_{{\rm \pi} }$ is defined as

(3.8)\begin{equation} f_{\rm \pi} = \frac{8}{9{\rm \pi}}\sin\theta \left\{1 + \frac{\cot(\theta/2)}{3{\rm \pi}} + {O}[\cot^2(\theta/2)]\right\}. \end{equation}

In figure 4 we compare our simulation results for a spherical particle with reported theoretical drag coefficient values from the literature, as a function of the contact angle. The drag coefficients are reported relative to the drag coefficient on a particle at a contact angle of $90^{\circ }$. The absolute value of $f$ at a contact angle of $90^{\circ }$ is $3{\rm \pi}$ corresponding to 50 % of the value of Stokes’ solution for a particle moving in an unbounded fluid. The more the particle sinks into the bulk phase upon decreasing its contact angle, or decreasing $\mathcal {H}$, the higher the drag coefficient. At a contact angle of about $\theta =25^{\circ }$, corresponding to an immersion level $\mathcal {H}\approx -0.9$, there is an increase of drag force up to 50 % compared with the particle at $\theta =90^{\circ }$. As is obvious from figure 4, our simulation results perfectly confirm the theoretical expressions (3.6)–(3.8) for a partially immersed sphere in the absence of surface stresses. We do not find evidence for a relevance of the mentioned higher-order terms.

Figure 4. Drag coefficient of a spherical particle versus contact angle or immersion length $\mathcal {H}\in [-1,0]$ at the clean interface between a viscous fluid and inviscid air. The drag coefficients are reported relative to the drag coefficient of a sphere at contact angle $90^{\circ }$, $f=3{\rm \pi}$. The diamonds and triangles represent data reported by Danov et al. (1995) and Zabarankin (2007), respectively. The red and blue curves are recovered from (3.6) and (3.7), respectively, while black circles represent our simulations.

Here we study flat interfaces and neglect possible deformations of the interface due to the presence of the particle. Based on previous studies of the drag coefficients obtained in the present work, we expect to be only marginally affected upon including deformation effects. Petkov et al. (1995) first studied experimentally deformation effects for the case of a spherical particle at the pure water–air interface. They showed that for small spheres, which do not create any substantial deformation of the fluid interface, the drag coefficient does not change significantly due to the deformation. For example, they showed that at $\theta = 82^{\circ }$ for a sphere with $R=222 \ \mathrm {mm}$, the drag coefficient $f/3{\rm \pi} \approx 1.08$. This value is very close to our simulation result at the same contact angle ($\,f/3{\rm \pi} \approx 1.07$). However, for a very large and heavy particle with a large deformation of the interface, the drag coefficient could be higher than Stokes’ drag coefficient ($\,f=6{\rm \pi}$). To the best of our knowledge there are no further experimental studies similar to (Petkov et al. 1995) for the effects of deformation on the drag coefficient of a spherical or non-spherical particle. A few experimental works had been dedicated to pairwise interactions of interfacial particles at liquid–liquid interfaces (Vassileva et al. 2005; Madivala, Fransaer & Vermant 2009). A recent numerical study, for a two-dimensional particle with $\theta =90^{\circ }$ at the distorted interface between two fluids (in a confined domain), showed that interfacial deformations do not seem to yield significant drag variations compared with the case of a flat interface and the simulation data indicated that interfacial distortions may increase or decrease the drag coefficient by no more than 10 % within the explored physical parameter space (Loudet et al. 2020).

3.2. Spherical particle at a viscous interface dominated by extra stresses

Next, we add material properties to the surface viscosity-dominated regime to see how they affect the drag coefficient of the sphere. While so far only the bulk force $\boldsymbol {F}^b$ gave rise to $f$, now $\boldsymbol {F}^s$ contributes as well. The Boussinesq numbers characterize the amount of extra surface stresses reflecting the material properties of the interface in the presence of a surface flow gradient. Except in some special cases, such as some biological membranes, in most systems the dilatational and shear viscosities are of the same order. Danov et al. (1995) reported the drag coefficient at different ${\textit {Bq}_ {1}}={\textit {Bq}_ {2}}$ (${\textit {Bq}_ {1,2}}$ for convenience) as a function of contact angle and showed that $f$ is almost independent of the contact angle at high interface viscosities ${\textit {Bq}_ {1,2}}$. Our simulation results show slightly smaller values for $f$; see figure 5. At small interface viscosities ${\textit {Bq}_ {1,2}}=1$, the particle dynamics is governed by the bulk phase which experiences the immersed particle volume. Therefore, $f$ decreases with increasing contact angle (decreasing immersion depth, increasing $\mathcal {H}$).

Figure 5. Drag coefficient of a spherical particle versus contact angle $\theta$, or equivalently, immersion length $\mathcal {H}\in [-1,0]$ for interfaces with different viscosities, i.e. different ${\textit {Bq}_ {1,2}}$. For ${\textit {Bq}_ {1,2}}=0$, the result is shown in figure 4.

Figure 6 shows the relative contribution $F^s/F$ of the interface on the total drag force $\boldsymbol {F} = \boldsymbol {F}^s + \boldsymbol {F}^b$ acting on the particle. At ${\textit {Bq}_ {1,2}}=1$, the interface contribution is up to 35 % for the symmetrically immersed sphere, $\theta = 90$. At ${\textit {Bq}_ {1,2}}=5$, bulk and interface have a comparable contribution to the total drag, the drag coefficient is therefore almost independent of the contact angle; see figure 5. At high ${\textit {Bq}_ {1,2}}=10$, the forces are mainly determined by the viscosities of the interface. For this reason, the drag on the particle increases by up to 50 % upon increasing the contact angle, until the immersion depth vanishes. For higher contact angles, $f$ becomes almost independent of $\theta$, because the particle is now exposed mainly to the inviscid air.

Figure 6. Interface contribution $F^s$ to the total drag force $F$ on a sphere as a function of contact angle $\theta$ or immersion length $\mathcal {H}\in [-1,0]$ for interfaces with different viscosities. At $\mathcal {H}=-1$, the particle is fully immersed in the viscous phase, and $F^s\rightarrow 0$, while the particle is symmetrically located at the interface at $\mathcal {H}=0$.

Fischer et al. (2006) solved the problem of translation of a spherical particle of radius $R$ embedded in an incompressible viscous monolayer, incorporating Marangoni effects by using a virtual image force source to impose surface incompressibility, with the surface shear viscosity $\eta ^s$, i.e. ${\textit {Bq}_ {1}}>0$, between two viscous phases with viscosities $\eta _a$ and $\eta$. They have obtained the following result for the translational drag coefficient $f$ as a series expansion of Boussinesq number ${\textit {Bq}} = \eta ^s/[R(\eta _a + \eta )]$ for $0<{\textit {Bq}} \ll 1$. For our set-up, ${\textit {Bq}}={\textit {Bq}_ {1}}$ and, thus,

(3.9)\begin{equation} f = f_0 + f_1({\textit{Bq}_{1}}) + {{O}}[({\textit{Bq}_{1}})^2]. \end{equation}

Fitted expressions (Fischer et al. 2006) for $f$ from numerical results gave the formulae for $f_0$ and $f_1$,

(3.10)\begin{equation} f_0 \approx 6{\rm \pi} \left\{\tanh \left[\frac{32(1-\mathcal{H})}{9{\rm \pi}^2}\right]\right\}^{1/2} \end{equation}

for any $\mathcal {H}$, and

(3.11)\begin{equation} f_1 \approx \left\{ \begin{array}{@{}ll} -4(2-\mathcal{H})^{{-}3/2}\ln \left[\dfrac{2}{\rm \pi} \arctan \left(\dfrac{2}{3}\right)\right] , & \mathcal{H} <{-}1, \\ -4\ln \left[\dfrac{2}{\rm \pi} \arctan \left(\dfrac{1-\mathcal{H}}{3}\right)\right] , & \mathcal{H} >{-}1, \end{array} \right. \end{equation}

where the original work used a different notation, $d=-(1+\mathcal {H})R$, the signed distance from the apex of the sphere to the plane of the interface. In other words, their negative $d$ coincides with the largest $y$-coordinate of the sphere (figure 1).

Figure 7 shows our simulation results for the drag coefficient of a spherical particle as a function of the contact angle at the incompressible interface. It was shown that an incompressible interface can be numerically realized with ${\textit {Bq}_ {2}}=1000$ at the limit of ${\textit {Bq}_ {1}} \to 0$ (Pourali et al. 2021). Recall that we have already explained (see the introduction) that incompressibility can be achieved by high ${\textit {Bq}_ {2}}$ or high ${\textit {Ma}}\,\textit {Pe}_s$ or a combination of both effects and the drag coefficient will be independent of the origin of the incompressibility (Pourali et al. 2021). The mechanism of high ${\textit {Ma}}\,\textit {Pe}_s$ is particularly important for colloidal or biological systems, where the relevant length (velocity) scales are often on the order of microns (per second). Here, ${\textit {Ma}}$ remains large even for trace surfactant concentrations well into the surface-gaseous regime (Bławzdziewicz, Cristini & Loewenberg 1999). It also appears that, under typical circumstances, surface diffusivity of the surfactant is insufficient to relax interfacial incompressibility (Chisholm & Stebe 2021).

Figure 7. Measured drag coefficient $f$ (black circles) of a spherical particle as a function of immersion level $\mathcal {H}\in [1,1]$ at an incompressible interface with ${\textit {Bq}_ {2}} =1000$ at the limit of ${\textit {Bq}_ {1}}\to 0$. Blue and red lines show the analytical expressions (3.9) and (3.12), respectively. The latter reproduces our data very well and is also compatible with the reported value for $f$ at $\mathcal {H}=0$ by Fischer et al. (2006), while $f$ vanishes at $\mathcal {H}=1$ in every case.

For a translational drag on a half-immersed sphere in a non-viscous monolayer, $f = 11.66$ which is in very good agreement with $f \approx 11.7$ or $f/3{\rm \pi} \approx 1.24$ reported by Fischer et al. (2006). It is higher than the drag coefficient on the particle at a free surface ($\,f/3{\rm \pi} =1$). The value of $f_0=11.7$ obtained by Fischer et al. (2006), as already noted by Pourali et al. (2021), is, however, not captured by their approximant (3.10). We thus fitted our result shown in figure 7, as it is also compatible with previous results for $\mathcal {H}=0$, to obtain

(3.12)\begin{equation} f_0 \approx 6{\rm \pi} \left\{\tanh \left[\frac{32(1-\mathcal{H})}{9{\rm \pi}^2}\right]\right\}^{5/11}. \end{equation}

The exponent $5/11$ is not physically motivated but represents the fitted value $0.455\pm 0.002$. This equation can be regarded as an improved version of (3.10). It captures our data by a maximum relative error of less then 1 % over the whole range $\mathcal {H} \in [-1,1]$, while the maximum relative error using the original formula exceeds 11 % for the largest $\mathcal {H}$. The reported value $f\approx 11.7$ is also recovered using this improved fitting formula, which differs in form by the original one only in the exponent, $5/11$ instead of $1/2$.

3.3. Spheroidal particle symmetrically located at a surfactant-free, fully compressible and inviscid interface

Another extremal case that serves to test analytical expressions is the non-spherical spheroid symmetrically ($\mathcal {H}=0$) located at the clean interface, in the absence of surfactant and surface viscosities. Happel & Brenner (1981) calculated the translational drag coefficient acting on a spheroidal particle translating in an unbounded domain at the $x$-direction with velocity $U$, permanently oriented such that its axis of rotational symmetry is aligned in the $x$-direction as well (so-called parallel motion). For the bulk drag coefficient, they obtained $f_{bulk} = 6{\rm \pi} K/\mathcal {D}$, so that $K/\mathcal {D}$ is the dimensionless friction-wise equivalent radius, a multiple of the length of the spheroidal particle. For prolate and oblate spheroids, their result for $K$, the ratio between equivalent sphere and the spheroidal particle radius, is given by

(3.13)\begin{equation} \frac{1}{K} = \frac{3}{4}\frac{\vartheta_p}{\mathcal{D}}\left[(\vartheta_p^2+1)\coth^{{-}1}(\vartheta_p) - \vartheta_p\right], \end{equation}

and

(3.14)\begin{equation} \frac{1}{K} = \frac{3}{4}\frac{\vartheta_o}{\mathcal{D}}\left[\vartheta_o+(1-\vartheta_o^2)\cot^{{-}1}(\vartheta_o)\right], \end{equation}

respectively. The dimensionless quantities $\vartheta _p$ and $\vartheta _o$ are defined by

(3.15a,b)\begin{equation} \vartheta_p = \frac{\mathcal{D}}{\sqrt{\mathcal{D}^2-1}}, \quad \vartheta_o =\frac{\mathcal{D}}{\sqrt{1-\mathcal{D}^2}}. \end{equation}

The corresponding expressions for the case of perpendicular motion for prolate and oblate spheroids are

(3.16)\begin{equation} \frac{1}{K^\perp} = \frac{3}{8} \frac{\vartheta_p}{\mathcal{D}} \left[\vartheta_p + (3-\vartheta_p^2)\ln(\mathcal{D}+\mathcal{D}/\vartheta_p)\right] \end{equation}

and

(3.17)\begin{equation} \frac{1}{K^\perp} = \frac{3}{8} \frac{\vartheta_o}{\mathcal{D}} \left[(3+\vartheta_o^2)\csc^{{-}1}(\sqrt{1+\vartheta_o^2})-\vartheta_o\right], \end{equation}

respectively. Note that our $\mathcal {D}$ is identical with the $\phi$ in Happel & Brenner (1981) and that we denoted by $a$ (our length unit) the length of the axis of rotational symmetry, while Happel & Brenner (1981) denoted by $a$ the length of the longest half-axis.

For a prolate particle symmetrically immersed at a clean interface between a viscous fluid and air, its drag is half of the value of a fully immersed particle, i.e.

(3.18)\begin{equation} f_c = 3{\rm \pi} \frac{K}{\mathcal{D}}. \end{equation}

Figure 8 shows the drag coefficient $f_c$ (subscript ‘$c$’ for clean) as a function of spheroid aspect ratio $\mathcal {D}$, which is reproduced using (3.18) with $K$ from (3.13). The drag coefficients from simulations for prolate particles with $\mathcal {D} = 1,2,3,4,5$ are also shown in the figure. The simulation results show very good agreement with the analytical values.

Figure 8. Dimensionless drag coefficients $f_c$ and $f^{\perp }_c$ versus aspect ratio $\mathcal {D}$ of a spheroidal particle symmetrically located at a clean, fully compressible, surfactant-free interface (${\textit {Bq}}_{1,2}=0$), moving in a direction of its symmetry axis, or perpendicular to it ($\perp$). The solid and dashed lines are (3.18) with $K$ from (3.13)–(3.17). Black squares are our data points for prolate spheroids in parallel motion.

3.4. Effect of shape for equal-volume or equal-area spheroids

In accord with our choice of units to make all quantities dimensionless, we throughout present a dimensionless drag coefficient $f$, whose dimensional counterpart is $a\eta f$, versus aspect ratio $\mathcal {D}$, or alternatively, the relative $f/f_c$ with respect to the clean interface. Because the volume of a spheroid is identical with the one of the equal-volume sphere of radius $\mathcal {D}^{-2/3}$ in units of $a$, friction coefficients as a function of aspect ratio at fixed particle volume are captured (here and below) upon multiplying the reported $f$ with $\mathcal {D}^{2/3}$. In several cases we are going to show drag coefficients or drag forces versus aspect ratio not only for spheroids with constant semi-axis $a$, but also for particles with constant volume (Appendix C). Considering particles with equal volume, the minimum relative frictional force (compared with the one of the equal-volume sphere) then occurs for movement parallel to the main axis of a slightly elongated prolate spheroid with aspect ratio $\mathcal {D}\approx 1.952$, and spheroids with $\mathcal {D}>3.813$ have more resistance than a sphere, in agreement with our results.

Rather than exploring the effect of aspect ratio $\mathcal {D}$ for spheroids with constant $a$ or constant volume, one could also choose to compare spheroids with equal surface area, using the transformation behaviour that follows from (A3) in Appendix A. There is no most appropriate choice, as the invariant quantity may depend on processing or biological conditions, but we found it appropriate to mention the transformation rules here.

3.5. Prolate particle at a viscous interface dominated by extra stresses

When a particle translates at the interface, the interface symmetrically compresses at the front of the particle and dilates at its rear. The interface compressibility has been quantified by calculating the local interface dilatation $\boldsymbol {\nabla }_s \boldsymbol {\cdot } \boldsymbol {u}_s$ in our previous work (Pourali et al. 2021). Interface compressibility also depends on the particle's contact angle. The interface compressibility for a prolate spheroidal particle with $\mathcal {D}=3$ at three different immersion levels $\mathcal {H}=-0.6,0,0.6$ is shown in the middle column of figure 9 for the choice ${\textit {Bq}}_{1,2}=1$. The flow field $\boldsymbol {u}_s$ and shear component of the surface velocity gradient, necessary to fully characterize the velocity gradient, are also shown in figure 9. These results demonstrate that with increasing particle immersion the interface is getting less compressible, as indicated by a decreasing $\boldsymbol {\nabla }_s \boldsymbol {\cdot } \boldsymbol {u}_s$. To express these findings quantitatively, we use the maximum dilatation at the rear of the interface, denoted as $(\boldsymbol {\nabla }_s \boldsymbol {\cdot } \boldsymbol {u}_s)_{max}$ as a measure for the maximum interface compressibility. The maximum dilatation for $\mathcal {D}=1,2,3$ as a function of $\mathcal {H}$ (figure 10) shows that the more the particle sinks in the viscous fluid the less the interface is compressible.

Figure 9. Interfacial flow field $\boldsymbol {u}_s$, interface velocity components, $u_{s,x}$ and $u_{s,z}$, interface compression $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}_s$ and shear component for a prolate spheroid in parallel motion (axis of rotational symmetry aligned in the vertical $x$-direction) with aspect ratio $\mathcal {D}=3$ at three different immersion lengths (a) $\mathcal {H} = -0.6$, (b) $\mathcal {H} = 0$ and (c) $\mathcal {H} = 0.6$, at the surface viscosity-dominated regime, partially compressible and viscous interface. In all simulations ${\textit {Bq}_ {1,2}}=1$. We used the symmetry of the system (the particle is moving in the $x$-direction) and show half of the interface. The corresponding fields for the case of perpendicular motion are shown in Appendix C.

Figure 10. Maximum interface expansion of the surface viscosity-dominated regime, partially compressible and viscous interface for the cases of a spherical ($\mathcal {D}=1$) and a prolate spheroid ($\mathcal {D}=3$) as a function of the particle's immersion level $\mathcal {H}\in [-0.75,0.75]$. In all simulation ${\textit {Bq}_ {1,2}}=1$ and the particle is moving in parallel direction of its axis of rotational symmetry.

Figure 11 highlights effects of particle shape and ${\textit {Bq}_ {1}}$ on the parallel and perpendicular drag coefficients of a particle partially immersed at interfaces with low ${\textit {Bq}_ {2}} = 1$ and high ${\textit {Bq}_ {2}} = 1000$. The drag coefficients are reported relative to the drag coefficient of a spherical particle (at the same condition). In figures 11(a) and 11(c) particles translate parallel to their principle axis, while the corresponding results for perpendicular translation are presented in figures 11(b) and 11(d). Depending on values of ${\textit {Bq}_ {1}}$, three different behaviours of $f$ with particle shape can be observed. At low ${\textit {Bq}_ {1}}$, the drag coefficient decreases with increasing aspect ratio qualitatively in agreement with the theoretical result for ${\textit {Bq}_ {1}}={\textit {Bq}_ {2}}=0$ (figure 8). With increasing ${\textit {Bq}_ {1}}$, the dependency of the drag coefficient on $\mathcal {D}$ diminishes and at high ${\textit {Bq}_ {1}} > 10$, it tends to increase linearly with $\mathcal {D}$. This behaviour can be related to the nature of the flow of a spheroidal particle which is a mixed flow with both shear and dilatational contributions. The effects of interface shear viscosity are more pronounced in particles with high aspect ratio. This can be seen in figures 12(a) and 12(b) where we compare the shear component of the interface flow gradient, which enters the stress tensor via the Boussinesq–Scriven constitutive law (2.3), for $\mathcal {D} = 2$ and $\mathcal {D}=5$. The results confirm that the shear component is higher for $\mathcal {D}=5$. At low ${\textit {Bq}_ {1}}$, upon increasing the particle aspect ratio, due to a corresponding decrease in particle volume, the drag coefficient decreases. At very high ${\textit {Bq}_ {1}}$, the shear effects determine the drag coefficient on the particle, and the drag coefficient therefore increases with $\mathcal {D}$. At intermediate ${\textit {Bq}_ {1}}$, these two effects (particle size and shear effect) cancel out each other, hence, the drag coefficient is independent of the particle's aspect ratio. The bulk and interface contributions to the drag force as a function of ${\mathcal {D}}$ for three different ${\textit {Bq}_ {1}} = 0.1, 1, 10$ are shown in figure 12(c). These results show that bulk and interface contributions exhibit opposing trends upon varying ${\mathcal {D}}$ and also upon varying ${\textit {Bq}_ {1}}$ so that for each ${\textit {Bq}_ {1}}$ (each ${\mathcal {D}}$) there seems to exist a critical ${\mathcal {D}}$ (critical ${\textit {Bq}_ {1}}$) that marks the transition between bulk- and interface-dominated drag. To appreciate the effect of particle shape at given particle volume, we show the scaled form of figure 12(c) in Appendix C. The parallel drag coefficient of prolate particles with $\mathcal {D}=2$, $3$ and $5$ versus immersion length at a clean interface is reported in figure 13(a). Corresponding results for spheroids at the incompressible interface with ${\textit {Bq}_ {2}}=1000$ are shown side-by-side in figure 13(b). The behaviour of $f$ for both types of interfaces is similar to a spherical particle.

Figure 11. (a,c) Parallel and (b,d) perpendicular drag coefficients of prolate particles symmetrically ($\mathcal {H}=0$) located at the surface viscosity-dominated regime as a function of particle shape $\mathcal {D}$ for different interface shear viscosities ${\textit {Bq}_ {1}}$. Particle located at (a,b) low ${\textit {Bq}_ {2}}=1$ and (c,d) high ${\textit {Bq}_ {2}}=1000$. Drag coefficients are reported relative to the drag coefficient of a spherical particle $f^{{sp}}$ at otherwise unchanged conditions. For spheres ($\mathcal {D}=1$), the reduced $f$ reaches unity in each case.

Figure 12. Shear component of the interface flow gradient for (a) $\mathcal {D} =2$ and (b) $\mathcal {D} =5$. The spheroid is symmetrically ($\mathcal {H}=0$) located at the interface, moving in the direction of its symmetry axis (parallel translation). The regime of the flow is surface viscosity-dominated and the interface is characterized by ${\textit {Bq}_ {1}}=0.01$ and ${\textit {Bq}_ {2}}=1$. The colourbar indicates a quantitative difference of the field in these two plots. (c) Bulk (solid line) and interface (dashed line) contributions to the drag force as a function of ${\mathcal {D}}$ at three different ${\textit {Bq}_ {1}} = 0.1, 1, 10$ and ${\textit {Bq}_ {2}}=1$, $\mathcal {H}=0$ as for (a,b), on a semilogarithmic scale. The data shown in (c) is plotted differently, as to eliminate the effect of particle volume, in Appendix C.

Figure 13. Parallel drag coefficient for three different prolate particles versus their immersion level $\mathcal {H}\in [-1,0.8]$ at (a) a surfactant-free, fully compressible, inviscid (${\textit {Bq}}_{1,2}=0$) and (b) surface viscosity-dominated regime, incompressible, non-viscous interface (${\textit {Bq}_ {1}} \to 0$, ${\textit {Bq}_ {2}} = 1000$). (a) Compressible. (b) Incompressible.

Figures 14(a) and 14(c) show drag coefficients for prolate particles $\mathcal {D}=2$ and $\mathcal {D}=5$ at different immersion in complex interfaces specified by ${\textit {Bq}_ {1}}$ and ${\textit {Bq}_ {2}}$. At very low ${\textit {Bq}_ {1,2}}$ for $\mathcal {H} > -0.8$, the drag coefficient linearly decreases with increasing $\mathcal {H}$. At high ${\textit {Bq}_ {1,2}}$ the drag coefficient is independent of $\mathcal {H}$. This trend is very similar in two prolate particles. This very weak dependency of $f$ on $\mathcal {H}$ can also be observed in figures 14(b) and 14(d), which show drag coefficient versus ${\textit {Bq}_ {1,2}}$ at $\mathcal {H}=-0.9,-0.7,-0.5,-0.3,0$. Except near full immersion ($\mathcal {H} = -0.9$), all curves coincide at high ${\textit {Bq}_ {1,2}}$.

Figure 14. Parallel drag coefficient $f$ relative to drag coefficient at a clean interface $f_c$ as a function of (a,c) immersion level $\mathcal {H}\in [-1,0]$ and (b,d) surface viscosities ${\textit {Bq}_ {1,2}}$ for prolate spheroids with (a,b) $\mathcal {D}=2$ and (c,d) $\mathcal {D}=5$ at a surfactant-free, complex (partially compressible and viscous) interface with ${\textit {Bq}_ {1}}={\textit {Bq}_ {2}}$.

3.6. Prolate particle at a surfactant-laden, partially compressible and viscous interface: Marangoni flow

When a particle translates at the interface covered with a surfactant, it changes the surfactant distribution at the interface by pushing the surfactant, resulting in a accumulation of surfactant in the front of the particle, and a depletion at its rear. This variation in the surfactant concentration causes a fluid flow from high to low concentration regions. This flow is called the Marangoni flow. The surfactant variation consequently results in a variation in interface tension causing a Marangoni force on the particle. The relaxation of the surfactant concentration variation depends on the Péclet number $\textit {Pe}_s$ (Pourali et al. 2021).

The Marangoni effects can also depend on the particle aspect ratio and immersion length. The Marangoni flow and surfactant concentration fields at three different immersion lengths $\mathcal {H}=-0.5,0,0.5$ for a prolate particle $\mathcal {D}=2$ are shown in figure 15. The figure shows Marangoni velocity, $\Delta \boldsymbol {u}_s=\boldsymbol {u}_s - \boldsymbol {u}_s^0$, surfactant concentration and contribution of the Marangoni flow in the interface stress tensor. Representative time evolutions of $\varGamma$ at the front of a spherical particle and Marangoni drag force components $F^M= \boldsymbol {F}^M \boldsymbol {\cdot } \boldsymbol {e}_x$ at various $\mathcal {H}$ are also shown in figure 16. In all simulations ${\textit {Ma}}=10$, $\textit {Pe}_s=0.5$ and ${\textit {Bq}}_{1,2}=0.1$. In the definition of the Marangoni velocity, $\boldsymbol {u}_s^0$ is the velocity at the surface viscosity-dominated regime with ${\textit {Bq}_ {1,2}}=0.1$. The surfactant concentration reaches a steady state at $t > 10$. The results also show that the more the particle sinks in the viscous phase, the less accumulation of surfactant at the front of the particle occurs. To illustrate this effect, the surfactant concentration profiles in the $x$-direction are shown in figure 17. The difference between the surfactant concentration at the front and rear of the particle, $\Delta \varGamma$, is shown in the inset plot. When the particle sinks more in the viscous fluid, the maximum accumulation and depletion occur further away from the particle surface, whereas for $\mathcal {H}>0$, it happens at the particle surface. Another aspect is a wider distribution of $\varGamma$ for $\mathcal {H} < 0$; see also the $\varGamma$ field in figure 15. For $\mathcal {H}>0$, there is a high accumulation (depletion) at the front (rear) of the particle which decays very fast with distance from the particle surface. To investigate this effect, we can use the Marangoni flow field, first column in figure 15. It shows a higher Marangoni velocity for $\mathcal {H}=0.5$ compared with $\mathcal {H}=-0.5$. This high Marangoni flow distributes surfactant easier at the interface at $\mathcal {H}=0.5$.

Figure 15. Marangoni interfacial flow field $\Delta \boldsymbol {u}_s=\boldsymbol {u}_s-\boldsymbol {u}_s^0$ (left column) and surfactant interfacial concentration field $\varGamma$ (right column) for a prolate particle in parallel motion (in vertical $x$-direction) with aspect ratio $\mathcal {D}=2$ at three different immersion levels: (a) $\mathcal {H} = -0.5$, (b) $\mathcal {H} = 0$ and (c) $\mathcal {H} = 0.5$. In all these simulations, ${\textit {Bq}_ {1,2}}=0.1$, ${\textit {Ma}}=10$ and $\textit {Pe}_s=0.5$. We used the symmetry of the system and show only half of the interface.

Figure 16. Representative time evolution of (a) surfactant concentration $\varGamma$ and (b) Marangoni force component $F^M=\boldsymbol {F}^M\boldsymbol {\cdot }\boldsymbol {e}_x$ at various $\mathcal {H}$ up to $t=30$ for a spherical particle moving in the positive $x$-direction. In all simulations ${\textit {Ma}}=10$, $\textit {Pe}_s=0.5$ and ${\textit {Bq}}_{1,2}=0.1$. At $t=30$, the stationary regime has been reached, while the particle has travelled only $3/20=15\,\%$ of the time needed to hit the box boundary.

Figure 17. Surfactant concentration profile $\varGamma$ vs relative distance $x-X$ from the sphere's centre, at various $\mathcal {H}$. The sphere is moving at constant speed in the positive $x$-direction. In all simulations ${\textit {Ma}}=10$, $\textit {Pe}_s=0.5$ and ${\textit {Bq}}_{1,2}=0.1$.

The Marangoni force on a spherical particle and a prolate particle with $\mathcal {D}=2$ is shown in figure 18(a). The results are for a simulation at ${\textit {Ma}}=10$, $\textit {Pe}_s=0.5$ and ${\textit {Bq}_ {1,2}}=0.1$. The corresponding bulk forces are presented in figure 18(b). The Marangoni force linearly decreases with the immersion length for $\mathcal {H}<0$. When the particle is more in the inviscid phase at $\mathcal {H} \approx 0.5$, the Marangoni force reaches its maximum value and decreases with a further increase in $\mathcal {H}$. The Marangoni force on a prolate particle is smaller than for a spherical particle. The parallel drag coefficient on a spherical particle is also higher than for the prolate particle (figure 19).

Figure 18. (a) Marangoni force component $F^M = \boldsymbol {F}^M \boldsymbol {\cdot } \boldsymbol {e}_x$ and (b) bulk drag force component $F^b = \boldsymbol {F}^b \boldsymbol {\cdot } \boldsymbol {e}_x$ as a function of immersion level $\mathcal {H}\in [-1,0.8]$ for a spherical ($\mathcal {D}=1$) and a prolate ($\mathcal {D}=2$) particle moving in the positive $x$-direction, coinciding with its symmetry axis, at a surfactant-laden and complex interface. In these simulations ${\textit {Bq}_ {1,2}}=0.1$, ${\textit {Ma}}=10$ and $\textit {Pe}_s=0.5$. (a) Marangoni force. (b) Bulk force.

Figure 19. Parallel drag coefficient $f$ as a function of immersion length $\mathcal {H}\in [-1,0.8]$ for a spherical ($\mathcal {D}=1$) and prolate ($\mathcal {D}=2$) particle at a surfactant-laden and complex interface. In these simulations ${\textit {Bq}_ {1,2}}=0.1$, ${\textit {Ma}}=10$ and $\textit {Pe}_s=0.5$, as for figure 18. To highlight the effect of shape at a given particle volume, the same data has also been plotted differently in Appendix C. While the particle volume is responsible for the overall shift of $f$, the aspect ratio at constant particle volume has a more pronounced influence at large $|\mathcal {H}|$.

4. Conclusion

We studied spherical and spheroidal particles of various aspect ratios at clean and surfactant-laden interfaces between a viscous fluid and air using the FEM. Prolate particles in contact with the interface are translating parallel and perpendicular to their principle axis within the interfacial plane. The immersion in the viscous fluid, specified by the contact angle, or alternatively, the dimensionless immersion length $\mathcal {H}\in [-1,1]$, had been varied. The interface is characterized by two Boussinesq numbers, ${\textit {Bq}_ {1}}$ and ${\textit {Bq}_ {2}}$. We investigated the effects of particle aspect ratio, orientation, contact angle and ${\textit {Bq}_ {1}}$ on the drag coefficient and the Marangoni surface flow field.

For a spherical particle, the drag coefficients at a compressible interface show good agreement with reported values by Danov et al. (1995), Dörr et al. (2016) and Zabarankin (2007). At an incompressible interface, we have compared our results with Fischer et al. (2006) and proposed a modified equation for the drag coefficient as a function of particle submergence. Results show that even a small immersion of the particle in the viscous fluid can alter the drag on the particle.

For both compressible and incompressible interfaces, the drag coefficient of a prolate particle, regardless of whether it is translating parallel or perpendicular to its principle axis, linearly decreases with increasing aspect ratio $\mathcal {D}$ at low ${\textit {Bq}_ {1}}$. When ${\textit {Bq}_ {1}}$ is comparable with the particle aspect ratio, the drag coefficient becomes independent of the particle size. At high ${\textit {Bq}_ {1}}$, the drag coefficient linearly increases with $\mathcal {D}$.

We also studied the Marangoni flow, at the interface covered with insoluble surfactant (Langmuir monolayer), and we believe that these results are the first to account for the particle submergence, into the viscous subphase, and aspect ratio on the Marangoni flow for the practically relevant case of an incompressible interface. For a spherical particle and a prolate particle with $\mathcal {D}=2$, we observed that when particles sink more in the bulk phase the Marangoni effects diminish. In a monolayer the drag coefficient of particles also increases with the immersion depth.

There is a wide range of data in the literature for the interface shear and dilatational viscosities. Some values have been collected in table 1 of our previous work (Pourali et al. 2021). The surface viscosity typically varies over the range $10^{-8}$ to $10^{-3}\ \textrm {Ns}\,\textrm {m}^{-1}$ (Dimova et al. 2000), it can be also as small as $10^{-10}\ \textrm {Ns}\,\textrm {m}^{-1}$ (Ortega et al. 2010) or as high as $2\ \textrm {Ns}\,\textrm {m}^{-1}$, for films stabilized by proteins (Dimova et al. 2000). If one considers the water–air interface, with the viscosity of water $\eta = 0.89 \times 10^{-3}\ \textrm {Ns}\,\textrm {m}^{-2}$, then the ratio ${\mathcal {L}}=\eta ^s / \eta$, which is the length of the spheroidal probe, $a$, times ${\textit {Bq}_ {1}}$ is typically in the range of $\sim 10^{-5}$–1 m. For a system with a typical interfacial shear viscosity of $\eta ^s\simeq 10^{-6}\ \textrm {Ns}\,\textrm {m}^{-1}$, the ${\textit {Bq}_ {1}} \simeq 10^6, 10^3$ and $1$ for spheroidal particles with $a=1\ \textrm {nm}$, $a=1\ \mathrm {\mu }\textrm {m}$ and $a=1\ \textrm {mm}$, respectively. According to Danov et al. (1995), for most practically important cases, the dilatational and shear surface viscosities exhibit the same order of magnitude. Only in extreme cases, such as some biological membranes, they can differ by several orders of magnitude. While typical values for $\eta ^s$ and $\kappa ^s$ are of the order of $10^{-6}\ \textrm {Ns}\,\textrm {m}^{-1}$, they may vary over multiple orders of magnitude as a function of surfactant concentration, or alternatively, the Marangoni number. In practice, $\eta ^s$, $\kappa ^s$ and ${\textit {Ma}}$ are therefore not independent parameters, but their relationship remains to be explored further for specific systems.

There are of course many more cases that can be explored in the seven-dimensional parameter space, spanned by ${\textit {Bq}_ {1}}$, ${\textit {Bq}_ {2}}$, $\textit {Pe}_s$, $\mathcal {D}$, $\mathcal {H}$, ${\textit {Ma}}$, and the orientation of the spheroid with respect to the particle velocity. The present selection served to provide trends and specific examples, while an approximate analytic expression for the measurable quantities remains to be developed.