Dynamics of mass polar spheroids during sedimentation

Abstract

The dynamics of sedimenting particles under gravity are surprisingly complex due to the presence of effective long-ranged forces. When the particles are polar with a well-defined symmetry axis and non-uniform density, recent theoretical predictions suggest that prolate objects will repel and oblate ones will weakly attract. We tested these predictions using mass polar prolate spheroids, which are composed of 2 mm spheres glued together. We probed different aspect ratios ($\kappa$) and centre of mass variations ($\chi$) by combining spheres of different densities. Experiments were done in both quasi-two-dimensional (2-D) and three-dimensional (3-D) chambers. By optically tracking the motion of single particles, we found that the dynamics were well described by a reduced mobility matrix model that could be solved analytically. Pairs of particles exhibited an effective repulsion, and their separation roughly scaled as $(\kappa - 1)/\chi ^{0.39}$, i.e. particles that were more prolate or had smaller mass asymmetry had stronger repulsion effects. In three dimensions, particles with $\chi >0$ were distributed more uniformly than $\chi =0$ particles, and the degree of uniformity increased with $\kappa$, indicating that the effective 2-body repulsion manifests for a large number of particles.

1. Introduction

Sedimentation is a longstanding and important problem in fluid dynamics. In its simplest form, particles far from equilibrium settle in a fluid through some external forcing, typically gravity, at low Reynolds number (Stokes 1851). Throughout its storied history, one can observe a microcosm of physics problems that span multiple fields. Starting from basic hydrodynamics, the long range velocity fields generated by sedimenting particles lead to several interesting phenomena (Stokes 1851; Brady & Bossis 1988; Xue et al. 1992; Ramaswamy 2001; Guazzelli, Morris & Pic 2011). Examples include unbounded velocity fluctuations (Caflisch & Luke 1985), chaotic behaviour (Brady & Bossis 1988; Jánosi et al. 1997) and periodic orbits (Claeys & Brady 1993; Ekiel-Jeżewska & Felderhof 2005; Jung et al. 2006; Chajwa, Menon & Ramaswamy 2019). Sedimentation is found throughout nature; from silt and sand in a river to biogenic particles in the ocean (Monroy et al. 2019). Most sedimentation work has been done on uniform particles or particles with simple symmetries. However, within nature, most particles are not uniform. They can be rough and polygonal, and they can be made of many different materials, causing their mass to be distributed non-uniformly (Domokos et al. 2020). For example, it has been found that some phytoplankton adjust their centre of mass to respond to external environmental flows for better survival in turbulent environments (Sengupta, Carrara & Stocker 2017).

Gravitational sedimentation at low Reynolds number (Stokes flow) is a special case of the Navier–Stokes equation where inertia is negligible. Because of this, Stokes flow is quasistatic and time reversible. For a single spherical particle of radius $R$ and density $\rho _p$ settling in an unbounded fluid of density $\rho _f$ and viscosity $\eta$, balancing the Stokes drag force with gravitational and buoyant forces leads to the following expression for the steady state terminal velocity:

(1.1)\begin{equation} U_T=\frac{2}{9}\frac{\rho_p-\rho_f}{\eta}g R^2.\end{equation}

Here, $g$ is the gravitational acceleration. The addition of many other particles in the fluid complicates this picture. To leading order, the fluid disturbance at a distance $r$ from a sedimenting sphere with velocity $U_{s}$ and radius $R$ scales as $U_{s}R/r$. In sedimenting suspensions of many particles, these long range hydrodynamic interactions complicate a local description of particle dynamics. Batchelor solved the problem of a diverging mean sedimentation velocity (Batchelor 1972), but Caflisch and Luke pointed out that the velocity fluctuations were still unbounded as the system size increases (Caflisch & Luke 1985).

To illustrate the Caflisch–Luke paradox, consider the variance of the sedimentation velocity of a group of $N$ particles contained in a volume of size $L$. The volume fraction $\phi$ of particles is $NV_{p}/L^{3}$, where $V_{p} = \frac {4}{3}{\rm \pi} R^{3}$ is the volume of a single particle. Within this region, if the particles are randomly and independently distributed, the fluctuation in particle number is simply $\sqrt {N}$. To find the velocity fluctuations, we can balance the total change in the Stokes’ drag force over the suspension with the change in gravitational and buoyant forces due to these number fluctuations: $6{\rm \pi} \eta L{\rm \Delta} v \approx (\rho _p-\rho _f)V_pg\sqrt {N}$. Solving for ${\rm \Delta} v$, we arrive at the fractional change in velocity, ${\rm \Delta} v/v_{o} = L^{1/2}\sqrt {\phi R^{2}/V_{p}}$. This would indicate that the velocity fluctuations depend on the system size, $L$. Simulations agree with these predictions in unbounded fluids (Koch 1993, 1994; Ladd 1996, 1997; Cunha et al. 2002; Mucha et al. 2004), while experiments generally observe a limit to the size of the fluctuations (Ham & Homsy 1988; Xue et al. 1992; Nicolai & Guazzelli 1995; Segrè, Herbolzheimer & Chaikin 1997).

To reconcile this paradox, several different physical mechanisms have been proposed. The long-ranged interactions must be screened out by some large length scale, or by changing the interactions themselves. For example, wall effects at the size of the experimental container (Brenner 1999), correlated particle positions arising from a pre-imposed structure factor (Koch & Shaqfeh 1989, 1991), polydispersity (Nguyen & Ladd 2005), stochasticity in the concentration (Levine et al. 1998), stratification (Mucha et al. 2004) or shape effects (Doi & Makino 2005; Krapf, Witten & Keim 2009; Goldfriend, Diamant & Witten 2017; Palusa et al. 2018; Chajwa et al. 2019, 2020; Witten & Diamant 2020). The latter example is of particular interest since it is a local change to particle interactions. Shape effects can be captured within Stokes flow using a response matrix that only depends on particle geometry and couples to external forces and torques.

We start by considering the Navier–Stokes equation for an incompressible fluid in the low Reynolds number regime:

(1.2)$$\begin{gather} \boldsymbol{\nabla} P = \eta\boldsymbol{\nabla}^2 \boldsymbol{v} + \boldsymbol{f}_{b}, \end{gather}$$

(1.3)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{v} = 0, \end{gather}$$

where $P$ is the pressure, $\eta$ is the dynamic viscosity, $\boldsymbol {v}$ is the velocity field and $\boldsymbol {f}_{b}$ are any body forces per unit volume on the fluid, such as gravity. The linearity of these equations allows us to write the equations of motion for a single particle suspended in the fluid and subjected to an external force or torque as

(1.4)$$\begin{gather} \boldsymbol{v}(t) = \boldsymbol{\mathsf{T}}_{vF}\boldsymbol{\cdot}\boldsymbol{F} + \boldsymbol{\mathsf{T}}_{v\tau}\boldsymbol{\cdot}\boldsymbol{\tau}, \end{gather}$$

(1.5)$$\begin{gather}\boldsymbol{\omega}(t) = \boldsymbol{\mathsf{T}}_{\omega F}\boldsymbol{\cdot}\boldsymbol{F} + \boldsymbol{\mathsf{T}}_{\omega \tau}\boldsymbol{\cdot}\boldsymbol{\tau}, \end{gather}$$

which can be written in matrix form as

(1.6)\begin{equation} \begin{pmatrix} \boldsymbol{v} \\ \boldsymbol{\omega} \end{pmatrix} = \begin{pmatrix} \boldsymbol{\mathsf{T}}_{vF} & \boldsymbol{\mathsf{T}}_{v\tau}\\ \boldsymbol{\mathsf{T}}_{\omega F} & \boldsymbol{\mathsf{T}}_{\omega\tau} \end{pmatrix} \begin{pmatrix} \boldsymbol{F} \\ \boldsymbol{\tau} \end{pmatrix}.\end{equation}

Here, $\boldsymbol {\omega }$ is the angular velocity of rotation about the centre of geometry, and $\boldsymbol {F}$ and $\boldsymbol {\tau }$ are the external forces and torques, respectively. The convention we use is the same as Witten & Diamant (2020). The shape dependent $\boldsymbol{\mathsf{T}}$ matrices couple the velocities of the particle to external forces and torques. In the fixed lab frame, the matrices depend on the particle's orientation to the imposed flow. We can also put restrictions on the matrices by physical insight. The dissipated power of the object, $\boldsymbol {F}\boldsymbol {{\cdot }}\boldsymbol {v} + \boldsymbol {\tau }\boldsymbol {{\cdot }}\boldsymbol {\omega }$, must be positive, which implies the diagonal blocks, $\boldsymbol{\mathsf{T}}_{vF}$ and $\boldsymbol{\mathsf{T}}_{\omega \tau }$ must be symmetric, and $\boldsymbol{\mathsf{T}}_{v\tau }$ and $\boldsymbol{\mathsf{T}}_{\omega F}$ must be transposes of each other but not necessarily positive or symmetric. Taken together, these matrices comprise the mobility matrix $\boldsymbol{\mathsf{T}}$ of an object. If you invert the relation, the matrix is called the resistance matrix. As an illustration, for a uniform sphere in an unbounded fluid, the mobility matrix is

(1.7)\begin{equation} \begin{pmatrix} \boldsymbol{v} \\ \boldsymbol{\omega} \end{pmatrix} = \begin{pmatrix} \dfrac{1}{6{\rm \pi}\eta R}\delta_{ij} & 0\\ 0 & \dfrac{1}{8{\rm \pi}\eta R^{3}}\delta_{ij} \end{pmatrix} \begin{pmatrix} \boldsymbol{F} \\ \boldsymbol{\tau} \end{pmatrix},\end{equation}

where $\delta _{ij}$ is the Kronecker delta.

The dynamics of a single particle are determined by the time evolution of $\boldsymbol{\mathsf{T}}$. As the particle moves through the fluid, its orientation can change with respect to the centre of mass velocity. The orientation of the particle relative to the force determines what $\boldsymbol{\mathsf{T}}$ looks like in the lab frame. Analogously, if we move to the body frame of the particle, $\boldsymbol{\mathsf{T}}$ becomes fixed and the force and torque become time dependent. The motion of the particle cannot change the magnitude of the force, so only the force's direction changes with time. Depending on the symmetries of $\boldsymbol{\mathsf{T}}$, different classes of trajectories can be found. For a comprehensive list of these trajectories and symmetries, refer to Doi & Makino (2005), Krapf et al. (2009) and Witten & Diamant (2020).

In the case of gravitational sedimentation, asymmetric particles with mass distribution polarity will undergo rotation in response to external forcing (Witten & Diamant 2020). This is because the total form and skin drag on the particle can apply a net torque when the centre of mass is in a different location than the geometric centre of the particle. Consequently, an external force leads to a net torque, and the particle will rotate so that the external force is parallel to an eigendirection of $\boldsymbol{\mathsf{T}}_{\omega F}$ (Witten & Diamant 2020). The response of a single particle can have important implications for the sedimentation dynamics of many particles. Recent work has theoretically explored the sedimentation of ‘mass polar’ prolate spheroids, whose centre of mass lies along the major axis away from the geometric centre (Goldfriend et al. 2017). These particles are defined by two parameters: the ratio of major to minor axes, $\kappa$, and the centre of mass offset from the geometric centre, $\chi$. Using a linear stability analysis of a uniform suspension of particles in Stokes flow, they predicted a repulsive interaction for $\kappa >1$ (prolate), and an attractive interaction for $\kappa <1$ (oblate). The effect is surprisingly enhanced for smaller values of $\chi$. These effects, over a large collection of particles, can either enhance particle clustering and velocity fluctuations ($\kappa <1$), or inhibit them ($\kappa >1$).

Inspired by Goldfriend, Diamant & Witten (2015, 2016) and Goldfriend et al. (2017), we experimentally tested these predictions by fabricating prolate, mass polar ‘dimers’ and ‘trimers’. The particles were composed of multiple spheres of varying materials bonded together. Our experiments tracked the position and rotation of pairs of particles in a quasi-2-D environment. First, we examined the motion of single particles to quantify the mobility matrix. Using the symmetry properties of prolate spheroidal particles, we derive an analytic solution for the particle dynamics that shows excellent agreement with the experimental data. Then, by sedimenting pairs of particles in the same quasi-2-D environment, we found that prolate particles experienced an effective repulsion that increased with $\kappa$ and decreased with $\chi$, in agreement with Goldfriend et al. (2017). Finally, we sedimented hundreds of particles in a 3-D container and analysed the distribution of their post-sedimented positions. The inherent repulsion manifested as wider spatial distributions of particles on the floor of the experimental apparatus. This shows local changes in particle interactions have a large effect on global sedimentation patterns.

2. Experimental methods and particle fabrication

Composite particles were fabricated by gluing together smooth ball bearings using a cyanoacrylate based glue. Each sphere had a diameter of 2 mm, and the material and mass density of each sphere were chosen to produce various numerical values of $\chi$. We used the minimal amount of glue possible to adhere the spheres by applying a low-viscosity glue instead of a viscous glue. The remaining thin layer of glue that extended away from the contact point possibly affected the motion of the sedimentation of the particles, but the repeatability of the experiments indicated that this has only a minimal effect. The materials used were aluminium, stainless steel, copper, tungsten carbide, zirconium dioxide and Delrin. Spheres were glued in either a dimer ($\kappa = 2$) or linear trimer ($\kappa = 3$) configuration. The accessible range of $\chi$ was 0.0–0.43. To analytically calculate $\chi$ for any linear chain of $n$ spherical particles, we assumed all particles were ‘light’ with density $\rho _l$ except for a single ‘heavy’ particle with density $\rho _h$ positioned at the end of the chain. The result is

(2.1)\begin{equation} \chi = \frac{1}{n}\frac{(n-1)|\rho_h - \rho_l|}{\rho_h + (n-1)\rho_l }, \quad n\geq 2. \end{equation}

The centre of mass is displaced by a distance $\kappa \chi R$, for a physical representation of $\kappa$ and $\chi$, see figure 1.

Figure 1. (a) Schematic diagram of our quasi-2-D experimental setup. The tank dimensions are $19\ {\rm cm} \times 15\ {\rm cm} \times 0.4\ {\rm cm}$. The top of the tank has a gating mechanism that allows us to drop multiple particles simultaneously. The mechanism consists of a slotted piece of acrylic and a metal rod in a U shape. By moving the prongs of the rod, the horizontal part can be rotated out of the plane, releasing the particles simultaneously. (b) Schematic depicting a $\kappa =2$ composite particle and coordinates in the lab frame. Here, $\theta$ is defined as the angle between the composite particle's major axis and the vertical direction. The orange sphere has a larger mass density in this case, so the centre of mass is shifted away from the centre of geometry (black cross) to the position indicated by the black cross. (c) Physical representation of $\kappa$ and $\chi$ shown on an example particle. The centre of mass of the particle is offset by an amount $\kappa \chi R$. For all our experiments, $R=1$ mm and the typical Reynolds number is $\sim 10^{-4}$.

Two sets of experiments used a quasi-2-D tank made out of cast acrylic (figure 1). We laser cut sheets of cast acrylic and used SCIGRIP 4 acrylic plastic cement to glue them together to create a tank of dimensions 19 cm high, 15 cm wide, with a gap of thickness 4 mm. The tank was filled with pure silicone oil of kinematic viscosity 10 000 cSt and density of 0.971 g cm$^{-3}$. A gating mechanism was placed at the top of the chamber consisting of a thin rectangle of acrylic with 2.5 mm holes spaced out evenly. The holes helped to align the particles so that the initial orientations are fixed before sedimentation. A thin metal rod held them in place and facilitated a simultaneous release of the particles at the beginning of an experimental run.

After the particles were released, we imaged their sedimentation using a CCD camera (Point Grey) at 6 frames per second with a spatial resolution of 12 pixels per mm. After recording, we processed the images using ImageJ (Schindelin et al. 2012) for easier detection of each sphere in a composite particle. Images were first binarized with a brightness threshold, then each sphere was separated with a watershedding algorithm. The resulting image was eroded, leaving us with easily trackable objects composed of white pixels. Particle tracking and linking between frames were done with TrackPy (Allan et al. 2021). The resulting trajectories of the individual spheres were used to calculate various quantities associated with the dynamics of the composite particles.

The second set of experiments were done in a cylindrical 3-D chamber of diameter of 12 cm and a height of 21 cm (see § 4). The chamber was fabricated from a cast acrylic tube with wall thickness of 12 mm. The chamber was also filled with silicone oil of the same viscosity (10 000 cSt). We placed 100 particles of a single $\kappa$ and $\chi$ combination in the fluid and sealed the chamber so that there were no trapped air bubbles. Particles were allowed to sediment under gravity to the bottom of the chamber, and the distribution of particles was imaged from above. We then flipped the chamber and repeated the experiment 50 times for each set of particles. Due to finite-size wall effects driving convection and particles resting on top of one another, identifying the individual spheres from each particle was not feasible, as was done in the 2-D experiments. Thus, images were cropped and binarized and the spatial distributions of black pixels were analysed.

The quasi-2-D geometry allows us to easily track the position and rotation of particles, but it also imposes a form of screening for the interactions between particles. The divergence of velocity fluctuations in suspensions arises from the $1/r$ decay of velocity around a sedimenting particle; however, in confined 2-D environments, the fluid flow decays as $1/r^2$. A detailed discussion of the differences can be found in Beatus, Bar-Ziv & Tlusty (2012). The faster decay allows convergence of the velocity fluctuations found in three dimensions, meaning that the majority of the screening is provided by the confining walls of our chamber. Although this is important for a statistically large number of particles, our results show that mass polarity strongly affects sedimentation dynamics in both 2-D and 3-D geometries.

3. Results and discussion

3.1. Single particle dynamics

After fabricating the composite, prolate particles, we observed the sedimentation of single, isolated particles to better understand their dynamics and to extract the terms in the mobility matrix (1.6). The response of a single particle to an external force or torque informs its effective interactions with neighbouring particles (Goldfriend et al. 2017; Witten & Diamant 2020). For example, a rod-shaped particle of uniform mass density will sediment without a change in its initial angle (Ramaswamy 2001; Witten & Diamant 2020). This results in a diagonal drift. However, the mass polarity of our objects causes them to align with the external gravitational field, meaning that a mass polar object will rotate until its centre of geometry lies directly above its centre of mass ($\theta = 0$). For our experiments, mass polar particles were released from an initial angle of $\theta = {\rm \pi}$, so that they rotated a total of ${\rm \pi}$ radians throughout the sedimentation process. A trajectory for a single $\kappa = 2$ particle composed of Cu$+$St (see table 1) is shown in figure 2(a). Particles with larger values of $\chi$ rotated much more rapidly due to the larger gravitational torque applied to the geometric centre of the particle. This can be compared with an St–St particle in figure 2(b), which shows no preference for rotations since it has no mass polarity ($\chi =0$). For particles with $\chi =0$, we occasionally observed ‘fluttering’, or oscillations of angular orientation during sedimentation. This was likely due to interactions with the walls of the experimental chamber during slight rotations out of the quasi-2-D plane of the experiment (Mitchell & Spagnolie 2014; D'Angelo et al. 2017).

Table 1. The different types of particles used in our experiments along with their corresponding $\kappa$ and $\chi$ values. Materials used are steel (St), aluminium (Al), copper (Cu), Delrin plastic (Pl), tungsten carbide (Tc) and zirconium dioxide (ZrO$_{2}$). Values of $\chi$ are kept to two significant digits.

Figure 2. Two representative examples of the particle trajectories in our single particle experiments. Here, $x=0$ is defined as the geometric centre of the particle at the earliest time. (a) Particle with $\chi > 0$ (Cu$+$St, see table 1). The left part of the panel is a composite image of the particle during the length of the experiment. The right graph shows the corresponding particle orientations, with the arrows pointing from the heavier sphere (Cu) to the lighter sphere (St). The colour bar represents time. Gravity points downward in all pictures. (b) Particle with $\chi = 0$ (St$+$St).

To quantitatively capture the coupling between the external force and dynamics of single particles, we applied the mobility matrix formalism (1.6). Because we are using a quasi-2-D geometry, the complexity of the problem is reduced since the particle can only rotate in the plane. However, the mobility coefficients will be different from those measured in an unbounded, 3-D fluid. With two planar walls, our experimental setup is most similar to a Hele-Shaw cell, where the mobility matrix formalism has already been successfully implemented (Bet et al. 2018) and tested (Georgiev et al. 2020). Because we are considering symmetric prolate particles, the mobility matrix in the body frame (indicated by superscript $b$) is reduced to

(3.1)\begin{equation} \begin{pmatrix} v_{x}^b\\ v_{y}^b\\ \omega_{z}^b \end{pmatrix} = \frac{1}{6{\rm \pi}\eta R} \begin{pmatrix} a_t & 0 & 0\\ 0 & b_t & 0\\ 0 & 0 & \dfrac{3a_r}{4R^2} \end{pmatrix} \begin{pmatrix} F_{x}^b\\ F_{y}^b\\ \tau_{z}^b \end{pmatrix},\end{equation}

where $v_{x}^b$ and $v_{y}^b$ are the translational velocities in the body frame and $\omega _{z}^b$ is the angular velocity perpendicular to the plane of motion. Here, $F_{x}^b$ and $F_y^b$ are the components of the gravitational force in the body frame, and $\tau _{z}^b$ is the external torque from gravity about the particle's centre of geometry (see figure 1). The dimensionless translational mobility coefficients $b_t$ and $a_t$ represent mobility along the major and minor axes of the particle ($b_t>a_t$). The dimensionless rotational mobility coefficient is $a_r$. These coefficients should be identical for all of our particles with the same $\kappa$ and $R$, regardless of the internal density distribution ($\chi$). They characterize the drag from the external flow, which applies stress on the surface of the particle.

Our experimental data, however, are collected in the lab frame. Thus, we first rotate all vectors and the mobility matrix by an angle $\theta$ (figure 1) to obtain the equations of motion in the lab frame:

(3.2)$$\begin{gather} \boldsymbol{\varOmega} = \begin{pmatrix} \cos(\theta) & \sin(\theta) & 0\\ -\sin(\theta) & \cos(\theta) & 0\\ 0 & 0 & 1 \end{pmatrix}, \end{gather}$$

(3.3)$$\begin{gather}\boldsymbol{\varOmega}\boldsymbol{\cdot} \begin{pmatrix} v_{x}^b\\ v_{y}^b\\ \omega_{z}^b \end{pmatrix} = \left(\boldsymbol{\varOmega}\boldsymbol{\cdot} \frac{1}{6{\rm \pi}\eta R} \begin{pmatrix} a_t & 0 & 0\\ 0 & b_t & 0\\ 0 & 0 & \dfrac{3a_r}{4R^2} \end{pmatrix} \boldsymbol{\cdot}\boldsymbol{\varOmega}^{{-}1}\right) \boldsymbol{\varOmega}\boldsymbol{\cdot} \begin{pmatrix} F_{x}^b\\ F_{y}^b\\ \tau_{z}^b \end{pmatrix}. \end{gather}$$

After multiplying and collecting terms, we use the substitutions $2c_1 = a_t + b_t$, $2c_2 = b_t - a_t$ and $c_3 = 3a_r/4$ to write the result in the following form:

(3.4)\begin{equation} \begin{pmatrix} v_{x}\\ v_{y}\\ \omega_{z} \end{pmatrix} =\frac{1}{6{\rm \pi}\eta R} \begin{pmatrix} c_{1} - c_{2}\cos(2\theta) & c_{2}\sin(2\theta) & 0\\ c_{2}\sin(2\theta) & c_{1} + c_{2}\cos(2\theta) & 0\\ 0 & 0 & \dfrac{c_{3}}{R^2} \end{pmatrix} \begin{pmatrix} 0\\ F_{y}\\ \tau_{z} \end{pmatrix}.\end{equation}

We have chosen this parametrization out of convenience. For example, in the case of a perfect sphere, $b_t=a_t$, thus $c_1=1$, $c_2=0$ and $c_3=3/4$ (1.7). We have dropped the superscript since we are referring to the lab frame where the gravitational force only points in the $y$-direction.

The matrix multiplication above gives us the following equations of motion for our particles in the lab frame:

(3.5)$$\begin{gather} v_{x}=\dot{x} = \frac{c_{2}\sin(2\theta)}{6{\rm \pi}\eta R}F_{y}, \end{gather}$$

(3.6)$$\begin{gather}v_{y}=\dot{y} = \frac{c_{1} + c_{2}\cos(2\theta)}{6{\rm \pi}\eta R}F_{y}, \end{gather}$$

(3.7)$$\begin{gather}\omega_{z}=\dot{\theta} = \frac{c_{3}}{6{\rm \pi}\eta R^{3}}\tau_{z}. \end{gather}$$

The dotted variables denote differentiation with respect to time. Similar simplified equations for single-particle dynamics in quasi-2-D geometries have been derived by Bet et al. (2018) and Ekiel-Jeżewska & Wajnryb (2009). In our experiments, the net force and torque on a particle will depend on the values of $\kappa$ and $\chi$. For $\kappa =2$ particles, the net gravitational force and torque about the centre of geometry are

(3.8)$$\begin{gather} F_{y} ={-}\tfrac{4}{3}{\rm \pi} R^{3}(\rho_{h} + \rho_{l} - 2\rho_{f})g, \end{gather}$$

(3.9)$$\begin{gather}\tau_{z} ={-}\tfrac{4}{3}{\rm \pi} R^{4}(\rho_{h} - \rho_{l})g\sin{\theta}. \end{gather}$$

Equations (3.5)–(3.7) are coupled through $\theta$, and can be solved analytically. However, the solution can be generalized by making the equations dimensionless. We used the sphere radius $R$ for a characteristic length scale and $\tau =R/U_T$ for the characteristic time scale, where $U_T$ is the terminal velocity of the lighter sphere (1.1). This non-dimensionalization results in the following equations of motion, where all variables are considered dimensionless for clarity of notation:

(3.10)$$\begin{gather} \dot{x} ={-}K_1c_{2}\sin(2\theta), \end{gather}$$

(3.11)$$\begin{gather}\dot{y} ={-}K_1(c_{1} + c_{2}\cos(2\theta)), \end{gather}$$

(3.12)$$\begin{gather}\dot{\theta} ={-}K_{2}c_{3}\sin(\theta), \end{gather}$$

(3.13)$$\begin{gather}K_1 = \frac{\rho_{h} + \rho_{l} - 2\rho_{f}}{\rho_{l} - \rho_{f}}, \end{gather}$$

(3.14)$$\begin{gather}K_{2} = \frac{\rho_{h} - \rho_{l}}{\rho_{l} - \rho_{f}}. \end{gather}$$

Equation (3.12) can be immediately solved since it is independent of the other equations. The result is

(3.15)\begin{equation} \cot\left(\frac{\theta(t)}{2}\right)= \cot\left(\frac{\theta_{0}}{2}\right) {\rm e}^{K_2c_{3}t},\end{equation}

where $\theta _0$ is the initial value of $\theta$ at $t=0$. Plugging this back into (3.10) and (3.11) and simplifying algebraically, we get

(3.16)$$\begin{gather} x(t) = x_0 + \frac{4c_{2}FK_1\cot(\theta_0/2)}{c_{3}K_{2} - c_{3}F^2 K_{2} \cot(\theta_0/2)^2} -\frac{2c_{2}K_1\sin(\theta_0)}{c_{3}K_{2}}, \end{gather}$$

(3.17)$$\begin{gather}y(t) = y_0 - \frac{K_1}{K_2 c_3}\left((c_1 + c_2)c_3K_2t + 2c_2\left(\cos(\theta_0)+\frac{1-F^2\cot^2 \left(\dfrac{\theta_0}{2}\right)}{1+F^2\cot^2 \left(\dfrac{\theta_0}{2}\right)} \right)\right), \end{gather}$$

where $F = {\rm e}^{K_{2}c_{3}t}$ is a function of time, and used here for compactness. In the limit of particles with uniform mass density ($K_2\rightarrow 0$, $\chi \rightarrow 0$), these functional forms simplify to

(3.18)$$\begin{gather} \theta(t)=\theta_0, \end{gather}$$

(3.19)$$\begin{gather}x(t) = x_0 - c_2 K_1 t \sin(2\theta_0), \end{gather}$$

(3.20)$$\begin{gather}y(t) = y_0 - K_1 t (c_1+c_2\cos(2\theta_0)). \end{gather}$$

Equations (3.18) and (3.19) verify the prediction that for polar particles of uniform density, the angle of inclination does not change, and the particle drifts laterally in the $x$-direction (Ramaswamy 2001).

After taking the inverse cotangent of (3.15) and using standard least-squares nonlinear regression, we can fit these analytic forms to the experimental data with very good agreement. Figure 3 shows five identical experiments and their corresponding fits. For $\theta (t)$, there are only two fitting parameters, $c_3$ and $\theta _0$. Once they are determined by the fit, then $x(t)$ can be fit for the parameters $c_2$ and $x_0$. Finally, $y(t)$ can then be fit for $c_1$ and $y_0$. The curves are compared to each other by assigning $t=0$ when the particles are completely horizontal, i.e. $\theta ={\rm \pi} /2$. We also moved the $x$ and $y$ origin to correspond to $t=0$. Open symbols represent data and curves are the fits to (3.15)–(3.17). The fits for the $x$-position show more systematic deviation from the data, yet the overall displacement is also much smaller. For example, as shown in the inset in $y$ versus $t$, the residuals of these fits are comparable to the variability in $x$ versus $t$, which is a fraction of a particle radius in displacement. Although the source of the systematic asymmetry is unclear, we suspect that when particles are released from the gating mechanism, they are not perfectly parallel with the walls of the quasi-2-D chamber. If a particle's alignment varies during the rotation from $\theta = {\rm \pi}$ to $\theta = 0$, we would expect variations in the mobility coefficients (i.e. $c_2$) due to wall effects (Brenner 1999; Mitchell & Spagnolie 2014), resulting in an asymmetry in $x(t)$ about $\theta = {\rm \pi}/2$. Additionally, we do not expect errors in particle tracking to lead to systematic asymmetry even though the $x$-motion is of the order of the particle size. Tracking errors would manifest more as random noise rather than systematic deviations from theory. The data for $\theta$, $x$ and $y$ can also be fit simultaneously using a global least squares regression for all parameters, since parameters appear in multiple equations. We found less than 5 % difference in the fitted parameter values using this method, so we have only chosen to report the results of the sequential fitting. Similar analytic solutions and quality of fits were recently found in the alignment of mirror-symmetric particles in a microfluidic device (Bet et al. 2018; Georgiev et al. 2020).

Figure 3. Data for five experiments with a single Al$+$Pl particle (table 1). Only 5 % of points are plotted for clarity. Initially, the heavy aluminium sphere begins above the lighter Delrin sphere. Open symbols represent data, and curves are model fits from (3.15)–(3.17). Different symbols and colours are separate experiments. Inset shows the residual difference between the model fit $y_m$ and the data $y$ for the vertical position of the particle.

One of the major assumptions of our model was that all coefficients are independent of $\chi$, and only depend on the shape of the composite, prolate particles. This is evident from (3.1), since $a_t$, $b_t$ and $a_r$ are dimensionless coefficients that only depend on the particle shape, not the density distribution. We confirmed this prediction using all fits of single particle experiments with $\kappa = 2$, as shown in figure 4(a–c). The coefficients $c_1$, $c_2$ and $c_3$ are computed directly from nonlinear least-squares regression of the data (3.10)–(3.12). For particles with $\chi = 0$ (uniform density), we used (3.18) and (3.19) to fit the data. In this form, there is no torque from gravity, so $c_3$ cannot be determined and is not shown. However, $c_1$ and $c_2$ can be determined, but are not very reliable because of experimental artefacts that affect the angle (and thus translational velocity) during sedimentation. These artefacts include small differences in the distribution of glue used between the particles, rotations out of the quasi-2-D plane and other 2-D confinement effects such as ‘fluttering’ (Brenner 1999; Mitchell & Spagnolie 2014; D'Angelo et al. 2017). For finite $\chi$, the particles rotate significantly due to gravitational torque, and $c_1$, $c_2$ and $c_3$ can be determined reliably. There appears to be some small systematic trend in $c_1$, but the overall variation is small and the data for all parameters is consistent with a constant value over the range $0<\chi <0.25$.

Figure 4. Best-fit parameters versus $\chi$ from (3.15)–(3.17) for single particle sedimentation experiments with $\kappa =2$. The material combinations used were circle, St$+$St; diamond, Cu$+$St; plus, St$+$ZrO$_{2}$; right triangle, Cu$+$ZrO$_{2}$; square, Al$+$Pl; star, Al$+$St (see table 1). Each data point is the weighted mean of five different trials with error bars representing the standard error of the mean. (a–c) Parameters $c_1$–$c_3$ directly computed from the nonlinear regression of data in the lab frame. (d–f) Body frame coefficients: $a_t$, $b_t$ and $a_r$. Here, St$+$St is missing from $c_3$ and $a_r$ because of the limiting form of $\theta (t)$ when $\chi = 0$ (3.18), which has no dependence on $c_3$.

Using $2c_1 = a_t + b_t$, $2c_2 = b_t - a_t$ and $c_3 = 3a_r/4$, we computed the shape-dependent drag coefficients of our symmetric particles in the body frame, as shown in figure 4(d,e). Here again, $a_r$ cannot be determined from the $\chi = 0$ data, and data with finite $\chi$ are most reliable. The average values of these mobility coefficients with $\chi >0$ are shown by the dashed lines: $\bar {a}_t=0.328\pm 0.018$, $\bar {b}_t=0.404\pm 0.012$ and $\bar {a}_r=0.221\pm 0.008$. We suggest that these experimental values for the mobility coefficients can be compared directly to simulations of particles composed of spheres (Garcia de la Torre & Carrasco 2002). In comparison to sedimentation in an unbounded, 3-D fluid, we expect our measured values of $a_t$, $b_t$ and $a_r$ to be somewhat smaller since the particles experience a larger drag due to the confining walls.

3.2. Sedimenting particle pairs

Goldfriend et al. (2017) theoretically examined a sedimenting suspension of mass polar spheroids using a continuum linear stability analysis. To briefly summarize their results, they considered a suspension of particles settling due to an external body force $F$ in the direction of gravity in a fluid of viscosity $\eta$. A sinusoidal concentration perturbation was applied in a direction perpendicular to the force with amplitude $c(x)$ and a characteristic wavelength $\lambda$. These fluctuations in the concentration create velocity fluctuations, $U(x)$. Balancing the change in gravitational force of the suspension versus the change in the drag force gives $c\lambda F \approx U \eta /\lambda$. By solving for the amplitude $U$, we see that $U \sim c\lambda ^2 F/\eta$. The indefinite scaling of $U$ with $\lambda$ is a demonstration of the Caflisch–Luke paradox described in the introduction. These slabs of particles will also experience vorticity of the magnitude $U/\lambda \sim c\lambda F/\eta$. For uniform spheres, this will cause a rotation of the sphere, but no drift. However, self-aligning objects will be tilted away from their preferred alignment. This causes a drift in the $x$-direction with velocity $\sim \gamma R c F/\eta$, where $\gamma =\gamma (\kappa,\chi )$ is a proportionality constant determined by the shape and mass distribution of an individual particle. For positive $\gamma$, which requires $\kappa >1$ (Goldfriend et al. 2017), the relative velocity of the particles is positive, meaning they drift away from each other. This is the screening mechanism that stabilizes the suspension. For negative $\gamma$, which requires $\kappa <1$, they drift towards each other, leading to unbounded growth of the instability.

In our experiments, we examined the particle-level interactions by measuring the relative separation of pairs of prolate ($\kappa >1$) particles as they repel each other during sedimentation. We placed two particles heavy-side down in adjacent slots of the plastic gate so that their initial separation was 3.3 mm. Each experiment was conducted five times for reproducibility. Figure 5 shows a representative selection of settling trajectories for various values of $\kappa$ and $\chi$. These are composite images of the particles during sedimentation, spaced 3.33 s apart. The arrows to the right of each panel show the orientations of each particle during sedimentation and the colour represents time. First, particles with $\chi = 0$ heavily influenced each other. Their dynamics were typically characterized by one of the particles rotating or flipping completely. This particle often lagged behind the other one, which did not flip, but followed a curved trajectory. This can be seen in both figures 5(a) and 5(d). The particles did not preferentially align to gravity, and instead produced a variety of dynamics. For example, the periodic variation in separation visible in figure 5(d) is reminiscent of Kepler orbits observed in sedimenting pairs of disks (Chajwa et al. 2019). In fact, a periodic variation in the relative position between adjacent, sedimenting prolate particles was theoretically predicted by Claeys & Brady (1993) (see their figure 4). On average, we did not observe a net repulsion or attraction between our particles with a uniform mass distribution ($\chi =0$).

Figure 5. Experimental trajectories of two-particle interaction experiments. In each alphabetic panel, the image shows a composite of the particles’ trajectories during sedimentation. The graph shows the orientation of the particles with arrows pointing from the heavier particle to the lighter particle(s). The colour is used to show when two arrows are at the same time during their transit. Here, $x=0$ is defined as the halfway point between the particle centres on the first frame. Panels (a–c) are $\kappa =2$ particles, panels (d–f) are $\kappa =3$ particles. Panels (a) and (d) show particles with $\chi =0$. Panels (b) and (e) show particles with the smallest $\chi$, Cu$+$St (see table 1). Panels (c) and (f) show particles with the largest $\chi$, Cu$+$Pl (see table 1).

For particles with $\chi > 0$, there was an immediate rotation and repulsion between the particles leading to a horizontal separation that grew with time. Eventually the particles would align with the external gravitational field, and the separation saturated. This is shown in figure 5(b,c) for $\kappa = 2$ and figure 5(e,f) for $\kappa = 3$. The finite width of the quasi-2-D chamber, 4$R$, introduced a length scale that could potentially set an upper limit on the range of the repulsive interaction. However, we observed that the final separation between the particles could be as much as 30$R$ (figure 5e) for smaller values of $\chi$. The repulsive effect was most prominent for particles composed of materials with closely matched densities (i.e. copper and steel). Although this may seem counter-intuitive at first, particles with $0 < \chi \ll 1$ can rotate away from vertical more easily, and thus experience a larger repulsion and horizontal drift. As $\chi \rightarrow 0$, we expect one of the particles to be able to flip entirely if they are close enough to interact strongly, leading to the periodic type of interactions observed for $\chi =0$ (figures 5a and 5d). In this limit, the eventual behaviour of the sedimenting particles should be determined both by $\chi$ and by the initial separation.

The inverse relationship between $\chi$ and the mutual repulsion was also predicted by Goldfriend et al. (2017). The authors found that the growth rate of the horizontal velocity fluctuations scaled as $\gamma = \kappa ^{2/3}/3\chi$ for highly prolate particles ($\kappa \gg 1$). To quantify this effect in our experiments, we chose to measure the total change in horizontal separation, ${\rm \Delta} H$, between the particles’ geometric centres in each experiment. This is plotted in figure 6(a) as a function of $\kappa$. Generally, the separation increased with $\kappa$. However, to compare between each set of experiments that corresponded to different values of $\chi$, we multiplied the final separation by $\chi ^\alpha$, where $\alpha$ was determined by simultaneous fitting of all the data to the following form:

(3.21)\begin{equation} \chi^\alpha\frac{{\rm \Delta} H}{2R}=A(\kappa-1),\end{equation}

where we have imposed the requirement that there be no repulsion for $\kappa =1$ (i.e. single spheres). The fit was performed by subtracting the left- and right-hand sides of (3.21), squaring the difference and summing over all data points. The best fit values for the parameters were $\alpha =0.39\pm 0.05$ and $A=1.28\pm 0.20$, where the errors represent one standard error. The fit shows very good agreement with the data, as plotted in figure 6(b).

Figure 6. Graphs of our experimental response parameter, ${\rm \Delta} H$, versus $\kappa$. Here, ${\rm \Delta} H$ is defined as the difference between the final and initial horizontal separation. The colours and shapes represent different material combinations of the composite particles. By order of increasing $\chi$, they are Cu$+$Pl (upright triangle), Al$+$St (upside down triangle), Tc$+$St (star) and Cu$+$St (circle) (see table 1). Panel (a) is the raw data, with each data point being an average over five runs and error bars representing one standard deviation. Panel (b) is the same data collapsed using the best fit parameters obtained from (3.21).

In general, the predictions from Goldfriend et al. (2017) are in excellent qualitative agreement with our experiments, yet the scaling, $\alpha \sim 0.39$, is quite different than that predicted by Goldfriend et al. (2017): $\alpha \sim 1$. There are a few reasons that can explain this discrepancy. (1) $\gamma$ represents an instantaneous response for an initially uniform concentration of particles. Here we are using the final separation, ${\rm \Delta} H$, which is essentially an integral of the repulsion between the particles in time. (2) The quasi-2-D environment should screen the long-ranged, $1/r$ hydrodynamic interactions (Beatus et al. 2012), so one may expect a different theoretical scaling between $\gamma$ and $\chi$ based purely on geometry. (3) Our quasi-2-D chamber may introduce other effects that depend on the thickness of the chamber, for example, it is well known that the net viscous drag force on a sedimenting particle can be dependent on the distance to a nearby wall (Brenner 1999; Mitchell & Spagnolie 2014). (4) Our particles are not perfect examples of the prolate and oblate ellipsoids discussed by Goldfriend et al. (2017). Despite these differences, the experimental data with different values of $\kappa$ and $\chi$ can be reasonably collapsed using the dependence listed in (3.21).

Lastly, we verified that this mutual repulsion led to more uniformly distributed suspensions of many particles. We filled our quasi-2-D chamber with 29 particles with $\kappa = 2$. The left column of figure 7 shows that for particles with $\chi =0$, there is no preferential alignment to gravity, resulting in a large spread of particle separations, both vertically and horizontally. Particles can flip very easily and often come into contact. Some of the particles experienced small rotations out of the plane as well. The right column of figure 7 illustrates that particles with $\chi =0.24$ followed a more uniform spatial distribution. All particles tended to align with gravity, resulting in a mutual repulsion. When particles are in close proximity, they tilted away from the vertical and drifted apart, similar to figure 5. Surprisingly, the particles with $\chi = 0.24$ did not spread as much in the vertical direction as $\chi = 0$, suggesting that vertical fluctuations in concentration may be suppressed for $\chi >0$. An intuitive explanation for this behaviour stems from the variations in vertical velocities of particles. For $\chi = 0$, particle rotations lead to a spread in vertical terminal velocities (3.17), whereas particles with $\chi >0$ are mostly aligned to gravity and sediment at the same rate. Although vertical fluctuations were not directly addressed by Goldfriend et al. (2017), we hypothesize that the mutual repulsion in mass polar particles also suppresses the ‘clumping instability’ observed in uniform suspensions (Chajwa et al. 2020) that leads to large vertical separations between particles.

Figure 7. Series of three sequential images of 29 particles sedimenting in our quasi-2-D chamber. (a) ($\chi = 0$) shows St$+$St particles and (b) ($\chi = 0.24$) shows St$+$Al particles (table 1). The top row shows images at the same vertical position near the top of the experimental chamber at early times, the middle row shows the same particles later in time and the bottom row shows the particles near the end of the experiment, at the bottom of the chamber.

4. 3-D particle suspensions

Although the effective repulsion between our prolate particles is apparent in a confined, quasi-2-D environment, it is possible that the dynamical evolution of these particles in three dimensions could hinder the repulsion since the particles have more motional degrees of freedom. Figure 8 illustrates the experimental procedure, described in § 2, where particles are sedimented repeatedly and their resulting spatial distribution is imaged after each repeated experiment. For $\chi =0$, particles tended to cluster during sedimentation, resulting in a rapid increase in their velocity due to mutual drag reduction at finite distances. For $\chi =0.18$, there is a visible alignment of particles to the direction of gravity (vertical), and a broader spatial distribution with less clustering.

Figure 8. (a) Experimental procedure for 3-D sedimentation. The sealed chamber was repeatedly flipped and imaged, as described in § 2. (b) Sample image from the side during a single sedimentation experiment of particles with $\chi =0$ and $\kappa =3$ (St$+$St$+$St). (c) Sample image during sedimentation of particles with $\chi =0.18$ and $\kappa =3$ (St$+$Al$+$Al).

To quantify the post-sedimentation spatial distribution of particles, one would ideally extract the centre of mass position of each particle and calculate the radial distribution function of their positions. However, after sedimentation, we found that particles often overlapped by stacking in the vertical direction, making identifying the centre of mass impossible. Instead, we choose to threshold the images so that particles became black pixels, and the background became white. Samples of these images are shown in figure 9. We then calculated the radial distribution function of the positions of the black pixels. This was done by first finding the centre of mass of all black pixels, corresponding to $r=0$ in each image, and then binning pixel positions radially along $r$. We divided the number of pixels in each bin range by the area of the annulus associated with the bin. To compare between $\kappa =2$ and $\kappa =3$, we also normalized the radial positions by the root mean squared radius of the data for $\chi =0$, calculated by

(4.1)\begin{equation} r_{rms} = \sqrt{\frac{1}{N} \sum_r r^2}.\end{equation}

Here, the sum runs over every black pixel in all 50 images associated with $\chi =0$, and $N$ is the sum total of all black pixels from all images with $\chi =0$.

Figure 9. (a) Sample images of particles resting on the bottom of the 3-D chamber after sedimentation. The left column shows $\kappa =2$ particles with $\chi =0$ (St$+$St; top) and $\chi =0.24$ (St$+$Al, bottom). The right column shows $\kappa =3$ particles with $\chi =0$ (St$+$St$+$St; top) and $\chi =0.18$ (St$+$Al$+$Al, bottom). We quantified the radial distribution of black pixels at a different radii $r$ from the centre of the sedimented pattern, as indicated by the magenta arrow. (b) Radial probability density function (PDF) of black pixels from the images. The radius has been normalized by $r_{rms}$ for $\chi = 0$ (4.1). The legend indicates the values of $\kappa$ and $\chi$ for the different PDFs. Each curve was produced using data from 50 post-sedimented images.

The resulting distribution functions are shown in figure 9. After normalizing the radial position, we see that all data for $\chi =0$ collapses onto the same distribution (dashed lines). As expected, when $\chi >0$ (solid lines), these distributions broaden due to the net repulsion between the particles. Furthermore, particles with larger $\kappa$ and smaller non-zero $\chi$ should experience a larger repulsion, as predicted by (3.21). This is consistent with our data, since the distribution for $\kappa =3$, $\chi =0.18$, is broader than for $\kappa =2$, $\chi =0.24$. We note that because the initial state of each round of sedimentation was set by the final state of the previous one, the sequential images of the final sedimentation pattern were not statistically independent. Nevertheless, we do not expect these effects to qualitatively change our results, and taken together, figures 8 and 9 confirm that the effective pairwise repulsion between mass polar particles also suppresses clumping in three dimensions.

5. Conclusion

Particles with mass polarity are forced to align with the direction of gravity during sedimentation. This alignment arises because the centre of mass is displaced from the geometric centre of each particle, resulting in a net torque imposed by the fluid flow. Our work examined the motion of single particles sedimenting in a viscous fluid, and we derived a simple analytic expression for the their motion in a quasi-2-D environment. Fitting trajectories of the individual particles allowed us to reconstruct the parameters of the mobility matrix. When two or more prolate particles are present, we showed that they experience a mutual repulsion, as first described by Goldfriend et al. (2017). Surprisingly, this repulsion is strongest for small values of $\chi$, i.e. when the centre of mass is only slightly displaced from the centre of geometry. The repulsion also increases as the particles become more asymmetric (more prolate, large $\kappa$). We also showed that this overall repulsion persists in 3-D experiments with hundreds of particles.

There are still many open questions facilitated by this work. First, Goldfriend et al. (2017) showed that there should be a mutual attraction between particles for $\kappa <1$. We found that our particle fabrication method, i.e. gluing individual spheres together, did not easily lend itself to making oblate particles with $\kappa <1$. Such particles would cluster rapidly during sedimentation, and may dramatically increase the overall sedimentation rate of a suspension of particles. Additionally, Goldfriend et al. (2017) predicted the existence of hyperuniformity in the density distribution of a sedimenting suspension. Our experimental results in three dimensions demonstrate a net repulsion and a more uniform concentration, yet we would need many more particles with accurate tracking in three dimensions to quantify hyperuniformity. One alternative route could be simulating many particles efficiently with a parametrized interaction based on our results. We hope our simplified mobility matrix may serve as a starting point for such idealized simulations of many interacting particles.