1. Introduction
There are numerous examples of non-spherical atmospheric particles from both natural and anthropogenic sources. Examples include volcanic (Rossi et al. Reference Rossi, Bagheri, Beckett and Bonadonna2021) and wildfire (Bodí et al. Reference Bodí, Martin, Balfour, Santín, Doerr, Pereira, Cerdà and Mataix-Solera2014) ashes, sand (Merrison Reference Merrison2012), dust (Mahowald et al. Reference Mahowald, Albani, Kok, Engelstaeder, Scanza, Ward and Flanner2014; van der Does et al. Reference van der, M., P., P. and Stuut2018) and microplastics (Brahney et al. Reference Brahney, Hallerud, Heim, Hahnenberger and Sukumaran2020; Tatsii et al. Reference Tatsii, Bucci, Bhowmick, Guettler, Bakels, Bagheri and Stohl2023; Reininger et al. Reference Reininger, Tatsii, Bhowmick, Bagheri and Stohl2025) (see figure 1
a). Those particles are known to have various impacts, from increasing health risks for humans (Horwell & Baxter Reference Horwell and Baxter2006; Manisalidis et al. Reference Manisalidis, Stavropoulou, Stavropoulos and Bezirtzoglou2020; Zhang et al. Reference Zhang, Xu, Li, Chen, Ma, Zeng and Shi2020a
) to the radiative properties of clouds (Baran Reference Baran2012; Matus & L’Ecuyer Reference Matus and L’Ecuyer2017), which influence the weather and climate (Ramanathan et al. Reference Ramanathan, Cess, Harrison, Minnis, Barkstrom, Ahmad and Hartmann1989; Norris et al. Reference Norris, Allen, Evan, Zelinka, O’Dell and Klein2016). However, our understanding of the transport properties of these particles is rather limited. This is due to the lack of theoretical models for heavy particles larger than 0.1 mm in air, and the existing experimental data come almost exclusively from studies conducted in liquids. The experiments in liquids (e.g. Willmarth, Hawk & Harvey Reference Willmarth, Hawk and Harvey1964; Jayaweera & Cottis Reference Jayaweera and Cottis1969; Field et al. Reference Field, Klaus, Moore and Nori1997; Belmonte, Eisenberg & Moses Reference Belmonte, Eisenberg and Moses1998; Lopez & Guazzelli Reference Lopez and Guazzelli2017; Roy et al. Reference Roy, Hamati, Tierney, Koch and Voth2019, Reference Roy, Kramel, Menon, Voth and Koch2023; Esteban, Shrimpton & Ganapathisubramani Reference Esteban, Shrimpton and Ganapathisubramani2020; Will et al. Reference Will, Mathai, Huisman, Lohse, Sun and Krug2021; Baker & Coletti Reference Baker and Coletti2022; Cabrera et al. Reference Cabrera, Sheikh, Mehlig, Plihon, Bourgoin, Pumir and Naso2022; Giurgiu et al. Reference Giurgiu, Caridi, De Paoli and Soldati2024), in which the dynamic similarity is observed using the particle Reynolds number
${\textit{Re}}_{\!p}$
, provide many insightful findings on the settling behaviour of non-spherical particles, with
\begin{equation} {\textit{Re}}_{\!p} = \frac {D_{\textit{eq}} v_T} {\nu }, \end{equation}
where
$D_{\textit{eq}}$
is the diameter of a volume-equivalent sphere,
$v_T$
is the terminal velocity of the particle, and
$\nu$
is the kinematic viscosity of the ambient fluid. Here,
${\textit{Re}}_{\!p}$
quantifies the relative importance of the inertial forces compared to the viscous forces acting on the particle (Happel & Brenner Reference Happel and Brenner1983; Khayat & Cox Reference Khayat and Cox1989; Gustavsson et al. Reference Gustavsson, Sheikh, Lopez, Naso, Pumir and Mehlig2019).
Figure 1. (a) Natural and anthropogenic sources of non-spherical solid particles that impact the weather, climate and ecosystems. (b) Atmospheric particles in the phase space of particle Reynolds number
${\textit{Re}}_{\!p}$
and particle to fluid density ratio
$R$
(Kajikawa Reference Kajikawa1972, Reference Kajikawa1992; van Hout & Katz Reference van Hout and Katz2004; Pruppacher & Klett Reference Pruppacher and Klett2010; Bagheri & Bonadonna Reference Bagheri and Bonadonna2016; Allen et al. Reference Allen, Allen, Phoenix, Jiménez, Simonneau, Binet and Galop2019; Zhang et al. Reference Zhang, Kang, Allen, Allen, Gao and Sillanpää2020b
). The relevant
${\textit{Re}}_{\!p}$
and
$R$
of this current study are shown with the black
$\diamond$
.
However, the ratio of particle to fluid density
\begin{equation} R = \frac {\rho _p}{\rho _f}, \end{equation}
which in air is typically of the order of 1000 for most solid particles, cannot exceed much beyond 10 in liquids, as summarised in figure 1(b). Recent experiments at low Reynolds numbers in air, particularly for axisymmetric shapes such as spheroids, have advanced our understanding considerably. Analytical and empirical approaches reveal that particle inertia strongly affects settling dynamics, often inducing pronounced transient oscillations (Bhowmick et al. Reference Bhowmick, Seesing, Gustavsson, Guettler, Wang, Pumir, Mehlig and Bagheri2024) or sustained orientation instabilities at Reynolds numbers much lower than those observed at low
$R$
(Tinklenberg, Guala & Coletti Reference Tinklenberg, Guala and Coletti2023). While these findings clarify the behaviour of axisymmetric particles or thin disks, the broader phase space remains unexplored, in particular regarding how non-axisymmetric features shape settling at Reynolds numbers applicable to a wide range of atmospheric particles, i.e.
$1 \lt {\textit{Re}}_{\!p} \lt 10$
at
$R = 1000$
as shown in figure 1(b). This raises a critical question: how do these non-axisymmetric features, which are key to generalising our insights across the diverse shapes in the atmosphere, affect settling dynamics?
This study targets the parameter range
$1\lt {\textit{Re}}_{\!p}\lt 10$
at
$R = 1000$
(
$\rho _p={1200}\,\textrm {kg}\boldsymbol{\,}\textrm {m}^{-3}$
and
$\rho _f = {1.2}\,\textrm {kg}\boldsymbol{\,}\textrm {m}^{-3}$
), which is relevant for many atmospheric particle sources yet remains largely unexplored experimentally or numerically. We examine the settling of heavy ellipsoids of various shapes in still air, leveraging their utility as surrogates for irregular particles (Bagheri & Bonadonna Reference Bagheri and Bonadonna2016). By combining novel experiments with particle-resolved direct numerical simulations (DNS), we uncover complex transient dynamics and highlight a largely uncharted regime of non-spherical particle behaviour, with important implications for understanding particle-rich atmospheric environments.
2. Methods
2.1. Experiments
The Göttingen Turret, the experimental set-up for this study, is developed at the Max Planck Institute for Dynamics and Self-Organization (MPI-DS). It features four synchronised high-speed cameras – TX, TY, BX and BY – imaging an air-filled glass chamber (GC) where particles are injected by a particle injector (PI) and fall under gravity in still air, as shown in figure 2(a,b). Two cameras capture each of two overlapping volumes from orthogonal directions (see figure 2
b), enabling full three-dimensional (3-D) tracking of particles with at least one dimension
${\geqslant}{100}\,{\unicode{x03BC} }\textrm {m}$
.
Figure 2. (a) A 3-D CAD drawing of the air-filled glass chamber (GC) and the particle injector (PI) for releasing particles inside the GC. (b) Schematic representation of the illumination of the GC and the positioning of the mirrors (M) for imaging the particles by the upper TX and TY cameras, and the lower BX and BY cameras. Axis
$Z$
is in the direction of gravity. (c) The shapes of ellipsoidal, non-spherical particles of this study are defined with two shape parameters – elongation
$\textit{EL}$
and flatness
$\textit{FL}$
– keeping the volume constant (equivalent to a sphere of
${140}\,{\unicode{x03BC} }\textrm {m}$
diameter
$D_{\textit{eq}}$
). The particle shape-space is divided into four categories: I – near spherical, II – prolate/elongated, III – highly non-spherical ellipsoidal, and IV – oblate/flattened shapes. (d) Laser microscope scans of two 3-D printed particles from categories I and III. The 3-D surface features are
${\lt}{1}\,{\unicode{x03BC} }\textrm {m}$
in size. (e) Experimental observation of oscillations in angular orientation and lateral drift when non-spherical particles are falling under gravity
$g$
. Images (a) and (b) are adjusted from figure S1 in the supplementary information of Bhowmick et al. (Reference Bhowmick, Seesing, Gustavsson, Guettler, Wang, Pumir, Mehlig and Bagheri2024).
All cameras are equipped with Nikon ED AF Micro Nikkor 200 mm lenses, operated at F-stop 22 and 1 : 1 magnification, resulting in an image pixel size of 6.75
$\unicode{x03BC}$
m px–1. The effective exposure time is 2.5
$\unicode{x03BC}$
s. The camera system is synchronised in a leader–follower configuration, with camera TX acting as the leader, and its strobe signal used to synchronise the remaining three follower cameras. Given the fixed data throughput of 9 Gpx s–1 (Phantom VEO-Series, Vision Research, AMETEK), maximising spatial resolution required a trade-off with temporal resolution. We employed the maximum available resolution 4096 px in the
$ Z$
(falling) direction, corresponding to a vertical observation window of approximately 28 mm per camera at 1 : 1 magnification, yielding a total vertical measurement distance of approximately 55 mm when the top and bottom camera fields of view are combined. The horizontal field of view is 720 px (4.86 mm). Images are captured at frame rate 2932 frames per second. At this temporal resolution, changes in the angular motion of the particles are sufficiently small to allow for accurate reconstruction of angular dynamics, typically below 100 Hz. A detailed description of the set-up is provided in Bhowmick et al. (Reference Bhowmick, Seesing, Gustavsson, Guettler, Wang, Pumir, Mehlig and Bagheri2024).
Ellipsoids are 3-D-printed using a NanoScribe system (Photonic Professional (GT) User Manual, Nanoscribe GmbH 2015), achieving surface resolution
$ {\leqslant}{0.3}\,{\unicode{x03BC} }\textrm {m}$
. The photoresist acrylic resin IP-S
$ (\rm CH_{1.72}N_{0.086}O_{0.37})$
is polymerised via two-photon polymerisation (Liu et al. Reference Liu, Campbell, Stein, Jiang, Hund and Lu2018). The particles are printed in stacked polymerised layers
$ {\leqslant} {0.6}\,{\unicode{x03BC} }\textrm {m}$
thick, with mass density
$ \rho _p = {1200}\,\textrm {kg}\boldsymbol{\,}\textrm {m}^{-3}$
, matching many atmospheric particles (see figure 1
b).
For imaging, particles are carefully transferred from the glass substrate to an injection needle for controlled release into the air-filled glass chamber, which is positioned in such a way that the cameras capture the transient settling motion of particles a few centimetres (adjustable) below the detachment point of the injection needle. All particles are initially placed on the tip of the injection needle so that their largest surface area is in contact with the needle. Such an arrangement results in their maximum projected area being oriented normal to gravity. However, high-speed imaging of the release process (performed in separate experiments) revealed that the initial orientation cannot be reliably controlled. In some instances, particles detach in arbitrary orientations with non-zero lateral velocity, despite consistent placement. The initial release vertical velocity, however, can be controlled to some extent by adjusting the vertical position from which the needle is released (Bhowmick et al. Reference Bhowmick, Seesing, Gustavsson, Guettler, Wang, Pumir, Mehlig and Bagheri2024). This height was adjusted such that for the majority of core experiments for a given shape, both the transient and terminal phases of the settling motion occurred entirely within the tracking volume.
The vertical distance between the detachment point of the particle from the injection needle and the start of the imaging field is adjustable and is selected to ensure that both transient dynamics and the terminal settling state are fully captured. The terminal state is defined as the point in time at which the vertical acceleration of the particle
$|a_z|$
falls below
${0.2}\,\textrm {m}\boldsymbol{\,}\textrm {s}^{-2}$
. This threshold was chosen based on an analysis of vertical acceleration profiles within the field of view, and serves as a conservative lower bound that accounts for all measurement uncertainties. We also calculated particle lateral drift
$\delta$
, which is defined as the total lateral path length traversed by the particle (normal to gravity) in the
$xy$
-plane over fixed settling distance 55 mm (corresponding to the full height of the combined observation volume). As mentioned above, there may be non-zero initial lateral velocities when particles are injected; however, these do not constitute systematic biases. Therefore, we focus on the general statistical analysis of lateral drifts based on the experimental measurements. Further details on the 3-D calibration procedure, sub-pixel spatial accuracy, and computation of velocity and acceleration can be found in Bhowmick et al. (Reference Bhowmick, Seesing, Gustavsson, Guettler, Wang, Pumir, Mehlig and Bagheri2024).
The particles of this study are of ellipsoidal shape, defined by two shape parameters: elongation
$\textit{EL}$
, and flatness
$\textit{FL}$
, as visualised in figure 2(c) and defined as
\begin{equation} \textit{EL} = \frac {I}{L}, \quad \textit{FL} = \frac {S}{I}, \end{equation}
where
$L$
,
$I$
and
$S$
are the longest, intermediate and shortest dimensions of the particle, which are orthogonal to each other. The main characteristics of the ellipsoids investigated here are summarised in table 1. The shapes considered in this study are divided into four categories, each containing an equal number of cases. The particle shape space is categorised as follows: category I – near-spherical, II – prolate/elongated, III – highly non-spherical ellipsoidal, and IV – oblate/flattened shapes (figure 2
c). The control parameter in our experiment is the particle mass (or volume), which is equal to that of a volume-equivalent sphere with diameter
$D_{\textit{eq}} = {140}\,{\unicode{x03BC} }\textrm {m}$
for all particles. The particles are produced using two-photon polymerisation 3-D printing, which enables us to achieve an explicit control over particle shape and size with sub-micrometre printing resolution; e.g. see microscope images of the manufactured particles in figure 2(d).
Table 1. Details of the particles, with category numbers, shape parameters elongation
$\textit{EL}$
and flatness
$\textit{FL}$
. Here,
${\textit{Re}}_{\!p} = D_{\textit{eq}} v_T/\nu$
is the particle Reynolds number based on the equivalent sphere diameter
$D_{\textit{eq}} = {140}\,{\unicode{x03BC} }\textrm {m}$
and on the observed terminal velocity
$v_T$
from the DNS. The Stokes time
$\tau _p = ({\rho _p}/{18\rho _f}) LI/\nu$
provides an estimate of the particle response time, where
$\rho _p = {1200}\,\textrm {kg}\boldsymbol\,{}\textrm {m}^{-3}$
and
$\rho _f = {1.2}\,\textrm {kg}\boldsymbol{\, }\textrm {m}^{-3}$
are the mass-densities of the particle and the air, respectively, and
$\nu = {1.5\times 10^{-5}}\,\textrm {m}^{2}\boldsymbol{\,}\textrm {s}^{-1}$
is the kinematic viscosity of air. For all the particles shown in this table,
$R=1000$
. A ✓ indicates the presence of data; otherwise, a
$\times$
is shown.
Higher non-sphericity causes increased drift from the observation volume, affecting imaging success rates. The sphere had the highest success rate (34.4 %), while the ellipsoid with most flatness and elongation (
$\textit{EL}=\textit{FL}=0.2$
) had the lowest (11.5 %). Data analysis follows a camera model (Tsai Reference Tsai1987) to generate 3-D world coordinates from calibration grid images. The centre of mass and longest particle length are determined via ellipsoid fitting, and the velocity and angular orientation are extracted accordingly (Bhowmick et al. Reference Bhowmick, Seesing, Gustavsson, Guettler, Wang, Pumir, Mehlig and Bagheri2024).
2.2. Simulations
We employed particle-resolved DNS based on the lattice Boltzmann method (LBM) (Succi Reference Succi2001) to model the settling dynamics of atmospheric particles in a quiescent fluid. The numerical solver Palabos Turret (Bhowmick et al. Reference Bhowmick, Latt, Wang and Bagheri2025a
,
Reference Bhowmick, Latt, Wang and Bagherib
), originally developed by the Palabos group (Latt et al. Reference Latt2021) at the University of Geneva, was further modified and verified at the MPI-DS for this study. It solves the quasi-incompressible Navier–Stokes equations via the LBM and Newton’s equations for particle motion, coupled through the immersed boundary method (IBM) (Ota, Suzuki & Inamuro Reference Ota, Suzuki and Inamuro2012). It is validated through comparisons with analytical predictions and experimental observations of the free-fall dynamics of particles of different shapes, including spheres, ellipsoids and volcanic fragments, across a wide range of
${\textit{Re}}_{\!p}$
(Bhowmick et al. Reference Bhowmick, Latt, Wang and Bagheri2025a
).
Palabos Turret utilises a regularised LBM with a BGK collision operator (Bhatnagar, Gross & Krook Reference Bhatnagar, Gross and Krook1954) to resolve the 3-D flow field with a uniform Eulerian mesh. The fluid is modelled as particle populations, which follow the evolution rule as
\begin{equation} f_i(\boldsymbol{x}+\boldsymbol{\xi }_i\,\delta t,t+\delta t) = f^{(0)}_i+\left (1-\frac {1}{\tau }\right )f^{(1)}_i. \end{equation}
Here
$f_i$
,
$f^{(0)}_i$
and
$f^{(1)}_i$
are the particle distribution function and its equilibrium and off-equilibrium populations, respectively (Bhowmick et al. Reference Bhowmick, Latt, Wang and Bagheri2025a
). The index
$i=0,\ldots , q-1$
identifies the elements of a discretised velocity space
$\{\boldsymbol{\xi }_i\}_{i=0}^{q-1}$
, which is represented by the D3Q27 lattice;
$\delta \boldsymbol{x}$
and
$\delta t$
are the spatial and temporal resolutions, respectively; and
$\tau$
is the relaxation time, which relates to the kinematic viscosity
$\nu$
through the speed of sound
$c_s$
as
\begin{equation} \nu =c_s^2\,\delta t\left (\tau -\frac {1}{2}\right )\!. \end{equation}
Macroscopic variables of the fluid, such as density
$\rho$
and velocity
$\boldsymbol{u}$
, can be calculated by taking a weighted sum of the particle populations:
\begin{equation} \rho =\sum _i f_i, \quad \boldsymbol{u}=\frac {\sum _i \boldsymbol{\xi }_i f_i}{\sum _i f_i}. \end{equation}
The hydrodynamic stress tensor
$\boldsymbol{\sigma }$
, which is used for the evaluation of the fluid forces on the solid surface, can be obtained directly from
$f_i$
as
\begin{equation} \boldsymbol{\sigma } = -p \boldsymbol{I} + \left (\frac {1}{2\tau }-1\right )\sum _i \boldsymbol{\xi }_i \boldsymbol{\xi }_i\big(f_i-f^0_i\big). \end{equation}
Here,
$p=c_s^2 \rho$
is the fluid pressure, and
$\boldsymbol{I}$
is the identity matrix. The first (second) term on the right-hand side of (2.5) corresponds to the pressure (friction) component of the stress tensor.
Solid particles, represented by non-deformable Lagrangian surface meshes, are immersed in the fluid and allowed to move and rotate freely (figure 3). Given a total force
$\boldsymbol F$
acting on a particle of mass
$m$
, the position of its centre of mass
$\boldsymbol x$
, translational velocity
$\boldsymbol v$
and angular velocity
$\boldsymbol \omega$
are governed by
\begin{equation} \frac {{\rm d}{\boldsymbol x}}{{\rm d}t} = {\boldsymbol v}, \quad \frac {{\rm d}{\boldsymbol v}}{{\rm d}t} = \frac {{\boldsymbol F}}{m}, \quad \frac {{\rm d}{\boldsymbol d}}{{\rm d}t} = {\boldsymbol \omega }, \end{equation}
respectively. Here,
$\boldsymbol d$
is an orthonormal triad attached to the particle and following its motion. With inertia tensor
$\boldsymbol J$
and torque
$\boldsymbol T$
, the evolution of the particle’s rotation is
\begin{equation} \frac {\rm d}{{\rm d}t}{\boldsymbol \omega } = {\boldsymbol J}^{-1}{\boldsymbol T}. \end{equation}
Forces on the particles, such as gravity and fluid-induced inertia, are computed using the IBM. In particular, the ‘multi-direct forcing’ approach (Ota et al. Reference Ota, Suzuki and Inamuro2012), with an iterative procedure, is used to couple translational and rotational dynamics between fluid and particle.
Figure 3. (a) The 3-D computational domain showing a particle (black ellipsoid) at an intermediate stage of settling under gravity in the
$ z$
-direction, after being released from rest at height
$ h_{\textit{initial}}$
. Periodic boundary conditions are applied laterally, while the top boundary is held at constant pressure, and the bottom boundary is a no-slip (wall) boundary. (b) An exemplary simulated particle with
$ \textit{EF} = 0.2$
and
$ {\textit{FL}} = 1.0$
, shown alongside a portion of the surrounding fluid domain. The cross-sectional colour map represents the local velocity field. (c) A close-up view of (b), showing the uniform Eulerian mesh in the fluid domain, and the Lagrangian surface mesh on the particle.
Similar to the experimental set-up, the computational domain is a rectangular chamber with constant pressure (i.e. density) imposed at the top, a no-slip condition at the bottom, and periodic boundaries on the lateral sides, as shown in figure 3(a). Based on mesh independence tests carried out on various shapes, we have found that the lateral width of the computational domain must be at least 10 times the largest particle dimension,
$ L$
, to minimise the wall effects. It is also found that at least 10 fluid points in the smallest particle dimension are needed to resolve the boundary layers, and to ensure reliable simulations of translational and rotational dynamics (figure 3
b). The spatial resolutions of the fluid and solid surface meshes are approximately matched, with resolution ratio close to unity (figure 3
c). These constraints make it computationally infeasible to simulate the full settling trajectory, from rest to terminal velocity, within a single domain. For a computational domain of size
$ 10L \times 10L \times 1000L$
, with mesh resolution
$ \Delta x = S/10$
, the total number of mesh points becomes
$ (10L / \Delta x)^2 \times (1000L / \Delta x) = 10^8 L^3 / S^3$
. For a typical case with
$ {\textit{EL}} = I/L = 0.3$
and
$ {\textit{FL}} = S/I = 0.3$
, this gives
$ S/L = {\textit{EL}} \times {\textit{FL}} = 0.09$
, thus the total number of mesh points is approximately
$ 10^8 / (0.09)^3 \approx 1.37 \times 10^{11}$
, which makes such a full-domain simulation currently computationally extremely challenging.
To address the computational challenges of simulating particle settling over long distances, we developed a truncated settling chamber, illustrated in figure 3(a), which is comprehensively discussed and validated in Bhowmick et al. (Reference Bhowmick, Latt, Wang and Bagheri2025a
). The domain vertically spans approximately three times the lateral width, amounting to approximately
$ 30L$
for the particles considered in this study, and is divided into three vertical sections: (i) a sponge zone at the top, occupying 10 % of the total height, designed to damp spurious pressure waves generated by numerical artefacts; (ii) a central region where the particle settles sufficiently at distance from all domain boundaries; and (iii) a near-ground zone at the bottom, also comprising 10 % of the total height. The bottom boundary, representing a no-slip wall, serves to prevent artificial air entrainment and downward flow that could arise if periodic boundary conditions were used on both the top and bottom boundaries.
The particles simulated in this study are released with initial velocity zero and oriented such that the longest dimension,
$ L$
, is aligned with gravity, resulting in the minimum projection area perpendicular to gravity at the initial release height
$ h_{\textit{initial}}$
(figure 3
a). Once the particle reaches a lower height
$ h_{\textit{final}}$
, both the particle and the surrounding flow field are duplicated and repositioned at the top of the computational domain, effectively resetting the particle to its initial height
$ h_{\textit{initial}}$
to continue its descent. This copy-and-reposition procedure is repeated until the particle reaches its terminal state, defined as the point at which the vertical acceleration satisfies
$ |a_z| \lt 0.01\,\mathrm{m\,s}^{-2}$
. This approach enabled efficient simulation of various particle configurations without requiring a domain long enough to capture the entire settling trajectory. Nonetheless, the simulations remained computationally intensive. For instance, the most demanding case, with
$ {\textit{EL}} = {\textit{FL}} = 0.3$
, was performed on the supercomputer Cobra at the Max Planck Computing and Data Facility, employing grid resolution
${4.2}\,{\unicode{x03BC} }\textrm {m}$
, a total of 4.11 billion mesh nodes, and 15 360 MPI tasks, amounting to 3.1 million computing hours.
As mentioned earlier, Bhowmick et al. (Reference Bhowmick, Latt, Wang and Bagheri2025a
) have found very close agreement between numerical results obtained from Palabos Turret and experimental measurements for various regular and irregular shapes in a wide range of particle Reynolds numbers. For the ellipsoids studied here, validations are also presented in § 3.1. The numerical simulations provided detailed information, including the particle’s position, velocity, orientation and surface force distribution, and the surrounding velocity field. This comprehensive dataset from numerical simulations complements the experimental measurements and allows the two approaches to be cross-validated. Quantities measurable by both approaches for cross-validation include the terminal velocity, the lateral drift
$\delta$
, and the oscillation frequency and decay rate of the oscillation amplitudes of the longest axis of the ellipsoid (
$L$
). In simulations,
$\delta$
is defined as the total lateral path length traversed by the particle (normal to gravity) in the
$xy$
-plane until it reaches the terminal state. One-to-one comparison of
$\delta$
between experiments and simulations is not possible since in experiments the initial conditions could not be controlled precisely enough (see § 2.1), although their statistical signatures can be compared. Additional quantities accessible only through simulations include the pressure and friction components of the drag force, as well as the total swept volumes in the lateral (
$\varPsi _\perp$
) and vertical (
$\varPsi _z$
) directions.
2.3. Normalisation
In the following, normalised parameters are denoted by
$^*$
. The spatial coordinates
$x^*, y^*, z^*$
and the drift
$\delta ^*$
are normalised by the equivalent sphere diameter
$D_{\textit{eq}}$
, and velocity
$v$
is normalised by the terminal settling speed
$v_T$
. Time
$t$
is normalised as
$T^* = t v_T / D_{\textit{eq}}$
. The pressure components of the drag force on a particle are expressed in dimensionless form as
$\boldsymbol{P}^* = \boldsymbol{P}/F_z$
, where
$F_z$
denotes the vertical component of the drag force. The vertical and lateral swept volumes,
$\varPsi _z^*$
and
$\varPsi _\perp ^*$
, are normalised by the particle volume
$\unicode{x03C0} D^3_{\textit{eq}} /6$
. In simulations,
$\varPsi _\perp$
is defined as the total volume swept by the particle (normal to gravity) in the horizontal directions until it reaches the terminal state, whereas
$\varPsi _z$
is the total volume swept by the particle in the vertical directions until it reaches the terminal state.
3. Results
3.1. Experimental observations
Figure 4. (a) Top view of the drifting trajectories of the non-spherical particles from experiments (zoomed in for better visibility) coloured by particle categories I–IV shown in figure 2(c). Here,
$x^*$
and
$y^*$
are the two lateral components normalised by
$D_{\textit{eq}}$
. (b) Histograms of the normalised lateral drifts (
$\delta ^* = \delta /D_{\textit{eq}}$
) with bin width
$0.25\,\delta ^*$
from experimental measurements presented here. Here,
$\delta$
is calculated as total lateral displacement in the
$ xy$
-plane over fixed settling distance 55 mm, corresponding to the full height of the combined observation volume in experiments. Also,
$N_d$
is the number of experimental runs for a specific amount of
$\delta ^*$
, and the median and standard deviation values are shown as med. and
$\sigma$
. (c) Top view of the drifting particle trajectories from the DNS (rotated for better visibility), where the particles are released with zero initial velocity and minimum projection area normal to the falling direction. Shape parameters (
${\textit{EL}}$
,
${\textit{FL}}$
) for the non-planar trajectories are annotated next to the tracks. (d) Comparison of the terminal velocities
$v_T$
and frequencies of the angular oscillations
$f$
of the particles from experiments (vertical axes) and DNS (horizontal axes). The error bars are equivalent to one standard deviation variation. The dashed diagonal line shows 1 : 1 correspondence between experiments and DNS. The maxima of relative differences between the experiments and DNS are less than 3 % and 10 %, respectively, for
$v_T$
and
$f$
. We use
$R = 1000$
for all the results presented in this figure.
We have conducted a total of 262 successful experimental runs for different particle shapes and density ratio
$R=1000$
, with each shape having from 9 to 22 runs (see example snapshots from experiments in figure 2
e). By covering a wide range of flattened and/or elongated ellipsoids, we comprehensively explore the phase space of diverse particle shapes. This approach also serves as an effective proxy for irregular particles with similar form (Bagheri & Bonadonna Reference Bagheri and Bonadonna2016), eliminating the need for infinite case-by-case investigations. The particle Reynolds number achieved here ranges from approximately
${\textit{Re}}_{\!p} =2.1$
to
$4.5$
, as also detailed in table 1. This is significantly lower than the threshold Reynolds number of approximately 100, beyond which an unstable falling pattern is observed due to hydrodynamic instabilities in the wake of the falling particles in air and in liquids (Willmarth et al. Reference Willmarth, Hawk and Harvey1964; Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012; Auguste, Magnaudet & Fabre Reference Auguste, Magnaudet and Fabre2013; Esteban et al. Reference Esteban, Shrimpton and Ganapathisubramani2020; Tinklenberg et al. Reference Tinklenberg, Guala and Coletti2023).
Our experiments show that ellipsoids of various shapes exhibit orientation oscillations and lateral drifts as they settle in still air at
$2.1\leqslant {\textit{Re}}_{\!p}\leqslant 4.5$
and
$R=1000$
. The more flattened (low
$\textit{FL}$
) or elongated (low
$\textit{EL}$
) the ellipsoids are, the more pronounced these effects become. The top view of the particle trajectories (figure 4
a) and the histograms of lateral drift (figure 4
b) show the intricate falling behaviour of the investigated particles observed in experiments. Across all shape categories, particularly in categories II–IV, there is a notable lateral drift exceeding
$10 D_{\textit{eq}}$
during experiments. Some particles exhibit wavy lateral trajectories, while some others also make multiple sharp turns (figure 4
a). The elongated particles in categories III and II have the highest median lateral drift, while the standard deviation of lateral drift in categories II–IV is four times higher than the near spherical particles in category I (figure 4
b). As will be shown below, these lateral movements have a significant effect on the swept volume of the particles as they fall.
3.2. Numerical trajectories and validation
To fully understand the dynamics behind the complex trajectory of these particles at
$R=1000$
, we have replicated the experiments with large-scale DNS. As mentioned earlier, all simulations are carried out with a fixed initial condition, i.e. release with zero initial velocity and minimum projection area normal to the falling direction. Top views of DNS trajectories also reveal the wavy lateral motion of the particles (figure 4
c), with sharp turns and lateral drifts up to
$6 D_{\textit{eq}}$
for category III, compared to only one
$D_{\textit{eq}}$
lateral drift for category I. Of the seven triaxial ellipsoids exhibiting non-planar motion, four belong to category III, two to category II, and one to category I, whereas none belongs to category IV. This suggests that non-planar motion is common in extremely flattened and elongated triaxial ellipsoids, but occurs more frequently in elongated particles than in flattened ones.
The excellent quantitative agreement between the experiments and DNS for the terminal velocity
$v_T$
and oscillation frequency
$f$
of the particles (figure 4
d) suggests that many features of the transient dynamics are independent of initial conditions and could be generalised to real-world scenarios. For each shape, the standard deviation of
$v_T$
is less than
$2\,\%$
of its magnitude, and for
$f$
it is less than
$16\,\%$
, despite a large variability expected in the initial conditions of the experimental trajectories. This shows that the effects of the transient dynamics captured in the experiments are not dependent on initial conditions.
3.3. Effect of density ratio
$R$
Having successfully validated our DNS against the experimental measurements, we turn our attention to distinguishing the role of particle to fluid density ratio
$R$
in the settling dynamics of non-spherical particles with the same
${\textit{Re}}_{\!p}$
. We have selected one particle shape (fixed
$\textit{EL}$
and
$\textit{FL}$
) from each of the particle categories I–IV, as shown in figure 5. The density ratios
$R$
considered in this subsection in addition to
$R=1000$
– i.e.
$R=100$
and 10 – are representative of typical conditions for volcanic ash/snow aggregates in the atmosphere and particles in water, respectively. As expected, particles need a longer settling distance before reaching terminal velocity as
$R$
increases. For instance, at
$R=1000$
, the required
$z^* = z/D_{\textit{eq}}$
is approximately an order of magnitude longer than at
$R=10$
and 100.
Figure 5. Settling dynamics of particles with four different shapes, each examined at density ratios
$R = 1000$
, 100 and 10. The particle Reynolds number
${\textit{Re}}_{\!p}$
is held constant within each shape category. Trajectories of particles for (a)
$\textit{EL}=0.2$
,
$\textit{FL}=1.0$
(category II), (b)
$\textit{EL}=\textit{FL}=0.8$
(category I), (c)
$\textit{EL}=\textit{FL}=0.3$
(category III), (d)
$\textit{EL}=1.0$
,
$\textit{FL}=0.2$
(category IV). The 3-D trajectories are coloured by the instantaneous acceleration
$a_z$
, and their two-dimensional projections are coloured grey. Distances are normalised by
$D_{\textit{eq}}$
.
Figure 6. Various
${\textit{EL}}=\textit{FL}=0.3$
particle quantities at different density ratios
$R=1000$
, 100 and 10: (a) normalised vertical velocity
$v_z^*=v_z/v_T$
, where
$v_T$
is the terminal velocity; (b) normalised magnitude of lateral velocity
$v_\perp ^* = (v_x^2+v_y^2)^{0.5}/v_T$
; (c) pitch angle
$\varphi _L$
(
$L$
-axis angle from horizon); (d) roll angle
$\varphi _I$
(
$I$
-axis angle from horizon); (e) normalised pressure component of the drag force in the vertical direction,
$P_z^*=P_z/F_z$
, where
$F_z$
is the vertical component of drag force; (f) normalised magnitude of pressure components of the drag force in the lateral direction,
$P_\perp ^*=(P_x^2+P_y^2)^{0.5}/F_z$
. The normalised time
$T^*$
is calculated as
$tv_T/D_{\textit{eq}}$
. The transient dynamics of
$R=10$
and 100 particles show a rapid damping of angular oscillations and a significantly lower
$P_\perp ^*$
compared to the
$R=1000$
particle.
A key new observation from figure 5 is that at
$R=1000$
, the trajectories exhibit marked twists and turns, unlike at lower density ratios, for all the cases shown. The sequence of accelerating and decelerating intervals, as represented by the colours of the trajectories, also differs significantly at
$R=1000$
, especially when compared to trajectories with
$R=10$
. Similar to figure 4(c), different particle shapes/categories show clear influence on the lateral drift of the 3-D trajectories in figure 5. Near spherical ellipsoids (
$\textit{EL}=\textit{FL}=0.8$
of category I in figure 5
b) show less lateral drift than the other categories. Spheroids (
$\textit{EL}=0.2$
,
$\textit{FL}=1.0$
of category II in figure 5(a), and
$\textit{EL}=1.0$
,
$\textit{FL}=0.2$
of category IV in figure 5(d)) exhibit lateral drift in a planar two-dimensional motion with amplitudes larger than the ellipsoid with
$\textit{EL}=\textit{FL}=0.8$
. For the most non-spherical ellipsoid simulated here with
$\textit{EL}=\textit{FL}=0.3$
, which is a category III shape in figure 5(c), significant lateral drifts up to five times
$D_{\textit{eq}}$
combined with a complex 3-D motion are observed.
We further compare the key particle quantities from the DNS of the
$\textit{EL}=\textit{FL}=0.3$
ellipsoids at
$R=1000$
, 100 and 10, but at the same
${\textit{Re}}_{\!p} = 3$
in figure 6. As observed in figure 5(c), particles need a longer settling distance before reaching terminal velocity as
$R$
increases – for instance, at
$R = 1000$
, the required
$z^* = z/D_{\textit{eq}}$
is approximately 6.3 (3.6) times longer than at
$R = 10$
(
$R = 100$
), which is an expected behaviour. However, during the transition to terminal state, orientation oscillations and force histories deviate significantly at
$R=1000$
compared with
$R = 10$
or 100 (see figure 6). For example, oscillations and the progression of the transient motion to the final state also differ significantly for vertical
$v_T$
and lateral
$v_\perp$
velocity components as shown in figure 6(a,b), especially the latter showing strong fluctuations at
$R=1000$
. Figure 6(c,d) show that the damping of the amplitudes of the angular oscillations for the pitch angle
$\varphi _L$
(
$L$
-axis angle from horizon) and roll angle
$\varphi _I$
(
$I$
-axis angle from horizon) weakens as the density ratio
$R$
increases, leading to change from rapidly damped to under-damped oscillations. Another significant difference can be seen in figure 6(e,f) with regard to the contribution of the pressure component of the drag force. In particular, a large lateral pressure component of the drag force
$P_\perp ^*$
development is recorded for
$R = 1000$
, which also correlates with the lateral drifts during each twist and turn. These simulations support the conclusion that the ratio of particle to liquid density
$R$
is a key parameter governing the settling pattern, lateral dispersion, and evolution of velocity and hydrodynamic forces for free-falling, non-spherical particles, even at low Reynolds numbers considered in this study.
3.4. Detailed settling dynamics at
$R=1000$
We now turn our attention to the particle with
$\textit{EL}=\textit{FL}=0.3$
and
$R=1000$
to understand the fundamental physical processes behind twists and turns (figure 7
a,b) in the particle trajectories. This particle serves as a representative case for all other shapes, as it exhibits a combination of all the dynamics observed in other cases studied here, allowing us to study settling dynamics in more detail without analysing each case individually. Figure 7(a) shows that the particle follows a non-planar, nonlinear 3-D trajectory characterised by alternating phases of acceleration and deceleration. As depicted in figure 7(b), the particle rotates about all three axes, with pronounced pitch (i.e.
$L$
-axis angle from horizon) and yaw (i.e.
$S$
-axis angle from horizon) motions clearly visible in the top-view projection. The transient motion of the
$R=1000$
particle occurs in three phases (figure 7): (i) an initial acceleration phase – the particle accelerates from zero initial velocity without changing the initial orientation; (ii) a damped angular oscillation phase – the particle undergoes angular motion that results in lateral drift; and (iii) a terminal state – the particle stabilises to a terminal orientation with which it reaches the terminal velocity.
Figure 7. Settling dynamics of
$\textit{EL}=\textit{FL}=0.3$
particle at density ratio
$R = 1000$
in air. (a) The 3-D particle trajectory (coloured by the instantaneous acceleration
$a_z$
) and its two-dimensional projections (grey) at
${\textit{Re}}_{\!p}=3$
. Different transient phases are marked by: circle – start of angular motion, star – a pitch angle
$\varphi _L$
(
$L$
-axis angle from horizon) of
${45}^{\circ }$
, and square – onset of terminal state. (b) Top view of the particle trajectory together with the superimposed particle images. Time
$t$
evolutions of: (c) particle velocity normalised by terminal velocity
$\boldsymbol{v}^* = \boldsymbol{v}/v_T$
(
$v_z^*$
is the vertical component, and
$v_\perp ^*$
is the magnitude of the lateral components); (d) particle orientation (
$\varphi _L$
is the pitch angle, and
$\varphi _I$
is the roll angle –
$I$
-axis angle from horizon); and (e) the particle’s pressure components of the drag force
$\boldsymbol{P}^* = \boldsymbol{P}/F_z$
normalised by the vertical component of the drag force
$F_z$
(
$P_z^*$
and
$P_\perp ^*$
are respectively the vertical component and the magnitude of the lateral components). At
$\varphi _L={45}^{\circ }$
(
$\varphi _I={-17}^{\circ }$
): (f) the velocity distribution and streamlines of airflow around the particle (coloured based on
$v_z^*$
), and the surface distribution of
$v_\perp ^*$
; and (g) the surface distributions of
$P_\perp ^*$
and
$P_z^*$
. (h) The surface distribution of
$P_z^*$
at
$\varphi _L={0}^{\circ }$
and
${90}^{\circ }$
.
During the initial phase – until 0.05 s of time in the transient motion, which is marked with the purple circles in figure 7 – the particle retains its orientation while accelerating, resulting in no lateral drift or lateral velocity
$v_\perp ^*$
(figure 7
c). Its pitch
$\varphi _L$
and roll
$\varphi _I$
angles remain at their initial values
${90}^{\circ }$
and
${0}^{\circ }$
, respectively (figure 7
d). At this stage, the vertical
$P_z^*$
and lateral
$P_\perp ^*$
pressure components of the drag force are minimal, indicating that skin friction dominates (figure 7
e,h).
Beyond 0.05 s of time (purple circles in figure 7), the particle enters an angular oscillation phase, marked by rapid changes in pitch and roll angles (figure 7
d), and a sharp rise in lateral velocity
$v_\perp ^*$
– up to 14 % of the terminal vertical speed (figure 7
c). The proportion of pressure components of the drag force
$\boldsymbol{P}^*$
also increases significantly (figure 7
e). As the pitch angle passes
$\varphi _L = {45}^{\circ }$
, the vertical velocity
$v_z^*$
overshoots by approximately 30 %. At this orientation (grey stars in figure 7), the lateral surface velocity becomes notably uneven (figure 7
f), with a lateral pressure component of the drag force contribution
$P_\perp ^*$
reaching
${\sim}50\,\%$
of its vertical counterpart
$P_z^*$
(figure 7
e). This ensuing asymmetric distribution of pressure components of the drag force (figure 7
g) induces torques in both lateral and vertical directions, causing pronounced twisting and turning as the particle settles towards its terminal state in a nonlinear manner.
The particle’s orientation oscillates with two distinct frequencies: 26.9 Hz for pitch
$\varphi _L$
, and 51.8 Hz for roll
$\varphi _I$
, decaying exponentially at rates 21.2 Hz and 40.2 Hz, respectively (figure 7
d). Oscillations in
$v_\perp ^*$
,
$P_z^*$
and
$P_\perp ^*$
match the pitch frequency, while the roll angle has minimal effect. Here,
$P_\perp ^*$
is in phase with
$\varphi _L$
, whereas
$v_\perp ^*$
and
$P_z^*$
lag by
${90}^{\circ }$
. By 0.24 s of time in the transient motion (green squares), the particle settles into its terminal orientation – maximising cross-sectional area against gravity – and all parameters remain constant:
$v_z^* = 1.0$
and
$v_\perp ^* = 0.0$
,
$\varphi _L = \varphi _I = {0}^{\circ }$
, and
$P_z^* = 0.82$
and
$P_\perp ^* = 0.0$
. Comparisons across density ratios (figures 5
c and 6) show that at
${\textit{Re}}_{\!p}\approx 1$
, these complex trajectories and recurrent angular oscillations only arise at high
$R$
, despite the absence of flow separation or recirculation as shown in figure 7(f).
3.5. Angular dynamics at
$R=1000$
Figures 8 and 9 show the evolutions of the pitch, roll and yaw angle deviations from the horizon as functions of time for four particle shapes corresponding to those presented in figure 5, but only from
$R = 1000$
simulations. These angles are defined as
$\varphi _\varGamma = \sin ^{-1}(n_z)$
, where
$\boldsymbol{n} = (n_x, n_y, n_z)$
is the unit vector along axis
$\varGamma \in \{L, I, S\}$
. For all these particles, the initial orientation is given by
$\varphi _L = 90^\circ$
and
$\varphi _I = \varphi _S = 0^\circ$
, while the terminal orientation is
$\varphi _S = 90^\circ$
and
$\varphi _L = \varphi _I = 0^\circ$
. As a result,
$\varphi _L$
and
$\varphi _S$
should definitely change their initial values, while for
$\varphi _I$
, it is not necessary to change value to reach the terminal orientation.
Figure 8. Evolution of pitch
$\varphi _L$
, roll
$\varphi _I$
and yaw
$\varphi _S$
angles at density ratio
$R=1000$
for four ellipsoids with (a)
$\textit{EL}=0.2$
,
$\textit{FL}=1.0$
(category II), (b)
$\textit{EL}=\textit{FL}=0.8$
(category I), (c)
$\textit{EL}=\textit{FL}=0.3$
(category III) and (d)
$\textit{EL}=1.0$
,
$\textit{FL}=0.2$
(category IV) from the simulations.
Figure 9. Evolution of the angular velocities of the pitch
$\dot {\varphi }_L$
, roll
$\dot {\varphi }_I$
and yaw
$\dot {\varphi }_S$
angles at density ratio
$R=1000$
for four ellipsoids with (a)
$\textit{EL}=0.2$
,
$\textit{FL}=1.0$
(category II), (b)
$\textit{EL}=\textit{FL}=0.8$
(category I), (c)
$\textit{EL}=\textit{FL}=0.3$
(category III) and (d)
$\textit{EL}=1.0$
,
$\textit{FL}=0.2$
(category IV) from the simulations.
The most straightforward cases are the spheroids shown in figures 8(a) and 8(d), which exhibit a damped oscillatory transient as
$\varphi _L$
(or
$\varphi _S$
) progresses from
$90^\circ$
(or
$0^\circ$
) in perfect synchrony towards
$0^\circ$
(or
$90^\circ$
), while
$\varphi _I$
remains unchanged. For the triaxial ellipsoid with
${\textit{EL}} = {\textit{FL}} = 0.3$
, shown in figure 8(c), all three angles show damped oscillations towards the terminal state, albeit each in a very different manner. While the behaviour of
$\varphi _L$
is similar to that of the spheroids,
$\varphi _S$
is not in sync with it. Furthermore,
$\varphi _S$
approaches its terminal value in a peculiar way, bouncing back below
$90^\circ$
twice before starting to oscillate around it. On the other hand,
$\varphi _I$
exhibits damped oscillations around
$0^\circ$
throughout the whole process. In contrast, the triaxial ellipsoid with
${\textit{EL}} = {\textit{FL}} = 0.8$
shown in figure 8(b) exhibits yet another type of angular dynamics. Remarkably for this shape,
$\varphi _I$
flips to
$90^\circ$
during the first half of the transient, then relaxes back to
$0^\circ$
during the second half, both times through damped oscillations. For this particle,
$\varphi _L$
changes only during the first half, while
$\varphi _S$
changes only during the second half.
The angular dynamics shown in figure 9, as expected, depict the same trends discussed above; however, here the magnitude of the angular velocities reveals another aspect of the settling dynamics. In particular, it is evident that the highest angular velocity occurs for the oblate spheroid with
${\textit{EL}} = 1.0$
,
${\textit{FL}} = 0.2$
, for which both
$\dot {\varphi _L}$
and
$\dot {\varphi _S}$
reach 227 rad s−1 (figure 9
d), while for the prolate spheroid they peak at approximately 140 rad s−1 (figure 9
a). For a triaxial ellipsoid with
${\textit{EL}} = {\textit{FL}} = 0.8$
,
$\dot {\varphi _L}$
and
$\dot {\varphi _I}$
reach 96 rad s−1 in the first phase, while in the second phase,
$\dot {\varphi _I}$
and
$\dot {\varphi _S}$
reach a higher peak at approximately 115 rad s−1 (figure 9
b). Finally, for the triaxial ellipsoid with
${\textit{EL}} = {\textit{FL}} = 0.3$
,
$\dot {\varphi _L}$
,
$\dot {\varphi _I}$
and
$\dot {\varphi _S}$
peak at approximately 142, 106 and 151 rad s−1, respectively, with
$\dot {\varphi _S}$
being the first to reach its maximum.
3.6. Phase-diagrams of settling dynamics of ellipsoids at
$R=1000$
In this subsection, we systematically explore a broad range of ellipsoid shapes in both simulations and experiments, parametrised by elongation
$\textit{EL}$
and
$\textit{FL}$
(figure 10). Comprehensive datasets of key particle quantities – such as terminal velocity
$v_T$
, oscillation frequency
$f$
, decay rate of the angular oscillation
$\mu$
, and other relevant quantities for each particle shape – are tabulated in tables 2–5. The particle to fluid density ratio
$R$
is 1000 for all cases shown in these tables, containing both simulation and experimental observations. Figure 10 contains the contour plots of these tabulated particle quantities in the phase space of particle elongation
$\textit{EL}$
and flatness
$\textit{FL}$
.
Figure 10. In the phase space of particle elongation
$\textit{EL}$
and flatness
$\textit{FL}$
, the distributions of (a) the terminal velocity
$v_T$
, (b) the oscillation frequency of the angular orientation
$f$
, (c) the decay rate of the oscillation amplitudes
$\mu$
, (d) the terminal state contribution of pressure component of the drag force in the vertical direction
$P_T^*$
, (e) the normalised lateral drift
$\delta ^*$
, and (f) the normalised lateral swept volume
$\varPsi _\perp ^*$
. The scattered data are linearly interpolated using grid data in Matlab (R2023b). In the transient motion of a sphere,
$f$
and
$\mu$
are not present, which are shown as white triangles in the upper right corners of (b) and (c). For ellipsoid
$\textit{EL}=\textit{FL}=0.2$
, only
$v_T$
was measured. The diagrams in (a) and (b) show values from the experiments, while (c–f) are from the DNS. See table 1 for the main characteristics of the particles. Note that
$R = 1000$
for all the particles shown in this figure. Detailed tabulated data used to produce parts of this figure are in tables 2–5.
Table 2. Data on particle quantities 1:
$v_T^{D}$
is the terminal velocity of the particle from the DNS (the magnitude of vertical velocity when the absolute magnitude of vertical acceleration
$|a_z^{D}|$
reaches
$\leqslant {0.01}\,\textrm {m}\boldsymbol{\,}\textrm {s}^{-2}$
),
$\langle v_T^{E}\rangle$
is the terminal velocity from the experiments (the average of
$v_T^{E}$
from different experimental runs of each shape), and
$v_T^{E}$
from an experimental run is calculated by taking the average of the vertical velocities of the particle for the last 80 recorded frames when the absolute magnitude of the average vertical acceleration during these 80 frames,
$|\langle a_z^{E}\rangle |$
, reaches
$\leqslant {0.2}\,\textrm {m}\boldsymbol{\,}\textrm {s}^{-2}$
. The standard deviation in the statistics for
$v_T^{E}$
is given as
$\sigma (v_T^{E})$
, and the number of experimental runs used for each particle shape is denoted
$\#^{E}$
.
Table 3. Data on particle quantities 2:
$f_L^{D}$
and
$f_I^{D}$
are the frequencies of the angular oscillations of the pitch angle
$\varphi _L$
and roll angle
$\varphi _I$
, respectively, from the DNS,
$\langle f_L^{E}\rangle$
is the frequency of the angular oscillations of the pitch angle
$\varphi _L$
from the experiments (the average of
$f_L^{E}$
from different experimental runs of each shape), and
$f_L^{E}$
from an experimental run is calculated by fitting an exponentially decaying model
$\exp {(-\mu _L^{E}\times t+C)}$
that fits the peaks of the experimental
$\varphi _L$
with an
$R^2$
value
$\geqslant 0.8$
when the relative error between the fitted model and
$\varphi _L$
peaks from the experiments is
$\leqslant 20\,\%$
. Here,
$\mu _L^{E}$
is the decay rate of the oscillation amplitudes of
$\varphi _L$
. The standard deviation in the statistics for
$f_L^{E}$
is given as
$\sigma (f_L^{E})$
.
Table 4. Data on particle quantities 3:
$\mu _L^{D}$
and
$\mu _I^{D}$
are the decay rates of the oscillation amplitudes of the pitch angle
$\varphi _L$
and roll angle
$\varphi _I$
, respectively, from the DNS,
$\langle \mu _L^{E}\rangle$
is the decay rate of the oscillation amplitudes of the pitch angle
$\varphi _L$
from the experiments (the average of
$\mu _L^{E}$
from different experimental runs of each shape), and
$\mu _L^{E}$
from an experimental run is calculated by fitting an exponentially decaying model
$\exp {(-\mu _L^{E}\times t+C)}$
that fits the peaks of the experimental
$\varphi _L$
with an
$R^2$
value
$\geqslant 0.8$
when the relative error between the fitted model and
$\varphi _L$
peaks from the experiments is
$\leqslant 20\,\%$
. The standard deviation in the statistics for
$\mu _L^{E}$
is given as
$\sigma (\mu _L^{E})$
.
Table 5. Data on particle quantities 4:
$P_T^*$
is the contribution of the pressure component of the drag force to the vertical component of the drag force in the terminal state, obtained from the DNS,
$\delta$
is the integral of the lateral displacements (the lateral drift, obtained from the DNS),
$\varPsi _z^*$
is the vertical swept volume also obtained from the DNS, normalised by the swept volume of the volume-equivalent sphere, and
$\varPsi _\perp ^*$
is the lateral swept volume obtained from the DNS, normalised by the particle volume.
Over the examined phase space of particle shape, the terminal velocity
$v_T$
varies by a factor of 2.5 (figure 10
a), having similar sensitivity to both elongation and flatness, resembling ellipsoids in the Stokes regime (Bagheri & Bonadonna Reference Bagheri and Bonadonna2016). This factor – the normalised terminal velocity – is inversely proportional to the ratio of the particle’s surface area to that of its volume-equivalent sphere. Elongated or flattened ellipsoids fall with increased surface area in the direction of gravity in the terminal state that leads to increased air resistance, thus reducing
$v_T$
compared to the volume-equivalent sphere.
In contrast, the angular oscillation frequency
$f$
peaks for flattened, non-elongated particles, showing up to a twofold increase (figure 10
b). Consistent with the trends observed for pitch angular velocities in the preceding subsection, oblate shapes (category IV) exhibit the highest oscillation frequencies, peaking at approximately 40.6 Hz for the oblate spheroid with
${\textit{EL}}=1.0$
and
${\textit{FL}}=0.2$
. Specifically, the category IV particles have median oscillation frequency 37.6 Hz and standard deviation
$\sigma = {1.9}\,\textrm {Hz}$
. In comparison, categories II and III display moderate frequencies, with median values 29.6 Hz and 29.5 Hz, and standard deviations
$\sigma = {2.8}\,\textrm {Hz}$
and
$\sigma = {4.5}\,\textrm {Hz}$
, respectively. Category I, representing the near-spherical shapes, exhibits the lowest oscillation frequencies, with median 25.1 Hz and standard deviation
$\sigma = {6.2}\,\textrm {Hz}$
, indicating greater variability among the particles in this group.
The highest exponential decay rates, denoted by
$\mu$
, are observed for the oblate spheroids. However, across all ellipsoidal shapes considered,
$\mu$
ranged from 18.9 Hz to 24.2 Hz with standard deviation 1.4 Hz, indicating that the decay dynamics is only weakly dependent on the particle shape. For the three experimentally and numerically measurable quantities – i.e. terminal velocity
$v_T$
, oscillation frequency
$f$
and exponential decay rate
$\mu$
– we find close agreement between experiments and simulations, as shown in tables 2–5. Although the experimental set-up did not allow precise control of the initial release velocity, and almost no control of the particle orientation, this agreement shows that the transient dynamics for the considered set of shapes and Reynolds numbers is essentially independent of the initial conditions. Furthermore, the observed agreement also validates the experimental design, the data post-processing procedures, and the configuration of the numerical simulations, e.g. mesh resolution and domain size.
The DNS allow us to consistently compare quantities such as lateral drift and swept volume under a fixed initial condition – zero initial velocity and minimum projection area normal to gravity – and to analyse force distributions and secondary oscillations around the intermediate axis. The terminal-state pressure component of the drag force contribution,
$P_T^*$
, expressed as the fraction of the total vertical drag, shows a strong dependence on particle flatness (figure 10
d). It increases from 43 % for a sphere (
$\textit{EL} = \textit{FL} = 1.0$
) to 82 % for a highly flattened ellipsoid with
$\textit{EL} = \textit{FL} = 0.3$
. For comparison, the pressure component of the drag force contribution for a sphere in the Stokes regime (
${\textit{Re}}_{\!p} \ll 1$
) is approximately one-third (
$P_T^* \approx 33\,\%$
), increasing to 43 % at the larger Reynolds number investigated here (
${\textit{Re}}_{\!p} = 4.5$
). This trend is expected: as
${\textit{Re}}_{\!p}$
increases, inertial effects become more pronounced and increasingly influence the pressure field around the particle, while the relative contribution of viscous stresses, important for frictional drag, diminishes. For all non-spherical shapes examined,
$P_T^*$
exceeds that of the sphere, indicating that the pressure component of the drag force becomes more dominant as the particle becomes more non-spherical.
To analyse transient quantities such as the swept volume and lateral drift from DNS, we consider the entire transient motion of the particles from the fixed initial condition to the terminal state. In contrast to
$P_T^*$
, lateral drift
$\delta ^* = \delta /D_{\textit{eq}}$
depends more on elongation (figure 10
e), echoing experimental trends (figure 4
b) despite unavoidable initial-condition variations in experiments. The DNS reveal no drift for a sphere (
$\textit{EL}=\textit{FL}=1.0$
) but up to
$8D_{\textit{eq}}$
for
$\textit{EL}=0.2$
and
$\textit{FL}=0.5$
. We also quantify the volume swept by the particle as it reaches terminal velocity via the lateral swept volume
$\varPsi _\perp ^*$
, normalised by particle volume, which reaches 3.5 for elongated ellipsoids, and is zero for a sphere (figure 10
f). Meanwhile, the normalised vertical swept volume
$\varPsi _z^*$
is proportional to the terminal-state orientation of the particle, and increases from 1 (sphere) to 3.7 for the shapes investigated (table 5).
4. Discussion and concluding remarks
Using a combination of experimental observations and numerical simulations, it is demonstrated that even at particle Reynolds numbers as low as 2.1–4.5, the settling of heavy particles (i.e.
$ R = 1000$
) exhibits complex angular dynamics and fully 3-D trajectories. The fact that these behaviours are consistently observed in both experiments and simulations – with excellent quantitative agreement – strongly supports their reliability, effectively ruling out the influence of experimental bias or numerical artefacts. Moreover, the results presented here clearly confirm that these dynamics arise from the large particle-to-fluid density ratio, as corroborated by additional simulations conducted at
$ R = 10$
and
$ R = 100$
.
Various aspects of the settling behaviour, including particle trajectories, angular velocities, oscillation frequencies, forces acting on the particle, and flow structures, are presented and discussed in detail for representative cases. All non-spherical particles examined at
$ R = 1000$
exhibit damped oscillatory motion during the transient phase of their settling motion. However, while the trajectories remain planar for spheroidal particles, they become fully 3-D for some triaxial (non-axisymmetric) ellipsoids. This transition arises from the asymmetric distribution of hydrodynamic forces acting both vertically and in the lateral direction, which in turn influences the angular oscillations and overall trajectory. Whereas spheres exhibit no lateral drift, triaxial ellipsoids in our study are observed to drift laterally by up to ten times their volume-equivalent diameter. This leads to a substantial increase in the swept volume, which is directly linked to collision and encounter rates (Arguedas-Leiva et al. Reference Arguedas-Leiva, Słomka, Lalescu, Stocker and Wilczek2022), thereby enhancing the potential for particle aggregation in particle-laden flows.
Furthermore, the systematic investigation of ellipsoids of various shapes presented in this study allows us to comprehensively demonstrate the complex and nuanced dependence of settling behaviour on particle shape. The relationship between particle shape and key dynamic quantities is shown to be rich and highly non-trivial. For instance, while the settling velocity is similarly sensitive to both elongation and flatness along the symmetry axis (i.e.
$ {\textit{EL}} = {\textit{FL}}$
), the oscillation frequency exhibits a different pattern, displaying symmetry about the orthogonal diagonal, and larger for the flattened particles. The contribution of pressure components of the drag force is also found to be strongly governed by flatness. In contrast, lateral drift and the resulting swept volume are primarily influenced by particle elongation. These findings highlight that the dynamics of heavy particles settling in air and their dependence on particle shape are remarkably complex, even in the relatively low Reynolds number regime considered here, where instabilities such as wakes and recirculation zones are absent.
Revisiting figure 1, our findings suggest that the dynamics observed in this study may have significant implications for the behaviour of naturally occurring particles. In nature, far more pronounced non-sphericity (Magono & Lee Reference Magono and Lee1966; Um et al. Reference Um, McFarquhar, Hong, Lee, Jung, Lawson and Mo2015; Bagheri & Bonadonna Reference Bagheri and Bonadonna2016; Zhang et al. Reference Zhang, Kang, Allen, Allen, Gao and Sillanpää2020b ; Xiao et al. Reference Xiao, Cui, Brahney, Mahowald and Li2023; Tatsii et al. Reference Tatsii, Bucci, Bhowmick, Guettler, Bakels, Bagheri and Stohl2023) suggests even stronger lateral transport, particle encounters and aggregation, which ultimately affect atmospheric residence times. For example, typical ice crystals (Auer & Veal Reference Auer and Veal1970) and microplastics (Zhang et al. Reference Zhang, Kang, Allen, Allen, Gao and Sillanpää2020b ; Tatsii et al. Reference Tatsii, Bucci, Bhowmick, Guettler, Bakels, Bagheri and Stohl2023) can have much smaller flatness or elongation than studied here, potentially leading to even larger swept volumes and lateral drift. This mechanism may help to explain large-scale transport and loose-aggregate formation in volcanic ash and snowflakes. Turbulence further complicates settling orientation (Pumir & Wilkinson Reference Pumir and Wilkinson2016; Voth & Soldati Reference Voth and Soldati2017; Mathai, Lohse & Sun Reference Mathai, Lohse and Sun2020; Tinklenberg, Guala & Coletti Reference Tinklenberg, Guala and Coletti2025) and can amplify or weaken settling instabilities, making the dynamics of atmospheric particles even more complex. The full implications of the findings presented here, and the extent to which they may influence the dynamics of solid atmospheric particles, remain to be explored under real-world conditions, where atmospheric turbulence and the naturally occurring diversity of particle shapes are expected to play a significant role.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2025.10674.
Acknowledgements
This work was supported by the Max Planck Society. Scientific activities were carried out at the Max Planck Institute for Dynamics and Self-Organization (MPI-DS). Computational resources from the Max Planck Computing and Data Facility and MPI-DS are gratefully acknowledged.
Funding
T.B. was funded by the German Research Foundation (DFG) Walter Benjamin Position (project no. 463393443).
Declaration of interests
The authors report no conflict of interest.
Author contributions
Conceptualisation and research design: T.B., Y.W. and G.B. Experiments: T.B. Post-processing of the experimental data: T.B., Y.W. and G.B. Simulations: T.B. Post-processing of the simulation data: T.B. and Y.W. Analysis of the results: T.B., Y.W. and G.B. Development of simulation software: T.B., Y.W. and J.L. Visualisation: T.B., Y.W. and G.B. Interpretation of results: T.B., Y.W. and G.B. Writing – original draft: T.B., Y.W., J.L. and G.B. Writing – review and editing: T.B., Y.W., J.L. and G.B.
Open access
