1. Introduction
Particle-laden flows arise in a wide range of natural and engineered systems, including sediment transport in rivers, atmospheric aerosols, biological suspensions and solid–liquid separation processes. In these systems, suspended particles interact with the surrounding fluid and with one another through long-range hydrodynamic disturbances, giving rise to collective dynamics across multiple spatial and temporal scales. Among the canonical manifestations of these interactions is sedimentation, in which particles settle under gravity or an imposed body force. Despite its apparent simplicity, sedimentation remains a fundamental benchmark problem in fluid mechanics because the emergent dynamics result from a multiscale coupling between particle geometry, suspension microstructure and disturbance–flow propagation (O’Melia Reference Nicolai, Herzhaft, Hinch, Oger and Guazzelli1998; van der Lans et al. Reference Uzair, Qasim, Bhatti, Zaidi, Molina, Yamamoto and Hamid2000; Ji et al. Reference Ji, Xiong, Wu, Chen and Wu2009).
For isolated particles settling in a viscous fluid, the motion approaches a terminal velocity determined by the balance between hydrodynamic drag and effective particle weight. In suspensions, however, disturbance flows generated by neighbouring particles modify the local fluid environment and reduce the mean settling velocity relative to the single-particle limit. This hindered-settling behaviour has been extensively characterised for spherical particles. Classical analyses by Batchelor (Reference Batchelor1972) and Batchelor & Wen (Reference Batchelor and Wen1982) established theoretical scaling laws for dilute Stokesian suspensions, while empirical relations such as the Richardson–Zaki correlation (Richardson & Zaki Reference Qiu, Lin, Zhou, Wang, Yuan, Wei, Dai, Luo and Chen1954) remain widely used to describe the concentration dependence of the mean settling velocity. Subsequent studies further connected hindered settling to suspension microstructure and weak inertial effects at low-to-intermediate volume fractions (Brenner Reference Brenner1999; Di Felice Reference Di Felice1999).
Beyond the mean settling rate, sedimenting suspensions develop velocity fluctuations and evolving microstructures that strongly influence macroscopic transport. Because hydrodynamic interactions are intrinsically long-ranged, early theoretical analyses predicted that fluctuation amplitudes could diverge with system size (Caflisch & Luke Reference Caflisch and Luke1985; Tee et al. Reference Stokes2002). Experiments and simulations demonstrated instead that these fluctuations saturate once the domain exceeds a characteristic correlation length associated with hydrodynamic screening (Nicolai et al. Reference Nakayama and Yamamoto1995; Ladd Reference Ladd1996, Reference Ladd1997; Segre et al. Reference Sanjeevi, Kuipers and Padding1997). In the Stokes regime, the resulting fluctuation amplitudes and relaxation times exhibit robust power-law scalings with particle volume fraction, providing a central framework for understanding collective sedimentation dynamics and hydrodynamic screening in particulate suspensions (Brenner Reference Brenner1999; Cunha et al. Reference Cunha, Abade, Sousa and Hinch2002; Padding & Louis Reference O’Melia2008; Hamid et al. Reference Hamid, Molina and Yamamoto2014; Shakeel et al. Reference Segre, Herbolzheimer and Chaikin2018).
While spherical particles provide a useful theoretical reference, most practical suspensions consist of particles whose geometry introduces orientation-dependent drag, rotational dynamics and anisotropic disturbance fields (Ketterhagen Reference Ketterhagen2011; Cheng et al. Reference Cheng, Cao, Li and Borthwick2023; Lv & Zhao Reference Van Der, Pieter, Pedersen, Jensen, Glarborg and Dam-Johansen2025). Studies of rods, disks, spheroids, ellipsoids and cylindrical particles have shown that deviations from spherical symmetry significantly modify settling velocity, velocity fluctuations and suspension microstructure through shape-dependent hydrodynamic coupling (Hölzer & Sommerfeld Reference Hussaini and Zang2008; Zastawny et al. Reference Zastawny, Mallouppas, Zhao and van Wachem2012; Hamid et al. Reference Hamid, Molina and Yamamoto2015; Hamid et al. Reference Hamid, Arshad, Mehdi, Qasim, Ullah, Molina and Yamamoto2020; Hamid et al. Reference Hamid, Shakaib, Molina, Qasim, Fareed, Uzair and Yamamoto2024; Gustavsson et al. Reference Gustavsson, Jucha, Naso, Lévêque, Pumir and Mehlig2017; Maxey Reference Madhav and Chhabra2017). These studies have established particle shape as a key control parameter in sedimentation dynamics.
Velocity fluctuations in sedimenting suspensions arise from long-range hydrodynamic interactions and from the evolving particle microstructure. Classical analyses predicted that fluctuation amplitudes could grow with system size; however, experiments, simulations and theoretical studies have shown that this growth is limited by hydrodynamic screening and by system-dependent effects. As reviewed in Guazzelli & Hinch (Reference Guazzelli and Hinch2011), sidewalls, bottom boundaries, finite container height and finite-Reynolds-number effects can all influence fluctuation saturation in experiments. Periodic numerical simulations, including those based on the smoothed profile method (Yamamoto et al. Reference Yamamoto, Molina and Nakayama2021), eliminate these wall-induced effects and provide an idealised setting to isolate bulk many-body hydrodynamic interactions. Therefore, fluctuation statistics obtained in periodic direct numerical simulations should be compared with bounded experimental systems with appropriate caution.
Far less understood, however, is the role of particle topology, distinct from conventional shape descriptors. Most previous studies have focused on convex particles that exclude fluid from their interior volume. In contrast, non-convex particles with internal flow pathways can fundamentally alter disturbance–flow connectivity and the resulting many-body hydrodynamic interactions. Ring-shaped particles provide a canonical example of this effect. Their toroidal topology introduces a central opening that permits fluid penetration through the particle interior, thereby modifying local recirculation, wake formation, hydrodynamic screening and short-range microstructural organisation. Such ring-like geometries arise in colloidal self-assembly, polymer loop dynamics, filtration media and engineered catalytic suspensions, yet the consequences of topology-mediated fluid penetration for collective sedimentation in the non-Brownian Stokes regime remain largely unexplored.
Here, we investigate the sedimentation of ring-shaped particles using direct numerical simulations based on the smoothed profile method in the Stokes-flow regime. Each particle is modelled as a rigid assembly of six spheres forming a non-convex toroidal geometry. By systematically varying the particle volume fraction, we quantify how topology-mediated fluid penetration affects mean settling velocity, velocity fluctuations, correlation dynamics and suspension microstructure. By isolating the role of non-convex topology, the present study establishes particle topology as an independent control parameter governing hydrodynamic screening and collective sedimentation in particulate suspensions. § 2 describes the numerical methodology. § 3 outlines the computational set-up. § 4 presents the main results, including single-particle drag, collective dynamics, fluctuation statistics, correlations, microstructure and transport. § 5 provides a mechanistic discussion of topology-driven hydrodynamic interactions. § 6 summarises the principal findings and their implications.
2. Numerical methodology
The smoothed profile method (SPM) (Hamid et al. Reference Hamid, Molina and Yamamoto2013; Hamid et al. Reference Hamid, Molina and Yamamoto2014; Hamid et al. Reference Hamid, Arshad, Mehdi, Qasim, Ullah, Molina and Yamamoto2020; Yamamoto et al. Reference Yamamoto, Molina and Nakayama2021) is employed to simulate the coupled dynamics of rigid particles suspended in a viscous fluid. In this approach, the particle–fluid interface is represented by a diffuse indicator function on a fixed Eulerian grid, avoiding body-fitted or dynamically remeshed grids. Particle rigidity is imposed through a distributed body force that couples the total velocity field to the rigid-body motion of each particle. The translational and rotational dynamics are advanced using the Newton–Euler equations, with hydrodynamic forces and torques evaluated from momentum exchange between the fluid and particle domains. Short-range repulsive interactions are included to prevent particle overlap, and the coupled system is advanced using a fractional-step time integration scheme. The numerical parameters, including grid resolution, particle size relative to the mesh, domain size, simulation duration, number of snapshots, and normalisation based on isolated-particle reference values, are specified below. Previous validations of SPM for sedimenting particles of different geometries (Yamamoto et al. Reference Yamamoto, Kim and Nakayama2006; Hamid & Yamamoto Reference Hamid and Yamamoto2013b ; Hamid et al. Reference Hamid, Molina and Yamamoto2024; Fareed et al. Reference Fareed, Nadeem, Ullah and Hamid2026a , Reference Fareed, Nadeem, Ullah, Molina, Yamamoto, Chamorro and Hamidb ) provide the basis for assessing the accuracy of our simulations.
Although the SPM represents the particle–fluid boundary using a diffuse interface, this smoothing is a prescribed regularisation rather than a numerical interface-diffusion error associated with under-resolved interface-capturing methods. The accuracy of hydrodynamic force and torque transfer depends on adequate smoothing-length resolution; insufficient resolution can degrade particle–fluid coupling. In these simulations, the interface thickness is resolved over multiple grid points, and the particle size is resolved consistently with validated SPM practice. Thus the smoothing is controlled by the prescribed interface thickness and grid resolution, rather than by uncontrolled numerical diffusion.
The flow field is discretised on a fixed Eulerian grid, which enables efficient resolution of many-body hydrodynamic interactions throughout the computational domain. The SPM has been extensively validated for particulate suspensions involving various particle geometries and flow regimes, yielding accurate predictions of drag, settling velocity and collective suspension dynamics (Hamid et al. Reference Hamid, Arshad, Mehdi, Qasim, Ullah, Molina and Yamamoto2020, Reference Hamid, Shakaib, Molina, Qasim, Fareed, Uzair and Yamamoto2024). Here, each ring-shaped particle is represented as a rigid bead assembly composed of spherical elements. Rigid-body motion is enforced through a Lagrangian description, while the translational and rotational dynamics are obtained by integrating the Newton–Euler equations (Yamamoto et al. Reference Yamamoto, Kim and Nakayama2006)
\begin{equation} \boldsymbol{M}_{\!I} \dot {\boldsymbol{V}_{\!I}} = \boldsymbol{F}_{\!I}^H + \boldsymbol{F}_{\!I}^c + \boldsymbol{F}_{\!I}^{\textit{ext}}, \end{equation}
\begin{equation} \frac {\text{d}}{\text{d}t}({I_I} \boldsymbol{\varOmega}_I) = \boldsymbol{N}_{\!I}^H + \boldsymbol{N}_{\!I}^c +\boldsymbol{N}_{\!I}^{\textit{ext}}, \end{equation}
\begin{equation} \dot {\boldsymbol{Q}}_I = \text{skew}(\boldsymbol{\varOmega }_{\!I})\,\boldsymbol{Q}_{\!I}. \end{equation}
Here,
$\boldsymbol{M}_{\!I}$
,
$\boldsymbol{I}_I$
and
$\boldsymbol{Q}_{\!I}$
denote the mass, inertia tensor and orientation matrix of the
$I$
th rigid particle, respectively. The translational position and velocity are
$\boldsymbol{R}_{\!I}$
and
$\boldsymbol{V}_{\!I}$
, with
$\dot {\boldsymbol{R}}_I = \boldsymbol{V}_{\!I}$
. Particle motion is driven by the combined action of hydrodynamic forces
$\boldsymbol{F}_{\!I}^{H}$
, short-range interparticle forces
$\boldsymbol{F}_{\!I}^{c}$
, and external forcing
$\boldsymbol{F}_{\!I}^{\textit{ext}}$
due to gravity. The corresponding torques are
$\boldsymbol{N}_{\!I}^{H}$
,
$\boldsymbol{N}_{\!I}^{c}$
and
$\boldsymbol{N}_{\!I}^{\textit{ext}}$
. The indices
$I$
and
$i$
refer to the rigid particle and its constituent spherical beads, respectively.
To prevent overlap between neighbouring beads, short-range interactions are modelled using a truncated Lennard-Jones potential (Ahmed & Sadus Reference Ahmed and Sadus2009)
\begin{equation} U_{\textit{LJ}}(r_{ij}) = \begin{cases} 4\epsilon \left [\left (\sigma /r_{ij}\right )^{36} - \left (\sigma /r_{ij}\right )^{18}\right ] + \epsilon , & r_{ij} \le 2^{1/18}\sigma , \\[5pt] 0, & r_{ij} \gt 2^{1/18}\sigma ,\end{cases} \end{equation}
where
$r_{ij} = |\boldsymbol{R}_{\!i}-\boldsymbol{R}_{\!j}|$
is the separation distance between beads,
$\epsilon$
is the interaction energy scale, and
$\sigma = 2a$
defines the characteristic length scale in terms of the bead radius
$a$
. The angular velocity of particle
$I$
is denoted by
$\boldsymbol{\varOmega }_{\!I}$
, and the skew-symmetric operator in (2.3) is given by
\begin{equation} \text{skew}(\boldsymbol{\varOmega }) = \begin{pmatrix} 0 & -\varOmega _z & \varOmega _{\!y} \\[3pt] \varOmega _z & 0 & -\varOmega _x \\[3pt] -\varOmega _{\!y} & \varOmega _x & 0 \end{pmatrix}\!. \end{equation}
The fluid phase is governed by the incompressible Navier–Stokes equations, solved in an Eulerian reference frame,
\begin{equation} \partial _t \boldsymbol{u}_{\!f} + (\boldsymbol{u}_{\!f} \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_{\!f} = \rho _{\!f}^{-1}\,\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }_{\!f}, \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_{\!f} = 0 , \end{equation}
where
$\boldsymbol{u}_{\!f}$
is the fluid velocity field, and
$\rho _{\!f}$
is the fluid density. The fluid stress tensor is expressed as
\begin{equation} \boldsymbol{\sigma }_{\!f} = - p\boldsymbol{I} + \eta _{\!f}\left [\boldsymbol{\nabla }\boldsymbol{u}_{\!f} + (\boldsymbol{\nabla }\boldsymbol{u}_{\!f})^{\rm T}\right ]\!, \end{equation}
where
$p$
denotes the hydrodynamic pressure,
$\eta _{\!f}$
is the dynamic viscosity, and
$\boldsymbol{I}$
is the identity tensor. Within the SPM framework, the fluid and particle phases are coupled through a single total velocity field
$\boldsymbol{u}$
that smoothly blends the fluid velocity
$\boldsymbol{u}_{\!f}$
and the rigid-body particle velocity
$\boldsymbol{u}_{\!p}$
,
\begin{equation} \boldsymbol{u} = (1-\varphi )\boldsymbol{u}_{\!f} + \varphi \boldsymbol{u}_{\!p} , \end{equation}
where
$\varphi (\boldsymbol{x},t)$
is a diffuse indicator function that distinguishes the fluid and particle domains. The rigid-particle velocity field is constructed as
\begin{equation} \varphi \boldsymbol{u}_{\!p} = \sum _{I} \varphi _I \left [ \boldsymbol{V}_{\!I} + \boldsymbol{\varOmega }_{\!I} \times \boldsymbol{r}_{\!I} \right ]\!, \end{equation}
where
$\boldsymbol{V}_{\!I}$
and
$\boldsymbol{\varOmega }_{\!I}$
are the translational and angular velocities of the
$I$
th rigid particle, and
$\boldsymbol{r}_{\!I}$
is the position vector relative to its centre of mass. The indicator function satisfies
\begin{equation} 0 \le \varphi (\boldsymbol{x},t) \le 1, \end{equation}
with
$\varphi =1$
inside the particle domain, and
$\varphi =0$
in the surrounding fluid. The smooth transition across the particle–fluid interface enables both phases to be represented consistently on a single Eulerian grid while preserving the rigid-body kinematics of the particle phase.
Momentum exchange between the fluid and particle phases, as well as the enforcement of rigid-body motion within the particle domain, are achieved through the body-force density
$\boldsymbol{f}_{\!p}$
. The term
$\varphi \boldsymbol{f}_{\!p}$
represents the force applied to impose the rigidity constraint and ensure consistency between the fluid and particle velocities. The governing equation for the total velocity field is therefore given by
\begin{equation} \partial _t \boldsymbol{u} + (\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u} = \rho ^{-1}\, \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma } + \varphi \boldsymbol{f}_{\!p}. \end{equation}
The coupled fluid–particle system is advanced in time using a fractional-step scheme adapted to account for the rotational dynamics of non-spherical rigid bodies. Let
$\boldsymbol{u}^n$
denote the velocity field at time
$t^n = nh$
, where
$h$
is the time step. In the first substep, the advective and viscous contributions are evaluated to obtain an intermediate velocity field:
\begin{equation} \boldsymbol{u}^* = \boldsymbol{u}^n + \int _{t_n}^{t_n + h} \mathrm{d}s \, \boldsymbol{\nabla }\boldsymbol{\cdot }\big [ \rho ^{-1}\,\boldsymbol{\sigma }(s) - \boldsymbol{u}(s)\,\boldsymbol{u}(s) \big ]. \end{equation}
The particle positions and orientations are then updated consistently with the rigid-body kinematics:
\begin{equation} \boldsymbol{R}_{\!I}^{n+1} = \boldsymbol{R}_{\!I}^{n} + \int _{t_n}^{t_n + h} \mathrm{d}s \, \boldsymbol{V}_{\!I}(s), \end{equation}
\begin{equation} \boldsymbol{Q}_{\!I}^{n+1} = \boldsymbol{Q}_{\!I}^{n} + \int _{t_n}^{t_n + h} \mathrm{d}s \, \text{skew}(\boldsymbol{\varOmega }_{\!I}(s)) \, \boldsymbol{Q}_{\!I}(s). \end{equation}
Since the profile function
$\varphi$
depends explicitly on particle position and orientation, it is updated consistently using (2.13) and (2.14). The corresponding intermediate rigid-particle velocity field is then constructed as
\begin{equation} \varphi \boldsymbol{u}_{\!p}^{*} = \sum _{I} \varphi _I^{n+1} \left [ \boldsymbol{V}_{\!I}^{n} + \boldsymbol{\varOmega }_{\!I}^{n} \times \boldsymbol{r}_{\!I}^{n+1} \right ]\!. \end{equation}
To ensure consistency between the momentum exchanged within the diffuse particle domain and the hydrodynamic loading acting on each rigid body, the hydrodynamic force and torque are obtained directly from momentum conservation. The time-integrated hydrodynamic force is
\begin{equation} \int _{t_n}^{t_n + h} \mathrm{d}s \, \boldsymbol{F}_{\!I}^{H} = \int \mathrm{d}\boldsymbol{x} \, \rho \, \varphi _I^{n+1} \big ( \boldsymbol{u}^{*} - \boldsymbol{u}_{\!p}^{*} \big ), \end{equation}
and the corresponding hydrodynamic torque is
\begin{equation} \int _{t_n}^{t_n + h} \mathrm{d}s \, \boldsymbol{N}_{\!I}^{H} = \int \mathrm{d}\boldsymbol{x} \, \boldsymbol{r}_{\!I}^{n+1} \times \rho \, \varphi _I^{n+1} \big ( \boldsymbol{u}^{*} - \boldsymbol{u}_{\!p}^{*} \big ). \end{equation}
The resulting hydrodynamic force and torque are then used to advance the translational and rotational velocities of each rigid particle according to
\begin{equation} \boldsymbol{V}_{\!I}^{n+1} = \boldsymbol{V}_{\!I}^{n} + \boldsymbol{M}_{\!I}^{-1} \!\int _{t_n}^{t_n+h}\!\! \left ( \boldsymbol{F}_{\!I}^{H} + \boldsymbol{F}_{\!I}^{c} + \boldsymbol{F}_{\!I}^{\textit{ext}} \right ) {\rm d}s, \end{equation}
\begin{equation} \boldsymbol{\varOmega }_{\!I}^{n+1} = \boldsymbol{\varOmega }_{\!I}^{n} + \boldsymbol{I}_I^{-1} \!\int _{t_n}^{t_n+h}\!\! \left ( \boldsymbol{N}_{\!I}^{H} + \boldsymbol{N}_{\!I}^{c} + \boldsymbol{N}_{\!I}^{\textit{ext}} \right ) {\rm d}s. \end{equation}
Following the rigid-body velocity update, particle rigidity is enforced through the body-force term
$\varphi \boldsymbol{f}_{\!p}$
, yielding the corrected velocity field
\begin{equation} \boldsymbol{u}^{n+1} = \boldsymbol{u}^{*} + \int _{t_n}^{t_n + h} {\rm d}s \, \varphi \boldsymbol{f}_{\!p}, \end{equation}
with
\begin{equation} \int _{t_n}^{t_n + h} \mathrm{d}s \, \varphi (s)\, \boldsymbol{f}_{\!p}(s) = \varphi ^{n+1} \big( \boldsymbol{u}_{\!p}^{n+1} - \boldsymbol{u}^{*} \big ) - \frac {h}{\rho } \boldsymbol{\nabla }p_p. \end{equation}
The pressure field
$p_p$
associated with the rigidity constraint is obtained by imposing incompressibility on the corrected velocity field,
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^{n+1}=0$
. Within the smoothed profile framework, additional kinematic constraints are imposed to preserve rigid-body motion of the bead assembly (Uzair et al. Reference Tee, Mucha, Cipelletti, Manley, Brenner, Segre and Weitz2023):
\begin{equation} \boldsymbol{V}_{\!i} = \boldsymbol{V}_{\!I} + \boldsymbol{\varOmega }_{\!I} \times \boldsymbol{G}_{\!i},\quad \boldsymbol{\varOmega }_i = \boldsymbol{\varOmega }_{\!I}, \end{equation}
\begin{equation} \tilde {\boldsymbol{G}}_i = \boldsymbol{Q}_{\!I}^{\rm T} \boldsymbol{G}_{\!i} = \text{constant}. \end{equation}
Here,
$\boldsymbol{G}_{\!i} = \boldsymbol{R}_{\!i} - \boldsymbol{R}_{\!I}$
is the position vector of the
$i$
th spherical bead relative to the centre of mass of the rigid particle. When expressed in the particle reference frame as
$\tilde {\boldsymbol{G}}_i$
, these relative vectors remain invariant in time. The instantaneous bead locations are therefore reconstructed from
$\boldsymbol{R}_{\!i} = \boldsymbol{R}_{\!I} + \boldsymbol{G}_{\!i}$
, with
$\boldsymbol{G}_{\!i} = \boldsymbol{Q}_{\!I} \tilde {\boldsymbol{G}}_i$
.
Computational domain and ring particle geometry. Simulations are performed in a domain of size
$2^7\varDelta \times 2^7\varDelta \times 2^8\varDelta$
, with gravity
$\boldsymbol{g}$
acting in the
$-z$
direction. The particle consists of six spherical beads forming a ring. The ring orientation is defined by the angle
$\theta$
between the normal vector
$\boldsymbol{n}$
and the vertical
$z$
axis.
Figure 1. Long description
Panel A: A three-dimensional rectangular computational domain with dimensions labeled as 2^7Δ in the x and y directions and 2^8Δ in the z direction. Gravity acts downward in the z direction. A particle composed of six spherical beads arranged in a ring is positioned within the domain. Panel B: A close-up view of the ring particle geometry. The ring is composed of six spherical beads, with the diameter of the ring labeled as D = 16Δ and the width of each bead labeled as w = 6Δ. The orientation of the ring is defined by the angle θ between the normal vector n and the vertical axis z.
3. Numerical set-up
We place our focus on the collective settling dynamics of randomly oriented, initially uniformly distributed ring-shaped particles in a Newtonian fluid under gravity. Each particle is constructed as a rigid assembly of six spheres of radius
$3\varDelta$
, where
$\varDelta$
denotes the lattice spacing. Throughout the simulations, the lattice spacing is taken as the unit grid length
$\varDelta = 1$
in lattice units (figure 1). This bead-based geometry differs from an ideal smooth torus with a circular tube because it retains small-scale surface discreteness near the bead contacts. Nevertheless, for the present low-
$\textit{Re}$
sedimentation problem, the principal features of interest are associated with the non-convex topology, namely the central opening, fluid penetration through the ring, and the resulting disturbance field. A smooth torus may therefore modify near-surface flow details, but the topology-driven mechanisms discussed here are expected to remain qualitatively similar. At
$\textit{Re}=0.1$
, strong inertial wake-memory effects are not expected; instead, the flow through the ring adjusts approximately quasi-steadily to the instantaneous particle motion. The resulting toroidal geometry has diameter
$D = 16\varDelta$
and width
$w = 6\varDelta$
, corresponding to aspect ratio
$D/w = 2.6$
. The simulations follow the standard SPM formulation of Yamamoto et al. (Reference Yamamoto, Molina and Nakayama2021), with a diffuse-interface smoothing length
$\xi =2\varDelta$
. The time step is selected to satisfy the stability requirements associated with viscous diffusion, advection, and particle–fluid coupling. Owing to the diffuse-interface formulation, the particle boundary is represented as a smoothed, or soft, toroidal geometry rather than an ideal sharp-interface torus. The interface thickness remains small compared with the particle size, and the chosen particle resolution,
$D/\varDelta=16$
, is consistent with previous validated SPM studies (Yamamoto et al. Reference Yamamoto, Kim and Nakayama2006; Hamid & Yamamoto Reference Hamid and Yamamoto2013a
; Hamid et al. Reference Hamid, Arshad, Mehdi, Qasim, Ullah, Molina and Yamamoto2020, Reference Hamid, Shakaib, Molina, Qasim, Fareed, Uzair and Yamamoto2024; Uzair et al. Reference Tee, Mucha, Cipelletti, Manley, Brenner, Segre and Weitz2023), where comparable resolutions were shown to yield reliable predictions of hydrodynamic forces, torques and sedimentation characteristics.
The six-bead ring is used as a minimal rigid assembly that preserves the essential non-convex topology and internal opening while maintaining computational feasibility for many-particle simulations. This bead-based model is intended as a controlled geometry for isolating topology effects, rather than as a direct representation of naturally occurring sedimenting particles, which are often irregular, rough and polydisperse. The main mechanisms examined here, including fluid penetration through the internal opening, modification of the disturbance field, and changes in hydrodynamic screening, are expected to remain relevant beyond the specific six-bead configuration. The model is particularly applicable to engineered suspensions containing bead-like or ring-shaped inclusions, whereas extensions to realistic natural particles would require additional geometric complexity.
The suspension volume fraction
$\phi$
is varied over the range
$10^{-3} \le \phi \le 10^{-1}$
, corresponding to particle counts between 10 and 617. This representation preserves rigid-body kinematics while resolving near-field hydrodynamic interactions at the particle scale.
The ring geometry represents a minimal non-convex configuration with an internal flow pathway, enabling direct examination of how particle topology influences drag, rotational dynamics, and collective sedimentation behaviour. The simulations remained stable over all investigated time intervals, with no observed numerical drift or unbounded growth of particle velocities. Although trajectory-level differences can accumulate over long times in many-body systems, the ensemble- and time-averaged quantities – including mean settling velocity, velocity fluctuations, correlation functions and diffusion coefficients – remained statistically stationary over the sampling intervals. These checks indicate that round-off errors do not affect the conclusions drawn from the statistics. The computational domain has dimensions
$L_X \times L_Y \times L_Z = 2^{7}\varDelta \times 2^{7}\varDelta \times 2^{8}\varDelta$
, corresponding to aspect ratio
$1:1:2$
, with periodic boundary conditions imposed in all directions. The elongated vertical extent ensures adequate representation of particle trajectories and wake interactions while mitigating finite-size effects, consistent with prior studies (Kobayashi et al. Reference Kobayashi, Jung, Matsuoka, Nakayama, Molina and Yamamoto2023; Uzair et al. Reference Tee, Mucha, Cipelletti, Manley, Brenner, Segre and Weitz2023). The particle diameter spans
$D/\varDelta = 16$
grid points, providing sufficient resolution of particle geometry and near-field hydrodynamics. The particle Reynolds number is defined as
$\textit{Re} = \rho V_t D / \eta$
, where
$V_t$
is the terminal settling velocity of an isolated particle,
$\rho$
is the fluid density, and
$\eta$
is the dynamic viscosity. The reference velocity
$V_t$
is obtained from single-particle simulations under identical conditions. All simulations are performed at
$\textit{Re}=0.1$
to ensure operation within the Stokes regime, where viscous forces dominate and inertial effects remain weak. The simulations are not performed in the strict zero-
$\textit{Re}$
limit, and both fluid and particle inertia can influence the orientational dynamics. Although
$\textit{Re}=0.1$
remains predominantly viscous, weak fluid inertia introduces Oseen-type fore–aft asymmetry in the disturbance flow. The Newton–Euler equations are solved with
$\rho _p/\rho _{\!f}=5$
, so translational and rotational particle inertia are explicitly retained. Thus unlike the ideal Stokes limit, in which an isolated falling toroidal particle would not select a preferred orientation, the present weakly inertial regime can permit irreversible orientational relaxation. These inertial effects are distinct from the additional orientational dynamics produced by many-body hydrodynamic interactions in concentrated suspensions.
Composite particle geometries formed by assemblies of spherical beads: (a) disk, (b) cube, (c) rod, (d) triangular, and (e) ring. The ring configuration is the primary case, while the remaining geometries serve as references for comparison of shape-dependent sedimentation dynamics.
Figure 2 illustrates the non-spherical particle geometries constructed from spherical bead assemblies used as reference cases. The configurations include disk-shaped (Hamid et al. Reference Hamid, Shakaib, Molina, Qasim, Fareed, Uzair and Yamamoto2024), cube-shaped (Fareed et al. Reference Fareed, Nadeem, Ullah, Molina, Yamamoto, Chamorro and Hamid2026b
), rod-shaped (Hamid et al. Reference Hamid, Arshad, Mehdi, Qasim, Ullah, Molina and Yamamoto2020), triangular-shaped (Fareed et al. Reference Fareed, Nadeem, Ullah and Hamid2026a
) and ring-shaped particles. The particle-to-fluid density ratio is fixed at
$\rho _p/\rho _{\!f} = 5$
, with gravity acting in the negative
$z$
direction. The chosen density ratio is not intended to represent a specific water–sediment system, but is used as a controlled numerical parameter to maintain stable sedimentation at fixed particle Reynolds number
$\textit{Re}=0.1$
. Under these conditions, the dynamics remains predominantly governed by viscous hydrodynamic interactions, with only weak inertial contributions. Comparable density ratios have been used in previous SPM studies of sedimenting suspensions (Hamid et al. Reference Hamid, Arshad, Mehdi, Qasim, Ullah, Molina and Yamamoto2020, Reference Hamid, Shakaib, Molina, Qasim, Fareed, Uzair and Yamamoto2024; Fareed et al. Reference Fareed, Nadeem, Ullah and Hamid2026a
,Reference Fareed, Nadeem, Ullah, Molina, Yamamoto, Chamorro and Hamid
b
). The density ratio affects particle inertia, and the present results are therefore interpreted within the low-Reynolds-number regime specified above.
Particle–particle interactions are mediated by hydrodynamic coupling and short-range repulsive forces, and all simulations are evolved until statistically stationary states are reached. Periodic-domain simulations use a Fourier pseudo-spectral formulation, with spatial derivatives and pressure projection evaluated in Fourier space (Hussaini & Zang Reference Hölzer and Sommerfeld1986; Nakayama & Yamamoto Reference Mucha, Tee, Weitz, Shraiman and Brenner2005). A uniform vertical pressure gradient is imposed to balance the net gravitational forcing and ensure that the total fluid–particle volumetric velocity remains zero (Ladd Reference Ladd1994). The zero-wavenumber mode is explicitly treated to prevent spurious mean-flow drift, consistent with standard spectral formulations of incompressible periodic flows (Hussaini & Zang Reference Hölzer and Sommerfeld1986; Canuto et al. Reference Canuto, Quarteroni, Hussaini, Zang and Thomas2007). All simulations are performed on a uniform grid containing approximately
$4\times 10^{6}$
lattice points. Each case is integrated over
$1000t_t{-}1500t_t$
, where the characteristic settling time scale is
$t_t = D/V_t$
. Statistical quantities are computed from approximately
$5\times 10^{4}$
snapshots collected over at least
$2.5\times 10^{7}$
time steps. Numerical stability was monitored throughout the simulations. No spurious mean-flow drift, unbounded particle velocities or loss of statistical stationarity were observed over the sampling intervals. Although trajectory-level differences may accumulate in long many-body simulations due to round-off error, the reported ensemble- and time-averaged quantities remained statistically stationary.
Brownian motion is neglected, as sedimentation at
$\textit{Re} = 0.1$
corresponds to an effectively infinite Péclet number. The particles are therefore treated as non-Brownian, with gravitational forces dominating over thermal fluctuations, consistent with prior studies on non-Brownian suspensions (Hamid & Yamamoto Reference Hamid and Yamamoto2013a
,
Reference Hamid and Yamamotob
). The domain size and aspect ratio are selected to minimise finite-size effects while maintaining computational efficiency (Qiu et al. Reference Padding and Louis2015; Shakeel et al. Reference Segre, Herbolzheimer and Chaikin2018). At higher volume fractions, many-body interactions reduce artefacts associated with periodic boundary conditions. At very low volume fractions, interactions with periodic images may still influence the dynamics (Zastawny et al. Reference Zastawny, Mallouppas, Zhao and van Wachem2012; Uzair et al. Reference Tee, Mucha, Cipelletti, Manley, Brenner, Segre and Weitz2023). To mitigate such effects, all quantities are normalised using single-particle reference values obtained under identical conditions. A summary of the numerical and physical parameters is provided in table 1.
Numerical and physical parameters used in the simulations.
Table 1. Long description
A table with three columns and twelve rows detailing parameters, symbols, and values used in simulations. The columns are labeled Parameter, Symbol, and Value. The table includes parameters such as Domain size, Total grid points, Total time steps, Number of snapshots, Density ratio, Reynolds number, Volume fraction, Number of particles, Ring diameter, Ring width, Lattice spacing, Ring aspect ratio, and Domain aspect ratio. Each row lists the corresponding symbol and value for each parameter.
4. Results
We first examine the hydrodynamic behaviour of an isolated ring-shaped particle. The drag coefficient is defined as
\begin{equation} C_D = \frac {2F_D}{\rho _{\!f} V_t^2 A_{\!p}}, \end{equation}
where
$F_D = Q_p g (\rho _p - \rho _{\!f})$
is the hydrodynamic drag force,
$Q_{{\kern-1pt}p}$
is the particle volume, and
$A_p$
is the projected area normal to the settling direction.
All parameters, except particle geometry, are held fixed to isolate topology-induced effects. Accurate prediction of
$C_D$
for non-spherical particles remains challenging due to its sensitivity to particle geometry, orientation and local flow structure (Sanjeevi et al. Reference Richardson and Zaki2018; Vergara et al. Reference Vergara, Wei and Fuentes2024). Existing correlations (Ganser Reference Ganser1993; Kurose & Komori Reference Kurose and Komori1999; Hölzer & Sommerfeld Reference Hussaini and Zang2008) are largely developed for convex particles and therefore do not explicitly account for internal flow pathways associated with non-convex geometries.
Drag coefficient of an isolated ring-shaped particle as a function of Reynolds number. Symbols denote present simulations; lines show reported drag correlations (Levenspiel & Haider Reference Levenspiel and Haider1989; Chien Reference Chien1994; Madhav & Chhabra Reference Lv and Zhao1995; Ceylan et al. Reference Ceylan, Altunbaş and Kelbaliyev2001; Xu et al. Reference Xu, Song, Li, Pang and Zhu2019).
Normalised average sedimentation velocity
$V_{\textit{sed}}/V_t$
as a function of particle volume fraction for ring-shaped particles and reference geometries: rods (Hamid et al. Reference Hamid, Arshad, Mehdi, Qasim, Ullah, Molina and Yamamoto2020), spheres (Hamid & Yamamoto Reference Hamid and Yamamoto2013c
), disks (Hamid et al. Reference Hamid, Shakaib, Molina, Qasim, Fareed, Uzair and Yamamoto2024), cubes (Fareed et al. Reference Fareed, Nadeem, Ullah, Molina, Yamamoto, Chamorro and Hamid2026b
) and triangular particles (Fareed et al. Reference Fareed, Nadeem, Ullah and Hamid2026a
). The shaded region denotes a narrow envelope about the
$(1-\phi )^7$
scaling.
Figure 3 shows the variation of
$C_D$
with Reynolds number
$10^{-2} \le \textit{Re} \le 1$
for the ring-shaped particle. In the Stokes limit (
$\textit{Re} \ll 1$
), the drag follows the expected scaling
$C_D \sim \textit{Re}^{-1}$
(Stokes Reference Shakeel, Hamid, Ullah, Molina and Yamamoto2007). However, the magnitude of
$C_D$
remains systematically higher than that predicted for convex particles. This trend is consistent with additional viscous dissipation associated with flow penetration through the central opening, although the present drag data do not uniquely isolate internal recirculation from other geometric effects. As
$\textit{Re}$
increases, deviations from Stokes scaling emerge due to weak inertial effects, leading to asymmetry in the disturbance field and gradual flattening of the
$C_D$
–
$\textit{Re}$
relation. Comparison with empirical correlations shows that the Levenspiel & Haider (Reference Levenspiel and Haider1989), Madhav & Chhabra (Reference Lv and Zhao1995) and Chien (Reference Chien1994) models provide closer agreement with our results, whereas the Ceylan et al. (Reference Ceylan, Altunbaş and Kelbaliyev2001) and Xu et al. (Reference Xu, Song, Li, Pang and Zhu2019) correlations underpredict the drag, particularly in the viscous regime.
We next consider the collective settling dynamics that emerges from hydrodynamic interactions among many ring-shaped particles. The mean settling velocity decreases with increasing concentration due to disturbance–flow interactions and the associated fluid backflow. This behaviour is commonly described by the Richardson–Zaki relation (Richardson & Zaki Reference Qiu, Lin, Zhou, Wang, Yuan, Wei, Dai, Luo and Chen1954)
\begin{equation} V_{\textit{sed}}/V_t = (1 - \phi )^n, \end{equation}
where
$V_{\textit{sed}}=\langle V_{\textit{Iz}}\rangle$
is the ensemble-averaged settling velocity. Figure 4 shows
$V_{\textit{sed}}/V_t$
as a function of
$\phi$
for the particle geometries considered. Here,
$V_t$
denotes the long-time mean settling velocity of an isolated particle after the initial transient has decayed, providing the reference velocity for the suspension results; the orientation dependence of the ring settling dynamics is considered separately through its evolving orientation distribution. At low volume fractions, isolated particles tend to orient horizontally, while increasing volume fraction promotes more vertical orientations and enhanced settling. However, deviations of
$V_{\textit{sed}}/V_t$
from unity may also be sensitive to finite-domain effects under periodic boundary conditions, particularly in the dilute regime where the number of particles is small and the disturbance flow is long-ranged. At higher volume fractions, particle interactions increasingly hinder settling and reduce the velocity. This sensitivity is consistent with the classical result that Stokes sedimentation in periodic systems depends on the periodic cell size (Hasimoto Reference Hasimoto1959). Since
$V_t$
is obtained from an isolated-particle simulation under the same periodic-domain conditions, the leading single-particle finite-domain effect is included in the normalisation. Residual deviations at low
$\phi$
should therefore be interpreted cautiously, as they may include finite-domain sensitivity of the dilute many-particle configuration in addition to direct particle–particle interactions. For ring-shaped particles, the settling velocity exhibits a pronounced reduction with concentration for
$\phi \gtrsim 0.005$
, and follows a Richardson–Zaki-type scaling with exponent
$n \approx 7$
over an intermediate range of
$\phi$
. The Richardson–Zaki exponent
$n \approx 7$
is used here only as an intermediate-concentration reference trend, and is not intended to describe the full range of
$\phi$
considered. At very low concentrations (
$\phi \lt 0.01$
), the statistics are more sensitive to finite-size effects and the limited number of particles, whereas at higher concentrations (
$\phi \gt 0.08$
), near-field constraints, geometric crowding and many-body hydrodynamic interactions lead to departures from a simple power-law hindered-settling description. The local deviation around
$\phi \approx 0.07{-}0.08$
is therefore not overinterpreted as a separate scaling regime. A statistical-convergence check for the mean settling velocity is provided in Appendix A, supporting the stationarity of the reported values over the sampling interval.
Particle motion in sedimenting suspensions exhibits significant velocity fluctuations driven by long-range hydrodynamic interactions. These fluctuations arise from spatial and temporal variations in the disturbance flow induced by neighbouring particles, and depend strongly on the volume fraction
$\phi$
. The instantaneous fluctuation of the vertical velocity of particle
$i$
is defined as
\begin{equation} \delta V_{\textit{Iz}}(t) = V_{\textit{Iz}}(t) - V_{\textit{sed}}. \end{equation}
Following the onset of cluster settling, the mean magnitude of these fluctuations across all
$N$
particles is evaluated as
\begin{equation} \varDelta V_z(t) = \frac {1}{N}\sum _{I=1}^{N} |\delta V_{\textit{Iz}}(t)|. \end{equation}
The temporal variability is then quantified by the standard deviation
\begin{equation} \sigma _{\varDelta V_z} = \left [ \frac {1}{T}\sum _{k=1}^{T} \left ( \varDelta V_z(t_k) - \langle \varDelta V_z \rangle \right )^2 \right ]^{1/2}, \end{equation}
where
$T$
is the number of temporal samples, and
$\langle \Delta V_z \rangle$
denotes the time-averaged value of
$\Delta V_z$
. The horizontal fluctuations
$\sigma _{\Delta V_x}$
are defined analogously.
(a) Horizontal,
$\sigma _{\Delta V_x}/V_t$
, and (b) vertical,
$\sigma _{\varDelta V_z}/V_t$
, velocity fluctuations as functions of particle volume fraction
$\phi$
for different particle geometries. Solid lines indicate reference scalings
$\phi ^{1/2}$
and
$\phi ^{1/3}$
in the dilute regime. The shaded region denotes a narrow envelope about the scalings.
Figure 5. Long description
Two scatter plots illustrate the relationship between velocity fluctuations and particle volume fraction for various particle geometries. Panel A: The scatter plot shows horizontal velocity fluctuations (σΔVx/Vt) on the vertical axis and particle volume fraction (ϕ) on the horizontal axis. Different symbols represent different particle geometries, including disk, rod, sphere, cube, triangle, and ring. Solid lines indicate reference scalings, and the shaded region denotes a narrow envelope about the scalings. Panel B: The scatter plot displays vertical velocity fluctuations (σΔVz/Vt) on the vertical axis and particle volume fraction (ϕ) on the horizontal axis. Similar to Panel A, different symbols represent different particle geometries, and solid lines indicate reference scalings with a shaded region around them.
Figure 5 shows the normalised horizontal,
$\sigma _{\Delta V_x}/V_t$
, and vertical,
$\sigma _{\Delta V_z}/V_t$
, velocity fluctuations as functions of the particle volume fraction
$\phi$
. In the dilute regime (
$\phi \lt 0.02$
), the fluctuations increase with
$\phi$
and exhibit pronounced anisotropy, with vertical fluctuations approximately
$1.5$
–
$2$
times larger than the horizontal component. This anisotropy reflects the preferential alignment of disturbance flows and wake interactions along the direction of gravity.
The vertical velocity fluctuations follow an approximate
$\phi ^{1/3}$
scaling, consistent with hydrodynamic screening models of sedimentation (Levine et al. Reference Levine, Ramaswamy, Frey and Bruinsma1998; Mucha et al. Reference Maxey2004; Guazzelli & Hinch Reference Guazzelli and Hinch2011), which predict a finite correlation length governing velocity correlations. The observed scaling indicates that hydrodynamic screening in ring suspensions is governed by anisotropic disturbance propagation, arising from topology-induced flow penetration and enhanced wake coupling. For orientations near
$\theta \approx 90^\circ$
, the central opening provides an additional pathway for fluid transport, which may contribute to the larger horizontal velocity fluctuations of the ring relative to the disk. The comparable vertical velocity fluctuations indicate that this internal flow pathway has a weaker influence on vertical particle motion.
As
$\phi$
increases, many-body interactions progressively modify this scaling, with the vertical component exhibiting stronger growth due to persistent gravity-aligned correlations. Related flow configurations arise in systems involving topological or toroidal colloids, porous filtration media, and engineered catalytic or active toroidal particles, where internal openings or connected pore spaces can modify drag, disturbance propagation, hydrodynamic screening and collective transport. These examples provide a broader context for the topology-driven mechanisms examined here. The magnitude and anisotropy of the fluctuations also depend strongly on particle geometry. Nearly isotropic particles, such as spheres, exhibit weaker directional dependence, and typically follow the classical
$\phi ^{1/2}$
scaling associated with isotropic microstructures. In contrast, anisotropic or non-convex particles, including rings, rods, disks and triangular particles, exhibit stronger anisotropy and follow
$\phi ^{1/3}$
scaling in the dilute regime (Hamid et al. Reference Hamid, Arshad, Mehdi, Qasim, Ullah, Molina and Yamamoto2020, Reference Hamid, Shakaib, Molina, Qasim, Fareed, Uzair and Yamamoto2024). Similar behaviour has been reported in sedimenting fibre suspensions, where anisotropic hydrodynamic interactions and particle orientation govern both fluctuation magnitude and spatial correlations (Chiarini et al. Reference Chiarini, Gallorini and Rosti2026).
The PDFs of the orientation angle
$\theta$
for ring-shaped particles using
$5^\circ$
angular bins: (a)
$\phi =0.001$
, comparing
$0\le t\le T_0/2$
and
$T_0/2\le t\le T_0$
, with
$T_0=5\times 10^4$
; (b)
$\phi =0.05$
for the same temporal intervals. Here,
$\theta$
is measured relative to the reference direction.
Figure 6 shows the probability density function (PDF) of the orientation angle
$\theta$
. Figure 6(a) compares the PDFs at
$\phi =0.001$
over two temporal windows,
$0\le t\le T_0/2$
and
$T_0/2\le t\le T_0$
. Residual differences between the two intervals, most pronounced over
$20^\circ \lesssim \theta \lesssim 80^\circ$
, indicate that the orientational statistics have not yet reached a stationary state in the dilute regime. At the higher concentration
$\phi =0.05$
(figure 6
b), the two distributions collapse, indicating convergence to a statistically stationary orientational state. Similar behaviour at
$\phi =0.03$
and
$0.1$
confirms that increasing concentration suppresses temporal variability and accelerates the relaxation of orientation statistics through collective hydrodynamic interactions. Across all concentrations, the distribution remains biased towards
$\theta \approx 90^\circ$
. Figure 7 illustrates the instantaneous flow field around the ring-shaped particle at two orientations, showing fluid motion through the central opening, and the formation of localised recirculating structures whose spatial organisation varies with particle orientation.
Instantaneous velocity field around an isolated settling ring-shaped particle at orientations (a)
$\theta =0^\circ$
and (b)
$\theta =90^\circ$
. Colour contours show the vertical velocity component, and arrows denote the in-plane velocity vectors. The circle marks the projected particle boundary. The flow fields illustrate orientation-dependent disturbance structures and fluid motion through the central opening of the ring.
Figure 7. Long description
Panel A: A heat map displays the instantaneous velocity field around an isolated settling ring-shaped particle at a specific orientation. The vertical axis represents the vertical velocity component (Vz), and the horizontal axis represents the horizontal velocity component (Vx). Color contours indicate the magnitude of the vertical velocity component, with red representing higher values and blue representing lower values. Black arrows denote the in-plane velocity vectors. A circle marks the projected particle boundary. The velocity field shows complex patterns and disturbances around the particle. Panel B: Another heat map illustrates the instantaneous velocity field around the same ring-shaped particle but at a different orientation. The vertical axis represents the vertical velocity component (Vz), and the horizontal axis represents the horizontal velocity component (Vx). Color contours indicate the magnitude of the vertical velocity component, with red representing higher values and blue representing lower values. Black arrows denote the in-plane velocity vectors. A circle marks the projected particle boundary. The velocity field shows different patterns and disturbances compared to Panel A, highlighting orientation-dependent fluid motion through the central opening of the ring.
To further characterise the collective dynamics, we examine the spatial structure of velocity fluctuations through particle velocity correlations, which quantify the range of hydrodynamic coupling. The correlation function of vertical velocity fluctuations is defined as
\begin{equation} C_z(r)= \left \langle \sum _{I\lt J} \delta V_{\textit{Iz}}\,\delta V_{\textit{Jz}} \,\delta (r-r_{\textit{IJ}}) \right \rangle , \end{equation}
where
$r_{\textit{IJ}}$
is the particle separation distance, and
$\delta V_{\textit{Iz}}=V_{\textit{Iz}}-V_{\textit{sed}}$
denotes the fluctuation relative to the mean settling velocity. Figures 8(a) and 8(b) show the spatial correlations of velocity fluctuations in the horizontal and vertical directions. The correlations are strongly anisotropic. The horizontal component
$C_z(x)$
decays rapidly and becomes negative beyond
$x/D \approx 3{-}4$
, indicating weak anti-correlated motion induced by neighbouring wakes. In contrast, the vertical component
$C_z(z)$
remains positive over larger separations, reflecting persistent streamwise coherence of disturbance flows. Increasing
$\phi$
primarily enhances the correlation magnitude, while the characteristic decay length remains of order a few particle diameters, consistent with hydrodynamic screening. Temporal correlations provide complementary information on the persistence of velocity fluctuations. The autocorrelation functions are defined as
$C_x(t)=\langle V_{Ix}(t)\,V_{Ix}(0)\rangle$
and
$C_z(t)=\langle \delta V_{\textit{Iz}}(t)\,\delta V_{\textit{Iz}}(0)\rangle$
. The normalised autocorrelations (figures 8
c,d) show that correlation times decrease with increasing
$\phi$
, indicating accelerated decorrelation due to many-body interactions. Dilute suspensions exhibit longer-lived correlations, particularly in the vertical direction, reflecting sustained propagation of gravity-aligned disturbances. Across all concentrations, vertical correlations remain more persistent than horizontal ones, consistent with the anisotropic structure of the sedimentation-induced flow field.
(a) Horizontal and (b) vertical spatial correlations; (c) horizontal and (d) vertical temporal autocorrelation functions.
As a complementary measure, the correlation time is obtained from the temporal autocorrelation function as
\begin{equation} \tau _\beta =\frac {1}{C_\beta (0)}\int _0^\infty C_\beta (t)\,\text{d}t, \end{equation}
where
$\beta$
denotes the horizontal or vertical component.
(a) Horizontal and (b) vertical correlation times. Included for comparison are rod-like (Hamid et al. Reference Hamid, Arshad, Mehdi, Qasim, Ullah, Molina and Yamamoto2020), spherical (Hamid & Yamamoto Reference Hamid and Yamamoto2013c
), disk-shaped (Hamid et al. Reference Hamid, Shakaib, Molina, Qasim, Fareed, Uzair and Yamamoto2024), cube-shaped (Fareed et al. Reference Fareed, Nadeem, Ullah, Molina, Yamamoto, Chamorro and Hamid2026b
) and triangular (Fareed et al. Reference Fareed, Nadeem, Ullah and Hamid2026a
) particles. The shaded region denotes a narrow envelope about the
$\phi ^{-1/2}$
scaling.
Figure 9. Long description
Two scatter plots depict correlation times for various particle shapes. Panel A shows horizontal correlation times, while Panel B shows vertical correlation times. The x-axis in both panels represents the variable phi, ranging from 10^-3 to 10^-1. The y-axis in both panels represents the ratio of correlation time to a reference time, ranging from 10^-1 to 10^2. Different symbols and colors represent different particle shapes: disks, rods, spheres, cubes, triangles, and rings. A shaded region denotes a narrow envelope about the scaling. The data points show how correlation times vary with phi for each particle shape, highlighting differences in horizontal and vertical components due to particle geometry and anisotropy.
Figure 9 shows the correlation times normalised by the terminal time
$t_T$
. For ring-shaped particles,
$\tau _x$
and
$\tau _z$
decrease with increasing
$\phi$
, and exhibit an approximate
$\phi ^{-1/2}$
scaling. This trend reflects the progressive shortening of correlation times due to enhanced many-body interactions. The correlations remain strongly anisotropic, with
$\tau _z$
significantly larger than
$\tau _x$
, indicating preferential persistence of velocity disturbances along the settling direction. This behaviour is consistent with the anisotropic structure of disturbance flows and wake-mediated interactions, which promote sustained streamwise coherence.
The suspension microstructure provides additional insight into the origin of the observed velocity correlations, and is quantified through the radial distribution function
\begin{equation} g(r)=\frac {2L^3}{N^2} \left \langle \sum _{I\lt J}\delta (r-r_{\textit{IJ}}) \right \rangle , \end{equation}
where
$r_{\textit{IJ}}=|\boldsymbol{R}_{\!I}-\boldsymbol{R}_{\!J}|$
is the particle separation distance.
Radial distribution function
$g(r)$
for different volume fractions. The separation distance
$r$
is normalised by the bead diameter
$\alpha$
of the spheres composing the ring particle.
Figure 10 shows
$g(r)$
for different volume fractions. A pronounced primary peak appears near
$r/\alpha \approx 3$
for
$\phi \lt 0.01$
, followed by a broader secondary maximum around
$r/\alpha \approx 4{-}5$
, indicating preferred separation distances imposed by the ring topology. As
$\phi$
increases, the peak amplitudes decrease, and
$g(r)$
approaches unity, reflecting the progressive weakening of spatial ordering as many-body interactions become dominant. The characteristic separation distances identified by
$g(r)$
are consistent with the correlation lengths inferred from the velocity correlation analysis, indicating that the microstructure and hydrodynamic correlations are governed by the same underlying topology-mediated interactions.
Finally, the long-time transport induced by velocity fluctuations is quantified through the hydrodynamic self-diffusion coefficient derived from the mean-square displacement (MSD):
\begin{equation} D_X(t)=\frac {1}{2t}\big \langle \big (R_{\textit{IX}}(t)-R_{\textit{IX}}(0)\big )^2 \big \rangle , \end{equation}
\begin{equation} D_Z(t)=\frac {1}{2t}\big \langle \big (R_{\textit{IZ}}(t)-R_{\textit{IZ}}(0)-V_{\textit{sed}}t\big )^2 \big \rangle , \end{equation}
where
$R_{I\alpha }$
denotes the particle position in direction
$\alpha$
, and
$V_{\textit{sed}}$
is the mean settling velocity.
Self-diffusivity in horizontal,
$D_x$
, and vertical,
$D_z$
, directions for
$\phi = 0.1$
, illustrating the approach to the long-time diffusive regime.
(a) Horizontal and (b) vertical self-diffusion coefficients, normalised by
$V_tD$
, as functions of volume fraction for ring-shaped particles. Included for comparison are rod-like (Hamid et al. Reference Hamid, Arshad, Mehdi, Qasim, Ullah, Molina and Yamamoto2020), spherical (Hamid & Yamamoto Reference Hamid and Yamamoto2013c
), disk-shaped (Hamid et al. Reference Hamid, Shakaib, Molina, Qasim, Fareed, Uzair and Yamamoto2024), cube-shaped (Fareed et al. Reference Fareed, Nadeem, Ullah, Molina, Yamamoto, Chamorro and Hamid2026b
) and triangular (Fareed et al. Reference Fareed, Nadeem, Ullah and Hamid2026a
) particles. The shaded region denotes a narrow envelope about the
$\phi ^{1/2}$
scaling.
Figure 11 shows the temporal evolution of the diffusion coefficients obtained from the MSD slopes, indicating convergence to a diffusive regime at long times. Figure 12 shows the resulting diffusion coefficients as functions of
$\phi$
. For ring-shaped particles, both horizontal and vertical diffusion increase with concentration, and follow an approximate
$\phi ^{1/2}$
scaling. The vertical component remains systematically larger than the horizontal one, reflecting anisotropic particle transport induced by gravity-aligned disturbance flows. The
$\phi ^{1/2}$
dependence shown in figure 12 is used only as an approximate low-concentration reference trend, rather than as a universal scaling law. The vertical diffusivity follows this trend mainly for
$\phi \lesssim 0.009$
, whereas the horizontal diffusivity remains closer to it over a broader range, up to
$\phi \lesssim 0.09$
. Deviations at larger concentrations reflect the increasing influence of many-body hydrodynamic interactions and concentration-dependent effects. Compared with other particle geometries, ring-shaped particles exhibit intermediate diffusivity levels: spheres and disks generally show larger diffusivities, whereas triangular particles exhibit weaker transport.
5. Discussion
The results consistently indicate that the internal flow pathway of ring-shaped particles plays an important role in modifying the disturbance field relative to convex particles. Unlike solid particles that displace fluid through an externally confined wake, the toroidal geometry permits fluid motion through the particle interior, producing orientation-dependent flow structures that interact with the surrounding flow. These features are consistent with modified hydrodynamic disturbances and enhanced viscous dissipation at the single-particle level, although internal recirculation is not uniquely isolated from other geometric effects in our study.
At the suspension scale, this topology-induced modification of the disturbance field provides a unifying mechanism for the observed collective dynamics. The altered flow structure promotes enhanced hydrodynamic coupling along the direction of gravity, resulting in anisotropic velocity fluctuations and extended streamwise correlations. Then hydrodynamic screening is not solely governed by particle concentration, but is also influenced by topology through its effect on disturbance propagation. The observed
$\phi ^{1/3}$
scaling of fluctuations and the finite correlation length can thus be interpreted as consequences of topology-modulated screening, in which internal flow pathways promote coupling over distances of a few particle diameters while limiting long-range divergence.
This mechanism may also contribute to the emergence of preferred microstructural arrangements and anisotropic transport. The characteristic separation distances observed in the radial distribution function reflect geometric constraints imposed by the ring topology in combination with wake-mediated interactions, while the enhanced vertical diffusion arises from the persistence of gravity-aligned disturbance structures. The convergence of orientation distributions at higher concentrations further indicates that many-body interactions stabilise configurations that enhance hydrodynamic coupling to the surrounding flow.
The observations suggest that particle topology constitutes a distinct control parameter in particulate suspensions, separate from conventional shape descriptors such as aspect ratio or sphericity. By enabling internal flow pathways, non-convex geometries modify the mechanisms of disturbance generation, propagation and screening, thereby strongly influencing macroscopic transport properties through topology-driven hydrodynamic interactions. However, the present results do not uniquely isolate internal recirculation from other geometric effects, including surface discreteness, projected area and orientation-dependent drag. Controlled comparisons between smooth tori and bead-based rings, open rings and closed non-convex particles without through-flow, and rings with systematically varied opening size and bead count, would be required to separate these contributions quantitatively.
6. Conclusions
Direct numerical simulations using the smoothed profile method were performed to investigate the sedimentation dynamics of suspensions of non-convex ring-shaped particles in a viscous Newtonian fluid at
$\textit{Re}=0.1$
. By varying the particle volume fraction over
$0.001 \le \phi \le 0.1$
, the role of particle topology in governing hydrodynamic interactions and collective transport was systematically quantified. The simulations were initialised from randomly oriented, spatially uniform particle distributions, and the reported statistics correspond to the subsequent statistically developed suspension dynamics.
For an isolated particle, the drag coefficient follows Stokesian scaling but remains systematically larger than that of convex particles, consistent with enhanced viscous dissipation associated with flow penetration through the ring opening. At the suspension scale, this leads to enhanced hindered settling, with the mean settling velocity following a Richardson–Zaki-type scaling with exponent
$n \approx 7$
. Velocity fluctuations exhibit pronounced anisotropy and follow an approximate
$\phi ^{1/3}$
scaling in the dilute regime, consistent with hydrodynamic screening. Spatial and temporal correlations reveal persistent streamwise coherence, with correlation times decreasing as concentration increases. The suspension microstructure exhibits multiple preferred separation distances imposed by the ring topology, and particle orientations stabilise near
$\theta \approx 90^\circ$
at higher concentrations.
These results demonstrate that particle topology, beyond conventional shape descriptors, constitutes an independent control parameter governing hydrodynamic interactions, microstructure and transport in particulate suspensions. By enabling internal flow pathways, non-convex geometries modify the propagation of disturbances, and lead to enhanced anisotropy and altered collective dynamics. These findings provide insight into sedimentation beyond convex particles, and are directly relevant to systems involving porous or toroidal inclusions, including filtration media, fluidised beds and engineered multiphase flows. Finally, we stress that our study is restricted to the low-Reynolds-number regime and to a specific monodisperse ring geometry; extensions to mixed particle sizes, shapes and topologies, as well as to different initial spatial distributions, may introduce additional segregation, orientational and hydrodynamic-screening effects, and remain important directions for future work.
Acknowledgements
The authors gratefully acknowledge the support of the National Research Programme for Universities (no. 20-15306/NRPU/R&D/HEC/2021), funded by the Higher Education Commission of Pakistan.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Convergence analysis
Figure 13 shows representative results for
$\phi =0.05$
and
$\phi =0.075$
. The blue curves correspond to the cumulative average of the normalised settling velocity, whereas the red curves represent block-averaged values obtained over finite time intervals.
For both volume fractions, the cumulative average gradually approaches an approximately constant value as the simulation proceeds, indicating that the system reaches a statistically stationary state. Although the block-averaged mean exhibits fluctuations due to the intrinsic temporal variability associated with many-body hydrodynamic interactions, these fluctuations remain bounded and do not display any systematic drift. The persistence of oscillations around a nearly constant mean value demonstrates that the reported averages are sufficiently converged. Therefore, the deviations from the Richardson–Zaki scaling observed at higher volume fractions are not attributable to inadequate statistical sampling, but instead reflect physical effects associated with enhanced many-body interactions and near-field hydrodynamic constraints.
Statistical convergence of the normalised mean settling velocity. The blue line denotes the cumulative average, and the red line denotes the block-averaged mean. Results are shown for (a)
$\phi =0.05$
and (b)
$\phi =0.075$
. In both cases, the statistics approach a stationary value and exhibit no systematic drift, indicating adequate temporal averaging.
27Δ×27Δ×28Δ
g
−z
θ
n
z
Vsed/Vt
(1−ϕ)7
σΔVx/Vt
σΔVz/Vt
ϕ
ϕ1/2
ϕ1/3
θ
5∘
ϕ=0.001
0≤t≤T0/2
T0/2≤t≤T0
T0=5×104
ϕ=0.05
θ
θ=0∘
θ=90∘
ϕ−1/2
g(r)
r
α
Dx
Dz
ϕ=0.1
VtD
ϕ1/2
ϕ=0.05
ϕ=0.075
Open access
