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 Reference Leal2007). 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 Reference Brenner1963). 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 Reference Chisholm and Stebe2021). 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 Reference Manikantan and Squires2020). 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. Reference Pourali, Kröger, Vermant, Anderson and Jaensson2021).
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 Reference Jaensson, Anderson and Vermant2021). The drag coefficient is an important quantity in predicting and analysing interfacial rheological properties or trajectories obtained from particle tracking experiments (Bonales et al. Reference Bonales, Ritacco, Rubio, Rubio, Monroy and Ortega2007; Maestro et al. Reference Maestro, Bonales, Ritacco, Fischer, Rubio and Ortega2011). Knowledge about the motion and diffusion of an anisotropic particle at the interface is crucial for the understanding of biological systems (Ding, Warriner & Zasadzinski Reference Ding, Warriner and Zasadzinski2001), micro-organism motion (Lauga et al. Reference Lauga, DiLuzio, Whitesides and Stone2006; Masoud & Stone Reference Masoud and Stone2014; Shaik & Ardekani Reference Shaik and Ardekani2017), micrometre sized ‘Marangoni surfers’ (Dietrich et al. Reference Dietrich, Jaensson, Buttinoni, Volpe and Isa2020), and inclusions such as proteins, or other membrane-bound particles, in biological membranes, and the design of artificial membranes (Ally & Amirfazli Reference Ally and Amirfazli2010). Small particle probes also give microrheological methods increased sensitivity compared with macroscopic methods (Samaniuk & Vermant Reference Samaniuk and Vermant2014).
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 (Reference Saffman and Delbrück1975). 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 Reference Saffman and Delbrück1975; Saffman Reference Saffman1976),
\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 Reference Saffman1976; Dimova et al. Reference Dimova, Danov, Pouligny and Ivanov2000):
\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. (Reference Gambin, Lopez-Esparza, Reffay, Sierecki, Gov, Genest, Hodges and Urbach2006) 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 (Reference Levine and MacKintosh2002) and Levine, Liverpool & MacKintosh (Reference Levine, Liverpool and MacKintosh2004). 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,
\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 Reference Kirkwood and Auer1951; Klopp, Stannarius & Eremin Reference Klopp, Stannarius and Eremin2017). 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 (Reference Fischer2004b) 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. (Reference Danov, Aust, Durst and Lange1995) 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. (Reference Dimova, Danov, Pouligny and Ivanov2000) 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 (Reference Fischer, Dhar and Heinig2006) used an approach different from Levine et al. (Reference Levine, Liverpool and MacKintosh2004), 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 (Reference Stone and Masoud2015) 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. Reference Pourali, Kröger, Vermant, Anderson and Jaensson2021) 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 Reference Klingler and McConnell1993; Steffen et al. Reference Steffen, Heinig, Wurlitzer, Khattari and Fischer2001; Wurlitzer, Schmiedel & Fischer Reference Wurlitzer, Schmiedel and Fischer2002; Fischer Reference Fischer2004a; Fischer et al. Reference Fischer, Dhar and Heinig2006; Manikantan & Squires Reference Manikantan and Squires2020), 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. Reference Pourali, Kröger, Vermant, Anderson and Jaensson2021). 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. Reference Pourali, Kröger, Vermant, Anderson and Jaensson2021). 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. Reference Peppicelli, Jaensson, Tregouet, Schroyen, Alicke, Tervoort, Monteux and Vermant2019). 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 Reference Samaniuk and Vermant2014), 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
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\rm \pi}=\boldsymbol{0}, \end{gather}
\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 Reference Brenner1991; Jaensson & Vermant Reference Jaensson and Vermant2018; Venerus & Öttinger Reference Venerus and Öttinger2018). It is assumed to be given by the Boussinesq–Scriven constitutive law (Boussinesq Reference Boussinesq1913; Scriven Reference Scriven1960)
\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 Reference Brenner1991).
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 Reference Venerus and Öttinger2018), 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,
\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 Reference Brenner and Leal1978; Stone Reference Stone1990; Brenner Reference Brenner1991)
\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
\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
\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
\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 Reference Lopez and Hirsa2000; Manikantan & Squires Reference Manikantan and Squires2020). 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 Reference Scriven1960; Ortega, Ritacco & Rubio Reference Ortega, Ritacco and Rubio2010). 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 Reference Earnshaw1981).
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.
\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.
\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
\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,
\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
\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. (Reference Pourali, Kröger, Vermant, Anderson and Jaensson2021), 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. Reference Pourali, Kröger, Vermant, Anderson and Jaensson2021). 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 Reference Hu, Patankar and Zhu2001). 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
\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 Reference Carrozza, Hulsen and Anderson2020; Balemans, Hulsen & Anderson Reference Balemans, Hulsen and Anderson2016; Jaensson, Hulsen & Anderson Reference Jaensson, Hulsen and Anderson2017)
\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}
\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 Reference Donea and Huerta2003).
The weak form of the SCD equation amounts to find 
$\varGamma$ such that
\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 Reference Taylor and Hood1973) are used for 
$\boldsymbol {u}$ and 
$p$ whereas triangular 
$P_2$ elements are used for 
$\varGamma$. Mesh generation is done using Gmsh (Geuzaine & Remacle Reference Geuzaine and Remacle2009), 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:
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} (K_e \boldsymbol{\nabla} (\Delta x_m)) = 0 \quad (\boldsymbol{x}_m \in \varOmega), \end{gather}
\begin{gather}\Delta x_m = 0 \quad (\boldsymbol{x}_m \in \partial \varOmega), \end{gather}
\begin{gather}\Delta x_m = \Delta X \quad (\boldsymbol{x}_m \in S_\text{p}), \end{gather}
\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. Reference Hu, Patankar and Zhu2001). 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,
\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$,
\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
\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}$,
$$\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 (Reference Zanini and Isa2016), 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. Reference Davies, Krüger, Coveney and Harting2014)
\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
\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
\begin{equation} \cos \theta ={-}\mathcal{H}. \end{equation}Using this definition of the contact angle, the energy of the equilibrium configuration is given by
\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
\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 Reference Ranger1978). With some additional assumptions, the drag coefficient of a particle at an interface has been quantitatively evaluated by hydrodynamic calculations. Zabarankin (Reference Zabarankin2007) 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. (Reference Dörr, Hardt, Masoud and Stone2016) and Villa et al. (Reference Villa, Boniello, Stocco and Nobili2020) for hydrophilic particles, 
$\theta < 90^{\circ }$ (
$\mathcal {H}\in [-1,0]$),
\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]$,
\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
\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. (Reference Danov, Aust, Durst and Lange1995) and Zabarankin (Reference Zabarankin2007), 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. (Reference Petkov, Denkov, Danov, Velev, Aust and Durst1995) 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. Reference Petkov, Denkov, Danov, Velev, Aust and Durst1995) 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. Reference Vassileva, van den Ende, Mugele and Mellema2005; Madivala, Fransaer & Vermant Reference Madivala, Fransaer and Vermant2009). 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. Reference Loudet, Qiu, Hemauer and Feng2020).
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. (Reference Danov, Aust, Durst and Lange1995) 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. (Reference Fischer, Dhar and Heinig2006) 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,
\begin{equation} f = f_0 + f_1({\textit{Bq}_{1}}) + {{O}}[({\textit{Bq}_{1}})^2]. \end{equation}
Fitted expressions (Fischer et al. Reference Fischer, Dhar and Heinig2006) for 
$f$ from numerical results gave the formulae for 
$f_0$ and 
$f_1$,
\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
\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. Reference Pourali, Kröger, Vermant, Anderson and Jaensson2021). 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. Reference Pourali, Kröger, Vermant, Anderson and Jaensson2021). 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 Reference Bławzdziewicz, Cristini and Loewenberg1999). It also appears that, under typical circumstances, surface diffusivity of the surfactant is insufficient to relax interfacial incompressibility (Chisholm & Stebe Reference Chisholm and Stebe2021).
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. (Reference Fischer, Dhar and Heinig2006), 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. (Reference Fischer, Dhar and Heinig2006). 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. (Reference Fischer, Dhar and Heinig2006), as already noted by Pourali et al. (Reference Pourali, Kröger, Vermant, Anderson and Jaensson2021), 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
\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 (Reference Happel and Brenner1981) 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
\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
\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
\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
\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
\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 (Reference Happel and Brenner1981) and that we denoted by 
$a$ (our length unit) the length of the axis of rotational symmetry, while Happel & Brenner (Reference Happel and Brenner1981) 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.
\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. Reference Pourali, Kröger, Vermant, Anderson and Jaensson2021). 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. Reference Pourali, Kröger, Vermant, Anderson and Jaensson2021).
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. (Reference Danov, Aust, Durst and Lange1995), Dörr et al. (Reference Dörr, Hardt, Masoud and Stone2016) and Zabarankin (Reference Zabarankin2007). At an incompressible interface, we have compared our results with Fischer et al. (Reference Fischer, Dhar and Heinig2006) 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. Reference Pourali, Kröger, Vermant, Anderson and Jaensson2021). The surface viscosity typically varies over the range 
$10^{-8}$ to 
$10^{-3}\ \textrm {Ns}\,\textrm {m}^{-1}$ (Dimova et al. Reference Dimova, Danov, Pouligny and Ivanov2000), it can be also as small as 
$10^{-10}\ \textrm {Ns}\,\textrm {m}^{-1}$ (Ortega et al. Reference Ortega, Ritacco and Rubio2010) or as high as 
$2\ \textrm {Ns}\,\textrm {m}^{-1}$, for films stabilized by proteins (Dimova et al. Reference Dimova, Danov, Pouligny and Ivanov2000). 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. (Reference Danov, Aust, Durst and Lange1995), 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.
Acknowledgements
The authors thank M.A. Hulsen at the Eindhoven University of Technology (TU/e) for access to the TFEM software libraries.
Funding
M.K. would like to acknowledge computing resources at the Swiss National Computing Center (CSCS project s987) and support from the Swiss National Science Foundation (grant number 200021L_185052).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available within this article.
Appendix A. Equilibrium position of the spheroid
To generalize the equations presented in § 3.1 to a spheroid with half-axes 
$a$ and 
$b$, we begin with the parametric representation 
$\boldsymbol {r}(\phi ,\xi )$ of the spheroidal surface, aligned and positioned as shown in figure 1. In dimensionless form (length unit 
$a$), and upon making use of our abbreviations 
$\mathcal {H}=h/b\in [-1,1]$ and 
$\mathcal {D}=a/b\in [0,\infty ]$ (2.10a,b), the parametric form of our spheroid surface centred at 
$y=h$ (dimensionless form 
$y=h/a=\mathcal {H}/\mathcal {D}$) is given by
\begin{equation} \boldsymbol{r} = \sqrt{1-\xi^2}\cos\phi \,\boldsymbol{e}_x + \frac{\mathcal{H} +\sqrt{1-\xi^2}\sin\phi}{\mathcal{D}}\,\boldsymbol{e}_y + \frac{\xi}{\mathcal{D}} \,\boldsymbol{e}_z, \end{equation}
where 
$\boldsymbol {e}_\mu$ denotes a Cartesian base unit vector, and 
$\phi$ and 
$\xi$ vary in the range 
$\phi \in [0,2{\rm \pi} ]$ and 
$\xi \in [-1,1]$ for a full spheroid. With the help of
\begin{equation} \left|\frac{\partial \boldsymbol{r}}{\partial \phi} \times \frac{\partial\boldsymbol{r}}{\partial \xi}\right| = \frac{1}{\mathcal{D}^2}\sqrt{(1-\xi^2)\cos^2\phi + \mathcal{D}^2 [\xi^2+(1-\xi^2)\sin^2\phi]}\end{equation}the dimensionless surface area of the full (prolate or also oblate) spheroid is
\begin{equation} S(\mathcal{D}) = \int_0^{2{\rm \pi}} \int_{{-}1}^1 \left|\frac{\partial \boldsymbol{r}}{\partial \phi} \times \frac{\partial\boldsymbol{r}}{\partial \xi}\right| \,\textrm{d}\xi\,\textrm{d}\phi = \frac{2{\rm \pi}}{\mathcal{D}^2}\left[1 + \frac{\mathcal{D} \sinh^{{-}1}\sqrt{\mathcal{D}^{{-}2}-1}}{\sqrt{\mathcal{D}^{{-}2}-1}}\right]. \end{equation}
For the sphere (
$\mathcal {D}=1$), this evaluates to 
$S(1)=4{\rm \pi}$. Faraudo & Bresme (Reference Faraudo and Bresme2003) wrote this differently starting from the normal vector form of the spheroid. For the energy considerations, we are interested in the parts of the spheroidal surface exposed to air and fluid, 
$S_a$ and 
$S_f$, respectively. They are given by
\begin{equation} S_a(\mathcal{D},\mathcal{H}) = \int_0^{2{\rm \pi}} \int_{{-}1}^1 \left|\frac{\partial \boldsymbol{r}}{\partial \phi} \times \frac{\partial\boldsymbol{r}}{\partial \xi}\right| \varTheta\left(\mathcal{H}+\sqrt{1-\xi^2}\sin\phi\right) \,\textrm{d}\xi\,\textrm{d}\phi \end{equation}
and 
$S_f=S-S_a$, where 
$\varTheta (y)$ denotes the Heaviside step function, 
$H(y)=1$ for 
$y>0$ and 
$H(y)=0$ otherwise. For a sphere, because the last term in (A2) is linear in 
$\xi$, the 
$S_a$ is linear in 
$\mathcal {H}$, more specifically, 
$S_a(1,\mathcal {H})=2{\rm \pi} (1+\mathcal {H})$ and, thus, 
$S_f(1,\mathcal {H})=4{\rm \pi} -S_a=2{\rm \pi} (1-\mathcal {H})$. For a spheroid, the integral (A4) cannot be evaluated analytically, but it is most conveniently evaluated to arbitrary precision using 
$N\gg 1$ pairs 
$(\zeta _1,\zeta _2$) of independent random numbers equally distributed on the interval 
$[0,1]$. Each pair is used to calculate 
$\phi =2{\rm \pi} \zeta _1$, 
$\xi =2\zeta _2-1$ and the corresponding value of the integrand in (A4) using (A2). The area 
$S_a$ is then just the arithmetic mean of the 
$N$ integrands, multiplied by 
$4{\rm \pi}$.
In addition, we need the circular area in the interfacial plane that is occupied by the spheroid. As the cross-section of the spheroid within the 
$y=0$ plane is an ellipse with half-axes 
$\sqrt {a^2-h^2}$ and 
$\sqrt {b^2-h^2}$, we obtain the dimensionless
\begin{equation} S_c(\mathcal{D},\mathcal{H}) = {\rm \pi}\frac{\sqrt{(\mathcal{D}^2-\mathcal{H}^2)(1-\mathcal{H}^2)}}{\mathcal{D}^2}. \end{equation}
The total dimensionless energy of this system can then be expressed, following our arguments from § 3.1, and with the help of the three surface tensions and the three surface areas 
$S_a$, 
$S_f$ and 
$S_c$, as
\begin{align} E(\mathcal{D},\mathcal{H}) &= S_a(\mathcal{D},\mathcal{H}) \gamma_{pa} + S_f(\mathcal{D},\mathcal{H}) \gamma_{pf} - S_c(\mathcal{D},\mathcal{H}) \gamma \nonumber\\ &= S_a(\mathcal{D},\mathcal{H}) (\gamma_{pa}-\gamma_{pf}) - S_c(\mathcal{D},\mathcal{H}) \gamma + S(\mathcal{D}) \gamma_{pf}. \end{align}For the sphere, this immediately evaluates to
\begin{equation} E(1,\mathcal{H}) = 2{\rm \pi} (1+\mathcal{H})(\gamma_{pa}-\gamma_{pf}) - {\rm \pi}(1-\mathcal{H}^2) \gamma + 4{\rm \pi} \gamma_{pf} \end{equation}
and, thus, exactly coincides with (3.1) divided by 
$R^2$, i.e. with the non-dimensionalized (3.1). For the spheroid, we have to use the above expressions for areas. The equilibrium position of the spheroid's centre of mass minimizes the total energy 
$E$ given by (A6). While 
$dS_c/d\mathcal {H}$ can be written down analytically, 
$S_a$ depends on 
$\mathcal {H}$ only through the Heaviside function 
$\varTheta (x)$. Because 
$\textrm {d}\varTheta (\mathcal {H}+\ldots )/\textrm {d}\mathcal {H}=\delta (\mathcal {H}+\ldots )$ is Dirac's delta distribution, the equilibrium position 
$\mathcal {H}$ minimizing the energy can be found at minor numerical effort. Moreover, according to (A6) and because 
$S(\mathcal {D})$ does not depend on 
$\mathcal {H}$, the reduced equilibrium position 
$\mathcal {H}$ (figure 20) is a function of two variables, aspect ratio 
$\mathcal {D}$ and the dimensionless 
$(\gamma _{pa}-\gamma _{pf})/\gamma$, not only for a sphere but for the spheroid as well.
Figure 20. (a) Reduced equilibrium position 
$\mathcal {H}$ of the spheroid as a function of the dimensionless difference 
$(\gamma _{pa}-\gamma _{pf})/\gamma$ of surface tensions for various aspect ratios 
$\mathcal {D}$. (b) To highlight the differences between the curves in (a), panel (b) shows 
$\mathcal {H}\gamma /(\gamma _{pa}-\gamma _{pf})$. For the sphere, 
$\mathcal {H}=-(\gamma _{pa}-\gamma _{pf})/\gamma$ in agreement with (3.2). Each integral has been estimated using 
$N=10^6$ evaluations.
Appendix B. Implicit scheme
We have used the explicit-ALE scheme introduced in § 2.6.2 to solve the weak form of the governing equations. An alternative approach is to use an implicit scheme. In the implicit time integration scheme, the weak forms of the governing equations, (2.15)–(2.17), are solved together in one system. The same time integration schemes for 
$\partial \varGamma /\partial t$ as for the explicit scheme are used. Due to the linearity of the equation of state used here, the only nonlinear terms appear in the SCD, which are handled by a Picard iteration,
$$\begin{gather} \cdots+[ (\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\\ \approx \cdots+[ (\boldsymbol{u}_s)_{n+1} - (\boldsymbol{u}_m)_{n+1} ] \boldsymbol{\cdot} \boldsymbol{\nabla}_s\varGamma^i_{n+1} + [\boldsymbol{\nabla}_s\boldsymbol{\cdot} (\boldsymbol{u}_s)_{n+1}]\varGamma^i_{n+1}, \end{gather}$$
where 
$\varGamma ^i_{n+1}$ is the surface concentration from the previous step in the Picard iteration, with the initial guess given by 
$\varGamma ^0_{n+1}=\varGamma _n$. The iteration is terminated when the normalized difference between iteration steps for the velocity magnitude and the surface concentration has become smaller than 
$10^{-6}$.
In figure 21(a) we have reproduced the case of figure 17 using the implicit scheme, with and without ALE formulation. While in figure 21(a) (without ALE) the mesh does not move with the particle, in figure 21(b) (with ALE), similar to the explicit-ALE method used in the manuscript, the mesh moves with the particle using the ALE formulation. This exemplary comparison shows that results are basically insensitive to the choice of scheme, apart from minor deviations for large 
$\mathcal {H}$ when the particle is predominantly exposed to air. In figure 22 we additionally compare the two implicit schemes at a different, higher 
${\textit {Ma}}=100$. Again, there are no notable differences between the two implicit methods.
Figure 21. (a) Implicit and (b) explicit-ALE simulations of the surfactant concentration profile 
$\varGamma$ at various 
$\mathcal {H}$ of a spherical particle. Figure 17 in the manuscript is the corresponding simulation result implemented with the explicit method. In all simulations 
${\textit {Ma}}=10$, 
$\textit {Pe}_s=0.5$ and 
${\textit {Bq}}_{1,2}=0.1$.
Figure 22. Same as figure 21 for a larger 
${\textit {Ma}}=100$.
Appendix C. Additional data
For the set-up of a spheroidal particle in parallel motion, for which fields have been presented in figure 9, we furthermore performed simulations in perpendicular motion. Results are shown in figure 23. As explained within § 3.4 and figures of the manuscript, our particle volumes change upon varying the aspect ratio 
$\mathcal {D}$, because the length of the 
$a$-axis is used to introduce dimensionless quantities. To appreciate the effect of particle shape on the dimensionless drag coefficient, or equivalently, the dimensionless drag force, at constant particle volume, we provide two graphs highlighting the transformation issue. Both are shown in figure 24.
Figure 23. Same as figure 9 for the case of perpendicular motion, i.e. interfacial flow field 
$\boldsymbol {u}_s$, interface velocity components, 
$u_{s,x}$ and 
$u_{s,z}$, interface compression 
$\boldsymbol {\nabla }_s\boldsymbol {\cdot }\boldsymbol {u}_s$ and shear component for a prolate spheroid in perpendicular motion (axis of rotational symmetry aligned in horizontal 
$z$-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 a surfactant-free, 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.
Open access
