On the dynamics and wakes of a freely settling Platonic polyhedron in a quiescent Newtonian fluid

Abstract

We investigate systematically the free settling of a single Platonic polyhedron in an unbounded domain filled with an otherwise quiescent Newtonian fluid. We consider a particle–fluid density mimicking a rock in water. Five Platonic polyhedrons of increasing sphericity are studied for a range of Galileo numbers $10 \leqslant \mathcal {G}a \leqslant 300$. We construct a regime map in the parameter space of Galileo number and particle volume fraction ($\mathcal {G}a, \phi$), highlighting how the angularity of the Platonic polyhedron impacts its settling path and the onset of instabilities. We find that the initial angular position solely affects the transient settling process. All the Platonic polyhedrons maintain a stable settling angular position at low $\mathcal {G}a$. Higher angularity leads to path unsteadiness at lower $\mathcal {G}a$. Path instability progresses from steady vertical to unsteady vertical, followed by oblique settling observed for highly spherical particles, but helical settling (HS) for more angular particles. The particle autorotation is found to be the pivotal factor influencing path instability and the regime transition of angular particles. Beginning in the unsteady oblique and helical regimes, particle autorotation becomes more prevalent, escalating further in the chaotic regime as $\mathcal {G}a$ increases. The particle angular velocity vector is shown to be predominantly situated in the horizontal plane. A thorough force balance in the horizontal plane reveals that the Magnus force is the primary driving force of the HS regime. Additionally, we establish two new empirical correlations to predict the particle settling velocity and the disturbed wake length that solely require the physical properties of the system ($\mathcal {G}a$ and $\phi$). Our numerical results suggest that an increase of the density ratio from $2$ to $3$ exerts only a marginal impact on the path instability of the most angular particle, the settling tetrahedron.

1. Introduction

Particle-laden multiphase flows are frequently encountered in various natural and industrial processes, including sediment transport in rivers, pollutants dispersion in oceans (Martin 1981), fluidized beds in chemical engineering (van der Hoef et al. 2008), bubbly flows (Magnaudet & Eames 2000) and fuel combustion (Fréret, Laurant & de Chaisemartin 2008) in mechanical engineering. The gravity-driven settling of a solid body within a quiescent Newtonian fluid exhibits a diverse range of dynamic features, as the dynamics of one phase influences the other phase in a myriad of ways. Path instabilities are commonly observed, with particle motion patterns ranging from steady vertical to fully chaotic depending on a combination of parameters such as particle size, shape, moment of inertia and the properties of the surrounding flow. Settling or rising bodies under gravity in a Newtonian fluid, which is otherwise at rest, are governed by three control parameters: the body-to-fluid density ratio ($m=\rho _p^*/\rho _f^*$); the Galileo number ($\mathcal {G}a$), based on a characteristic length scale (typically diameter or equivalent diameter) of the particle and the gravitational velocity scale; and a geometric parameter characterizing the shape of the rigid body such as aspect ratio $\chi$ or sphericity $\phi$ (Ern et al. 2012). Path instability in such physical cases can be attributed to two primary causes: (i) the evolution of hydrodynamic loads (forces and torques) in response to a fluid disturbance applied to the body, and (ii) wake instability, which occurs beyond a critical Reynolds number associated with the first bifurcation of the wake (Ern et al. 2012). Consequently, analysing the motion of settling or rising objects necessitates investigating the features of particle wakes.

Many studies choose the sphere as their focus due to its omnidirectional symmetry (Johnson & Patel 1999; Auguste & Magnaudet 2018). Jenny, Bouchet & Dušek (2003), Jenny, Dušek & Bouchet (2004) and Jenny & Dušek (2004) conducted comprehensive numerical research on path instabilities and regime transitions for a freely moving sphere. They explored a wide range of Galileo numbers ($150 \leqslant \mathcal {G}a \leqslant 350$) and density ratios ($0 \leqslant m \leqslant 10$), finding that the critical $\mathcal {G}a$ for the first regular bifurcation is slightly lower for a freely moving particle compared with a fixed one. With a $\mathcal {G}a \geqslant 155$, a steady oblique (SO) path emerges, replacing the initial steady vertical (SV) rectilinear path. In the SO regime, studies using Schlieren visualization (Veldhuis & Biesheuvel 2007) and particle image velocimetry (PIV) measurements (Horowitz & Williamson 2010) identified two counter-rotating trailing vortices. These vortices were found to connect and form a single-sided chain of vortex rings, exerting a lift force that guides the sphere's oblique descent. The critical $\mathcal {G}a$ for the onset of unsteadiness was determined to be $\mathcal {G}a_{cr} = 167$ for massless spheres ($m=0$) and $\mathcal {G}a_{cr} = 196$ for fixed spheres ($m \rightarrow \infty$). In this unstable regime, two oscillatory oblique (OO) path modes coexist, showing low-frequency oscillations for $m < 2.5$ and high-frequency oscillations for $m \geqslant 2.5$.

In addition to solid spheres, axially symmetric particles such as infinitely long cylinders (Namkoong, Yoo & Choi 2008), short cylinders (disks) (Field et al. 1997; Fernandes et al. 2007, 2008; Auguste, Magnaudet & Fabre 2013) and annular disks (Bi et al. 2022), have long been a focal point in investigations of settling dynamics. These bodies, which also include ellipsoids, are among the simplest three-dimensional objects for which the interplay between geometric anisotropy and wake-induced loads can be examined in detail. Such studies provide critical insights into why a freely falling body may deviate from a vertical path. These investigations typically concentrate on non-rectilinear and non-vertical paths, as well as on path instabilities resulting from vortex generation on the disk surface. When viscous effects are prominent, the disk descends along a vertical path with a stable orientation (Mercier et al. 2020). As inertia increases, we start to observe planar path types, which may exhibit fluttering motion or inclined paths due to the disk autorotation (Auguste et al. 2013). At higher Galileo numbers, three-dimensional regimes emerge. Wake instability has been noted to strongly depend on the aspect ratio ($\chi$) and density ratio ($m$) of the particle. Importantly, the angle formed between the body symmetry axis and the settling velocity can modify vortex formation. Beyond axially symmetric particles, recent studies have also explored the free settling of porous particles (Ahmerkamp et al. 2022) and thin curved particles within a quiescent fluid (Chan et al. 2021).

Among all settling or rising paths, the dynamics and mechanisms of zigzag or helical motions of axisymmetric bodies – like spheres, spheroids and disks – have sparked considerable interest (Fernandes et al. 2005, 2008). These motions, driven by lift force and associated torque, are especially prevalent in the case of rising spheroidal bubbles (Mougin & Magnaudet 2001, 2002, 2006). Studies reveal that force and torque balances enable the existence of zigzag and helical motions, where added-mass effects play a crucial role in counterbalancing the bubble wake effect. Similar helical rising patterns have been noted in lightweight rising particles, where high horizontal velocity fluctuations appear concurrent with the emergence of the helical regime (Zhou & Dušek 2015; Ardekani et al. 2016; Auguste & Magnaudet 2018). In the case of freely settling disks, Zhong, Chen & Lee (2011) drew a connection between the helical settling (HS) path and helicoidal vortex shedding in the particle wake. The disk rotation around its symmetry axis was found to have an angular velocity about half that of the rotation around the transverse direction. Later, Auguste et al. further added two non-planar regimes to the lexicon: hula-hoop and helical autorotation (Auguste et al. 2013). In Kim, Chang & Kim (2018) the freely settling motion of two identical, rigidly connected disks was examined, with consistent helical motion across different disk separations. In a recent study, Bi et al. studied freely falling annular disks in quiescent water at relatively low Reynolds numbers ($10 \leqslant {\mathcal {R}e} \leqslant 500$) (Bi et al. 2022). They also observed a helical motion characterized by periodic oscillations in the horizontal plane and a constant inclination angle relative to the vertical direction. The key role of asymmetrical pressure distribution in disk dynamics is highlighted, which exerts fluid forces and torques on the particle.

Now we know that wake instability in non-spherical three-dimensional bodies is predominantly influenced by particle shape and vortex structures. Despite the extensive studies on axially symmetric particles, research on free settling of angular bodies, such as cubes, is scarcely available in the existing literature. Some research has delved into the interaction between flows and a stationary angular particle. In 2004, Saha et al. conducted numerical simulations of laminar flows past a fixed cube for $20 \leqslant {\mathcal {R}e}_{edge} \leqslant 300$, where ${\mathcal {R}e}_{edge}$ represents the Reynolds number based on the edge length of the cube (Saha 2004, 2006). The transition sequence and flow structures were found to closely resemble those of a fixed sphere. Subsequently, Klotz et al. (2014) carried out a comprehensive experimental investigation of the wake behind a fixed cube for $20 \leqslant {\mathcal {R}e}_{edge} \leqslant 400$ using PIV and laser-induced fluorescence visualization tools. They identified the onset of two bifurcations related to regime transitions from steady to unsteady flows. More recently, Meng et al. (2021) performed high-fidelity body-fitting direct numerical simulations of flows past a fixed cube at $1\leqslant {\mathcal {R}e}_{edge} \leqslant 400$ and proposed accurate critical Reynolds numbers for regime transitions using weakly nonlinear instability analysis. In our previous studies, we conducted a systematic investigation of flows past a fixed Platonic polyhedron at Reynolds numbers ranging from $100 \leqslant {\mathcal {R}e} \leqslant 500$ (Gai & Wachs 2023c,a,d). Our findings demonstrate that particle angularity and angular position significantly influence wake structures in both steady and unsteady regimes. In the steady regime, the particle front inclined faces determine the number and distribution of vortex pairs. The shape of the particle determines the recirculation region form, leading to various vortex-shedding patterns with different shedding frequencies. In turn, the dynamic features within the particle wake contribute to the modification of hydrodynamic forces and torques (Gai & Wachs 2023b).

In contrast to the freely settling sphere or spheroid, angular isometric particles display some distinct and unique dynamic characteristics. The helical regime for settling angular particles was first reported in Rahmani & Wachs (2014), replacing the transitional oblique regime for spheres. Both the cube and tetrahedron transitioned through vertical and helical/spiral paths before transitioning to chaotic motions. Force analysis suggests that Magnus forces drive the lateral movement along the helical/spiral path for cubes, while local vortical forces act as counterbalancing forces. However, for a tetrahedron, these roles interchange during the helical motion, revealing a strong dependence of motion dynamics on particle shape. Particles with high angularity, such as tetrahedrons or cubes, exhibit unsteady behaviour at relatively low $\mathcal {G}a$. Seyed-Ahmadi & Wachs (2019) expanded the parameter map for freely settling cubes. A notable feature of cubes moving in helical trajectories is the appearance of a new vortex-shedding mode, where four hairpin (4H) vortices are shed per motion cycle. The drag coefficient $C_d$ increases abruptly when the cube begins to move helically. This increased drag force was considered to originate from two contributing factors: (i) the concept of induced drag, resulting from the misalignment of the main vorticity threads with the particle velocity vector, and (ii) the orientation of the cube with respect to the direction of motion, which alters the drag force.

Despite these significant insights, the dynamics, path instability and the applicability of the current physical understanding to a broader range of angular particles continue to pose unanswered questions and challenges in the field. In this study we investigate the free settling of five Platonic polyhedrons with an increasing number of surface faces in an unbounded quiescent fluid. We aim to provide an improved understanding within a numerical framework of four aspects: (i) a comprehensive regime transition map for freely settling Platonic polyhedrons; (ii) an empirical model of the settling velocity and disturbed wake length for angular particles; (iii) an in-depth comprehension of the HS regime; and (iv) a clearer link between fixed and freely settling particles in fluid flow.

The paper is structured as follows. In § 2 we discuss the governing equations, numerical method, dimensionless numbers and simulation set-up. In § 3 we present numerical validations. In § 4 we delve into various regimes of path instability and transitions, as a function of the Galileo number $\mathcal {G}a$, particle sphericity $\phi$ and density ratio $m$. This section discusses intricate wake structures, rotation characteristics and offers a thorough analysis of drag and drift coefficients. Furthermore, we establish models to comprehend the settling velocity and disturbed wake length. We specifically study the HS regime, gaining insights from a fixed particle in the flow. In § 5 we present the major conclusions and contemplate potential prospectives for future research.

2. Numerical method and dimensionless numbers

The distributed Lagrange multiplier/fictitious domain (DLM/FD) method describes a rigid body immersed in a fluid by creating a fictitious domain $P$, where the rigid-body constraints are enforced so that the fictitious domain behaves like a rigid particle. In the following, we use the $^*$ superscript to denote dimensional quantities throughout the text, while non-dimensional variables are presented without the superscript. We consider a cubic computational domain $\mathcal {D}$ with boundary $\varGamma$ containing a solid sub-domain $P$, which is a rigid body immersed in a Newtonian fluid of constant density $\rho _f^*$ and viscosity $\mu _f^*$. To maintain a strict separation between the fluid and solid domains, the fluid sub-domain is denoted as $\mathcal {D}\setminus P$, such that $\mathcal {D}\setminus P \bigcap P = 0$. The DLM/FD method employed in this paper has been extensively discussed in previous works (Glowinski et al. 1999; Wachs 2009, 2011; Wachs et al. 2015; Selçuk et al. 2021). For the sake of brevity, we provide only a brief summary of the governing equations for the fluid motion and the combined equation of motion for the particle–fluid mixture.

2.1. The DLM/FD method

The equations governing the conservation of mass and momentum of an incompressible Newtonian fluid are expressed as

(2.1)$$\begin{gather} \rho_f^*\left(\frac{\partial \boldsymbol{u}^*}{\partial t^*} + \boldsymbol{u}^*\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}^*\right) ={-}\boldsymbol{\nabla} p^* + \mu_f^* \Delta\boldsymbol{u}^* \quad \mbox{in} \ \mathcal{D}\setminus P , \end{gather}$$

(2.2)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^*=0 \quad \mbox{in} \ \mathcal{D}\setminus P, \end{gather}$$

where $\boldsymbol {u}^*$ represents the velocity vector of the fluid and $p^*$ denotes the fluid pressure; $\mathcal {D}$ is the computational domain and $P$ refers to the solid particle.

The equations governing the motion of the particle are

(2.3)$$\begin{gather} M^*\frac{{\rm d}\boldsymbol{U^*}}{{\rm d}t^*} = M^*\boldsymbol{g^*} + \boldsymbol{F^*}, \end{gather}$$

(2.4)$$\begin{gather}\frac{{\rm d}\boldsymbol{I^*\boldsymbol{\varOmega^*}}}{{\rm d}t^*} = \boldsymbol{T^*}, \end{gather}$$

where $M^*$ denotes the particle mass, $\boldsymbol {I^*}$ the moment of inertia tensor of the particle and $\boldsymbol {g^*}$ the gravity acceleration; $\boldsymbol {F^*}$ and $\boldsymbol {T^*}$ are two vectors representing the hydrodynamic force and torque (with respect to the centre of mass) on the particle. Expressions for these two quantities are as follows:

(2.5)$$\begin{gather} \boldsymbol{F^*} = \int_{\partial P} \boldsymbol{\sigma^*}\boldsymbol{\cdot}\boldsymbol{\hat{n}}\,{\rm d}S^*, \end{gather}$$

(2.6)$$\begin{gather}\boldsymbol{T^*} = \int_{\partial P} \boldsymbol{r_i^*} \times\left(\boldsymbol{\sigma^*}\boldsymbol{\cdot}\boldsymbol{\hat{n}}\right){\rm d}S^*. \end{gather}$$

Here $\boldsymbol {\hat {n}}$ denotes the outward-oriented unit normal vector to the surface of the particle $\partial P$; $\boldsymbol {\sigma ^*}= -p^*\boldsymbol {I_d}+2\mu _f^*\boldsymbol {D}[\boldsymbol {u}]^{\boldsymbol {*}}$ is the stress tensor for a Newtonian fluid.

The DLM/FD method employs a combination of the weak formulation of the fluid motion equation and that of the rigid particle to obtain an equation of motion for the fluid–particle mixture. To begin with, the weak formulation for the fluid sub-domain $\mathcal {D} \setminus P$ is obtained by imposing rigid-body motion constraints on the particle surface $\partial P$. The formulation is then extended to the entire computational domain $\mathcal {D}$ by imposing the rigid-body motion constraints on the entire particle $P$, leading to the derivation of the combined equation of motion. The rigid-body motion constraint in the particle sub-domain is relaxed using Lagrange multipliers $\boldsymbol {\lambda }^*$. We denote $\mathcal {W}_\varGamma$, $\mathcal {W}_0$, $\mathcal {L}_0$ and $\varLambda$ the functional spaces of the solution satisfying the boundary conditions (Selçuk et al. 2021). Subsequently, we solve the following constrained optimization problem in the case of a single moving particle $P$: find $\boldsymbol {u^*} \in \mathcal {W}_\varGamma$, $p^*\in \mathcal {L}^2_0$ and $\boldsymbol {\lambda }^*\in \varLambda$ such that

1. (i) combined equations of motion:

   (2.7)\begin{align} &\int_{\mathcal{D}} \rho_f^* \left(\frac{\partial \boldsymbol{u}^*}{\partial t^*} + \boldsymbol{u}^*\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}^*\right)\boldsymbol{\cdot}\boldsymbol{v}^*\, {\rm d}\kern0.7pt \boldsymbol{x} - \int_{\mathcal{D}}p^*\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}^*\,{\rm d}\kern0.7pt \boldsymbol{x} \nonumber\\ &\quad+ \int_{\mathcal{D}}\mu_f^*\boldsymbol{\nabla}\boldsymbol{u}^*:\boldsymbol{\nabla} \boldsymbol{v}^*\,{\rm d}\kern0.7pt \boldsymbol{x} ={-}\int_{P}\boldsymbol{\lambda}^*\boldsymbol{\cdot} \boldsymbol{v}^*\,{\rm d}\kern0.7pt \boldsymbol{x}, \end{align}

2. (ii) continuity equation:

   (2.8)\begin{equation} \int_{\mathcal{D}} -q^*\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}^*\,{\rm d}\kern0.7pt \boldsymbol{x} = 0, \end{equation}

3. (iii) for the particle:

   (2.9)$$\begin{gather} \left(1-\frac{\rho_f^*}{\rho_p^*}\right)M^*\left(\frac{{\rm d}\boldsymbol{U^*}}{{\rm d}t^*} - \boldsymbol{g^*}\right)\boldsymbol{\cdot}\boldsymbol{V^*} = \int_{P}\boldsymbol{\lambda}^*\boldsymbol{\cdot}\boldsymbol{V}^*\,{\rm d}\kern0.7pt \boldsymbol{x} \end{gather}$$
   (2.10)$$\begin{gather}\left(1-\frac{\rho_f^*}{\rho_p^*}\right)\frac{{\rm d}\boldsymbol{I^*\varOmega^*}}{{\rm d}t^*}\boldsymbol{\cdot} \boldsymbol{\xi^*} =\int_{P}\boldsymbol{\lambda}^* \boldsymbol{\cdot}\left(\boldsymbol{\xi^*}\times\boldsymbol{r^*}\right){\rm d}\kern0.7pt \boldsymbol{x} \end{gather}$$
   (2.11)$$\begin{gather}\int_{P}\boldsymbol{\alpha}^*\boldsymbol{\cdot}\left(\boldsymbol{u^*} - (\boldsymbol{U^*} + \boldsymbol{\varOmega^*}\times\boldsymbol{r^*})\right){\rm d}\kern0.7pt \boldsymbol{x} = 0 \end{gather}$$
   or all $\boldsymbol {v}^*\in \mathcal {W}_0$, $(\boldsymbol {V}, \boldsymbol {\xi })\in \mathbb {R}^3\times \mathbb {R}^3$ $q^*\in \mathcal {L}^2(\mathcal {D})$ and $\boldsymbol {\alpha }^*\in \varLambda$. The equations presented in (2.7)–(2.8) can be modified to a non-variational form to make them more suitable for spatial discretization using finite volume or finite difference methods (Wachs et al. 2015; Selçuk et al. 2021).

We employ the DLM/FD solver developed and implemented by Selçuk et al. (2021) in Basilisk, an open-source code known for its ability to solve Navier–Stokes equations with octree adaptively refined Cartesian grids and the efficient convergence of its multigrid solvers for Poisson/Helmholtz-type problems (Popinet 2015). This method offers the advantage of unconditional stability with regard to the particle–fluid density ratio, improved fluid–solid coupling through the implicit hydrodynamic interaction computation and robust convergence properties of the Uzawa algorithm (Yu 2005; Wachs et al. 2015; Xia et al. 2022). For problems involving freely moving particles, we use our in-house Grains3D solver (Wachs et al. 2012; Rakotonirina et al. 2019), which is well suited for handling collisions between rigid particles of complex shapes. Additional details on the implementation and validation of the DLM/FD solver can be found in previous works by our group (Wachs 2011; Wachs et al. 2015; Selçuk et al. 2021; Gai & Wachs 2023a).

2.2. Dimensionless numbers

The system of (2.1)–(2.4) can be made dimensionaless by introducing the following scales: $L_{ref}^*=D_{sph}^*$ for length, where $D_{sph}^*$ denotes the diameter of the volume-equivalent sphere; $U_{ref}^* = \sqrt {|m-1|g^*D_{sph}^*}$ for velocity, where $m = \rho _p^* / \rho _f^*$ denotes the density ratio between the solid particle and the fluid; $L_{ref}^*/U_{ref}^*$ for time; $\mu _f^* U_{ref}^*/L_{ref}^*$ for stress; $\rho _f^* U_{ref}^{*,2}$ for pressure; $\rho _s^* L_{ref}^{*,2}$ for particle mass and $\rho _s^* L_{ref}^{*,4}$ for particle moment of inertia tensor. Consequently, the relevant dimensionless quantities in this study include: (i) the particle sphericity $\phi$, Galileo number $\mathcal {G}a$ and particle–fluid density ratio $m$ as input parameters; and (ii) the drag coefficient $C_d$, particle settling Reynolds number ${\mathcal {R}e}$, particle angular velocity $\boldsymbol {\varOmega }$ and dimensionless vorticity $\boldsymbol {\omega }$, among other relevant quantities, as output parameters.

The Galileo number $\mathcal {G}a$ is a measure of the ratio of gravity force to viscous force in a fluid flow, defined as

(2.12)\begin{equation} \mathcal{G}a = \frac{\rho_f^* D^*_{sph} U_{ref}^*}{\mu_f^*}. \end{equation}

Throughout this study, we refer to the sphere with the same volume as that of the angular particle as the volume-equivalent sphere, and therefore, $D_{sph}^*$ is defined as

(2.13)\begin{equation} D_{sph}^* = \left( \frac{6}{\rm \pi} V_{Platonic}^* \right)^{1/3}, \end{equation}

where $V_{Platonic}^*$ is the volume of the Platonic polyhedron.

The particle sphericity $\phi$ characterizes how much an angular particle looks like a perfect sphere. To be precise, $\phi$ is defined as the ratio of the surface area of a sphere $S_{sph}^*$ to the surface area $S_p^*$ of an angular particle with the same volume, as expressed in the equation below:

(2.14)\begin{equation} \phi = \frac{S_{sph}^*}{S_p^*}. \end{equation}

The closer $\phi$ is to $1$, the more closely the angular particle resembles a sphere.

The particle Reynolds number ${\mathcal {R}e}$ is defined as

(2.15)\begin{equation} {\mathcal{R}e} = \frac{\rho_f^* \overline{U_m^*}D_{sph}^*}{\mu_f^*}, \end{equation}

where $\overline {U_m^*}$ is the time-averaged settling velocity of the particle, calculated over a time period when the settling regime is fully established.

We want to remind the reader that the particle Reynolds number ${\mathcal {R}e}$ defined here is different from the edge Reynolds number ${\mathcal {R}e}_{edge}$ used in previous studies of the flow past a fixed cube (Saha 2006; Meng et al. 2021). The fact that Platonic polyhedrons with the same volume have different edge lengths justify the choice of the diameter of the volume-equivalent sphere $D_{sph}^*$ as the characteristic length throughout our analysis.

The drag coefficient $C_d$ of a freely settling particle is defined as

(2.16)\begin{equation} C_d = \frac{8|F_y^*|}{{\rm \pi} \rho_f^*\overline{U_m^*}^{2}D_{sph}^{*,2}}, \end{equation}

where $F_y^*$ denotes the streamwise component of the hydrodynamic force $\boldsymbol {F}^*$ exerted on the Platonic polyhedron which can be computed as follows:

(2.17)\begin{equation} \boldsymbol{F}^*= \int_{P} \boldsymbol{\lambda}^*\,{\rm d}\kern0.7pt \boldsymbol{x} + \frac{M^*}{m}\frac{{\rm d}\boldsymbol{U^*}}{{\rm d}t^*}. \end{equation}

We define the dimensionless hydrodynamic forces using the reference velocity $U_{ref}^*$ as

(2.18a–c)\begin{equation} F_x = \frac{8F_x^*}{{\rm \pi} \rho_f^* U_{ref}^{*,2} D_{sph}^{*,2}} ,\quad F_y = \frac{8F_y^*}{{\rm \pi} \rho_f^* U_{ref}^{*,2} D_{sph}^{*,2}},\quad F_z = \frac{8F_z^*}{{\rm \pi} \rho_f^* U_{ref}^{*,2} D_{sph}^{*,2}}. \end{equation}

Similarly, the dimensionless torques can be calculated as

(2.19a–c)\begin{equation} T_x = \frac{8T_x^*}{{\rm \pi} \rho_f^* U_{ref}^{*,2} D_{sph}^{*,3}},\quad T_y = \frac{8T_y^*}{{\rm \pi} \rho_f^* U_{ref}^{*,2} D_{sph}^{*,3}},\quad T_z = \frac{8T_z^*}{{\rm \pi} \rho_f^* U_{ref}^{*,2} D_{sph}^{*,3}}, \end{equation}

where $T_i^*$ denotes the $i\in (x,y,z)$ component of the hydrodynamic torque $\boldsymbol {T^*}$:

(2.20)\begin{equation} \boldsymbol{T^*}= \int_{P} \boldsymbol{r}\times\boldsymbol{\lambda}^*\,{\rm d}\kern0.7pt \boldsymbol{x} + \frac{1}{m}\frac{{\rm d}\boldsymbol{\boldsymbol{I^*}\varOmega^*}}{{\rm d}t^*}. \end{equation}

Please note that in the context of the DLM/FD method (Seyed-Ahmadi & Wachs 2020), the integral of Lagrange multipliers over $P$ is used to calculate the hydrodynamic force $\boldsymbol {F}^*$ and torque $\boldsymbol {T^*}$ in (2.17) and (2.20).

To retain consistency when analysing the rotation of different particles, we define the dimensionless particle angular velocity as

(2.21)\begin{equation} \boldsymbol{\varOmega} = \frac{\boldsymbol{\varOmega^*} D_{sph}^*}{2U_{ref}^*}, \end{equation}

where $\boldsymbol {\varOmega ^*}$ denotes the particle angular velocity.

The dimensionless vorticity, a measurement of fluid rotation, is defined as

(2.22a–c)\begin{equation} \omega_x = \frac{\omega^*_xD_{sph}^*}{U_{ref}^*}, \quad \omega_y = \frac{\omega^*_yD_{sph}^*}{U_{ref}^*}, \quad \omega_z = \frac{\omega^*_zD_{sph}^*}{U_{ref}^*}, \end{equation}

with $\omega ^*_x$ ($\omega ^*_y, \omega ^*_z$) the corresponding component of the dimensional vorticity in the $x$ ($y, z$) direction.

2.3. Numerical set-up

Platonic polyhedrons are five isometric polyhedrons with increasing sphericity from the tetrahedron ($\phi =0.67$) to the icosahedron ($\phi = 0.94$). These polyhedrons have an increasing number of faces from the tetrahedron ($4$ faces) to the icosahedron ($20$ faces), as illustrated in figure 1. They are good candidates to investigate the effects of particle angularity on the flow, particularly the transition from angular particles to a perfect sphere ($\phi =1$). While the free settling of a cube was systematically investigated by Seyed-Ahmadi & Wachs (2019), we aim to expand the range of sphericity of the freely settling particles to all the Platonic polyhedrons. This will allow us to gain more insights into the role of particle sphericity in the particle settling regime transitions.

Figure 1. Platonic polyhedrons in order of increasing number of faces and sphericity $\phi$.

Figure 2 presents the computational domain of the simulations performed in this study. The Platonic polyhedron is initially positioned close to the top centre of the cubic domain of edge length $L$ and released to freely settle in the Newtonian fluid. In our simulations we consider particles of density ratios $m=2$ and $m=3$. Note that the size of the settling particle and vortex structure in figure 2 is for illustrative purposes only. For a clearer representation of the particle size and the computational domain, please refer to figure 3.

Figure 2. Numerical set-up of a freely settling Platonic polyhedron in a cubic computational domain of edge length $L=700$ filled with a Newtonian fluid initially at rest.

Figure 3. Illustration of the octree grid in the three-dimensional cubic computational domain of side length $L=700$; red box: zoomed wake structures identified by the isosurface of $\lambda _2 = -1$ and coloured by $u_y$ for a freely settling cube at $\mathcal {G}a=180$ with adaptive grid; blue box: position of the cube with adaptive grid in the $x\unicode{x2013}y$ cut plane at $z=350$.

We refer to the cubic computational domain boundaries (i.e. the cube faces) as left and right in the $x$ direction, top and bottom in the $y$ direction and front and behind in the $z$ direction such that $\varGamma = left\cup right\cup top \cup bottom \cup front\cup behind$. Here $\boldsymbol {u}^*$ satisfies homogeneous Dirichlet boundary conditions on the left, right, bottom, front and behind boundaries as well as no-slip on the particle surface $\partial P^*$, and homogeneous Neumann boundary conditions on the top boundary. The particle is initially at rest and the complete set of boundary and initial conditions reads as

(2.23)$$\begin{gather} \boldsymbol{u}^*(\boldsymbol{x}^*,t^*) = (0,0,0) \quad \mbox{on} \ \varGamma\setminus top, \end{gather}$$

(2.24)$$\begin{gather}\boldsymbol{u}^*(\boldsymbol{x}^*,t^*) = \boldsymbol{U^*} + \boldsymbol{\varOmega^*}\times\boldsymbol{r^*} \quad \mbox{on} \ \partial P , \end{gather}$$

(2.25)$$\begin{gather}\frac{\partial{\boldsymbol{u}^*}}{\partial{x^*}}(\boldsymbol{x}^*,t^*) = (0,0,0) \quad \mbox{on} \ top , \end{gather}$$

(2.26)$$\begin{gather}\boldsymbol{u}^*(\boldsymbol{x}^*,0) = (0,0,0) \quad \mbox{in} \ \mathcal{D} , \end{gather}$$

(2.27)$$\begin{gather}\boldsymbol{U}^*(0) = \boldsymbol{\varOmega}^*(0) = (0,0,0) \quad \mbox{in} \ P, \end{gather}$$

where $\boldsymbol {x}^*=(x^*,y^*,z^*)$ is the position vector. We also assign an arbitrary zero reference value to the pressure $p^*$ at the right boundary.

In the Cartesian octree adaptive grid strategy implemented in Basilisk, a parent cube cell is partitioned into $8$ sub-cubes for local mesh refinement of specific regions of interest (Popinet 2015; Selçuk et al. 2021). At each time step, the Cartesian octree grid undergoes dynamic adaptation wherein the grid is refined in regions exhibiting strong gradient variations in any field of interest, while the regions with weak gradient variations are coarsened. The flow velocity is our primary field of interest in this study. To ensure accurate modelling of the fluid–particle interactions, a phase indicator field is employed, which takes a value of $0$ in the fluid and $1$ in the solid particle. This ensures that the region in the vicinity of the particle surface is always resolved with the finest grid resolution, guaranteeing that the flow structures in both the boundary layer and wake region are well captured. The hierarchical grid is constructed such that the cell size between two successive levels differs by exactly a factor of $2$. Consequently, the smallest cell size ought to be $\Delta x = L/2^{n_l}$, where $n_l$ denotes the maximum refinement level of the octree grid. In figure 3 the vortex structures in the wake of a freely settling cube at $\mathcal {G}a=180$ and $m=2$ are shown within a red box, where they are identified by the isosurface of $\lambda _2=-1$. The locally refined grid around the vortex structures are highlighted. In the blue box, a cut plane at $z=350$ reveals the instantaneous particle position in the grid. Notably, the difference in size between the particle and the computational domain is apparent. Figure 3 additionally displays the entire octree grid within the cubic three-dimensional computational domain, which underscores the benefits and advantages conferred by the adaptive mesh refinement to model a quasi unbounded domain.

In our study we implement a collocation point method (Glowinski et al. 1999; Yu et al. 2002; Wachs et al. 2015) within the DLM/FD solver to impose the rigid-body motion constraints on the Platonic polyhedron at the discrete level. The complete set of points is comprised of interior points chosen as grid nodes located inside the particle and surface Lagrangian points. The surface Lagrangian points are distributed uniformly parallel to the edges on the particle surface for the five Platonic polyhedrons, as shown in figure 4. To improve the shape resolution of the Platonic polyhedrons, we place points on all the edges and vertices of the particles. This distribution strategy, known as the parallel point set, has been shown to deliver a computed solution of satisfactory spatial accuracy (Wachs et al. 2015) and can be successfully extended to the Platonic polyhedrons (Gai & Wachs 2023c,a).

Figure 4. Lagrange points distribution (white dots) on the surface of the five Platonic polyhedrons, with the edge length $l_{edge}$ and the number of faces $n_{face}$.

In figure 4 we see that the surfaces of the tetrahedron (4 faces), octahedron (8 faces) and icosahedron (20 faces) are composed of triangles, and the parallel point set forms the vertices of a regular equilateral triangular distribution of grid points on each face of the particle. On the other hand, the surfaces of the cube (6 faces) and the dodecahedron (12 faces) are composed of squares and pentagons, respectively. We add a point on the arithmetic centre of the equilateral pentagon of the dodecahedron and divide the pentagon into five triangles by connecting the centre point with the face vertices. Then, we distribute points on the centre-vertex edges and inside each of the five small triangles using the same method as for tetrahedron, octahedron and icosahedron.

To ensure a reliable hydrodynamic force computation, we choose a uniform point-to-point distance of $l_{pp} = 2\Delta x$, which is compatible with the size of the grid cells surrounding the particle. As a result, the number of surface points on the particle increases with the refinement level $n_l$, leading to a more precise description of the particle. From figure 4, we observe that the edge length of the dodecahedron is significantly smaller (only $1/4$) compared with that of the tetrahedron. Therefore, for particles with high sphericity, it is crucial to choose an appropriate grid resolution to ensure satisfactory spatial accuracy and a proper geometric description.

Regarding the time resolution, we select the advective time scale $t_{ref}^* =D_{sph}^*/U_{ref}^*$ as the characteristic time scale. In all our simulations the dimensionless time step $\Delta t = \Delta t^*/t_{ref}^*$ is less than $10^{-3}$ to reduce the numerical error due to the use of an operator splitting solution algorithm. Furthermore, the Courant–Friedrich–Levy condition is fulfilled to ensure the numerical stability of the explicit treatment of the advection term in the momentum conservation equation. Mathematically, the time step $\Delta t$ is dynamically updated as follows:

(2.28)\begin{equation} \Delta t=\min \left\{ 10^{{-}3},\min_{i} \frac{0.8\Delta x_i^*}{\left| \boldsymbol{u}_i^* t_{ref}^*\right|} \right\}. \end{equation}

Here, $i$ represents an index covering all grid cells. By simulating over an adequately extended physical time, we ensure that the flow regimes are fully established in our simulations.

3. Validations

We begin by validating our numerical approach through simulations of a settling sphere in an unbounded domain. We consider a range of $\mathcal {G}a$ and particle density ratios $m$ to investigate the settling and drift velocities of the particle, ensuring that they agree with the established literature. Next, we perform a series of simulations to study the settling behaviour of a cube. Our focus is on the evolution of the settling velocity $U_m$ and the drag coefficient $C_d$ as a function of time.

3.1. Settling of a sphere

We consider four cases from the experiments of Mordant & Pinton (2000) that investigated the settling dynamics of solid spherical beads of various materials in water. These cases have also been used as references in previous numerical works (Uhlmann 2005; Uhlmann & Dušek 2014). To simulate the settling sphere, we use a cubic computational domain of size $L = 32$ with the particle released from the top and both the fluid and particle initially at rest. We choose the parameters listed in table 1 to match the experimental conditions, ensuring the same particle density ratio $m$, Froude number $\mathcal {F}r$ and particle Reynolds number ${\mathcal {R}e}$ across all cases.

Table 1. Experimental and numerical parameters for freely settling spheres in water (Mordant & Pinton 2000).

We use different reference values for velocity and time depending on the various case considered. For the $S1$, $S2$ and $S4$ cases, we adopt the same reference values as in the experiments (Mordant & Pinton 2000), namely $U_{ref}^{M,*} = \sqrt {g^*D_{sph}^*}$ and $t_{ref}^{M,*} = \sqrt {D_{sph}^*/g^*}$ for velocity and time, respectively. By contrast, for the $S5$ case, we use the reference value for velocity, namely $U_{ref}^* = \sqrt {|m-1|g^*D_{sph}^*}$, to enable comparison with the results of Uhlmann & Dušek (2014). The number of smallest grid cell per volume-equivalent diameter $1/\Delta x$ is chosen to be $42$.

Figure 5 presents the temporal evolution of the settling velocity $U_y$ and the drift velocity $U_d$ of a sphere settling freely under the flow conditions listed in table 1. Our numerical results for cases $S1$, $S2$ and $S4$ (in solid lines) agree well with the experimental results of Mordant & Pinton (2000) (in scatter points), as shown in figure 5(a), in terms of both the transient evolution and the steady settling velocities of $U_y$. We observe that in these three cases, the $U_y$ reaches a steady state after $t>20$ and remains constant throughout the settling process.

Figure 5. Temporal evolution of $U_y$ and $U_d$ of a freely settling sphere, in comparison with experimental results (Mordant & Pinton 2000) and numerical results (Uhlmann & Dušek 2014) in the literature. (a) Settling velocity. (b) Settling velocity. (c) Drift velocity.

In order to validate our DLM/FD solver, we performed a cross-validation with an in-house lattice Boltzmann method (LBM) solver developed in our group for case S5 (Cheng & Wachs 2022). The temporal evolution of the settling and drift velocities obtained from the DLM/FD and LBM solvers are shown in figures 5(b) and 5(c), respectively. In addition, we included the reference values proposed by Uhlmann & Dušek (2014) using the immersed boundary method (IBM) and the spectral-element method (SEM) simulations for comparison. As depicted in figure 5(b), we observe a satisfactory match between the DLM/FD and LBM solvers for the temporal evolution of the settling velocity. Both solvers converge to the same value as the SEM method, which has been claimed to be more accurate than the IBM method.

On the other hand, figure 5(c) shows that the temporal evolution of the drift velocity differs between the DLM/FD and LBM methods. However, both methods exhibit an increase in the drift velocity that reaches a steady value after a few oscillations during the transient period. The reference values of the IBM and SEM methods exhibit a difference of $0.025$ in figure 5(c) as reported in Uhlmann & Dušek (2014). We observe that the converged value of $U_d$ obtained by our DLM/FD solver in the steady state lies between the IBM and SEM reference values, while the LBM solver gives a value closer to the IBM reference value. As the reference value using the SEM method was claimed to be more accurate than the values obtained by the IBM method (Uhlmann & Dušek 2014), the steady $U_d$ obtained from our DLM/FD solver is deemed to be satisfactory. Consequently, we conclude that our DLM/FD solver can provide reliable results about the dynamic features and the regime determination of the settling sphere under the numerical configurations investigated.

3.2. Settling of a cube

We further validate our solver using a freely settling cube in an unbounded domain, similar to the configuration investigated by Seyed-Ahmadi & Wachs (2019). We consider a cubic computational domain of size $L = 700$ as illustrated in figure 3, which is large enough to be considered unbounded. The validation case has a Galileo number of $\mathcal {G}a = 70$ and a particle–fluid density ratio of $m=2$. The cube is initially released without any velocity from the top of the simulation domain, where the fluid is at rest. As gravity acts on the cube with density higher than the fluid ($m=2$), the cube accelerates until it reaches a steady state.

To thoroughly validate our solver and investigate the effect of mesh refinement, we compare the results computed by the present DLM/FD numerical method implemented on octree grids to results computed by our legacy code PeliGRIFF (parallel efficient library for grains in fluid flow). PeliGRIFF features a DLM/FD method implemented on regular fixed Cartesian grids (Wachs et al. 2015). Figure 6 shows the temporal evolution of $U_y$ of a freely settling cube in an unbounded domain, simulated by both solvers. A mesh refinement study is performed, where we increase the spatial resolution to $1/\Delta x = 93$ to investigate the effect of mesh refinement on the accuracy of the $U_y$, as shown in figure 6. We observe that using a mesh resolution of $1/\Delta x = 23$ in our adaptive DLM/FD solver captures the $U_y$ evolution in terms of both transient and steady state, and agrees well with the results of PeliGRIFF. As the spatial resolution increases, the evolutions of $U_y$ obtained by our DLM/FD solver tend to converge to the same value. Considering the balance between the computing cost and the spatial accuracy, we choose the spatial resolution $1/\Delta x = 30 \sim 40$ in our numerical simulations.

Figure 6. Temporal evolution of $U_y$ for a freely settling cube at $\mathcal {G}a=70$ and $m=2$ with different mesh resolution $1/\Delta x$; in comparison with the numerical results obtained using PeliGRIFF.

4. Results

4.1. Initial angular position

Figure 7(b) illustrates the settling paths in the $x\unicode{x2013}y$ plane of cubes at different initial angular positions, denoted by an angle $\theta$ in figure 7(a). Despite initial differences, all cubes eventually adopt a vertical path after settling over a sufficiently long distance, for instance, when $y < -30$. The steady settling velocity converges to $U_y = - 0.92$ for all settling cubes. For cubes at symmetric inclined angles, such as those with $\theta = 10^\circ$ and $\theta = 80^\circ$, their settling paths remain symmetric with respect to the central vertical line during the initial settling stage, specifically at $y \geqslant -20$. As the settling continues, the paths of these dual-angle cubes become less symmetric due to the wake instability. The initial angular positions primarily affect the early transient dynamics of particle settling. Given enough time, the particles invariably reach a fully developed state of vertical settling. In this study we classify the settling regimes according to their fully developed paths. Our discussions focus on the force balance and wake structures of settling Platonic polyhedrons within these distinctive regimes.

Figure 7. (a) Schematic representation of a cube inclined at an angle $\theta$ within the $x\unicode{x2013}y$ plane; (b) settling path of a cube in the $x\unicode{x2013}y$ plane at $\mathcal {G}a=70$ and $m = 2$, showcasing the influence of various initial angular positions.

4.2. Regimes of motion

Figure 8 exhibits the settling regime map for the five Platonic polyhedrons, characterized by sphericity spanning $0.67 \leqslant \phi \leqslant 0.94$ and Galileo numbers spanning $10 \leqslant \mathcal {G}a \leqslant 300$. Cases belonging to the same regime are marked by the same colour, with corresponding upper and lower bounds in $\mathcal {G}a$. We proceed to examine and discuss the following settling regimes.

1. (i) Steady vertical regime (SV, $\sqsubset \!\sqsupset$, filled black): in which the particle settles through a straight vertical path, with multi-planar symmetric wake structures identified by streamwise vorticity $\omega _y$.

2. (ii) Unsteady vertical regime (UV, $\sqsubset \!\sqsupset$, filled orange): in which the particle exhibits a vertical path with small horizontal fluctuations and unsteady wake structures in $\omega _y$.

3. (iii) Steady oblique regime (SO, $\sqsubset \!\sqsupset$, filled brown): in which the particle maintains a steady straight oblique settling path.

4. (iv) Unsteady oblique regime (UO, $\sqsubset \!\sqsupset$, filled maroon): in which the particle exhibits an oblique but fluctuating settling path.

5. (v) Helical settling regime (HS, $\sqsubset \!\sqsupset$, filled red): in which the particle settles following a helical (spiral) trajectory, as a result of horizontal drift and rotational movement.

6. (vi) Chaotic settling regime (CS, $\sqsubset \!\sqsupset$, filled blue): in which the particle settling path is unpredictable and irregular.

Figure 8. Particle settling regimes in function of $\mathcal {G}a$ for Platonic polyhedrons at $m=2$; steady vertical (SV, $\sqsubset \!\sqsupset$, filled black), unsteady vertical (UV, $\sqsubset \!\sqsupset$, filled orange), steady oblique (SO, $\sqsubset \!\sqsupset$, filled brown), unsteady oblique (UO, $\sqsubset \!\sqsupset$, filled maroon), helical settling (HS, $\sqsubset \!\sqsupset$, filled red) and chaotic settling (CS, $\sqsubset \!\sqsupset$, filled blue) regime; onset of vortex shedding for Platonic polyhedrons: tetra ($\blacktriangle$), cube ($\blacksquare$), octa ($\blacklozenge$), dodeca (, filled black), icosa (, filled black) and sphere ($\bullet$).

At low $\mathcal {G}a$, all Platonic polyhedrons follow a SV settling path. As $\phi$ increases, the SV regime persists at higher values of $\mathcal {G}a$. The icosahedron, with $\phi = 0.94$, maintains a SV settling path up to $\mathcal {G}a = 100$. Beyond the SV regime, an unsteady vertical (UV) settling regime emerges. In this regime the particle drift velocity becomes non-negligible, although the horizontal displacement during settling typically remains less than $2$ with a settling distance up to $600$. In the UV regime the frequency of horizontal fluctuations in the particle settling path differ among angular particles.

Subsequently, for particles with high angularity (low $\phi$), such as tetrahedrons and cubes, a HS regime emerges. This regime is observed for tetrahedrons at $60 \leqslant \mathcal {G}a \leqslant 150$ and for cubes at $180 \leqslant \mathcal {G}a \leqslant 200$. Please note that with a higher spatial resolution, the range of $\mathcal {G}a$ for cubes in this study is slightly higher but reasonably in agreement with the values reported in Seyed-Ahmadi & Wachs (2019). In the HS regime the Magnus force (due to particle rotation), the added mass (resulting from particle acceleration) and the vortex shedding significantly influence particle motion, leading to a HS path. For particles with high sphericity (octahedron, dodecahedron and icosahedron), two additional settling regimes are identified: SO and unsteady oblique (UO). In the SO regime the particle follows a stable, straight oblique path, whereas in the UO regime, small oscillations in the path occur at higher $\mathcal {G}a$. In the SO and UO settling regimes, the horizontal range of particle drift can span up to $10 \sim 20$ at a settling distance around $600$. At $\mathcal {G}a \geqslant 240$, the particle settling path becomes irregular and chaotic. The transition to the chaotic settling (CS) regime occurs at lower $\mathcal {G}a$ for highly angular particles, specifically at $\mathcal {G}a = 160$ for the tetrahedron and $\mathcal {G}a = 210$ for the cube. To better clarify the effects of particle angularity on the regime map, we incorporate the regimes of a settling sphere in figure 8. Based on the results of Jenny et al. (Jenny et al. 2004), we perform additional simulations to distinguish the SV and UV regimes and identify the critical $\mathcal {G}a$ for the onset of vortex shedding. The regime transitions observed for the sphere closely resemble that of the icosahedron, with the critical $\mathcal {G}a$ being slightly higher for the sphere. Furthermore, the icosahedron, owing to its high sphericity, also exhibits comparable settling dynamics to the sphere. Additional discussion on the critical ${\mathcal {R}e}$ for the onset of vortex shedding and a comparison between freely settling and fixed particles can be found in Appendix A.

A notable trend in figure 8 is that the oblique (SO, UO) and HS regimes are situated at the transition between the vertical (SV, UV) and chaotic (CS) regimes, with their upper and lower bounds in $\mathcal {G}a$ increasing with $\phi$. These two transitional regimes characterize the establishment of path unsteadiness. As we demonstrate, they share similar dynamics but exhibit distinct settling trajectories due to the autorotation of the angular particle.

4.2.1. Vertical settling: from steady to unsteady

At low $\mathcal {G}a$, the fluid viscous effects are significant, causing the particle to settle at a relatively low velocity. In the absence of hairpin vortex structures, the settling path transitions from SV to UV with the increase of $\mathcal {G}a$, a behaviour observed across particles of various shapes (Auguste et al. 2013; Bi et al. 2022). Figure 9 presents the paths of a freely settling tetrahedron and icosahedron at different $\mathcal {G}a$ in the SV and UV regimes in the $x\unicode{x2013}y$, $z\unicode{x2013}y$ and $x\unicode{x2013}z$ planes. The vertical settling regime is present at low $\mathcal {G}a$ values for all Platonic polyhedrons; here, we illustrate only two particles with the lowest (tetrahedron, $\phi =0.67$) and the highest sphericity (icosahedron, $\phi =0.94$). The transient lateral displacement of the tetrahedron observed in figures 9(a)–9(c) originates from the fact that the particle rotates from its initial angular position (illustrated at $t_{ini}$), to the stable face-down angular position (illustrated at $t_{fin}$) in both the SV and UV regimes. As reported in Gai & Wachs (2023d), angular positions with the maximum number of streamwise vortex pairs $\omega _y$ are preferred in the low $\mathcal {G}a$ cases. For a tetrahedron, both the vertex-down and face-down angular positions exhibit three pairs of vortices, more than the edge-down angular position. Considering that the mass centre of the tetrahedron is closer to the face than the vertex, the face-down angular position is more stable than the vertex-down position. Once the stable angular position is established, the tetrahedron follows a SV settling path at $\mathcal {G}a = 10$ in the SV regime. As $\mathcal {G}a$ increases, the settling path becomes unsteady, displaying minor fluctuations in the horizontal plane. We categorize the settling path as representative of the UV regime if the particle fluctuating horizontal displacement is around or less than $1$ after a settling distance of $300$. As $\mathcal {G}a$ increases, we note that the time required to establish the stable angular position decreases.

Figure 9. Path instability of a freely settling Platonic polyhedron in the SV and UV regime: (a–c) tetra, (d–f) icosa.

The icosahedron is initially released with an edge-down angular position at $t_{ini}$, as depicted in figures 9(d)–9( f). In the SV regime, at $\mathcal {G}a \leqslant 50$, the settling icosahedron maintains the edge-down angular position and exhibits a straight vertical trajectory. Since the initial angular position is the same as its stable angular position in the SV regime, there is no horizontal drift observed in the case of $\mathcal {G}a=50$. However, in the UV regime the stable angular position becomes vertex down (illustrated at $t_{fin}$). Consequently, we can clearly observe the horizontal drift caused by the change of angular position. Figure 9( f) presents the bottom-up view in the $x\unicode{x2013}z$ plane, highlighting the shift in stable angular position from the SV regime to the UV regime. It is important to highlight that, while a vertex generally faces downward, the angular position of the particle in the UV regime undergoes unsteady dynamic evolution. Once the stable angular position is established, the icosahedron settles following a straight vertical path with minor lateral fluctuations.

The transition from the SV regime to the UV regime stems from the interaction between the particle and vortex pairs in the particle wake. For flow past an angular particle, a pair of opposite-signed vortices is generated from each inclined face of the particle (Gai & Wachs 2023d). When the particle settling velocity is small, we observe minuscule vortical structures on the streamwise vorticity $\omega _y$ isosurfaces in the SV regime in figure 10. At the stable angular position, these vortical petals are relatively small compared with the particle size, exerting negligible effects on the particle motion, particularly in the horizontal plane, hence leading to a SV path. As the $\mathcal {G}a$ increases, the size of these induced vortex pairs expands, extending their tails into the particle wake region in the UV regime. The length of these vorticity isosurface tails, reminiscent of anemone tentacles, fluctuates over time, introducing unsteady drift perturbations to the particle motion.

Figure 10. Streamwise vorticity isosurfaces $\omega _y = \pm 0.03$ for freely settling Platonic polyhedrons in the SV and the UV regime.

4.2.2. Oblique settling: horizontal drift

As $\mathcal {G}a$ increases beyond the UV regime, we observe two distinct regimes for the settling Platonic polyhedrons: oblique regimes (SO and UO) for particles with high $\phi$ (octahedron, dodecahedron and icosahedron) and a HS regime for particles with low $\phi$ (tetrahedron, cube). We first examine the two oblique regimes (SO and UO), which have also been observed in freely settling sphere and disk cases, as reported in Jenny et al. (2004) and Bi et al. (2022).

Figure 11 illustrates the settling paths of the Platonic polyhedrons (octahedron, dodecahedron, icosahedron) in the SO and UO regimes. In this paper we categorize a particle as being in the SO regime when it exhibits an inclined straight settling path and its horizontal displacement is at least $1$ at $y = -300$. As $\mathcal {G}a$ escalates, low-frequency oscillations begin to appear in the horizontal plane, marking the transition to the UO regime.

Figure 11. Settling paths for Platonic polyhedrons in the SO and UO regimes: (a–c) octa, (d–f) dodeca and (g–i) icosa. The particle initial angular position is illustrated at $t_{ini}$, with the final sable angular position shown at $t_{fin}$ in the respective views.

Figures 11(a)–11(c) illustrate the settling path of an octahedron in both the SO and UO regimes in the $x\unicode{x2013}y$, $z\unicode{x2013}y$ and $x\unicode{x2013}z$ planes. At $\mathcal {G}a=120$, the horizontal displacement of the octahedron exceeds $15$ at $y=-600$ in the $x\unicode{x2013}y$ plane. Initially, the octahedron is released at a vertex-down angular position, as depicted at $t_{ini}$. In both the SO and UO regimes, the octahedron eventually turns to a face-down stable angular position, as illustrated at $t_{fin}$. Notably, at $\mathcal {G}a=90$, this change in angular position leads to a transient shift of the path to the right ($x>0$) as shown in figure 11(a). Upon the establishment of a stable face-down angular position, the particle undergoes a drift towards the left ($x<0$), a motion attributable to the distribution of vortex pairs on the face-down octahedron. Indeed, the very existence of the SO regime underscores the presence of a stable angular position, given the significant influence of the angular position on the particle lateral motion. The transition from the SO regime to the UO regime is observable in the range from $\mathcal {G}a=120$ to $\mathcal {G}a=130$. Considering that the stable angular positions in the UO regime remain face down but dynamically evolve over time, the low-frequency fluctuation in settling path is attributed to a combination of the dynamic evolution of the particle angular position and the ensuing instability of the vortex structures.

In the case of a settling dodecahedron, the horizontal displacement is observed to be around $5$ at $y=-400$ in the $z-y$ view for a particle at $\mathcal {G}a=150$, but it increases to $25$ at $\mathcal {G}a=170$. The inclination angle of the oblique path in the $z-y$ view sharply increases with $\mathcal {G}a$, and is found to be larger than that in the octahedron cases. For the icosahedron, the oblique settling path is curvilinear, as can be observed from figures 11(g)–11(h). Despite this, we still categorize these paths within the SO regime due to their substantial horizontal displacement and steady path. In both the dodecahedron and icosahedron cases, the particle is initially released with an edge-down angular position and subsequently rotates to a vertex-down position, of which the effect can be noticed in the initial stages of the particle settling path. The onset of the UO regime can be identified between $170 \leqslant \mathcal {G}a \leq 180$ for the dodecahedron and $180 \leqslant \mathcal {G}a \leqslant 190$ for the icosahedron, although it exhibits a very low fluctuation frequency at $\mathcal {G}a = 180$ in the icosahedron case.

The oblique settling observed in both the SO and UO regimes is driven by the emergence of steady/unsteady lateral forces, which stem from the vortex pairs generated on the particle surface. Given that low-pressure zones typically exist within the core of these vortices (Jeong & Hussain 1995), the symmetry of vortex pairs on the particle surface can significantly influence the pressure distribution, thereby directly affecting the total lateral forces exerted on the angular particle. Figure 12 provides a depiction of the generation and distribution of the streamwise vorticity, $\omega _y$, on the surface of an octahedron in both the SO and UO regimes. To better illuminate its evolution with increasing $\mathcal {G}a$, a higher value for $\omega _y$ isosurface is chosen in figure 12 than that in figure 10. The low-pressure zone on a horizontal plane ($x\unicode{x2013}z$) cutting the octahedron settling at $\mathcal {G}a=100$ is illustrated in yellow with isocontours on the left half of figure 12. The following two key observations can be discerned.

1. (i) The low-pressure zone (yellow and black) coincides with the region of high vorticity and exhibits planar symmetry.

2. (ii) The high pressure (white) on the right-top face drives the octahedron towards the left-bottom low-pressure region, elucidating the direction of the horizontal drift force $F_h$ and explaining the oblique settling path.

Figure 12. Streamwise vorticity $\omega _y$ and distribution of the low-pressure zone on an octahedron in the $x\unicode{x2013}z$ plane in the SO and the UO regime.

Especially noteworthy are the vortex pair distributions in the cases where $100 \leqslant \mathcal {G}a \leqslant 180$ for the staggered SO/UO regimes in figure 8 for a freely settling octahedron, as depicted on the right half of figure 12. The vortex pairs grow with $\mathcal {G}a$ and extend into the particle wake. When compared with $\mathcal {G}a=100$, the disparity in intensity and dynamic movement of the vortex pairs instigates the onset of unsteadiness in the oblique settling observed at $\mathcal {G}a=150$. At $\mathcal {G}a=170$, the merging of vortices results in the formation of two long vortex tails that counterbalance the horizontal fluctuations and reinstate the SO settling. We will also see the appearance of two vortical tails visualized using the $\lambda _2$ criterion in the following discussion. At $\mathcal {G}a=180$, the onset of vortex shedding induces UO settling, leading to the reappearance of the UO regime.

Further elucidation of the vorticity patterns in the wake of an oblique settling particle is provided in figure 13. First in the octahedron case, we observe the unsteadiness introduced by vortex merging in the case of the octahedron between $\mathcal {G}a=120$ and $\mathcal {G}a=130$ (figure 13a,b), and the subsequent re-establishment of steadiness by two long vortex tails at $\mathcal {G}a=170$, as shown in figure 13(c). Figures 13(e)–13(g) illustrate that the wake vorticity patterns for freely settling Platonic polyhedrons in the SO regime exhibit planar symmetry. Though the wake vorticity pattern of a sphere presents centrosymmetry, there appears to be a slight imbalance in vorticity intensity; the top-right (left-bottom) region displays a higher intensity than the top-left (bottom-right) region. In essence, the oblique settling path and its stability can be directly attributed to the non-zero lateral force, as a result of the planar symmetric vortex pairs generated on the particle surface.

Figure 13. Planar symmetry of the wake region in streamwise vorticity ($\omega _y >0$ in red and $\omega _y < 0$ in blue) in the $x\unicode{x2013}z$ plane located at $y = y_p + 1.5$ in the SO and the UO regime; (a) octa, $\mathcal {G}a = 120$ (orange shade), (b) octa, $\mathcal {G}a = 130$ (brown shade), (c) octa, $\mathcal {G}a = 170$ (orange shade) and (d) octa, $\mathcal {G}a = 180$ (brown shade), (e) dodeca, $\mathcal {G}a = 150$ (orange shade), ( f) icosa, $\mathcal {G}a = 130$ (orange shade), (g) sphere, $\mathcal {G}a = 180$ (orange shade).

4.2.3. Helical settling: spiraling drift

The helical motion of freely settling particles has been investigated for particles of various shapes, including spheres (Zhou & Dušek 2015), spheroids (Mougin & Magnaudet 2001), oblate ellipsoidal bubbles (Mougin & Magnaudet 2006), short cylinders and disks (Ern et al. 2012). Figure 14 showcases the settling path of a tetrahedron and a cube in the HS regime in the $x\unicode{x2013}y$, $z\unicode{x2013}y$ and $x\unicode{x2013}z$ planes. The tetrahedron at $\mathcal {G}a=80$ exhibits a circular path in the $x\unicode{x2013}z$ plane, and the helix pitch progressively grows with $\mathcal {G}a$. This circular helical path, however, is only discernible for $80 \leqslant \mathcal {G}a \leqslant 90$. As $\mathcal {G}a$ further increases, the helical path stretches horizontally, adopting an oval shape in the $x\unicode{x2013}z$ plane. Concurrently, the axis of the helix deviates from its parallel alignment with the $y$ axis. As depicted in figure 14(c), the tetrahedron at $\mathcal {G}a=120$ exhibits a markedly less ordered settling path compared with the case at $\mathcal {G}a = 90$. In accordance with the SV and UV regimes, the settling tetrahedron predominantly maintains a face-down angular position. Yet the leading face increasingly inclines relative to the horizontal $x\unicode{x2013}z$ plane as $\mathcal {G}a$ rises. For a better clarity, the temporal evolution of the tetrahedron angular position in the HS regime is depicted in figure 15(a). At $\mathcal {G}a=60$, the shift from the UV regime to the HS regime is noticeable, marked by the minimal helix radius that imparts the impression of near-vertical settling. Meanwhile at $\mathcal {G}a=90$, a circular helical path is distinctly noticeable in the horizontal $x\unicode{x2013}z$ plane. Figure 15(a) reveals the counter-clockwise spiral motion of the tetrahedron at $\mathcal {G}a=90 \sim 140$, observed from the bottom-up view in the $x\unicode{x2013}z$ plane.

Figure 14. Settling paths for Platonic polyhedrons in the HS regime: (a–c) tetrahedron, (d–f) cube. The particle's initial angular position is illustrated at $t_{ini}$, with the final angular position shown at $t_{fin}$ in the respective views.

Figure 15. Temporal evolution of the particle angular position (from purple to beige) in the HS regime; (a) tetrahedron, (b) cube.

In figures 14(d)–14(e) we present the paths of a freely settling cube in the HS regime at $\mathcal {G}a = 170, 180, 190$. The range of $\mathcal {G}a$ for the HS cube is narrower than that of the tetrahedron. At $\mathcal {G}a=170$, the cube lateral motion is significantly less pronounced than at $\mathcal {G}a=180$. With the transition from the UV to the HS regime, the cube at $\mathcal {G}a=170$ maintains a face-down angular position, yet it displays a helical path with a minor magnitude of horizontal spiral motion. In the HS regime $\mathcal {G}a=180 \sim 190$, a shift in stable angular position occurs. A vertex-down angular position replaces the face-down position in the HS regime, as illustrated in figure 15(b). The settling path initially exhibits a zigzag pattern, followed by the emergence of a circular helical path at $\mathcal {G}a=180$, as illustrated in figure 14(e). At $\mathcal {G}a=190$, the helix axis laterally drifts up to a distance of $7.5$ at $y=-600$.

Similar to figure 12, the streamwise vorticity generation and the low-pressure zone surrounding a settling tetrahedron is illustrated in figure 16. In the left part of figure 12, the generation of $\omega _y$ on the tetrahedron at $\mathcal {G}a=60, 80, 140$ is depicted. Notably, the face-down angular position in the HS regime tilts with respect to the horizontal plane (highlighted by the black line), leading to asymmetric vortex pair generation and the resulting heterogeneous pressure distribution. The right part of figure 16 provides a side view of the settling tetrahedron at $\mathcal {G}a=140$ to highlight the disparity in the vortex structures in the particle wake region. The merging of two blue vortices in the wake forms a larger blue rudder vortex in the wake, creating a more intense low-pressure region. This rudder vortex alters the face-down angular position by elevating one vertex of the bottom face, effectively steering the autorotation of the tetrahedron.

Figure 16. Streamwise vorticity $\omega _y$ and distribution of the low-pressure zone on a tetrahedron in the HS regime.

A closer inspection and visualization of wake structures of angular particles in the HS regime reveals two distinct types of wake vortices. The tetrahedron demonstrates a twisted double-helix wake structure, while the cube exhibits a vortex-shedding pattern comprising several hairpin vortices in each cycle.

The high angularity of the tetrahedron results in a helical motion starting from $\mathcal {G}a = 60$. At this value of $\mathcal {G}a$, only a few vortical fragments can be discerned in the wake region of the particle, as identified by the $\lambda _2$ criterion ($\lambda _2=-0.1$) in figure 17. Similar to the SO regime, an increase in $\mathcal {G}a$ results in the emergence of two elongated vortical tails behind the particle. The autorotation of the tetrahedron twirls the vortical tails to form two intertwined helices at $\mathcal {G}a = 80$ as depicted in figure 17. This double-threaded wake continues to expand with increasing $\mathcal {G}a$ and the elongated vortex ring structures begin to manifest at $\mathcal {G}a = 130$. Here, it is evident that a segment of the double-threaded tail has snapped off and starts to morph into a hairpin structure. By contrast, we depict the wake structure of a settling cube at $\mathcal {G}a=190$ in the last panel of figure 17. Clearly, even in the presence of vortex shedding, the settling cube can sustain a HS path.

Figure 17. Wake structures identified by $\lambda _2 = -0.1$ for the settling tetrahedron and cube in the HS regime on the $x\unicode{x2013}y$ (left) and $z\unicode{x2013}y$ plane (right).

These two unique wake structures stem from the pronounced particle streamwise autorotation induced by the merging of wake vortices, which in turn affect the particle motion. The distinction lies in the fact that the double-helix structure in the case of the tetrahedron primarily results from the influence of shape, and manifests before the onset of vortex shedding. By contrast, the HS of the cube is propelled by the shedding of hairpin vortices. Compared with the other three particles (octahedron, dodecahedron, icosahedron), the particle angularity emerges as a pivotal factor amplifying the autorotation. We delve deeper into the statistics of particle rotation in the following sections.

4.2.4. Chaotic settling: disordered motion

In the previously discussed regimes, once the flow regime stabilizes after the transition period, the settling motion of the particle becomes predictable. As the $\mathcal {G}a$ increases above a certain threshold, the settling path for all solid particles in the fluid evolves in a chaotic and unpredictable manner, characterized by random fluctuations in the particle settling path. Figure 18 illustrates the CS paths for a tetrahedron and a dodecahedron at high $\mathcal {G}a$. For brevity, we focus on these two Platonic polyhedrons as examples of the particle in the helical and oblique settling regimes, respectively. When compared with the HS regime (figure 14b), the chaotic regime exhibits amplified lateral displacements in the case of the tetrahedron, as shown in figure 18(a). The unpredictability also manifests in the wake structures, where no stable vortex structures can be identified. At $\mathcal {G}a=160$, the settling path of the tetrahedron becomes unorganized due to the intermittent shedding of vortex rings from the particle wake. The tetrahedron at $\mathcal {G}a = 200$ exhibits sharp direction changes in its settling path as shown in figure 18(a). Similar CS paths are observed for the dodecahedron in figure 18(c), although the frequency of fluctuation is reduced, likely due to its high $\phi$. Interestingly, the horizontal drift in the CS regime (figure 18c) is less pronounced than in the SO and UO regimes (depicted in figure 11e). A prominent characteristic of the chaotic regime is the significant particle angular velocity, which disrupts the ability of the particle to maintain a regular drift force.

Figure 18. Path instability for settling of a Platonic polyhedron in the CS regime; (a,b) tetra, (c,d) dodeca.

4.3. Wake vortex structures

In this section we explore the wake vortex dynamics during regime transitions for each of the five Platonic polyhedrons. Figure 19 displays vortex structures in the wake region for a freely settling Platonic polyhedron at various $\mathcal {G}a$. Using the $\lambda _2=-0.1$ criterion for isosurface identification, we present a snapshot taken once the settling regime is fully established. Figure 19(a) displays the mesh grid on the $x\unicode{x2013}y$ plane, intersecting the centre of the settling cubes, highlighting well-resolved vortex structures. For a better clarity, we omit the mesh grid in the panels of the other Platonic polyhedrons in figure 19.

Figure 19. Instantaneous vortical structures for a freely settling Platonic polyhedron at increasing $\mathcal {G}a$ in the $x\unicode{x2013}y$ plane; isosurface identified by $\lambda _2 = - 0.1$ and coloured by the streamwise velocity $u_y$. (a) Cube, $\phi = 0.806$. (b) Tetra, $\phi = 0.67$. (c) Octa, $\phi = 0.846$. (d) Dodeca, $\phi = 0.91$. (e) Icosa, $\phi = 0.94$.

At low $\mathcal {G}a$, all particles adopt a SV settling path, yielding no discernible vortex structures in the wake region. Some vortical fragments, originating from the particle rear surface, can be observed in the UV (e.g. $\mathcal {G}a=100$ in the cube case in figure 19a) and SO regimes (e.g. $\mathcal {G}a=120$ in the dodecahedron case in figure 19d). In these low $\mathcal {G}a$ regimes the isosurfaces of the streamwise vorticity $\omega _y$ offer a more suitable approach to visualizing wake structures as discussed in the previous sections.

As $\mathcal {G}a$ increases, we observe two vortex tails extending to the wake for all five Platonic polyhedrons, an observation consistent with prior studies on a fixed or settling sphere (Johnson & Patel 1999; Jenny et al. 2004). For the Platonic polyhedrons with higher $\phi$ (dodecahedron and icosahedron), two vortical tails remain straight and elongated ($\mathcal {G}a=160$ in figure 19(d) and $\mathcal {G}a=170$ in figure 19e), akin to the case of a settling sphere. In the case of the octahedron, we identify a third, comparatively faint tail in the wake region. The third tail breaks the planar symmetry and contributes to the wide range of Galileo numbers ($130 \leqslant \mathcal {G}a \leqslant 200$) observed in the UO regime. As shown in figure 19(c), the asymmetry of the vortex structures propels the particle off its original centreline, and low-frequency vortex shedding commences at $\mathcal {G}a = 160$ and $\mathcal {G}a = 190$. By contrast, the rotation of highly angular particles significantly influences the particle settling path. For settling tetrahedrons, the two vortical tails are no longer straight, instead intertwining to form a double-helix structure as seen in $\mathcal {G}a=120$ (figure 19b). This structure visibly alters the particle settling path in the HS regime.

As $\mathcal {G}a$ continues to rise, the double vortical tail structure begins to destabilize, leading to the onset of vortex shedding. The critical $\mathcal {G}a$ for the onset of vortex shedding tends to increase in accordance with $\phi$. Given that all Platonic polyhedrons rotate during their descent, there is no fixed plane for vortex shedding. Both the orientation and magnitude of the angular velocity vector vary over time for a freely settling Platonic polyhedron. Combined with the effects of vortex shedding, this typically results in irregular particle motion. As can be seen from figure 8, the critical $\mathcal {G}a$ for transitioning to CS paths correlates positively with $\phi$. As $\mathcal {G}a$ increases, the frequency of vortex shedding also rises. For example, instead of two rings shed during one period (the duration required for the particle to complete a single circular trajectory in the horizontal plane) of vortex shedding, we observe five vortex rings within a single period for a freely settling cube at $\mathcal {G}a = 180$ in the HS regime. At even higher $\mathcal {G}a$, such as $\mathcal {G}a = 300$ in figure 19(a), the number of rings shed per period can further increase.

4.4. Vertical settling characteristics

In the following we offer an in-depth discussion of particle motion along both vertical and horizontal directions. We aim to provide a comprehensive understanding of the particle dynamics, force/torque balance and highlight the intrinsic differences between the oblique regime and helical regime. In the vertical direction, two empirical models are drawn for the settling ${\mathcal {R}e}$ and the disturbed wake length $L_w$.

4.4.1. Settling drag $\overline {C_d}$ and Reynolds number ${\mathcal {R}e}$

Figure 20 depicts the evolution of the time-averaged drag coefficient $\overline {C_d}$ (along the $y$ direction) as a function of ${\mathcal {R}e}$, in the range $10\leqslant {\mathcal {R}e} \leqslant 300$, for freely falling Platonic polyhedrons with $m=2$. Please note that ${\mathcal {R}e}$ is computed with the time-averaged settling velocity $\overline {U_{m}}$. The empirical correlation of Haider & Levenspiel (1989) is also plotted in red as a reference curve. In general, our numerical simulation results align well with Haider's correlation for highly spherical particles such as the octahedrons in figure 20(c), dodecahedrons in figure 20(d) and icosahedrons in figure 20(e). As ${\mathcal {R}e}$ increases, discrepancies between the two curves become slightly more marked, but remain relatively small. For particles with higher angularity (lower $\phi$), the discrepancy becomes more visible. This is particularly noticeable for freely settling cubes, where an enhancement of $\overline {C_d}$ in the HS regime (denoted by red squares) is observed in figure 20(b), a phenomenon that deviates from the empirical correlation prediction. For a settling angular particle, the drag force inherently relies on the angular position of the particle with respect to the flow. In the previous work of our group (Seyed-Ahmadi & Wachs 2019), the drag enhancement was noted and the surge of the $\overline {C_d}$ was attributed to the collective effect of the alteration in its angular position and the misalignment between the threads of vorticity and the instantaneous direction of motion. As illustrated in figure 15(b), the settling cube maintains a face-down angular position at $\mathcal {G}a=170$ (${\mathcal {R}e}=210$) before transitioning to a vertex-down position at $\mathcal {G}a=180$, which coincides with an increase in $\overline {C_d}$ and a decrease in the Reynolds number to ${\mathcal {R}e}=182$. By contrast, the freely settling tetrahedron maintains a face-down angular position throughout the whole HS regime. However, despite the tetrahedron instantaneous direction of motion not being aligned with the vortex threads, as shown in figures 16 and 17, there is no abrupt increase in the $\overline {C_d}$ curve in figure 20(a). This indicates that the enhancement of $\overline {C_d}$ is not directly linked to the misalignment of wake vortices in the HS regime. Consequently, it can be stated more precisely that the shift in the stable angular position, specifically in the case of a settling cube, is the primary cause of the observed $\overline {C_d}$ enhancement.

Figure 20. Time-averaged drag coefficient $\overline {C_d}$ of a freely settling Platonic polyhedron with density ratio $m=2$ as a function of ${\mathcal {R}e}$; comparison between the numerical results of DLM/FD (- - - blue dash line) and the correlation of Haider et al. (——, orange solid line) (Haider & Levenspiel 1989); marker colour in the SV ($\sqsubset \!\sqsupset$, filled black), UV ($\sqsubset \!\sqsupset$, filled orange), SO ($\sqsubset \!\sqsupset$, filled brown), UO ($\sqsubset \!\sqsupset$, filled maroon), HS ($\sqsubset \!\sqsupset$, filled red) and CS ($\sqsubset \!\sqsupset$, filled blue) regimes. (a) Tetra. (b) Cube. (c) Octa. (d) Dodeca. (e) Icosa.

Additionally, we present the evolution of $\overline {C_d}$ and ${\mathcal {R}e}$ of the particle, based on the time-averaged settling velocity $\overline {U_m}$, as a function of $\mathcal {G}a$ in figure 21. The Platonic polyhedrons in different regimes are depicted with markers of respective shapes and colours. For any given particle shape, $\overline {C_d}$ tends to decrease with increasing $\mathcal {G}a$. The tetrahedron, with its high angularity, exhibits the greatest $\overline {C_d}$, whereas the icosahedron displays the lowest $\overline {C_d}$, across the range $10 \leqslant \mathcal {G}a \leqslant 300$. Interestingly, from the UV regime to the HS regime, we notice a change of slope of $\overline {C_d}$ as a function of $\mathcal {G}a$, as highlighted by the yellow and red dashed lines. This slope change is mainly attributed to the augmentation of the particle rotation as a result of the Magnus effect and can also be noticed for the other particles during the transition from the SO regime to the UO regime. Typically, $\overline {C_d}$ diminishes with an increase in $\phi$. However, there are exceptions: within the range of $70 \leqslant \mathcal {G}a \leqslant 170$, the cube, due to its face-down stable angular position, has a lower $\overline {C_d}$ than the octahedron. A notable enhancement in drag for the cube begins at $\mathcal {G}a = 180$ in the HS regime. For the other Platonic polyhedrons, the changes in $\overline {C_d}$ are more gradual, either decreasing very slowly or showing a minor increase in the CS regime, as indicated by the sky-blue markers.

Figure 21. (a) Time-averaged drag coefficient $\overline {C_d}$ of a freely settling Platonic polyhedron as a function of Galileo number in the range $10 \leqslant \mathcal {G}a \leqslant 300$; (b) settling ${\mathcal {R}e}$ as a function of the $\mathcal {G}a$ for Platonic polyhedrons: tetra ($\blacktriangle$), cube ($\blacksquare$), octa ($\blacklozenge$), dodeca (, filled black) and icosa (, filled black); colour at regimes SV ($\sqsubset \!\sqsupset$, filled black), UV ($\sqsubset \!\sqsupset$, filled orange), SO ($\sqsubset \!\sqsupset$, filled brown), UO ($\sqsubset \!\sqsupset$, filled maroon), HS ($\sqsubset \!\sqsupset$, filled red) and CS ($\sqsubset \!\sqsupset$, filled blue) regimes; (c) correlation prediction $\widetilde {{\mathcal {R}e}}$ versus numerical results ${\mathcal {R}e}$.

In figure 21(b) we observe that for a specific particle shape, ${\mathcal {R}e}$ typically increases linearly with $\mathcal {G}a$, with the notable exception being the significant drag enhancement of the cube in the HS regime. Analysing this linear relationship between ${\mathcal {R}e}$ and $\mathcal {G}a$ reveals that the slope of the ${\mathcal {R}e}-\mathcal {G}a$ curve is also a linear function of $\phi$. Based on our numerical observations, we can formulate the following empirical correlation in the range of $30 \leqslant \mathcal {G}a \leqslant 300$:

(4.1)\begin{equation} \widetilde{{\mathcal{R}e}} = 2.57\phi \mathcal{G}a - 0.92\mathcal{G}a - 0.032 \exp(7.50\phi). \end{equation}

Here $\widetilde {{\mathcal {R}e}}$ is an estimate of the settling Reynolds number. With an average relative error under $3\,\%$, the correlation solely relies on $\phi$ and physical properties of the particle–fluid system, represented by $\mathcal {G}a$. We plot the predicted values of $\widetilde {{\mathcal {R}e}}$ as red lines in figure 21(b). A strong correspondence is observed for all Platonic polyhedrons, with a minor exception being the settling cube case. In addition, the accuracy of this model is also depicted in figure 21(c), where the predicted $\widetilde {{\mathcal {R}e}}$ from (4.1) is plotted against the observed ${\mathcal {R}e}$ in our numerical results. Impressively, the correlation yields a determination coefficient of $r^2=0.9943$, providing a satisfactory fit to the numerical data. Please note that due to limited data, the current form of the correlation is not applicable for $\mathcal {G}a < 30$, particularly for particles with high $\phi$. In this specific range of $\mathcal {G}a$, the effect of particle shape $\phi$ is less dominant. Thus, the relation between $\mathcal {G}a$ and ${\mathcal {R}e}$ is expected to deviate from the proposed correlation in (4.1).

4.4.2. Disturbed wake length: decay of flow velocity

In particle-laden flows, interactions between dispersed particles are significant, except in some extreme dilute systems. When a particle settles in a fluid, it entrains the fluid in the wake region, which influences the motion of trailing particles. In order to quantify the effects of particle shape on wake fluid dynamics, and potentially on other particles around it, we compute the particle disturbed wake length $L_w$ as a function of the disturbed velocity $u_d$, defined as

(4.2)\begin{equation} u_d = \frac{u_{crit}}{\overline{U_{m}}}, \end{equation}

where $u_{crit}$ is the threshold fluid velocity. Fluid motion with velocities below this threshold is deemed outside the particle wake and can be neglected. The evolution of $L_w$ as a function of $u_d$ is illustrated in figure 22(a). Please note that a reduced $u_d$ corresponds to an extended disturbed wake length $L_w$.

Figure 22. (a) Disturbed wake length $L_w$ as a function of $u_d$; (b) disturbed wake volume $V_w$ as a function of $u_d$ for settling Platonic polyhedrons: tetra ($\blacktriangle$), cube ($\blacksquare$, filled dark blue), octa ($\blacklozenge$, filled light blue), dodeca (, filled yellow) and icosa (, filled orange).

Our analysis shows that the computed $L_w$ aligns remarkably well with the power law

(4.3)\begin{equation} L_w = 1.65 u_d ^{{-}1.205} \end{equation}

with a determination coefficient $r^2 = 0.96$. In figure 22(a) the master line, represented in red, underscores that the decay of velocity in the wake region adheres to a power law and remains unaffected by the specific angular shape of the particle. This behaviour aligns with expectations, given that fluid motion diffusion primarily hinges on intrinsic fluid properties, including density and viscosity. Please note that the existence of the master line does not diminish the influence of particle shape on $L_w$. Instead, the settling velocity $\overline {U_m}$, which factors into the definition of $u_d$, is influenced by the particle sphericity $\phi$ as indicated in (4.1).

In addition, we compute the volume of the disturbed wake region $V_w$ by integrating the cell volumes with velocity larger than $u_d$ and plot its evolution as a function of $u_d$ in figure 22(b). Similar to the $L_w$, the $V_w$ decreases as $u_d$ rises. The influence of particle shape is distinctly more pronounced on $V_w$ than on $L_w$. Particles with higher angularity (tetrahedron, cube, octahedron) tend to have a larger value of $V_w$ compared with a particle with high sphericity (dodecahedron, icosahedron). Through multi-variable regression, we establish the following empirical law for $V_w$:

(4.4)\begin{equation} V_w = 1.02 u_d ^{{-}1.93} \end{equation}

with the determination coefficient $r^2=0.9484$. From figure 22(b), we see that the influence of $\mathcal {G}a$ and $\phi$ on $V_w$ becomes more evident at a larger disturbed velocity ($u_d > 0.6$). Given that the disturbed wake region typically resembles a long spindle, with its transverse dimension being significantly smaller than its streamwise length $L_w$, it is not surprising that both (4.4) and (4.3) follow a similar power law relation.

In summary, knowing the physical properties of the particle–fluid system, we can estimate the settling velocity $\overline {U_m}$ (derived from the settling $\widetilde {{\mathcal {R}e}}$) and subsequently $L_w$ (given a targeted $u_d$) by combining (4.1) and (4.3). Knowledge of the approximative values of $\overline {U_m}$ and $L_w$ can significantly aid in the large-scale numerical modelling of dilute angular particle suspensions.

4.5. Horizontal drift characteristics

As discussed above, the lateral force contributes to the emergence of two transitional regimes: oblique (SO, UO) and HS, as illustrated in figure 8. In this section we try to understand the underlying similarities and differences between the oblique regime and the helical regime.

4.5.1. A comparison between fixed and freely settling particles

Figure 23 exhibits the regime transition map as a function of ${\mathcal {R}e}$ for the flow past a fixed Platonic polyhedron (Gai & Wachs 2023b). Overlaying this map, we indicate the stable angular position and regime of the freely settling Platonic polyhedrons corresponding to the same settling ${\mathcal {R}e}$ using distinct coloured markers. This juxtaposition aims to shed light on the interplay between fixed-particle and freely moving particle scenarios, especially in the oblique and HS regimes.

Figure 23. Regime transition map for flow past a fixed Platonic polyhedron in three angular positions: an edge (E), a face (F) and a vertex (V) facing the flow; regimes in fixed particle case: multi-planar symmetry ($\sqsubset \!\sqsupset$, filled peach), planar symmetry ($\sqsubset \!\sqsupset$, filled light blue), hairpin vortex shedding ($\sqsubset \!\sqsupset$, filled light peach) and chaotic vortex shedding ($\sqsubset \!\sqsupset$, filled dark blue); double-symmetric vortex shedding ($\sqsubset \!\sqsupset$, filled very light blue); overlaying the map, angular position of the freely settling Platonic polyhedron at the same settling ${\mathcal {R}e}$: tetra ($\blacktriangle$), cube ($\blacksquare$), octa ($\blacklozenge$), dodeca (, filled black) and icosa (, filled black); regimes in the freely settling case: SV ($\sqsubset \!\sqsupset$, filled black), UV ($\sqsubset \!\sqsupset$, filled orange), SO ($\sqsubset \!\sqsupset$, filled brown), UO ($\sqsubset \!\sqsupset$, filled maroon) and HS ($\sqsubset \!\sqsupset$, filled red).

The tetrahedron, when freely settling, consistently displays a single stable face-down angular position across regimes from the SV regime to the HS regime. On the other hand, the cube settles in a face-down position in the SV and UV regimes, but transitions to a vertex-down position in the HS regime. Figure 23 shows that the other three particles (octahedron, dodecahedron and icosahedron) likewise demonstrate shifts in their stable angular positions during regime transitions. Please note that in the freely settling particle case, the stable angular position is not perfect and may fluctuate in time.

Now we take a look at the interconnection between the regimes of a freely settling particle and its fixed counterpart. In the SO and UO regimes the settling icosahedron has a range of Reynolds number $206 \leqslant {\mathcal {R}e} \leqslant 238$. This coincides with the planar symmetric regime ($\sqsubset \!\sqsupset$, filled very light blue) observed in the fixed vertex icosahedron subject to the flow at the same ${\mathcal {R}e}$. This planar symmetric vortex structure results in an augmented lift force, driving the lateral motion of the particle if it is freely settling. Conversely, the oblique settling octahedron and dodecahedron exhibit Reynolds number ranges of $77 \leqslant {\mathcal {R}e} \leqslant 132$ and $123 \leqslant {\mathcal {R}e} \leqslant 213$, respectively. These values align with the late multi-planar symmetric regime ($\sqsubset \!\sqsupset$, filled peach) seen in the case of flow past a fixed particle. In the tetrahedron case the ${\mathcal {R}e}$ varies from ${\mathcal {R}e}=42$ to ${\mathcal {R}e}=117$ in the HS regime, while for the settling cube, it ranges from $182 \leqslant {\mathcal {R}e} \leqslant 207$. In the fixed particle case the face tetrahedron remains in the multi-planar symmetric regime at $42 \leqslant {\mathcal {R}e} \leqslant 117$. In the cube case the settling regime transition from the UV regime to the HS regime is accompanied by a shift from a face cube (multi-planar symmetric regime, $\sqsubset \!\sqsupset$, filled light peach) to a vertex cube (hairpin vortex-shedding regime, $\sqsubset \!\sqsupset$, filled very light peach) in the fixed particle case. This again underlines that the drag enhancement and wake structure evolution of a freely settling cube in the HS regime is mainly due to its change of stable angular position.

In summary, the oblique and helical regimes exhibit ranges of ${\mathcal {R}e}$ in the late multi-planar symmetry, planar symmetry and hairpin vortex-shedding regimes in the fixed particle case. The critical ${\mathcal {R}e}$ for the symmetry breakup in the freely settling particle wake usually has a smaller value than in the fixed particle case, primarily due to the fluid fluctuations. In fact, the lateral motion of the freely settling particle induces extra lift force components in the horizontal plane and an extra drag force in the vertical direction. It is hard to separate the lift force simply induced by the particle vertical motion from the lateral drag force originating from the horizontal motion of the particle. Yet, the results of the flow past a fixed angular particle in figure 23 can isolate the particle streamwise motion and subsequently provide helpful insights that aid understanding of the path instability. We present the critical ${\mathcal {R}e}$ corresponding to the onset of vortex shedding in the case of settling and fixed Platonic polyhedrons in Appendix A. A closer look at the lift force of a fixed angular particle in a flow is presented in Appendix B.

4.5.2. Drift hydrodynamic force

Compared with vertical regimes (SV, UV), the lateral motion of particles in the SO and HS regimes significantly influences the particle dynamics. The key distinction between oblique and HS regimes lies in their respective trajectories on the horizontal plane: a straight path characterizes the SO regime, while a spiral or circular path is indicative of the HS regime. In the SO regime, horizontal drift force guides the particle laterally, as demonstrated in figure 24(a). By contrast, the streamwise component of the particle torque, or the autorotation around the vertical axis, takes precedence in the HS regime. Consequently, the combined effects of the streamwise torque $T_y$ and horizontal drift force $F_h$ result in the particle following a HS path, as illustrated in figure 24(b). Therefore, the presence of two factors is crucial to ensure HS: substantial torque around the settling direction and a considerable, relatively steady lift force.

Figure 24. Schematics of the velocity $U$, vortical torque $T_y$ and horizontal force $F_h$; (a) oblique settling of an icosahedron, (b) HS of a tetrahedron.

Figure 25 displays the evolution of horizontal force $F_h$ during the free settling of the Platonic polyhedron in the SO regime (octahedron, dodecahedron, icosahedron) and the HS regimes (tetrahedron, cube) in the $x\unicode{x2013}z$ plane. The colour change from blue to red illustrates the temporal evolution of horizontal forces on the particle. In the oblique regimes the horizontal force exhibits a preferred direction during the transient acceleration phase. The elliptical shape of the blue dots in figure 25(a) suggests that the $z$ component of the hydrodynamic force, $F_z$, is more pronounced than its $x$-axis counterpart. As the particle settles, the yellow dots, exhibiting a smaller magnitude than the blue dots, imply a gradual balance of the horizontal force. Ultimately when the settling regime is fully established, we observe an even smaller magnitude of the clustering of the red dots at the centre of the plot. The horizontal drag force dynamically counterbalances the lift force induced by the vertical settling of the icosahedron. A similar process is observable for the dodecahedron in figure 25(b) and the octahedron in figure 25(c). The lateral force on the horizontal plane primarily aligns with the diagonal direction in figure 25(c) and the force fluctuations reduce once the lateral motion stabilizes.

Figure 25. Horizontal hydrodynamic force evolution in function of time, points coloured from blue to red as time going; (a–c) SO regime, (d–e) HS regime, ( f) settling icosahedron at $\mathcal {G}a=300$ in the CS regime, (g–i) vortex structures of a settling cube, tetrahedron and icosahedron at $\mathcal {G}a$ corresponding to (d–f). (a) Icosa, $\mathcal {G}a = 170$ (brown shade). (b) Dodeca, $\mathcal {G}a = 170$ (brown shade). (c) Octa, $\mathcal {G}a = 100$ (brown shade). (d) Cube, $\mathcal {G}a = 180$ (red shade). (e) Tetra, $\mathcal {G}a = 80$ (red shade). ( f) Icosa, $\mathcal {G}a = 300$ (blue shade). (g) Cube, $\mathcal {G}a = 180$ (red shade). (h) Tetra, $\mathcal {G}a = 80$ (red shade). (i) Icosa, $\mathcal {G}a = 300$ (blue shade).

In the HS regime $F_h$ serves as a centripetal force, enabling the particle to maintain circular motion in the horizontal plane. This results in an annular force distribution, as seen in figures 25(d) and 25(e), which is markedly different from the oblique settling regime. The horizontal force stems from the vortex structures in the particle wake: a double straight vortical tail in the SO and UO regimes, and an intertwined double-helix tail in the HS regime of the tetrahedron as depicted in figure 25(h). In the case of the settling cube, the lateral force distribution in figure 25(g) is determined by vortex shedding and autorotation of the cube around the vertical axis. In the last period (dark red dots) of shedding in figure 25(d), we can see five clear deviations of direction in the $F_x-F_z$ plot. Figure 25(g) reveals that five hairpin vortices are shed from the particle wake in each cycle. Therefore, the vortex force, originating from vortex shedding, primarily drives the particle lateral motion.

It is worth noting that during the CS regime of the icosahedron, a pseudo-helical path at $\mathcal {G}a=300$ is observed, as illustrated in figure 25( f). Correspondingly in figure 25(i), we see that the settling icosahedron at $\mathcal {G}a=300$ exhibits hairpin structures reminiscent of those seen in the settling cube in figure 25(g), albeit in a less ordered arrangement. Vortex shedding in the wake of the settling icosahedron provides the horizontal drift force $F_h$ and substantial streamwise torque $T_y$. A plausible explanation might be that the frequency of vortex shedding aligns with that of the particle spiral motion in the horizontal plane, resulting in a continuous HS path. We should point out that the helical motion of the icosahedron at $\mathcal {G}a=300$ does not contradict the oblique regime at lower ${\mathcal {R}e}$. In figure 8 we observe that in the SO regime, the vortex structure in the icosahedron wake has two long tails, and no vortex shedding occurs. The surge in streamwise torque $T_y$ at high $\mathcal {G}a$ accounts for this late emergence of HS motion in the CS regime.

4.5.3. Streamwise torque: from oblique to helical

Apart from the horizontal drift, a second key factor in the HS regime is the particle streamwise rotation around the vertical $y$ axis. To understand this, we calculate the time-averaged streamwise torque $\overline {T_y}$ and plot its evolution as a function of $\mathcal {G}a$ for all Platonic polyhedrons in figure 26. In the SV, UV and SO regimes, $\overline {T_y}$ of all angular particles exhibit values close to zero. With the increase of $\mathcal {G}a$, the streamwise torque takes on high values starting from the HS regime ($\sqsubset \!\sqsupset$, filled red) for both the tetrahedron and the cube. However, the dynamics differ between these two shapes: the tetrahedron shows a comparatively high streamwise torque with a minimal lateral drift force, whereas the cube demonstrates a substantial drift force but a relatively small streamwise torque. This distinction results in the tetrahedron exhibiting a HS path at lower $\mathcal {G}a$ ($60 \sim 150$) without any vortex shedding. The cube, on the other hand, exhibits a HS path at higher $\mathcal {G}a$ ($180 \sim 200$), with pronounced vortex shedding in the particle wake region. Given the small $\overline {T_y}$, the circular path radius in the horizontal plane in the cube case is larger than that in the tetrahedron case as shown in figures 25(g) and 25(h). At even higher $\mathcal {G}a$, in the UO and CS regimes, even though $\overline {T_y}$ increases to a non-zero value, the settling path cannot follow a helical pattern since the drift force is no longer stable.

Figure 26. Time-average streamwise torque $\overline {T_y}$ as a function of the $\mathcal {G}a$ for freely settling Platonic polyhedrons: tetra ($\blacktriangle$), cube ($\blacksquare$), octa ($\blacklozenge$), dodeca (, filled black) and icosa (, filled black); colour at regimes SV ($\sqsubset \!\sqsupset$, filled black), UV ($\sqsubset \!\sqsupset$, filled orange), SO ($\sqsubset \!\sqsupset$, filled brown), UO ($\sqsubset \!\sqsupset$, filled maroon), HS ($\sqsubset \!\sqsupset$, filled red) and CS ($\sqsubset \!\sqsupset$, filled blue).

Consequently, for an angular particle to maintain a HS regime, it requires a combination of: (i) a steady lateral force and (ii) a significant streamwise torque. In the SV and UV regimes, the horizontal drift and streamwise torque are both of minor importance. At higher $\mathcal {G}a$, the absence of either one of the two factors results in either the SO regime, characterized by a stable drift force and negligible streamwise torque, or the UO regime, marked by unstable drift force alongside substantial streamwise torque. Thus, both the oblique and helical regimes act as a transitional phase that bridges the vertical regime and the CS regime.

4.5.4. Horizontal force balance in the HS regime

In scenarios where a particle is fixed in a steady flow, or when a particle moves at a constant velocity, the hydrodynamic force arises solely from vortex forces. However, for freely moving angular particles, the impacts of particle acceleration and rotation become significant. Under acceleration, the force due to the inertia of the surrounding fluid increases, commonly referred to as the added mass of the particle. We adopt the force decomposition strategy proposed by Govardhan and Williamson (Govardhan & Williamson 2005; Seyed-Ahmadi & Wachs 2019), which separates the total hydrodynamic force into the sum of the potential force (added mass) $F_{AM}$ and a vortex force $F_{v}$. This strategy enables accurate estimation of the effects of vortical structures:

(4.5)\begin{equation} \boldsymbol{F_h} = \displaystyle{\sum \boldsymbol{F}} = \boldsymbol{F_{AM}} + \boldsymbol{F_{v}} ={-}\frac{M}{m}\frac{{\rm d}\boldsymbol{U}}{{\rm d}t}C_{AM} + \boldsymbol{F_{v}}. \end{equation}

Here ${{\rm d}\boldsymbol {U}}/{{\rm d}t}$ denotes the particle acceleration and $C_{AM}$ is the added-mass coefficient. In an axisymmetric flow around a sphere, one can deduce the apparent added-mass coefficient $C_{AM} = 0.5$ from the unsteady Bernoulli equation. In the existing literature, $C_{AM}$ values for angular particles are scarcely available. Blevins (2017) provided several values for the cube, including $C_{AM} = 0.64$, $0.67$ and $0.7$. However, to our best knowledge, there are no reported values of $C_{AM}$ pertaining to other Platonic polyhedrons.

When the angular particle is subject to high rotation and translation, especially when the rotation axis is perpendicular to the translation direction (settling along the $y$ axis), the Magnus force $\boldsymbol {F_M}$ makes an important contribution to the hydrodynamic force:

(4.6)\begin{equation} \boldsymbol{F_{v}} = \boldsymbol{F_M} + \boldsymbol{F'_{v}} = \frac{M}{m}(\boldsymbol{\varOmega} \times \boldsymbol{U})C_\varOmega + \boldsymbol{F'_{v}}. \end{equation}

Here $\boldsymbol {F'_{v}}$ is the remaining vortex force contribution when $\boldsymbol {\varOmega } \rightarrow 0$ and $C_\varOmega$ denotes the coefficient of the Magnus effect. Regarding the choice of $C_\varOmega$, Clanet (2015) recommends a Magnus coefficient $C_\varOmega = 0.4$ for spherical particles with $\varOmega < 0.5$, a value that holds true across various sports balls. However, literature on the Magnus coefficient for angular particles is limited. To obtain a clearer idea of the range of $C_\varOmega$ for the angular particle, we separately conducted numerical simulations of a rotating cube subject to an imposed angular velocity in an inertial flow. For a cube rotating around its opposite-face-centre axis at $\varOmega _{0} = 0.08$ and ${\mathcal {R}e} = 200$ (typical values for a settling cube in the HS regime), we found an approximate Magnus coefficient of $\widetilde {C_{\varOmega }} \approx 0.8$. Preliminary data analysis suggests that $\widetilde {C_{\varOmega }}$ varies significantly with ${\mathcal {R}e}$, the imposed angular velocity $\varOmega _{0}$ and the rotation axis. However, for a settling particle in the HS regime, these key factors determining $\widetilde {C_{\varOmega }}$ all dynamically evolve over time. Thus, for the settling cube in the HS regime, $\widetilde {C_{\varOmega }}=0.8$ cannot perfectly describe the Magnus force in (4.6).

Considering that the Magnus coefficient of a rotating angular Platonic particle has not been systematically studied, we use approximate coefficients to understand the qualitative impacts of particle rotation and added-mass forces on particle dynamics, particularly in the HS and CS regimes. For consistency, we choose the coefficients: $C_{AM} = 0.5$ and $C_\varOmega = 0.4$, to analyse the horizontal force balance of both the cube and tetrahedron in the HS and CS regimes. This choice is guided by two primary considerations: (i) the particle shape effect can be distinctly isolated and represented in the $\boldsymbol {F'_{v}}$ term in (4.6), applicable for both the tetrahedron and the cube; and (ii) use of estimated values for $C_\varOmega = 0.4 \sim 0.8$ and $C_{AM} = 0.5 \sim 0.7$ maintains the broad applicability and general nature of our discussion. Please note that this method is slightly different from that of our previous works on the settling of an angular particle of the same shape (Rahmani & Wachs 2014; Seyed-Ahmadi & Wachs 2019).

When the settling regime is fully established, as the particle transient acceleration ceases, the drag force counterbalances gravity in the vertical direction, with high $\mathcal {G}a$ fluctuations resulting from the evolution of vortical structures. Side drift becomes apparent in the SO regime, although the particle autorotation remains at a low level. As $\mathcal {G}a$ continues to increase, the particle rotation starts to intensify in the HS regime. The vortical tail remains straight in the SO and UO regimes, but transforms into a double helix in the HS regime due to the particle autorotation. To better understand the vortical forces, particularly the lateral forces exerted on the settling cube, we exhibit the force balance in the horizontal $x\unicode{x2013}z$ plane in figure 27. We apply a filter to the force data for a better visualization of force evolution as a function of time in the unsteady regimes.

Figure 27. (a–c) Force balance in the horizontal $x\unicode{x2013}z$ plane for the settling cube in the HS and CS regimes. (d–f) Settling path with lateral forces projected in the horizontal $x\unicode{x2013}z$ plane. Results are shown for (a) $\mathcal {G}a=180$ (red shade); (b) $\mathcal {G}a=190$ (red shade); (c) $\mathcal {G}a=300$ (blue shade); (d) $\mathcal {G}a=180$ (red shade); (e) $\mathcal {G}a=190$ (red shade); ( f) $\mathcal {G}a=300$ (blue shade).

At $\mathcal {G}a=180$, the magnitude of $\boldsymbol {F_M}$ is comparable to $\boldsymbol {F'_{v}}$ and it aligns in phase with the total force $\sum \boldsymbol {F}$ as shown in figure 27(a), indicating that the cube rotation has an influence on par with the shape effect. Its oscillation period expands when $t > 200$ $t_{ref}$, correlating with a shift in the settling path from two-dimensional zigzag to the three-dimensional HS. In the horizontal plane the cube follows a circular path, as denoted in figure 27(d). In this case, $\boldsymbol {F_M}$ provides the centripetal force for the helical path and remains the sole determinant of the settling direction. The added-mass force opposes the direction of the total force $\sum \boldsymbol {F}$, inhibiting particle acceleration and acting as the centrifugal force. There are still two settling stages observable at $\mathcal {G}a = 190$, as shown in figure 27(b). The helical axis moves to the $x^+$ direction in figure 27(e), where the angle between the Magnus force $\boldsymbol {F_M}$ and the particle motion direction diminishes. As a result, in addition to the centripetal role, $\boldsymbol {F_M}$ also acts as the driving force for the particle drift motion. At very high $\mathcal {G}a$, the Magnus effect becomes weaker than the shape effect and the vortex shedding. In figure 27(c) the remaining vortex force $\boldsymbol {F'_v}$ has a higher magnitude than that of $\boldsymbol {F_M}$, counterbalancing the driving force and introducing instability to the settling path.

Figure 28 showcases the force balance in the horizontal plane $x\unicode{x2013}z$ for the settling tetrahedron in the HS and CS regimes. Different from the settling cube, the magnitude of the Magnus force $\boldsymbol {F_M}$ is smaller than that of the remaining vortex force $\boldsymbol {F'_v}$ in the HS regime, as illustrated in figures 28(a)–28(d). In the HS regime, the oscillation period remains similar during the settling process. However, the time evolution of the forces loses its periodicity after the transition to the CS regime, as seen in figure 28(d). The radius of the circular path in the $x\unicode{x2013}z$ plane is relatively small in comparison to that of the cube, as shown in figures 28(e)–28( f). In the HS regime, $\boldsymbol {F_M}$ propels the particle horizontal motion, whereas the remaining vortex force $\boldsymbol {F'_v}$ mostly points to the centre of the circular path. This implies that the shape of the tetrahedron primarily contributes to the centripetal force. The horizontal lift force of the face-down tetrahedron gradually changes direction with respect to the streamwise autorotation, sustaining its helical motion. The pronounced $\boldsymbol {F'_v}$ can be attributed to the high angularity of the tetrahedron. As $\mathcal {G}a$ increases, the magnitude of $\boldsymbol {F_M}$ and $\boldsymbol {F'_v}$ shows greater variability in figure 28(c), the settling path can no longer remain circular, resulting in a variable radius for the helical path, as displayed in figure 28(g). As illustrated in figures 8 and 28(h), the tetrahedron settling becomes chaotic at $\mathcal {G}a \geqslant 160$, where the icosahedron can still sustain SO motion. Eventually, in the CS regime, $\boldsymbol {F_M}$ becomes almost tangential to the particle path. The random orientation of the total force $\sum \boldsymbol {F}$ reflects the chaotic motion of the particle in figure 28(h). In short, for the HS of a cube, the balance between $\boldsymbol {F_M}$ and $\boldsymbol {F'_v}$ is key, while in the case of the tetrahedron, the shape effect has a more pronounced role in contributing to the centripetal force. In the horizontal plane, the Magnus force has a component driving the direction of the particle motion in both the tetrahedron and the cube cases across the range of $\mathcal {G}a$ in this study.

Figure 28. (a–d) Force balance in the horizontal plane $x\unicode{x2013}z$ for the settling tetrahedron in the HS and CS regime. (e–h) Settling path with lateral forces in the horizontal $x\unicode{x2013}z$ plane; (a,e) $\mathcal {G}a=80$ (red shade), (b, f) $\mathcal {G}a=90$ (red shade), (c,g) $\mathcal {G}a=120$ (red shade) and (d,h) $\mathcal {G}a=200$ (blue shade).

4.6. Particle rotation statistics

Apart from the settling regime, various settling styles emerge due to the particle rotation such as fluttering or tumbling (Ern et al. 2012). In this section we explore and analyse the evolution of the particle angular velocity and the distinctive rotation styles of the settling Platonic polyhedrons at moderate $\mathcal {G}a$.

4.6.1. Angular velocity

We exhibit both the time-averaged particle angular velocity $\bar {\varOmega }$ and the maximum angular velocity $\varOmega _{max}$ as functions of $\mathcal {G}a$ in figure 29. As viscous effects diminish, the five Platonic polyhedrons exhibit a larger angular velocity when $\mathcal {G}a$ increases from $10$ to $300$. In this range of $\mathcal {G}a$, $\varOmega _{max}$ and $\bar {\varOmega }$ remain below $0.22$ and $0.11$, respectively, for all Platonic polyhedrons. Two distinct phases can be observed in figure 29(a): (i) when $\mathcal {G}a$ is small, $\bar {\varOmega }$ exhibits a plateau of values close to zero, the length of which increases with $\phi$; and (ii) as the transition to the HS and UO regimes occurs, $\bar {\varOmega }$ markedly increases.

Figure 29. (a) Time-averaged angular velocity $\bar {\varOmega }$, (b) maximum angular velocity $\varOmega _{max}$ in function of $\mathcal {G}a$ for all Platonic polyhedrons: tetra ($\blacktriangle$), cube ($\blacksquare$), octa ($\blacklozenge$), dodeca (, filled black) and icosa (, filled black); SV ($\sqsubset \!\sqsupset$, filled black), UV ($\sqsubset \!\sqsupset$, filled orange), SO ($\sqsubset \!\sqsupset$, filled brown), UO ($\sqsubset \!\sqsupset$, filled maroon), HS ($\sqsubset \!\sqsupset$, filled red) and CS ($\sqsubset \!\sqsupset$, filled blue) regimes.

For the tetrahedron, we observe an abrupt change in $\bar {\varOmega }$ at $\mathcal {G}a = 60$, marking the regime transition from the UV regime to the HS regime. Within the range $100 \leqslant \mathcal {G}a \leqslant 150$, the rate of increase for $\bar {\varOmega }$ reduces and $\bar {\varOmega }$ eventually diminishes from $\mathcal {G}a = 160$ to $\mathcal {G}a = 190$. The slowdown in the particle rotation is suggestive of the onset of vortex shedding. In comparison, $\bar {\varOmega }$ of a settling cube starts its sharp ascent at $\mathcal {G}a = 160$ in the HS regime. The rapid increase of the cube rotation is relatively short and ends at $\mathcal {G}a=180$. This range is significantly shorter than that seen in the tetrahedron case, aligning with the occurrence of the HS regime, as illustrated in figure 8. Once again, we observe a minor decrease of $\bar {\varOmega }$ when vortex shedding begins, followed by an increase in the CS regime. Here $\varOmega _{max}$ increases accordingly in the HS regime but we do not observe any minor decrease due to vortex shedding in figure 29(b). For particles with higher $\phi$, $\bar {\varOmega }$ increases in the range $140 \leqslant \mathcal {G}a \leqslant 180$ (octahedron, dodecahedron, icosahedron), corresponding to the transition from the SO regime to the UO regime. Furthermore, $\varOmega _{max}$ experiences an increase at $\mathcal {G}a \approx 100\sim 150$, which aligns with the transition from the UV regime to the SO regime.

In general, particles with lower $\phi$ tend to exhibit higher $\bar {\varOmega }$ and $\varOmega _{max}$ within the range of $10 \leqslant \mathcal {G}a \leqslant 300$, although this is not universally true. For instance, in figure 29(a), at $\mathcal {G}a > 180$, $\bar {\varOmega }$ of the cube either exceeds or closely matches that of the tetrahedron. Moreover, at $\mathcal {G}a >220$, the icosahedron displays higher values of both $\bar {\varOmega }$ and $\varOmega _{max}$ compared with those of the dodecahedron. It is worth noting that the moment of inertia $I$ that opposes particle rotation, decreases as $\phi$ of the Platonic polyhedrons rises. The fact that the icosahedron possesses a smaller $I$ than the dodecahedron might explain its higher angular velocities. Despite this, the role of $\phi$ in angular rotation is significant. Within the parameter scope of this study, particle shape seems to exert a greater influence than the moment of inertia.

4.6.2. Rotation style

Following the discussion of the maximum and average angular velocity, we are faced with the following two additional questions concerning the rotation style of angular particles. (i) Do Platonic polyhedrons exhibit a preferential direction of rotation? (ii) How does this preferential direction evolve with $\mathcal {G}a$? To shed light on these questions, we need to visualize the orientation of the particle angular velocity vector during free settling.

As shown in figure 29, the particle angular velocity is close to zero in the SV and UV regimes. It escalates to a significant value with the emergence of the HS and UO regimes. We present in figure 30 the probability density function (p.d.f.) of the particle angular velocity vector orientation of the tetrahedron and the dodecahedron at several $\mathcal {G}a$, projected onto a spherical surface. We select the tetrahedron and dodecahedron as they have the highest and lowest angular velocities among the Platonic polyhedrons across most of the investigated $\mathcal {G}a$ values in figure 29. Our results indicate that $\varOmega _y$ is relatively small when compared with the markedly larger magnitudes of the horizontal components ($\varOmega _x$, $\varOmega _z$), corroborating previous studies on freely settling cubes (Seyed-Ahmadi & Wachs 2019). As seen in figure 30, regions with a high probability of the angular velocity vector orientation are primarily located close to the global equator.

Figure 30. The p.d.f. map of the angular velocity vector orientation in the spherical coordinates projected on a globe surface; $\boldsymbol {\varOmega }$ orientation for a tetrahedron: (a) $\mathcal {G}a = 20$ (orange shade), (b) $\mathcal {G}a = 50$ (orange shade), (c) $\mathcal {G}a = 60$ (red shade), (d) $\mathcal {G}a = 300$ (blue shade); $\boldsymbol {\varOmega }$ orientation for a dodecahedron: (e) $\mathcal {G}a = 100$ (brown shade), ( f) $\mathcal {G}a = 170$ (brown shade), (g) $\mathcal {G}a = 180$ (dark brown shade) and (h) $\mathcal {G}a = 300$ (blue shade).

At $\mathcal {G}a=20$ for a tetrahedron and $\mathcal {G}a = 100$ for a dodecahedron, the magnitude of $\boldsymbol {\varOmega }$ is relatively small and its orientation exhibits substantial fluctuations on the sphere surface. In the UV regime in figure 30(a), the rotation orientation paths are blurred. However, in the SO regime, as depicted in figure 30(e), the path of the angular velocity vector orientation is discernible, which means that particle rotation exhibits a more stable orientation. Yet, it is worth noting that points on the north/south pole are scarcely observed in both cases, suggesting a rare occurrence of autorotation around the $y$ axis. As $\mathcal {G}a$ increases, the orientation of the angular velocity vector becomes more stable, as seen in figures 30(b) and 30( f). When the UV regime transitions to the HS regime at $\mathcal {G}a = 50 \sim 60$, the concentration of red dots around the globe equator intensifies, as illustrated by comparing figures 30(a) and 30(c). Similarly, a clear transition is observed between figures 30( f) and 30(g), from the SO regime to the UO regime.

To better visualize the evolution of the angular velocity vector, we plot the trajectory of the particle angular orientation on a two-dimensional plot as a function of azimuthal ($\varphi$) and elevation angle ($\theta$), as depicted in figure 31. The same $\mathcal {G}a$ values as in figure 30 are selected for easier comparison. Each scatter point represents a angular velocity vector orientation at a single time step. Joint p.d.f.s with respect to $\varphi$ and $\theta$ are displayed on the top and right axes.

Figure 31. Angular velocity vector orientation trajectory as a function of azimuth $\varphi$ and elevation angle $\theta$ during the freely settling of the Platonic polyhedron (scatter points); with joint p.d.f. distributions plotted above and to the right of each panel. (a) Tetra, $\mathcal {G}a= {20}-{50}$ (orange shade). (b) Tetra, $\mathcal {G}a={50}-{60}$ (orange and red shade). (c) Dodeca, $\mathcal {G}a={100}-{170}$ (brown shade). (d) Dodeca, $\mathcal {G}a={170}-{180}$ (brown and dark brown shade).

In the same settling regime, the scatter points become denser within the range $-50^\circ \leqslant \theta \leqslant 50^\circ$ as $\mathcal {G}a$ increases, yet they retain similar distribution patterns, as depicted in figures 31(a) and 31(c). However, transitioning from the UV regime to the HS regime reveals the emergence of curly striped structures around $\theta = 0$ in figure 31(b). The loop-d-loop tracks on these stripes suggest the HS path of the tetrahedron at $\mathcal {G}a = 60$. Similarly, for the dodecahedron, the shift from a sparse scatter distribution to the thin, wavy striped structures can be noted in figure 31(d). Correspondingly, the width of the p.d.f. curve with respect to the elevation $\theta$ becomes significantly narrower and the peak at $\theta = 0$ becomes more pronounced for both the tetrahedron and dodecahedron.

In figure 31(b), two distinct peaks are observed in the joint p.d.f. distribution with respect to the azimuth $\varphi$ (figure top axis). The difference between these peaks is $\Delta \varphi \approx 180^\circ$, indicating frequent switching of the angular velocity vector to its opposite direction. When the angular velocity vector orientation is close to the equator but evenly distributed on both sides, the rotation style of the Platonic polyhedrons is likely to be fluttering with a small rotation magnitude, which is the case for the freely falling tetrahedron at $\mathcal {G}a = 60$. Comparable peaks are also observed for the dodecahedron as displayed in figure 31. When the angular velocity vector orientation consistently remains on the same side of the equator, the particle is presumed to be tumbling with a larger magnitude or even commencing a rolling motion. For example, as depicted in figure 30(d), the tetrahedron is observed to be rolling at $\mathcal {G}a = 300$.

The fact that the angular velocity vector possesses a significant horizontal component, orthogonal to the particle settling direction ($y$ axis), underscores the importance of the Magnus force. We now demonstrate that the Magnus force also plays an important role in the vertical direction to amplify the settling drag. For simplicity and without compromising generality, we assume that the angular velocity vector $\boldsymbol {\varOmega }$ points in the $z^+$ direction. Since the particle settles along the $y^-$ direction, the Magnus force due to vertical motion can be dimensionally calculated as

(4.7)\begin{equation} \boldsymbol{F_{M,v}^*} = \frac{M^*}{m}C_\varOmega(\boldsymbol{\varOmega}^* \times \boldsymbol{U}^*) ={-}\frac{M^*}{m}C_\varOmega\varOmega_z^* U_y^* \boldsymbol{e_z}\times \boldsymbol{e_y} = \frac{M^*}{m}C_\varOmega\varOmega_z^* U_y^* \boldsymbol{e_x}, \end{equation}

where $e_y$ and $e_z$ are two unit vectors along the $y$ and $z$ axis. We see that the $\boldsymbol {F_{M,v}^*}$ aligns with the $x^+$ direction. As discussed in § 4.5.4, the Magnus force is the main driving force of the particle horizontal motion. Consequently, the particle velocity $U_x^*$ receives a significant contribution from $\boldsymbol {F_{M,v}^*}$, denoted as $U_{M,x}^*$. This rotation-induced horizontal drift can subsequently generate another Magnus force:

(4.8)\begin{equation} \boldsymbol{F_{M,h}^*} = \frac{M^*}{m}C_\varOmega(\varOmega_z^* \boldsymbol{e_z} \times U_{M,x}^*\boldsymbol{e_x}) = \frac{M^*}{m}C_\varOmega\varOmega_z^* U_{M,x}^*\boldsymbol{e_y} \end{equation}

that is oriented towards the $y^+$ direction. Hence, $\boldsymbol {F_{M,h}^*}$ can increase the drag force originating from the vertical settling. This sheds light on the change in decrease rate of $C_d$ from the UV regime to the HS regime in the tetrahedron case in figure 21(a), as the average angular velocity $\overline {\varOmega ^*}$ starts to rise (see figure 29a).

4.7. Particle–fluid density ratio

In numerous prior studies, the pair ($\mathcal {G}a, m$) has been deemed sufficient to determine the regime of freely rising or settling particles such as spheres and disks (Jenny et al. 2003; Auguste & Magnaudet 2018). However, the effects of the density ratio $m$ are less significant than $\mathcal {G}a$ for the settling of heavy particles in a Newtonian fluid (Seyed-Ahmadi & Wachs 2019). In this section we probe the sensitivity of regime transitions to variations in $m$ of the settling tetrahedron, the most angular particle among the Platonic polyhedrons. In natural environments the particle-to-fluid density ratio ranges from $m=2$ to $m=3$ for river sedimentation (Seyed-Ahmadi & Wachs 2019). Thus, we perform an additional series of simulations for freely settling tetrahedrons with $m=3$ in the range of $10 \leqslant \mathcal {G}a \leqslant 300$. The path instabilities of the settling tetrahedron are investigated to shed light on the influence of density ratio on regime transitions.

Figure 32(a) provides a comparison of regime transitions for freely settling tetrahedrons with density ratios $m=2$ and $m=3$. As $m$ increases, the settling velocity of the particle likewise increases, causing an overall shift to the left in the regime map. For instance, at $\mathcal {G}a=10$, the particle with $m=3$ settles in the UV regime, while the particle with $m=2$ remains in the SV regime. Following the UV regime, the HS regime commences at $\mathcal {G}a=60$, as in the case when $m=2$. However, the range of the HS regime is marginally narrower, ending at $\mathcal {G}a=140$ for a tetrahedron with $m=3$.

Figure 32. Comparison of freely settling tetrahedron at $m=2$ ($\blacktriangle$) and $m=3$ ($\blacktriangledown$). (a) Time-averaged drag coefficient $\overline {C_d}$. (b) Time-averaged angular velocity $\bar {\varOmega }$. Results are shown for the SV ($\sqsubset \!\sqsupset$, filled black), UV ($\sqsubset \!\sqsupset$, filled orange), SO ($\sqsubset \!\sqsupset$ filled brown), UO ($\sqsubset \!\sqsupset$, filled maroon), HS ($\sqsubset \!\sqsupset$, filled red) and CS ($\sqsubset \!\sqsupset$, filled blue) regimes.

The impacts of $m$ on $\overline {C_d}$ and $\bar {\varOmega }$ are depicted in figures 32(b)–32(c). At lower $\mathcal {G}a$ values, the $\overline {C_d}$ for particles with $m=2$ and $m=3$ are nearly identical. The difference becomes more pronounced when $\mathcal {G}a>200$. As for $\bar {\varOmega }$, the HS regime commences at $\mathcal {G}a=60$ for both density ratios. This initiates a concurrent increase in $\bar {\varOmega }$ for both particles. In conclusion, the flow dynamics and regime transitions of freely settling angular particles, even for the most angular particle, the tetrahedron, do not show significant sensitivity to moderate changes in the density ratio.

5. Conclusions and perspectives

We carried out a comprehensive investigation of the dynamics of a single Platonic polyhedron freely settling in an unbounded domain filled with quiescent Newtonian fluid. With five Platonic polyhedrons of increasing sphericity $\phi$, we classified settling regimes according to the path instability and wake structures over the range of Galileo numbers $10 \leqslant \mathcal {G}a \leqslant 300$. We provided new insights into the mechanism of HS of angular particles, elucidating the commonalities and differences in the particle dynamics between the oblique and HS regimes. Taking a step further, we introduced two novel empirical correlations based on $\mathcal {G}a$ and $\phi$, allowing for the accurate prediction of particle settling ${\mathcal {R}e}$ and disturbed wake length.

A regime map in the ($\mathcal {G}a, \phi$) parameter space is proposed, where we highlight the profound influence of the particle angularity on both its settling path and the onset of instabilities. We confirm that the initial angular position of the angular particle only affects the initial transient settling. Across all polyhedrons, a stable angular position emerges and may change during regime transitions. As the particle angularity amplifies, the settling path manifests unsteadiness at lower $\mathcal {G}a$. The progression of regime transitions is thoroughly investigated: starting from a SV descent, the particle evolves into UV settling due to the unsteadiness of the wake structures. With the increase of $\mathcal {G}a$, the particles with high $\phi$ (octahedron, dodecahedron and icosahedron) turn towards oblique settling, controlled by the drift force and asymmetry of the wake structures. By contrast, the tetrahedron and the cube exhibit HS trajectories, due to the simultaneous appearance of the stable horizontal force and the substantial streamwise torque. The HS in the tetrahedron case is shown to be a consequence of the high particle angularity, exhibiting a unique double-helix vortical structure in the wake, whereas in the cube case, the HS is propelled by the hairpin vortex shedding. Our detailed force balance in the horizontal plane reveals that the Magnus force acts as the driving force in the HS motion. Additionally, the Magnus force contributes to the increase of the particle vertical drag, once the angular velocity becomes important, such as in the helical and UO regimes. Novel empirical correlations are derived for ${\mathcal {R}e}$ and $L_w$ based on our high-fidelity numerical results. The particle angular velocity is observed to increase with $\mathcal {G}a$, with the particle angular velocity vector predominantly oriented in the horizontal plane. We show that even for the settling of the most angular particle, the tetrahedron, the change of density ratio $m$ from $2$ to $3$ induces minimal perturbation on path instability and regime transitions.

In summary, our study bridges critical gaps in understanding Platonic polyhedron settling dynamics and offers valuable benchmarks and predictive tools for future investigations and applications. In future research it would be highly beneficial to delve into the dynamics of multiple angular particle suspensions. The interactions between individual particles in such a system can provide significant insights into the collective motion such as sedimentation patterns and the potential turbulence modulation, to optimize industrial processes and understand natural environmental phenomena.