2. Experiments

The sedimentation experiments were performed in a quasi-two-dimensional slab geometry (see figure 2) with height $45$ cm, width $90$ cm and thickness $5$ cm filled with silicone oil with density $0.98$ g cm⁻ $^{3}$ and kinematic viscosity $5000$ cSt. Discs of diameter $2a = 8$ mm and thickness $2b=1$ mm were 3D-printed with resin of density $1.164$ g cm⁻ $^{3}$ . These discs are modelled in our theory as oblate spheroids of aspect ratio $b/a=0.125$ . The Reynolds number using the length scale set by the particle size was measured to be $\sim 10^{-4}$ , ensuring that we are in a Stokesian regime. The quasi-two-dimensional geometry of the container has two effects: it leads to a modest reduction of the single-particle sedimentation speed, and it cuts off hydrodynamic interactions with neighbours beyond a length scale set by the thickness.

Releasing multiple particles in a viscous fluid presents a challenge in synchronising release times and controlling initial orientations, as pointed out by Jung et al. (Reference Jung, Spagnolie, Parikh, Shelley and Tornberg2006) and Chajwa et al. (Reference Chajwa, Menon, Ramaswamy and Govindarajan2020), and described in Chajwa (Reference Chajwa2021). Our release mechanism consists of an array of slots where discs can be selectively placed to define the initial (dimensional) lattice spacing $d$ and dimensional perturbation wavenumber $2 \pi /\lambda$ (see figure 2). The release mechanism is centred along the thickness of the slab and immersed in the fluid to remove bubbles attached to the discs. The discs are ejected from the slots simultaneously using a ‘comb’ that inserts into the slots. This set-up allows us to study the dependence of the particle dynamics on $d$ and $\lambda$ . A D-SLR camera is used to capture images of the discs every 3 s. The images are thresholded, and an ellipse is fit to each disc to determine its orientation and the location of its centroid. In the later stages, the thresholded images had to be manually segmented to separate closely clustered discs.

Figure 2.

Experimental set-up: (a) shows the quasi-two-dimensional geometry of the container with gravity axis pointing along $-\hat {z}$ . The mechanism simultaneously releases discs with controlled lattice spacing and perturbation wavenumber. (b) The bottom view of the release mechanism shows the array of discs moments before they were released into the fluid. The initial spacing is chosen to control the dimensional lattice spacing $d$ and the dimensional perturbation wavenumber $2\pi /\lambda$ .

3. Governing equations

Our system consists of hydrodynamically interacting oblate spheroids (discs) labelled by $\alpha \in \{1, \ldots, \, N_d\}$ sedimenting under gravity in an otherwise quiescent unbounded fluid. The $\alpha$ th disc is described by semi-major and semi-minor axes $a$ and $b$ , respectively, a centre position $\boldsymbol x_\alpha$ , orientation unit vector $\boldsymbol p_\alpha$ which points along the symmetry axis of the disc and instantaneous velocity $\boldsymbol{V}_\alpha$ . For particle Reynolds number ${\rm Re}_p\ll 1$ , every disc generates a disturbance field that satisfies the Stokes equation (Kim Reference Kim1985; Kim & Karrila Reference Kim and Karrila2013). Moreover, when the particle relaxation time scale $\tau _p$ is much smaller than the time it takes to fall through its length, the Stokes number $St \equiv (\tau _p V_\alpha /a) \ll 1$ for each disc, so that its buoyancy-corrected weight $\boldsymbol{F}^g$ is instantaneously balanced by hydrodynamic drag $\boldsymbol{F}^h_\alpha$ for each disc. The simplest approximation for the hydrodynamic interaction experienced by a given disc would be to sum the flows due to the other discs, treating each disc as a force monopole placed at its centroid. The resulting velocity, vorticity and rate of strain at the location of the disc in question translate, rotate and align it, leading to the equations of motion

(3.1) \begin{align}& \boldsymbol V_\alpha ^{{pt}} = \frac {\boldsymbol{F}^g}{6\pi \mu a}\boldsymbol{\cdot }\left [\frac {\boldsymbol{p}_\alpha \boldsymbol{p}_\alpha }{X^A} + \frac {(\boldsymbol{\delta } - \boldsymbol{p}_\alpha \boldsymbol{p}_\alpha )}{Y^A}\right ] + \sum _{\beta \neq \alpha }^{N_d} \frac {\boldsymbol F^g}{8\pi \mu }\boldsymbol{\cdot }\unicode{x1D642}(\boldsymbol x_{\alpha \beta }), \end{align}

(3.2) \begin{align} \frac {d}{{\rm d}t}\boldsymbol p_\alpha &= \sum _{\beta \neq \alpha }^{N_d} \left ( \frac {\boldsymbol F^g \times \boldsymbol x_{\alpha \beta }}{8\pi \mu |\boldsymbol x_{\alpha \beta }|^3} \right ) \times \boldsymbol p_\alpha \nonumber \\ &\quad + \frac {3e^2}{2-e^2} \sum _{\beta \neq \alpha }^{N_d}\left (\frac {\boldsymbol F^g \boldsymbol{\cdot } \boldsymbol x_{\alpha \beta }}{8\pi \mu } \right ) \left (\frac {\boldsymbol p_\alpha \boldsymbol{\cdot }\boldsymbol x_{\alpha \beta }}{|\boldsymbol x_{\alpha \beta }|^5} \right ) [\boldsymbol x_{\alpha \beta } - \boldsymbol p_\alpha (\boldsymbol p_\alpha \boldsymbol{\cdot } \boldsymbol x_{\alpha \beta })], \end{align}

where

(3.3) \begin{equation} \unicode{x1D642}(\boldsymbol x) = \frac {\boldsymbol \delta }{|\boldsymbol x|} + \frac {\boldsymbol x \boldsymbol x}{|\boldsymbol x|^3} \end{equation}

is the Green’s function of the steady Stokes equation. In (3.1) and (3.2), $\boldsymbol x_{\alpha \beta } \equiv \boldsymbol x_\alpha - \boldsymbol x_\beta$ , $\mu$ is the dynamic viscosity of the fluid, $e = \sqrt {1-b^2/a^2}$ is the eccentricity and

(3.4) \begin{align}X^A &\equiv \frac {4}{3}e^3\big[\big(2e^2-1\big)K+e\sqrt {1-e^2}\big]^{-1},\end{align}

(3.5) \begin{align}Y^A &\equiv \frac {8}{3}e^3\big[\big(2e^2+1\big)K-e\sqrt {1-e^2}\big]^{-1},\end{align}

with

(3.6) \begin{equation}K \equiv \cot ^{-1} \frac {\sqrt {1-e^2}}{e}, \\[-6pt]\nonumber\end{equation}

are resistance functions for an oblate spheroid (Kim & Karrila Reference Kim and Karrila2013), which approximates the experimental discs with an aspect ratio of $b/a = 0.125$ , corresponding to $e = 0.992$ . The first term in (3.1) is the sedimentation velocity of an isolated disc falling under gravity, and the second is the sum of disturbance fields produced by other discs, in the point-force approximation. The two terms in (3.2) respectively are the contribution from the vorticity and the strain rate of the disturbance field to the orientational dynamics of the disc.

The point-force approximation captures the essential features of the periodic and unbounded dynamics of a pair of discs, and sharp transitions between these states, as the initial condition is varied (Chajwa et al. Reference Chajwa, Menon and Ramaswamy2019). It also provides a basic framework to understand the linear instability and wave-like modes observed in a one-dimensional array of sedimenting discs (Chajwa et al. Reference Chajwa, Menon, Ramaswamy and Govindarajan2020). However, at late times, the initial lattice state is disrupted, giving rise to the formation of configurations with small inter-particle separations, for which the point-force description is expected to be inadequate, since this approximation is accurate only up to $O(a/r)$ where $r=|\boldsymbol x_{\alpha \beta }|$ . Therefore, for nearly touching disc configurations observed during late stages of instability in our experiments, we need to account for finite size effects and to satisfy the no-slip boundary conditions on each particle. We use the method of reflections to study the close-range dynamics. This perturbative method is needed to iteratively satisfy the boundary conditions on the surface of each disc, with each iteration progressively correcting the errors from the previous one (Kim Reference Kim1985; Kim & Karrila Reference Kim and Karrila2013). At the $n^{}$ th reflection, the velocity at the location of particle $\alpha$ satisfies the no-slip condition on it due to the velocities induced by all other particles calculated at the $(n-1)^{}$ th reflection, and this corrected velocity of particle $\alpha$ will induce a small slip at the surfaces of all the other particles. As $n$ increases, the error goes to higher and higher orders in the interparticle spacing. We find that the first reflection, combined with a simple lubrication correction and near-contact repulsion, is sufficient to predict the structural arrangement of discs observed in the experiments. The first reflection is accurate up to $O(a^3/r^3)$ in linear velocity and $O(a^4/r^4)$ in the angular velocity. The linear velocity $\boldsymbol{V}_\alpha$ and angular velocity $\boldsymbol \varOmega _\alpha$ of a disc $\alpha$ can be expressed in terms of the hydrodynamic force $\boldsymbol{F}^h_\alpha$ and hydrodynamic torque $\boldsymbol{T}^h_\alpha$ acting on it, using the Faxén laws as (Kim & Karrila Reference Kim and Karrila2013)

(3.7) \begin{equation} \boldsymbol{V}_\alpha = \frac {-\boldsymbol{F}^h_\alpha }{6\pi \mu a}\boldsymbol{\cdot }\bigg[\frac {\boldsymbol{p}_\alpha \boldsymbol{p}_\alpha }{X^A} + \frac {\left(\boldsymbol{\delta } - \boldsymbol{p}_\alpha \boldsymbol{p}_\alpha \right)}{Y^A}\bigg] + \frac {1}{2c}\int _{-c}^{c} {\rm d}\xi \bigg\{1+ \frac {\left(c^2-\xi ^2\right)}{4e^2}\nabla^2 \bigg \} \boldsymbol{u}^\infty _\alpha \bigg |_{\boldsymbol x_\alpha (\xi )}, \end{equation}

and

(3.8) \begin{align} \boldsymbol{\varOmega }_\alpha &= \frac {-\boldsymbol{T}^h_\alpha }{8\pi \mu a^3}\boldsymbol{\cdot }\bigg[\frac {\boldsymbol{p}_\alpha \boldsymbol{p}_\alpha }{X^C} + \frac {\left(\boldsymbol{\delta } - \boldsymbol{p}_\alpha \boldsymbol{p}_\alpha \right)}{Y^C}\bigg] + \frac {3}{8c^3}\int _{-c}^{c} {\rm d}\xi \big(c^2-\xi ^2\big) \boldsymbol{\nabla } \times \boldsymbol{u}^\infty _\alpha \bigg |_{\boldsymbol x_\alpha (\xi )} \nonumber \\ &\quad - \frac {3}{4c^3} \frac {e^2}{2-e^2} \int _{-c}^{c} {\rm d}\xi \left \{ \big(c^2-\xi ^2\big) \left [ 1 + \big(c^2-\xi ^2\big)\frac {1}{8e^2}\nabla^2 \right ] \right \} \boldsymbol{p}_\alpha \times \bigg [ \boldsymbol{E}^\infty _\alpha \bigg |_{\boldsymbol x_\alpha (\xi )}\boldsymbol{\cdot } \boldsymbol{p}_\alpha \bigg ]. \end{align}

Here, $c \equiv a e$ , $\boldsymbol{E}^\infty _\alpha$ is the strain rate associated with $u^\infty _\alpha$ and

(3.9) \begin{align}&\qquad\qquad X^C \equiv \frac {2}{3}e^3\big[K - e\sqrt {1-e^2}\big]^{-1}, \end{align}

(3.10) \begin{align}& Y^C \equiv \frac {2}{3}e^3\big(2-e^2\big)\big[e\sqrt {1-e^2}-\big(1-2e^2\big)K\big]^{-1}. \end{align}

In writing (3.7) and (3.8), we have used the fact that the singularity distribution of an oblate spheroid (disc) can be written in terms of a line distribution placed along an imaginary focal length (Shatz Reference Shatz2004). This requires us to evaluate the background incidence field on disc $\alpha$ , $\boldsymbol{u}^\infty _\alpha$ , at $\boldsymbol{x}_\alpha (\xi ) \equiv \boldsymbol{x}_\alpha + i \xi \boldsymbol{p}_\alpha$ . Up to the first-reflection correction, the background incidence field on disc $\alpha$ , $\boldsymbol{u}_\alpha ^\infty$ , is given by the superposition of the disturbance fields generated by all the other discs, treating them as isolated and satisfying the no-slip boundary condition on their respective surfaces. Thus, we have

(3.11) \begin{align} \boldsymbol{u}_\alpha ^\infty ({\boldsymbol{x}}) &\approx \sum _{\beta \neq \alpha }^{N_d} \left[\frac {-\boldsymbol{F}^h_\beta }{8\pi \mu } \boldsymbol{\cdot } \frac {1}{2c}\int _{-c}^{c} \left \{ 1 + \frac {\big(c^2-\xi ^2\big)}{4e^2}\nabla^2 \right \} \unicode{x1D642}\left({\boldsymbol{x}} - {\boldsymbol{x}}_\beta - i\xi \boldsymbol{p}_\beta \right){\rm d}\xi \right .\nonumber \\ &\quad +\frac {1}{2}\frac {\boldsymbol{T}^h_\beta \times \boldsymbol{\nabla }}{8\pi \mu }\boldsymbol{\cdot } \frac {3}{4c^3}\int _{-c}^c \big(c^2-\xi ^2\big) \unicode{x1D642}\left({\boldsymbol{x}} - {\boldsymbol{x}}_\beta - i\xi \boldsymbol{p}_\beta \right){\rm d}\xi \nonumber \\ & \quad - \boldsymbol{\unicode{x1D64E}}_\beta ^h\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\cdot }\left . \frac {3}{4c^3}\int _{-c}^c\big(c^2-\xi ^2\big)\left \{1+\frac {c^2-\xi ^2}{8e^2}\nabla^2 \right \} \unicode{x1D642}\left({\boldsymbol{x}} - {\boldsymbol{x}}_\beta - i\xi \boldsymbol{p}_\beta \right){\rm d}\xi \vphantom{\frac {-\boldsymbol{F}^h_\beta }{8\pi \mu }} \right], \end{align}

where $N_d$ is the total number of discs and

(3.12) \begin{equation} \boldsymbol{\unicode{x1D64E}}_\beta ^h \equiv \frac {1}{2}\left (\frac {e^2}{2-e^2}\right )\Big[\Big(\boldsymbol{T}_\beta ^h\times \boldsymbol{p}_\beta \Big)\boldsymbol{p}_\beta + \boldsymbol{p}_\beta \Big(\boldsymbol{T}_\beta ^h\times \boldsymbol{p}_\beta \Big)\Big]. \end{equation}

Substituting (3.11) into (3.7) and (3.8), we get the far-field mobility matrix $\boldsymbol{\mathcal{M}}^\infty$ which relates the generalised velocities $\boldsymbol{\mathcal{V}}$ to generalised forces $\boldsymbol{\mathcal{F}}^h$ as

(3.13) \begin{align} \boldsymbol{\mathcal{V}} \approx -\boldsymbol{\mathcal{M}}^\infty \boldsymbol{\cdot } \boldsymbol{\mathcal{F}}^h; \, \, \boldsymbol{\mathcal{V}} \equiv \left[\boldsymbol{V}_1 \ldots \boldsymbol{V}_{N_d} \boldsymbol{\varOmega }_1 \ldots \boldsymbol{\varOmega }_{N_d}\right]^T \!, \,\,\,\, \boldsymbol{\mathcal{F}}^h \equiv \left[\!\boldsymbol{F}^h_1 \ldots \boldsymbol{F}^h_{N_d} \boldsymbol{T}^h_1 \ldots \boldsymbol{T}^h_{N_d}\right]^T .\\[-18pt]\nonumber\end{align}

As the discs approach each other, the far-field mobility matrix has to be supplemented with lubrication effects, which are dominant in the near contact configurations of the discs. The lubrication effects are preserved in the resistance formulation, which can be incorporated by inverting $\boldsymbol{\mathcal{M}}^\infty$ and adding the near-contact resistance matrix $\boldsymbol{\mathcal{R}}^L$ (Durlofsky, Brady & Bossis Reference Durlofsky, Brady and Bossis1987; Claeys & Brady Reference Claeys and Brady1993). We use a simple model to compute the localised lubrication interaction between the nearly touching discs (Berne & Pechukas Reference Berne and Pechukas1972; Ladd & Verberg Reference Ladd and Verberg2001; Günther et al. Reference Günther, Frijters and Harting2014). The hydrodynamic lubrication force on disc $\alpha$ , nearly touching another disc $\beta$ , is given by

(3.14) \begin{equation} \boldsymbol{F}^L_{\alpha \beta } = -\frac {6\pi \mu a \big[\big(\boldsymbol{V}_\alpha - \boldsymbol{V}_\beta \big)\boldsymbol{\cdot } \hat {\boldsymbol{\epsilon }}_{\alpha \beta }\big]}{8\sqrt {1-\big(\varGamma \boldsymbol{p}_\alpha \boldsymbol{\cdot }\boldsymbol{p}_\beta \big)^2}}\left (\frac {a}{|\boldsymbol{\epsilon }_{\alpha \beta }|} - \frac {1}{\varDelta _c} \right )\hat {\boldsymbol{\epsilon }}_{\alpha \beta }, \quad \varGamma \equiv \frac {a^2-b^2}{a^2+b^2}, \end{equation}

where $\varDelta _c = 2/3$ is the critical value below which the lubrication forces (3.14) are computed (Ladd & Verberg Reference Ladd and Verberg2001; Günther et al. Reference Günther, Frijters and Harting2014). Here, $\hat {\boldsymbol{\epsilon }}_{\alpha \beta } = \boldsymbol{\epsilon }_{\alpha \beta }/|\boldsymbol{\epsilon }_{\alpha \beta }|$ and $\boldsymbol{\epsilon }_{\alpha \beta }$ is the minimum separation vector between disc $\alpha$ and $\beta$ which can be computed by a simple extension of the method by (Claeys & Brady Reference Claeys and Brady1993), explained in Appendix A. The near-contact resistance matrix $\boldsymbol{\mathcal{R}}^L$ is computed by summing over the lubrication forces between the discs for which the minimum separation is less than $\varDelta _c$ . Therefore, the mobility matrix $\boldsymbol{\mathcal{M}}$ constructed by taking into account the far-field hydrodynamic interactions and near contact lubrication effects is given by

(3.15) \begin{equation} \boldsymbol{\mathcal{V}} \approx \boldsymbol{\mathcal{M}}\boldsymbol{\cdot } \boldsymbol{\mathcal{F}}, \quad \boldsymbol{\mathcal{M}} \equiv \big[\big(\boldsymbol{\mathcal{M}}^\infty \big)^{-1} + \boldsymbol{\mathcal{R}}^L\big]^{-1}. \end{equation}

Here, $\boldsymbol{\mathcal{F}} = -\boldsymbol{\mathcal{F}}^h$ represents the generalised external force, which in our case consists solely of buoyancy-corrected weights, $\boldsymbol{F}^g$ , for each disc, with zero torques. Although the lubrication force mitigates violent particle approaches, numerical time stepping can still cause unphysical particle overlap, which is generally prevented by incorporating repulsive near-contact forces (Butler & Shaqfeh Reference Butler and Shaqfeh2002; Drazer et al. Reference Drazer, Koplik, Khusid and Acrivos2002; Günther et al. Reference Günther, Frijters and Harting2014; Mari et al. Reference Mari, Seto, Morris and Denn2014; Townsend Reference Townsend2017). We therefore use a generalised external force $\boldsymbol{\mathcal{F}}$ obtained by supplementing the buoyancy-corrected weight $\boldsymbol{F}^g$ of each disc with a repulsive near-contact force (Mari et al. Reference Mari, Seto, Morris and Denn2014; Townsend Reference Townsend2017)

(3.16) \begin{equation} \boldsymbol{F}^r_{\alpha \beta } = f_* e^{-\tau |\boldsymbol{\epsilon }_{\alpha \beta }|}\hat {\boldsymbol{\epsilon }}_{\alpha \beta }. \end{equation}

Based on test runs, we choose $\tau =80/a$ and $f_*=50$ in the simulations. The choice of a large coefficient in $\tau$ ensures that the repulsive force is only active when two particles get extremely close, and, together with the choice of a large $f_*$ , is operational only to prevent overlap of particles.

4. Validation of the dynamics

The mobility matrix $\boldsymbol{\mathcal{M}}$ accounts for both far-field hydrodynamic interactions and near-contact lubrication effects, but omits higher-order effects coming from the second reflection and logarithmic corrections in the lubrication forces, which may not always be negligible. We also emphasise that, while we use the Greens function for an unbounded flow, we are comparing results from experiments performed in a container. As mentioned before, the presence of walls is modelled by cutting off the interaction range. The details of the cutoff are mentioned in the subsequent sections. The current approximation adequately captures the qualitative features of the cluster formation observed in our experiments. The rationale for approximating the wall effects in a minimal way but accounting for the method of reflections and lubrication comes from the fact that we are interested in the nature of cluster formations, which involve the close-range dynamics of discs. One could account for wall effects by placing image singularities beneath the walls, but those effects would come in at $\sim O(a/D)$ , where $D\approx 12.5 a$ is the thickness of the container. Therefore, as long as a given disc has neighbours within a distance of $D/a$ , the dominant contribution to the dynamics will be from the neighbouring discs, rather than the wall effects. We therefore expect that higher-order corrections and wall effects would primarily improve quantitative aspects, such as sedimentation velocities, but approaching such a precision is beyond the scope of this work. Our simulation model does not completely capture all the late-time dynamics of cluster evolution seen in our experiments. Specifically, while pairs of discs that settle into a ‘ $\perp$ ’ structure largely maintain this configuration thereafter in simulations, experimental observations show that only a few discs retain the ‘ $\perp$ ’ arrangement for significantly longer after pairing. However, in experiments involving isolated disc pairs, the long-term stability of the ‘ $\perp$ ’ structure has been confirmed (see figure 3 b). In any case, our analysis focuses on cluster formation, and we refrain from drawing conclusions about their late-time evolution based on our simulations.

Figure 3.

Comparison between experiments and simulations of pairs of discs sedimenting close to each other. (a) A bounded periodic oscillation of discs, categorised under the ‘rocking dynamics’ in Chajwa et al. (Reference Chajwa, Menon and Ramaswamy2019). (b) ‘Hydrodynamic screening’ where one falling disc enters the ‘hydrodynamic shadow’ of the other, forming a ‘ $\perp$ ’ structure. This ‘ $\perp$ ’ configuration is often observed, both experimentally and through simulations, in a one-dimensional lattice of sedimenting discs. (c) Another case of a ‘rocking dynamics’, arising from the linear instability of an equilibrium configuration forming an ‘ $=$ ’ shape.

Figure 3(a) demonstrates that our model successfully captures near-contact periodic oscillations, along with the phenomenon of ‘hydrodynamic screening’ shown in figure 3(b), where one disc enters the ‘hydrodynamic shadow’ of another, forming a ‘ $\perp$ ’ structure. The initial configuration for figure 3(a) is: ${\boldsymbol{x}}_{21}= \{0.68, 0, 1.18 \}$ , $\boldsymbol{p}_1 = \{-0.04, 0, 0.99\}$ and $\boldsymbol{p}_2 = \{0.946, 0, 0.324\}$ . Figure 3(b) has the initial configuration: ${\boldsymbol{x}}_{21}= \{1.5, 0, 1.813 \}$ , $\boldsymbol{p}_1 = \{0.879, 0, 0.476\}$ and $\boldsymbol{p}_2 = \{0.965, 0, -0.26\}$ . It may be noticed from the figure that quantitative details, such as the time it takes for the discs to fall through the same distance in experiments and simulations, differ. An interesting case is that of the ‘rocking dynamics’ (Chajwa et al. Reference Chajwa, Menon and Ramaswamy2019) shown in figure 3(c). This dynamics can be viewed as the evolution of perturbations from a steady-state configuration in which the discs are vertically stacked, with their symmetry axes aligned along gravity, forming an ‘ $=$ ’ shape. This enables a linear stability analysis of such a structural arrangement, which we shall perform in the following subsection, revealing that at least the first reflection is essential to capture the periodic oscillations characteristic of the rocking dynamics, a phenomenon that the point-force approximation ((3.1) and (3.2)) fails to predict at small separations.

4.1. Stability analysis of the ‘ $=$ ’, or horizontal pair, configuration

The configuration where both discs are horizontal, and one lies above the other at a vertical distance of $z^*$ times the size of the disc $a$ , is studied in this subsection in some detail, because it is a canonical example of pair interactions, and because it analytically demonstrates the inability of the point-force approximation to accurately capture the close-range dynamics. The rocking dynamics shown in figure 3(c) exhibits periodic oscillations in the orientations and positions of the pair of discs about the ‘ $=$ ’ configuration represented by $(|\boldsymbol{p}_1\boldsymbol{\cdot } \boldsymbol{g}|,\, |\boldsymbol{p}_2\boldsymbol{\cdot } \boldsymbol{g}|)=(1,1)$ . Here, $\boldsymbol{g}$ is a unit vector along gravity (along the $-z$ direction). The ‘ $=$ ’ configuration is a fixed point in the frame falling with the pair of discs. The system is non-dimensionalised using the semi-major axis $a$ of the discs as the length scale and $\tau _p={6\pi \mu a^2}/F$ as the time scale, where $F$ is the buoyancy-corrected weight of the spheroids. We shall henceforth work with dimensionless quantities, unless stated otherwise. Let the position of the disc 2 (upper disc) relative to disc 1 (lower disc) be given by $(\delta x,z^* + \delta z)$ and their corresponding orientation vectors be given by: $\boldsymbol{p}_1=(\cos (\pi /2+\delta \theta _1), \sin (\pi /2+\delta \theta _1))$ and $\boldsymbol{p}_2=(\cos (\pi /2+\delta \theta _2), \sin (\pi /2+\delta \theta _2))$ ; with their fixed point position and orientations being $(0,z^*)$ and $\boldsymbol{p}_1^*=\boldsymbol{p}_2^*=(0,1)$ , respectively. The dynamics of the pair, using the first reflection, leads to the following linear stability analysis:

(4.1) \begin{equation} \frac {\rm d}{{\rm d}t}\begin{pmatrix} \delta x \\ \delta \theta _1 \\ \delta \theta _2 \end{pmatrix} = \unicode{x1D63D}_{R_1} \begin{pmatrix} \delta x \\ \delta \theta _1 \\ \delta \theta _2 \end{pmatrix};\quad \unicode{x1D63D}_{R_1}=\begin{pmatrix} 0 & \kappa _0(e) & -\kappa _0(e) \\ J_1 & J_2 & J_3 \\ -J_1 & -J_3 & -J_2 \end{pmatrix}, \end{equation}

where

(4.2) \begin{align}&\qquad\qquad\qquad\qquad\qquad\quad \kappa _0(e) \equiv \left (\frac {1}{Y^A}-\frac {1}{X^A}\right ), \\[-8pt]\nonumber\end{align}

(4.3) \begin{align}&\qquad J_1 \equiv \frac {-3}{4} \int _{-1}^1\frac {3}{8}{\rm d}\xi _1\int _{-1}^1 {\rm d}\xi _2 \frac {1-\xi _1^2}{\tilde \xi _0^3} \left \{ 1-\frac {e^2}{2-e^2} \left [ -3 +\frac {24\kappa _1}{\tilde \xi _0^2} \right ] \right \},\\[-8pt]\nonumber \end{align}

(4.4) \begin{align} J_2 &\equiv \frac {-3}{4} \int _{-1}^1\frac {3}{8}{\rm d}\xi _1\int _{-1}^1 {\rm d}\xi _2 \frac {1-\xi _1^2}{\tilde \xi _0^3} \Bigg \{ i e\xi _1 - \frac {e^2}{2-e^2} \Bigg [ (z^*+i e\xi _2)\left ( -3 +\frac {30\kappa _1}{\tilde \xi _0^2} \right ) \nonumber \\ &\quad - \frac {6\kappa _1 (2z^*-ie\xi _1+2ie\xi _2)}{\tilde \xi _0^2} \Bigg ] \Bigg \}, \end{align}

and

(4.5) \begin{equation} J_3 \equiv \frac {-3}{4} \int _{-1}^1\frac {3}{8}{\rm d}\xi _1\int _{-1}^1 {\rm d}\xi _2 \frac {1-\xi _1^2}{\tilde \xi _0^3} \left \{-1+\frac {e^2}{2-e^2} \left [ -3 +\frac {24\kappa _1}{\tilde \xi _0^2} \right ] \right \}i e\xi _2. \end{equation}

For shorthand notation, we define the following quantities appearing in integrands:

(4.6) \begin{equation} \tilde \xi _0\equiv \sqrt {(z^*+ie(\xi _2-\xi _1))^2} \quad \textrm {and} \quad \kappa _1 \equiv \frac {3/2-\xi _1^2-\xi _2^2}{4}. \end{equation}

The vertical perturbation $\delta z$ appears only at second order in the algebra, and therefore does not feature in the linear stability analysis. The corresponding Jacobian matrix of the perturbation evolution in the point-force limit, obtained from (3.1) and (3.2), is given by

(4.7) \begin{equation} \unicode{x1D63D}_{pt} = \begin{pmatrix} 0 & \kappa _0(e) & -\kappa _0(e) \\[2pt] \dfrac {-3}{4z^3}\left [1+\dfrac {3e^2}{2-e^2}\right ] & \dfrac {-3}{4z^2}\left [\dfrac {3e^2}{2-e^2}\right ] & 0 \\[10pt] \dfrac {3}{4z^3}\left [1+\dfrac {3e^2}{2-e^2}\right ] & 0 & \dfrac {3}{4z^2}\left [\dfrac {3e^2}{2-e^2}\right ] \end{pmatrix} .\end{equation}

Figure 4.

The frequency of oscillations in the rocking dynamics about the vertically stacked, or ‘ $=$ ’ configuration, as a function of vertical separation $z^*$ at the fixed point. The aspect ratio here is $b/a=0.125$ , which is close to our experimental discs. The black vertical line denotes $z=z_c$ below which the point-force approximation fails to predict the rocking dynamics. The inset shows the dependence of the $z_c$ on the aspect ratio.

Under the first-reflection approximation, the eigenvalues of $\unicode{x1D63D}_{R_1}$ are purely imaginary for any $z^*$ and hence the linear dynamics is periodic about the ‘ $=$ ’ configuration. The oscillation frequency of this rocking dynamics is shown in figure 4. The point-force approximation, on the other hand, leads us to the incorrect conclusion that the fixed point ‘ $=$ ’ is a saddle below a certain value of $z=z_c$ , given by

(4.8) \begin{equation} z_c \equiv \left [\frac {16}{27}\left ( \frac {2+e^2}{e^4} - 1\right )\kappa _0(e) \right ]^{-1}. \end{equation}

A linear stability analysis using the point-force approximation thus cannot predict the rocking dynamics seen in the experiment for small disc separations, which shows that the first-reflection correction is crucial in this regime. Note that $z^*$ is bounded below by $2b/a$ , which corresponds to the case where the discs just touch.

5. Planar dynamics of a one-dimensional lattice of sedimenting discs

In this section, we study the sedimentation dynamics of a regular one-dimensional lattice of identical spheroids with their symmetry axes aligned with the horizontal. The corresponding problem for spheres was studied by Crowley (Reference Crowley1971), who established that any arbitrary perturbation in the initial positions of the spheres will always grow exponentially, i.e. the system is always unstable. Crowley also proposed a mechanism to explain the growth of perturbations, which leads to the formation of clumps of spheres. Chajwa et al. (Reference Chajwa, Menon, Ramaswamy and Govindarajan2020) investigated the planar sedimentation dynamics of a regular one-dimensional lattice of identical spheroids, both experimentally as well as theoretically, and showed, in contrast to spheres, that a system of spheroids is not always exponentially unstable: there is a regime where linear perturbations display algebraic rather than exponential growth. In the algebraically growing regime, Chajwa et al. (Reference Chajwa, Menon, Ramaswamy and Govindarajan2020) found that the modal perturbations are neutrally stable.

The analysis presented in Chajwa et al. (Reference Chajwa, Menon, Ramaswamy and Govindarajan2020) was based on the point-force approximation (see (3.1) and (3.2)) with nearest-neighbour interactions. In this section, we show that this picture is crucially modified when we take into account more accurate hydrodynamic interactions through the first reflection. We find that the qualitative behaviour of modes changes from pure oscillations to weak exponential growth – the sharp boundary between these two regions is lifted, and the growth rate changes smoothly over the entire regime rather than becoming strictly zero.

Figure 5.

Experiments (a) and (c), and simulations (b) and (d) showing clumping due to the Crowley mechanism. Panels (a) and (c) are aligned vertically, as are (b) and (d), to highlight that the locations where the clumps form correspond to the red ovals in (a) and (b). The lattice spacing and the dominant perturbation wavenumber are $\tilde {d}=1.875$ and $q=\pi /6$ , respectively. Panels (a) and (b) show early evolution ( $t\sim 5$ ). The valleys, marked by red ovals, represent regions where the discs begin to clump together. Panels (c) and (d) show late-time evolution ( $t\sim 15$ ). Clumps, occurring at the valleys, consist of three or more discs.

Figure 6.

Experiments (a) and (c), and simulations (b) and (d) showing pairing due to the drift mechanism. Panels (a) and (c) are aligned vertically, as are (b) and (d), to highlight that the locations where the pairs form correspond to the red boxes in (a) and (b). The lattice spacing and the perturbation wavenumber are $\tilde {d}=3.75$ and $q=\pi /2$ , respectively. Panels (a) and (b) show initial evolution ( $t\sim 20$ ). The valleys are highlighted with green ovals, while the discs that ultimately pair up are enclosed within red rectangular boxes, and are seen to be located away from the valleys. Panels (c) and (d) show late-time evolution ( $t\sim 45$ ). There is pairing but no clumping, and the pairing dynamics is characterised by two discs coming together in the form of a ‘ $\perp$ ’, or inverted ‘T’.

For ease of distinguishing, we refer to clusters of three or more spheroids as clumps, and clusters of two spheroids as pairs. Visuals of disc clusters reveal that, when clustering occurs, clumping takes place in the valleys when the Crowley mechanism dominates over the drift mechanism (figure 5). However, when the drift mechanism is stronger than the Crowley mechanism, there are no clumps, but pairing occurs away from the valley. Pairs typically take on a ‘ $\perp$ ’ shape, see figure 6.

To numerically investigate the distinction between clumps and pairs, as observed in the experiments, we use two different lattice spacings and wavenumbers, as indicated by the two coloured dots in figure 7. The simulations are conducted by solving (3.15) as a system of ordinary differential equations using the Runge–Kutta–Fehlberg method for time integration (Fehlberg Reference Fehlberg1969). Adaptive time stepping is employed with an error tolerance of $10^{-4}$ and a maximum step size of $0.1$ . In the Crowley-mechanism-dominated regime with lattice spacing $\tilde {d}=1.875$ and dominant perturbation wavenumber $q=\pi /6$ (red dot in figure 7), the resulting clumps formed at the valleys are shown in figure 5. In the drift-mechanism-dominated regime with lattice spacing $\tilde {d}=3.75$ and dominant perturbation wavenumber $q=\pi /2$ (green dot in figure 7), pairwise clumping is observed, as shown in figure 6.

Figure 7.

Contour plot of the log of the growth rate $\sigma$ (the only positive real part among the eigenvalues of $\unicode{x1D63C}_{R_1}$ ) for aspect ratio $b/a = 0.125$ . The red curve denotes the critical lattice spacing $\tilde {d}_c$ above which the point-force approximation (Chajwa et al. Reference Chajwa, Menon, Ramaswamy and Govindarajan2020) predicts neutral stability. The nonlinear evolution corresponding to the red dot, where $(q,\tilde {d})=(\pi /6,1.875)$ and the green dot, where $(q,\tilde {d})=(\pi /2,3.75)$ , were shown to support clustering and pairing, respectively (figures 5 and 6).

There are three degrees of freedom per disc: the horizontal position $x_\alpha$ of disc $\alpha$ , its vertical position $z_\alpha$ and its orientation angle $\theta _\alpha$ with respect to the horizontal ( $x$ ) axis. About the fixed point $\tilde {d}(\mathbb{Z},0,0)$ (this fixed point is in the frame of reference falling with the average terminal velocity of the undisturbed lattice), $\tilde {d}$ being the dimensionless lattice spacing, non-dimensionalised using the disc radius $a$ , an infinitesimal perturbation of the disc $\alpha \in \mathbb{Z}$ is given by $(\hat {x}_q (t),\hat {z}_q (t),\hat {\theta }_q (t)) e^{iq\alpha }$ , where $q$ is the dimensionless perturbation wavenumber. Under the first reflection and nearest-neighbour interactions, this perturbation evolves as

(5.1) \begin{equation} \frac {d}{{\rm d}t}\begin{bmatrix} \hat {x}_q \\[6pt] \hat {z}_q \\[6pt] \hat {\theta }_q \\[6pt] \end{bmatrix} = \unicode{x1D63C}_{R_1} \begin{bmatrix} \hat {x}_q \\[6pt] \hat {z}_q \\[6pt] \hat {\theta }_q \\[6pt] \end{bmatrix}; \quad \unicode{x1D63C}_{R_1} \equiv \begin{bmatrix} 0 & I_{12} & \kappa _0(e) + I_{13} \\[6pt] I_{21} & 0 & 0 \\[6pt] I_{31} & 0 & 0 \end{bmatrix}, \end{equation}

where $\kappa _0$ is given by (4.2). The terms in the matrix $\unicode{x1D63C}_{R_1}$ are given by

(5.2) \begin{align}&\,\,\,\quad I_{12} \equiv \dfrac {3}{2} \int _{-1}^1\frac {1}{2}{\rm d}\xi _1\int _{-1}^1\frac {1}{2}{\rm d}\xi _2 i\sin q \frac {i e \xi _{12}-\tilde {d}}{\xi _{0}^3}\left [1-\frac {6\kappa _2}{\xi _0^2} \right ], \\[-8pt]\nonumber\end{align}

(5.3) \begin{align}&\quad\, I_{21} \equiv -\dfrac {3}{2} \int _{-1}^1\frac {1}{2}{\rm d}\xi _1\int _{-1}^1\frac {1}{2} {\rm d}\xi _2 i\sin q \frac {i e \xi _{12}-\tilde {d}}{\xi _{0}^3}\left [1+\frac {6\kappa _2}{\xi _0^2} \right ], \\[-8pt]\nonumber\end{align}

(5.4) \begin{align} I_{13} &\equiv -\dfrac {3}{2} \int _{-1}^1\frac {1}{2}{\rm d}\xi _1\int _{-1}^1 \frac {1}{2}{\rm d}\xi _2 i e[\xi _1-\xi _2\cos (q)] \frac {i e \xi _{12}-\tilde {d}}{\xi _{0}^3}\left [1-\frac {6\kappa _2}{\xi _0^2} \right ], \end{align}

(5.5) \begin{align} I_{31} &\equiv \frac {3}{2} \int _{-1}^1\frac {3}{8}{\rm d}\xi _1\int _{-1}^1 {\rm d}\xi _2 \frac {1-\xi _1^2}{\xi _0^3} \left \{[1-\cos (q)]\left (-2+\frac {6e^2 \kappa _1}{(2-e^2)\xi _0^2} \right )\right . \nonumber \\ &\quad\left . + \frac {30e^2\kappa _1}{\left(2-e^2\right)\xi _0^2} \frac {i\sin q}{\left(i e\xi _{12}-\tilde {d}\right)} \right \}, \end{align}

where $\xi _{12} \equiv \xi _1-\xi _2$ , $\xi _0\equiv \sqrt {(i e\xi _{12}-\tilde {d})^2}$ and $\kappa _2 \equiv (2-\xi _1^2-\xi _2^2)/4$ . The eigenvalues of $\unicode{x1D63C}_{R_1}$ are $\{0, \sigma +i\omega , -\sigma -i\omega \}$ , where $\sigma$ represents the growth rate of perturbations, and $\omega$ denotes the oscillation frequency. The structure of the eigenvalues ensures that the system is always unstable provided $\sigma \neq 0$ .

On the other hand, under the point-force approximation, the corresponding Jacobian matrix $\unicode{x1D63C}_{pt}$ for the evolution of perturbations is given by (Chajwa et al. Reference Chajwa, Menon, Ramaswamy and Govindarajan2020)

(5.6) \begin{equation} \unicode{x1D63C}_{pt} \equiv \begin{bmatrix} 0 & \dfrac {-3 i \sin q}{2\tilde {d}^2} & \kappa _0(e) \\[10pt] \dfrac {3 i \sin q}{2\tilde {d}^2} & 0 & 0 \\[10pt] \dfrac {-6\sin ^2 q/2}{\tilde {d}^3} & 0 & 0 \end{bmatrix}. \end{equation}

For a given wavelength, the point-force approximation predicts a neutrally stable regime characterised by a zero growth rate beyond a critical lattice spacing $\tilde {d}_c$ , given by

(5.7) \begin{equation} \tilde {d}_c = \frac {3}{8\kappa _0(e)} \frac {\sin ^2 q}{\sin ^2 q/2}. \end{equation}

Figure 7 shows the logarithm of the growth rate of perturbations as a function of lattice spacing $\tilde {d}$ and perturbation wavenumber $q$ . While the point-force approximation predicts zero growth rate above the red curve, the first reflection has a small but finite growth rate at high lattice spacings. The growth rate shows a power-law behaviour $d^{-\gamma }$ for large lattice spacing, with the dependence of the power on the wavenumber shown in figure 8. One must interpret this power law with caution, as wall effects may alter the exponent. This requires the lattice spacings in the power-law regime to be smaller than the container thickness. Numerical data show that the scaling exponent in the large- $d$ limit is best fit by $\gamma \approx -4.5$ across perturbation wavenumbers, for a large range of spheroids with aspect ratios ranging from 0.19 to 0.8.

Figure 8.

(a) Perturbation growth rate $\sigma$ as a function of the lattice spacing for a wavenumber $q=\pi /6$ and aspect ratio $b/a=0.125$ . For large lattice spacing $\tilde {d}$ , the growth rate shows a power-law behaviour $\sigma \sim \tilde {d}^\gamma$ with $\gamma \approx -4.5$ . (b) The variation of the exponent $\gamma$ with wavenumber $q$ . Note that for a lattice of spheres the growth rate decays much more slowly, as $\sigma \sim \tilde {d}^{-2}$ , for all wavenumbers.

Figure 9.

The $\varepsilon$ -pseudospectrum $\sigma _{\varepsilon }(\unicode{x1D63C}_{pt})$ of $\unicode{x1D63C}_{pt}$ for the aspect ratio $b/a = 0.125$ . The blue dots denote the eigenvalues of $\unicode{x1D63C}_{pt}$ and the black dots denote the eigenvalues of $\unicode{x1D63C}_{R_1}$ . The matrix perturbation size $\varepsilon$ leads to the spreading of eigenvalues inside the yellow regions. The values of $\varepsilon$ are chosen so that the black dots lie at the boundary of the yellow regions, indicating that the first-reflection corrections to $\unicode{x1D63C}_{pt}$ are of the order of $\varepsilon$ : (a) $\varepsilon =4.8\times 10^{-3}$ for the set of parameters $(q,\tilde {d})=(\pi /2,3.75)$ (green dot in figure 7); (b) $\varepsilon =1.2\times 10^{-2}$ for the set of parameters $(q,\tilde {d})=(\pi /6,1.875)$ (red dot in figure 7). Note that the pseudospectra are much larger than $\epsilon$ .

The first reflection thus shows a qualitative change in the stability of the system as compared with the point-force approximation: the system is now linearly unstable at any lattice spacing, although with a growth rate that approaches zero with increasing spacing. An examination of the point-force approximation shows that this is to be expected. The non-normality of $\unicode{x1D63C}_{pt}$ provides for a possible large change in eigenvalues for a small perturbation in the stability operator ( $\unicode{x1D63C}_{pt}$ becoming $\unicode{x1D63C}_{R_1}$ ). Since in the neutral regime the eigenvalues of $\unicode{x1D63C}_{pt}$ lie on the imaginary axis, its pseudospectrum must protrude into the positive real axis, and $\unicode{x1D63C}_{R_1}$ is a member of the pseudospectrum of $\unicode{x1D63C}_{pt}$ which produces instability. The $\varepsilon$ -pseudospectrum $\sigma _{\varepsilon }(\unicode{x1D63C}_{pt})$ of $\unicode{x1D63C}_{pt}$ is defined as (Trefethen Reference Trefethen1997)

(5.8) \begin{equation} \sigma _{\varepsilon }(\unicode{x1D63C}_{pt}) \equiv \left\{z \in \mathbb{C}\, :\, \big|\big|\left(z-\unicode{x1D63C}_{pt}\right)^{-1}\big|\big|\gt \varepsilon ^{-1} \right\}. \end{equation}

Figures 9(a) and 9(b) show the pseudospectra of $\unicode{x1D63C}_{pt}$ corresponding to the conditions of figures 5 and 6, respectively. In each case, $\varepsilon$ is chosen such that the largest eigenvalue of $\unicode{x1D63C}_{R_1}$ lies on the boundary of the pseudospectrum. While the oscillation frequency remains close to its original value, a non-zero growth rate appears, which changes the dynamics qualitatively. Therefore, in the drift regime (green dot in figure 7), the first-reflection correction has a gentle destabilising effect. In contrast to this, the first reflection has a significant stabilising effect in the Crowley regime, and lies on a larger $\varepsilon$ -pseudospectrum in figure 9(b). The net growth rate remains positive (red dot in figure 7).

With the first reflection bringing a qualitative change in the nature of stability, one may surmise that modelling wall effects beyond a simple nearest-neighbour cutoff may lead to further qualitative change. This is not the case because the qualitative change occurred due to the structural instability of the system described by $\unicode{x1D63C}_{pt}$ : the eigenvalues lie exactly on the imaginary axis. With the first reflection, these eigenvalues are nudged off the imaginary axis, leading to a qualitative change. Inclusion of wall effects using image singularities might further nudge the eigenvalues, but they are not primarily responsible for the qualitative change already brought in by the more dominant first-reflection effects.

If only eigenvalue-based growth were in play, the extremely slow growth rate for large lattice spacing would contribute negligibly to perturbation growth during the time of our simulations. However, a large algebraic transient growth of the perturbations happens in this regime, due to the non-normal nature of $\unicode{x1D63C}_{R_1}$ . To highlight this, we look at the matrix norm of $e^{\unicode{x1D63C}_{R_1}t}$ for time $t\in [0,2\pi /\omega ]$ , where $\omega$ is the oscillation frequency given by the imaginary parts of non-zero eigenvalues of $\unicode{x1D63C}_{R_1}$ , see figure 10. The maximum non-modal growth can be quantified by the quantity $G_r$ defined as

(5.9) \begin{equation} G_r \equiv \frac { \big|\big|e^{\unicode{x1D63C}_{R_1} t^*} \big|\big|^2}{e^{2\sigma t^*}}; \quad t^* \equiv \mathop {\mathrm{arg\,max}}\limits _{t\in [0,2\pi /\omega ]}\big|\big|e^{\unicode{x1D63C}_{R_1}t}\big|\big| , \end{equation}

where $\sigma$ is the growth rate given by the real parts of eigenvalues of $\unicode{x1D63C}_{R_1}$ and the matrix norm can be obtained using the singular value decomposition. Figure 10 shows the logarithm of $G_r$ as a function of lattice spacing $\tilde {d}$ and perturbation wavenumber $q$ . By definition $\log (G_r) \geqslant 0$ . In the drift regime, the lattice is disrupted on the time scale $t^*$ .

Figure 10.

(a) Time evolution of the matrix norm of $e^{\unicode{x1D63C}_{R_1}t}$ and the exponential growth $e^{2\sigma t}$ for the aspect ratio $b/a=0.125$ . The non-normal nature of $\unicode{x1D63C}_{R_1}$ leads to transient growth of non-modal perturbations, which can be much larger than the exponential growth. The parameters $(q,\tilde {d})$ used to calculate $\unicode{x1D63C}_{R_1}$ correspond to the green dot in figure 7. The non-modal growth reaches its maximum at time $t=t^*$ . (b) Contour plot of the $\log _{10}$ of $G_r$ defined in (5.9). In regimes of small $\log _{10}(G_r)$ , the growth is exponential but very slow. Large values of $\log _{10}(G_r)$ indicate that transient algebraic growth will dominate.

5.1. Distinct clumping mechanisms in the exponential and algebraic growth regimes

As explained in Chajwa et al. (Reference Chajwa, Menon, Ramaswamy and Govindarajan2020), there are two mechanisms at play that decide the nature of clustering. The Crowley mechanism acts to form clumps at the valleys while the drift mechanism counteracts it (see figure 1). The fast decay of growth rate with increased lattice spacing (see figure 8) is because of the drift mechanism competing better against the Crowley mechanism. For a lattice of spheres the drift mechanism is absent, and the growth rate follows $\sigma \sim \tilde {d}^{-2}$ for large lattice spacing $\tilde {d}$ . In contrast, for increasingly flatter oblate ellipsoids, the growth rate is best fit by $\sigma \sim \tilde {d}^{-4.5}$ , as shown in figure 8.

Our objective is to develop a statistical measure that can differentiate between clustering in the Crowley regime and the drift regime. The drift mechanism leads to aggregation in pairs, whereas in the Crowley regime, a larger number of discs participate in forming an aggregate, thus suggesting $P_s(t)$ , the fraction of particles participating in aggregates of size $s$ , as a useful quantitative measure. Discs that are within a distance of $0.5a$ of each other are defined to be part of an aggregate. We conduct 500 simulations in both the Crowley and the drift regimes, each involving a lattice of 60 discs with periodic boundary conditions along the horizontal ( $x$ ) direction. Each disc interacts with eight neighbours on either side. Changing the number of neighbours about this cutoff scale does not affect cluster sizes and types, provided there are a sufficient number of discs interacting with their neighbouring ones upon spatial rearrangement as the system evolves. Although the effect of a distant disc (with the distance $r \gtrsim D/a$ , $D\approx 12.5 a$ being the thickness of the container) is weaker than that of the walls, the dynamics is primarily governed by the nearest neighbours. Thus, including the eighth disc while neglecting the wall effects has little impact, as both contribute marginally compared with the dominant influence of nearby discs. Each simulation is initialised with sinusoidal perturbations only in the $x$ -direction, at a chosen dominant wavenumber, with amplitude $0.625a$ , and additional random perturbations of amplitude $0.075 a$ , where $a$ is the semi-major axis of the discs. The initial orientation angle is chosen from a uniform distribution of width $8$ degrees. The amplitude of the perturbations and the noise in initial conditions emulate the perturbations in the experiments.

Figure 11.

Probability distribution $P_s(t)$ as a function of non-dimensional time for different aggregate size $s$ , in the Crowley regime $(q,\tilde {d})=(\pi /6,1.875)$ . Here, $P_s(t)$ is calculated as the ensemble average of the fraction of discs forming an aggregate of size $s$ relative to the total number of discs in each ensemble. Discs that are within the distance of $0.5a$ are considered to be aggregated or clustered. The error bars represent the standard deviation.

Figure 12.

Probability distribution $P_s(t)$ as a function of non-dimensional time for different aggregate size $s$ , in the drift regime $(q,\tilde {d})=(\pi /2,3.75)$ . Here, $P_s(t)$ is calculated as the ensemble average of the fraction of discs forming an aggregate of size $s$ relative to the total number of discs in each ensemble. The pairs of discs separate out in a ‘ $\perp$ ’ configuration from the rest and keep on settling together, maintaining their orientations.

Figure 11 shows the time evolution in the Crowley regime of the fraction of discs participating in clusters of different sizes. The 500 samples are divided into 20 sub-samples, over which the mean and standard deviation are calculated. The most frequently occurring clump consists of four discs, while the largest clump contains five discs. Upon examination, we find no preferred relative orientations among the discs within these clumps. A key characteristic of the drift regime is that pairing is by far the most frequent type of clustering observed. Each pair adopts a ‘ $\perp$ ’ shape (see figure 6). This is an instance of hydrodynamic screening, as discussed in § 4 and figure 3. With time, the paired discs separate out from the rest of the discs because they settle faster than isolated discs, usually without ever separating from each other. As a result, the fraction of pairs (or $s=2$ aggregates) only increases with time, as shown in figure 12. The relative orientation of a pair of discs is quantified by $\boldsymbol{p}_1\boldsymbol{\cdot }\boldsymbol{g}$ and $\boldsymbol{p}_2\boldsymbol{\cdot }\boldsymbol{g}$ . Here, $\boldsymbol{g}$ is the gravity unit vector and $\boldsymbol{p}_1$ and $\boldsymbol{p}_2$ are respectively the orientations of the lower and upper discs of a pair. As previously mentioned, the discs come together in a ‘ $\perp$ ’ configuration, i.e. $(|\boldsymbol{p}_1\boldsymbol{\cdot } \boldsymbol{g}|,\, |\boldsymbol{p}_2\boldsymbol{\cdot } \boldsymbol{g}|)=(1,0)$ as they begin to become a pair, as evident in figure 13. These ‘ $\perp$ ’ configurations, also observed in our experiments, are solely governed by two-body hydrodynamic interactions and are found to remain stable over time in numerical simulations. Although the experiments show the ‘ $\perp$ ’ signature at the time of pair formations, the pair orientations exhibit significant variability as the discs approach the base of the container. Given the uncertainty in the factors influencing the persistence of ‘ $\perp$ ’ structures in the experiments, we refrain from making conclusions about the late-time ( $t \gtrsim 60$ ) dynamics and evolution of clustered pairs based on our simulations.

Figure 13.

Probability density plot of orientations of pairs of discs in the drift regime $(q,\tilde {d})=(\pi /2,3.75)$ . Here, $\boldsymbol p_1$ refers to the orientation of disc 1, which is below disc $2$ , and $\boldsymbol{g}$ is the gravity unit vector. The black-green dots are experimentally observed orientations, and the background colour indicates the frequency of occurrence of a given pair of orientations in the simulations. A predominance of the ‘ $\perp$ ’ formation is seen in both experiments and simulations. The orientation statistics are acquired when the pairs start forming in the simulations, at a non-dimensional time of $t\sim 45$ . This time is later than $t^*$ , the time at which the non-normal growth reaches its maximum, as shown in figure 10. The black-green dots, indicating experimental samples, are taken when pair formation is observed.

6. Formation of the ‘ $\perp$ ’ shape

The phenomenon of pairing through drift merits attention as a distinctive property of the settling of non-spherical objects, ruled out for the Stokesian settling of two spheres. As two spheroids approach each other, they invariably rearrange themselves into a $\perp$ structure, both in our experiments and in our simulations. This process can be explained by a simple analysis similar to that of Koch & Shaqfeh (Reference Koch and Shaqfeh1989) (see figure 1 in their paper). According to their description, the neighbours of a disc tend to align their ‘thin’ sides parallel to the direction of the extensional component of the flow generated by that disc. However, we point to a key difference coming from the vorticity disturbance due to the rotation of the spheroid. Consider a force monopole at the origin and a disc of radius $a$ settling in its flow field at a radial distance $r$ , with position $(x,0,z)$ . With contributions from both the vorticity and strain rate of the flow produced by the monopole, it follows from (3.2) that the orientation of the disc $\boldsymbol{p}=(\cos \theta ,0,\sin \theta )$ evolves, to leading order in $a/r$ as

(6.1) \begin{equation} \frac {{\rm d}\theta }{{\rm d}t} = \frac {3}{4r^2} [\cos \phi -\zeta \sin \phi \sin 2(\phi -\theta ) ]; \quad \zeta \equiv \frac {3e^2}{4-2e^2}. \end{equation}

Here, $x=r\cos \phi$ and $z=r\sin \phi$ . Equation (6.1) is non-dimensionalised using the same scales as described in § 4. Note that the first term on the right-hand side of (6.1) represents the vorticity contribution, and the second term represents the contribution of the strain rate to the disc’s rotation.

Now, without the vorticity contribution, one gets the stable orientations seen in the schematic figure 1 of Koch & Shaqfeh (Reference Koch and Shaqfeh1989) and thereby an inward particle flux from all directions. In our study, both the vorticity and strain-rate contributions are crucial to the disc’s rotation. To find the stable orientations while incorporating the vorticity contribution, we proceed by defining the angle

(6.2) \begin{equation} \beta \equiv \frac {1}{2} \arcsin \frac {\cot \phi }{\zeta } . \end{equation}

For any $(r,\phi )$ , we set ${\rm d}\theta /{\rm d}t=0$ in (6.1) to get the fixed point

(6.3) \begin{equation} \theta ^* = \phi - \frac {n\pi }{2}-(-1)^n\beta , \quad n\in \mathbb{Z} ,\end{equation}

for the orientation angle $\theta$ . We can immediately see from (6.3) that solutions do not exist for all values of $\phi$ , which is a qualitatively distinct feature from the picture of Koch & Shaqfeh (Reference Koch and Shaqfeh1989). For $\varOmega \equiv \arctan (1/\zeta )$ , the solution for ${\rm d}\theta /{\rm d}t =0$ is defined only for the values of $\phi$ in the intervals $(\varOmega , \pi - \varOmega )$ and $(\pi + \varOmega , 2\pi - \varOmega )$ . Even where the solution exists, we show below that the stable orientations are different from what one gets with just the strain-rate contribution.

We can study the stability by considering just the solutions $\theta ^*_1 = \phi - \beta$ and $\theta ^*_2 = \phi + \beta - \pi /2$ , given the apolarity of discs. Perturbing the fixed points by small angles $\delta \theta _1$ and $\delta \theta _2$ respectively, substituting in (6.1) and Taylor expanding about the respective fixed points gives

(6.4) \begin{equation} \dot {\delta \theta }_k = (-1)^{k-1} \cos 2 \beta \, \frac { 3 \zeta \sin \phi }{2 \, {r^{2}}}\, \, \delta \theta _k , \quad k\in \{1,2\}. \end{equation}

Since $\cos 2\beta \gt 0$ , the fixed point $\theta ^*_1$ is unstable in the upper half-plane and stable in the lower half-plane, and vice versa for the other fixed point $\theta ^*_2$ .

Figure 14.

Stable orientations with both vorticity and strain rate contributing to a disc’s rotation: with a force monopole at the origin. Blue double-ended arrows and red-dashed lines are stable and unstable alignments of the thin side of the disc, respectively. The angle between stable and unstable orientations is $\pi /2 - 2\beta$ . The shaded grey region represents the range of angles $\phi$ where there are no stable orientations, with $\varOmega = \arctan (1/\zeta )$ .

To summarise, in the upper half-plane (when the disk is above the force monopole), the disc is stable with its orientation along the angle $\theta _2^* = \phi - \beta$ . In the lower half-plane, the stable angle of the orientation is $\theta _1^* = \phi + \beta - \pi /2$ . These stable orientations are shown in figure 14 with blue arrows indicating the stable alignment of the thin side of the disc. (Note that the thin side of the disc is perpendicular to its orientation vector.) In both half-planes, solutions exist only with the azimuthal angle $\phi$ lying in the intervals $(\varOmega , \pi - \varOmega )$ and $(\pi + \varOmega , 2\pi - \varOmega )$ , leaving out regimes shown as grey shaded region in figure 14, where there are no stable orientations.

With this picture in mind, we can explain the formation of ‘ $\perp$ ’ shapes as follows. For a pair of vertically separated settling discs, the force monopole flow created by the lower disc aligns the upper disc vertically, and that generated by the upper disc aligns the lower disc horizontally. Since a vertically aligned disc sediments more rapidly than a horizontal one, the upper disc catches up with the lower one, forming the ‘ $\perp$ ’ shape, frequently observed in disc pairs. In the drift regime, the lattice spacing is large enough, i.e. the system is dilute enough, for the above approximate analysis based on pair interactions and point-force flows to suffice. The absence of this ‘ $\perp$ ’ shape in the Crowley regime indicates that clusters in this regime are influenced by multi-particle interactions.