# Inertial settling of an arbitrarily oriented cylinder in a quiescent flow :
from short-time to quasi-steady motion

Jean-Lou Pierson<EMAIL_ADDRESS>IFP Energies Nouvelles, Rond-point
de l’échangeur de Solaize, 69360 Solaize, France

###### Abstract

In this article, we investigate the inertial settling of an arbitrarily
oriented cylinder settling under gravity. We focus on two regimes: the very
short-time and long-time dynamic. By using the generalized Kirchhoff equations
to describe the particle motion, we demonstrate that during the very short
dynamic regime, a cylinder starting from rest behaves with sedimenting
velocities and angular velocity proportional to $t$ and $t^{3}$, respectively.
We then explore the long-time behaviour and evaluate the validity of the
quasi-steady assumption under which the fluid unsteady term can be neglected.
Using a dimensional analysis, we establish that the quasi-steady assumption is
only applicable to Reynolds numbers much smaller than one. However, by
comparing the results of quasi-steady models to recent experiments and direct
numerical simulations, we demonstrate that this assumption is valid for a
broader range of Reynolds numbers, particularly for long fibres. We also
analyze the effect of particle inertia. We show particle inertia plays no
significant role in the magnitude of the sedimenting velocities and angular
velocity. However, for sufficiently large inertia we reveal that the quasi-
steady model takes the form of a damped oscillator when the particle
approaches its equilibrium position, which is broadside on to its direction of
motion. We discuss the relevance of this solution in light of direct numerical
simulations.

##  I Introduction

The settling of anisotropic particles is a common occurrence in various
environmental flows, such as the fall of microplastics in the ocean [1] or the
precipitation of ice crystals in the atmosphere [2, 3]. Despite its practical
significance, the accurate modeling of anisotropic particle settling in
turbulent or quiescent environments is challenging due to the coupling between
the particle motion and the surrounding fluid flow. Unlike spherical
particles, the orientation of the body has a significant impact on the rate of
sedimentation. For example, in Stokes flow and for slender particles, the
sedimentation velocity of a particle with its axis aligned with gravity is
twice that of a particle with its axis perpendicular [4]. Additionally, for
inertial flows, the orientation of the body is coupled to the translational
equation of motion due to a non-zero hydrodynamic torque [5]. This results in
an unsteady problem, as the orientation of the body can change over time in
response to torques. In this article, we examine the settling of a cylindrical
particle in a quiescent flow as a first step in understanding the effect of an
anisotropic shape on particle motion.

The most general equations for studying the gravitational settling of a single
body in a quiescent fluid are the generalized Kirchhoff equations originally
derived by Howe [6]. Under this framework, added mass and vorticity
contributions to the hydrodynamic forces and torque are non-ambiguously
separated. However, for most configurations of practical interest, the
vorticity contributions cannot be expressed in closed form as they depend on
fluid motion history [7]. In the limit of negligible inertia and for a
spherical particle translating and rotating, the force and torque can be
decomposed into a quasi-steady component and a history term that takes the
form of integro-differential equations [8]. To date, there is no equivalent
analytical formula for an arbitrary axisymmetric particles, particularly
cylinders. The explanation lies in the complexity of the history term for a
non-spherical body, whose expression in the frequency domain is often too
complicated to allow a closed-form expression in the time domain [9, 10].
Moreover, unlike spherical particles, non-spherical particles have distinct
high and low frequency expressions for the history term Lawrence and Weinbaum
[11].

The situation is even worst for finite Reynolds numbers for which very few
results exist for non-spherical particles both for the quasi-steady and
history loads. In contrast to Stokesian flow, for which, due to the
reversibility of the Stokes equation, an axisymmetric particle with fore-aft
symmetry embedded in a uniform flow of velocity $U$ experiences no torque, an
inertial torque appears for finite Reynolds number which scales as $U^{2}$ [5,
12]. This torque naturally induces a coupling between translation and rotation
for a sedimenting cylinder. Hence as a cylinder sediment in a fluid with non-
negligible inertia, it rotates toward its equilibrium orientation, which is
broad-side on to its direction of motion Khayat and Cox [12]. There is another
nonlinear coupling term in the force balance for a rotating and translating
axisymmetric body with fore-aft symmetry which scales as $\Omega U$ where
$\Omega$ is the angular velocity Cox [5]. This term is at the origin of the
lift force on a spinning sphere translating perpendicularly to its rotation
axis [13]. This coupling term has not been studied so far in the context of a
rotating cylinder settling perpendicular to its rotation axis. The history
loads for non-spherical particles in the inertial regime have also received
limited attention, with only a few studies providing force expressions for
arbitrarily shaped particles in the long-time limit [14]. These expressions
require knowledge of the steady velocity field created by the particle in
Stokes flow, which is unknown for a moderately long cylinder. Nevertheless
scaling arguments indicate a $t^{-2}$ long-time decay of the history force in
the finite-inertia regime and a slower $t^{-1/2}$ decay in the Stokes regime
[14].

Based on the previous literature review, it is challenging to make analytical
advancements without additional assumptions, primarily due to the absence of
closed-form expressions for the history terms and $\Omega U$ load
contribution. The quasi-steady assumption, introduced by Cox [5] in his
investigation of the settling of a small-eccentricity spheroid, posits that
unsteady terms related to fluid motion can be ignored. Through a scaling
analysis, Cox [5] demonstrated that this assumption is appropriate as long as
the Reynolds number based on the body length is much smaller than unity. This
assumption has been widely utilized in various practical configurations,
including the settling of fibrous aerosols in quiescent air Newsom and Bruce
[15] and fibres in liquids Roy  _et al._ [16], yielding satisfactory results
when compared to experiments. Also, in those applications, the Reynolds number
was not necessarily small. Shin  _et al._ [17] have even shown that the quasi-
steady theory remains accurate for Reynolds number based on the body length
close to one. Most prior studies [15, 16] have made use of the quasi-steady
loads derived by Khayat and Cox [12] for slender fibers and the leading-order
hydrodynamic torque resisting rotation provided by the slender body theory
Batchelor [4]. Recent research has shown that the lift force and inertial
torque provided by Khayat and Cox [12] and the leading-order expression for
the hydrodynamic torque resisting rotation are not accurate for moderately
long rods [18, 19], in line with the qualitative but non-quantitative
agreement reported by Cabrera  _et al._ [20] between the Khayat and Cox [12]
theory and their experiments for moderately long rods. Hence, the validity of
the quasi-steady assumption must be re-evaluated with accurate formulas for
the loads, and the inertial correction to the loads proportional to $\Omega U$
must be considered. Scaling analysis by Cox [5] and Pierson  _et al._ [18] has
shown that this term may be small compared to the $U^{2}$ contribution when
the inertia effect is not significant. However, the magnitude of this term for
moderate inertia remains a topic of debate [18].

There is another limit where analytical progress is possible. In the limit of
time shorter than the viscous time scale, added mass effects dominate over
viscous contribution [21]. In this limit, the equations of motion can be
solved provided that the added mass coefficients are known. Loewenberg [9] has
established the added mass forces through the application of potential flow
solutions. However, the added mass torque coefficient for a rotating cylinder
has no known solution at present.

The primary objective of this article is to address two key issues in the
study of an arbitrarily oriented cylinder settling under gravity. Firstly,
analytical solutions are provided for the problem in the short-time limit
where the added mass effects are dominant. Secondly, the validity range of the
quasi-steady assumption and the neglect of the $\Omega U$ load contribution
assumptions are investigated as a function of the relevant dimensionless
parameters. These parameters include the cylinder aspect ratio, the Archimedes
number which is a Reynolds number based on gravitational velocity and the
density ratio. The analysis is conducted using the generalized Kirchhoff
equations as proposed by Howe [6] and Mougin and Magnaudet [21], and the
quasi-steady load expressions derived in Fintzi  _et al._ [19], Pierson  _et
al._ [18] for moderately long rods and Khayat and Cox [12] for very elongated
fibers. The proposed model is validated through scaling analysis, experimental
measurements from Roy  _et al._ [16] and Cabrera  _et al._ [20] and direct
numerical simulations. The article is structured as follows. The governing
equations are presented in Section II. In Section III, the analytical
solutions for the short-time limit are described. The quasi-steady models are
derived in Section IV. The comparison between the quasi-steady models and
experimental measurements from Roy  _et al._ [16] and Cabrera  _et al._ [20]
as well as direct numerical simulations are discussed in Section V. Section VI
contains a discussion on the validity of the quasi-steady model as well as our
conclusions.

##  II Governing equations

$\mathbf{p}$$\mathbf{q}$$\mathbf{r}$$\mathbf{U}$$\theta$$\phi$$\mathbf{g}$$z$$y$
Figure 1: Finite-length cylinder submitted to the gravity acceleration
$\mathbf{g}$.

We consider a finite-length cylinder of length $L$ and diameter $D$ settling
under gravity with velocity $\mathbf{U}$ and angular velocity
$\mathbf{\Omega}$ (Figure 1). The main difficulty in studying this problem is
the coupling between the body motion and the surrounding flow field which
satisfies Navier-Stokes equations [7]. The equations of motion expressed in a
reference frame having its origin fixed with respect to the laboratory but
axes rotating with the body read [6, 21]

$(m\mathbb{I}+\mathbb{A})\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}+\boldsymbol{\Omega}\times((m\mathbb{I}+\mathbb{A})\mathbf{U})=\mathbf{F}^{\omega}+(m-\rho
V)\mathbf{g},$ (1)

$(\mathbb{J}+\mathbb{D})\frac{\mathrm{d}\boldsymbol{\Omega}}{\mathrm{d}t}+\boldsymbol{\Omega}\times((\mathbb{J}+\mathbb{D})\boldsymbol{\Omega})+\mathbf{U}\times(\mathbb{A}\times\mathbf{U})=\mathbf{T}^{\omega},$
(2)

where $m$ and $V$ are respectively the cylinder mass and volume, $\mathbb{I}$
is the identity matrix, $\mathbb{J}$ is the inertia tensor, $\mathbb{A}$ and
$\mathbb{D}$ are the second-order added inertia tensors. Using indicial
notation all those tensors, can be written as
$J_{ij}=J_{p}p_{i}p_{j}+J_{q}(\delta_{ij}-p_{i}p_{j})$,
$A_{ij}=A_{p}p_{i}p_{j}+A_{q}(\delta_{ij}-p_{i}p_{j})$ and
$D_{ij}=D_{q}(\delta_{ij}-p_{i}p_{j})$ where $\mathbf{p}$ is the unit
orientation vector and $\mathbf{q}$ is the unit vector perpendicular to
$\mathbf{p}$ in the plane $(y,z)$ (Figure 1). If the components of the inertia
tensor can be readily obtained in a closed form for a finite-length cylinder
($J_{p}=mD^{2}/8$, $J_{q}=m(3D^{2}/4+L^{2})/12$) this is not the case for the
components of the added mass tensors $A_{p},A_{q}$ and $D_{q}$. If $A_{p}$ and
$A_{q}$ have been already studied in the litterature [9], to the best of the
author knowledge, no expression for $D_{q}$ has not been published yet. Based
on the direct numerical simulation results of Kharrouba [22] and the potential
flow results of Loewenberg [9] we derive correlations for $A_{p}$, $A_{q}$ and
$D_{q}$ in Appendix A. In equations 1 and 2, $\mathbf{F}^{\omega}$ and
$\mathbf{T}^{\omega}$ are the force and torque due to the vorticity in the
flow. Except in inertia dominated regimes the motion of the cylinder is planar
in the $(p,q)$ plane [20] and the equation of motion simplifies to

$\displaystyle(m+A_{p})\frac{dU_{p}}{dt}-(m+A_{q})\Omega_{r}U_{q}$
$\displaystyle=F_{p}^{\omega}+(m-\rho V)g\cos\phi,$ (3)
$\displaystyle(m+A_{q})\frac{dU_{q}}{dt}+(m+A_{p})\Omega_{r}U_{p}$
$\displaystyle=F_{q}^{\omega}-(m-\rho V)g\sin\phi,$ (4)
$\displaystyle(J_{q}+D_{q})\frac{d\Omega_{r}}{dt}$
$\displaystyle=-U_{p}U_{q}(A_{q}-A_{p})+T_{r}^{\omega}.$ (5)

The physical origin of the torque can be exemplified by looking more closely
at equation 5. In steady potential flow, a torque manifests on a cylindrical
particle as indicated by the first term on the right-hand side of equation 5
[23]. As the aspect ratio, $\chi=L/D$, exceeds 1, $A_{q}>A_{p}$ (see Appendix
A), and the torque is positive, orienting the body broad-side on. As a result,
to leading order the torque expression provided by Khayat and Cox [12] for
small Reynolds number and large aspect ratio ($\chi\gg 1$) $T_{r}=-5\pi/24\rho
U_{p}U_{q}L^{3}/\ln^{2}(\chi)$ is the sum of two positive contribution due to
potential flow and vorticity. For $\chi\gg 1$ the potential contribution
$-U_{p}U_{q}(A_{q}-A_{p})$ scales as $-\rho\pi D^{2}U_{p}U_{q}L$ and is thus
negligible to leading order in comparison to the vorticity contribution. The
situation is less obvious for moderately large aspect ratio ($\chi\approx 2$)
for which both contributions may have the same order of magnitude since the
total torque on the body do not scales as $\rho U_{p}U_{q}L^{3}/\ln^{2}(\chi)$
[19].

##  III Short-time dynamics

In the high-frequency limit or equivalently for time shorter than the
diffusive scale added mass effect dominates over viscous contribution. This
may be proved by deriving the unsteady loads in the Stokes regime Kabarowski
and Khair [10], but also by performing a short time-analysis of the Navier-
Stokes equation Mougin and Magnaudet [21]. The latter is more general as it is
not limited to the Stokes flow regime. The proper length scale to be used in
the diffusive scale for the problem at hand is unknown a priori. However,
Kabarowski and Khair [10] have shown that for both transverse and longitudinal
oscillations, the added mass contributions dominate over the history load if
$t\ll D^{2}/\nu$. Hence, the proper length scale in the short dynamic is the
cylinder diameter, at least in the Stokes flow regime. Assuming $t\ll
D^{2}/\nu$, equations 3 \- 5 simplify to

$\displaystyle(m+A_{p})\frac{dU_{p}}{dt}-(m+A_{q})\Omega_{r}U_{q}$
$\displaystyle=(m-\rho V)g\cos\phi,$ (6)
$\displaystyle(m+A_{q})\frac{dU_{q}}{dt}+(m+A_{p})\Omega_{r}U_{p}$
$\displaystyle=-(m-\rho V)g\sin\phi,$ (7)
$\displaystyle(J_{q}+D_{q})\frac{d\Omega_{r}}{dt}$
$\displaystyle=-U_{p}U_{q}(A_{q}-A_{p}),$ (8) $\displaystyle\frac{d\phi}{dt}$
$\displaystyle=\Omega_{r},$ (9)

where we made explicit the equation ruling the dependency of $\phi$ with time.
To the best of our knowledge, there is no closed form analytical solution to
this non-linearly coupled system of equations, although there is a
straightforward analytical treatment in the case of zero gravity Lamb [24]. We
have to solve this system numerically but one may get insightful estimates
using asymptotic analysis. By balancing the acceleration of gravity with the
particle acceleration in equations 6 and 7 one get $U\sim gT$ where $U$ is a
characteristic velocity scale, and $T$ a characteristic timescale. Injecting
this scaling in equation 8, and since $(A_{q}-A_{p})/(I_{q}+D_{q})\sim
1/L^{2}$ for $\chi\gg 1$ one get the characteristic angular velocity
$\Omega\sim g^{2}T^{3}/L^{2}$. By using this estimates in equations 6, 7, 8
and 9 and defining the dimensionless (starred) quantities as
$U_{p}=UU_{p}^{*}$, $U_{q}=UU_{q}^{*}$ and $\Omega_{r}=\Omega^{*}\Omega_{r}$
we obtain

$\displaystyle\frac{dU_{p}^{*}}{dt^{*}}-\epsilon\frac{\mathcal{A}}{\mathcal{B}}\Omega_{r}^{*}U_{q}^{*}$
$\displaystyle=\mathcal{A}\cos\phi,$ (10)
$\displaystyle\frac{dU_{q}^{*}}{dt^{*}}+\epsilon\frac{\mathcal{B}}{\mathcal{A}}\Omega_{r}^{*}U_{p}^{*}$
$\displaystyle=-\mathcal{B}\sin\phi,$ (11)
$\displaystyle\frac{d\Omega_{r}^{*}}{dt^{*}}$
$\displaystyle=-\mathcal{C}U_{p}^{*}U_{q}^{*},$ (12)
$\displaystyle\frac{d\phi}{dt^{*}}$ $\displaystyle=\epsilon\Omega_{r}^{*},$
(13)

where $\mathcal{A}=(\bar{\rho}-1)/(\bar{\rho}+A_{p}^{*})$,
$\mathcal{B}=(\bar{\rho}-1)/(\bar{\rho}+A_{q}^{*})$,
$\mathcal{C}=(A_{q}^{*}-A_{p}^{*})/(\bar{\rho}J_{q}^{*}+D_{q}^{*})$,
$A_{p}^{*}=A_{p}/(\rho V),A_{q}^{*}=A_{q}/(\rho V)$, $J_{q}^{*}=J_{q}/(\rho
VL^{2})$, $D_{q}^{*}=D_{q}/(\rho VL^{2})$, $\bar{\rho}=\rho_{p}/\rho$ is the
density ratio and $\epsilon=g^{2}T^{4}/L^{2}$ can be understood as the ratio
of the characteristic timescale over a gravity time scale. In the following,
we consider the small time limit $T\ll(L/g)^{1/2}$ or equivalently
$\epsilon\ll 1$. We seek for solutions of equations 10 \- 13 in the form of
asymptotic expansions in powers of the small parameter
$U_{p}^{*}=U_{p}^{*(0)}+\epsilon U_{p}^{*(1)}+...$,
$U_{q}^{*}=U_{q}^{*(0)}+\epsilon U_{q}^{*(1)}+...$,
$\Omega_{r}^{*}=\Omega_{r}^{*(0)}+\epsilon\Omega_{r}^{*(1)}+...$ and
$\phi=\phi^{(0)}+\epsilon\phi^{(1)}+...$. The calculations are straightforward
and are detailed in appendix B. The solutions up to the order 1 with respect
to the small parameter $\epsilon$ are

$\displaystyle U_{p}^{*}$
$\displaystyle=\mathcal{A}t^{*}\cos\phi^{(0)}-\frac{\epsilon}{12}\mathcal{A}^{2}\mathcal{B}\mathcal{C}t^{*5}\cos\phi^{(0)}\sin^{2}\phi^{(0)}+\mathcal{O}(\epsilon^{2}),$
(14) $\displaystyle U_{q}^{*}$
$\displaystyle=-\mathcal{B}t^{*}\sin\phi^{(0)}-\frac{\epsilon}{12}\mathcal{A}\mathcal{B}^{2}\mathcal{C}t^{*5}\cos^{2}\phi^{(0)}\sin\phi^{(0)}+\mathcal{O}(\epsilon^{2}),$
(15) $\displaystyle\Omega_{r}^{*}$
$\displaystyle=\frac{\mathcal{A}\mathcal{B}\mathcal{C}}{3}t^{*3}\cos\phi^{(0)}\sin\phi^{(0)}+\mathcal{O}(\epsilon^{2}),$
(16) $\displaystyle\phi$
$\displaystyle=\phi^{(0)}+\epsilon\frac{\mathcal{A}\mathcal{B}\mathcal{C}}{12}t^{*4}\cos\phi^{(0)}\sin\phi^{(0)}+\mathcal{O}(\epsilon^{2}).$
(17)

(a) (b) (c)

(d) (e) (f)

(g) (h) (i)

Figure 2: Dimensionless velocities and angular velocity as function of time
for a $\chi=10$ cylinder starting from rest with $\phi=30^{\circ}$ and
$\epsilon=1$. (a), (b), (c) : $\bar{\rho}=1.5$. (d), (e), (f) :
$\bar{\rho}=10$. (g), (h), (i) : $\bar{\rho}=1000$. $-$ : numerical solution
of equations 10, 11, 12 and 13, $--$ : zero-th order asymptotic expansion,
$-\cdot-$ : first order asymptotic expansion.

Figure 2 displays the numerical and analytical solutions for a cylinder of
aspect ratio $10$ released with an orientation angle $\phi=30^{\circ}$. We
impose $\epsilon=1$ such that $T=(L/g)^{1/2}$. The numerical solution is
obtained by using a Runge-Kutta 4 algorithm. For $\bar{\rho}=1.5$ a good
agreement is observed between the numerical and solution and the zeroth order
analytical solution up to $t^{*}\approx 1.5$ and up to $t^{*}\approx 2.5$ for
the first order solution (Figures 2 (a), (b) and (c)). For moderately large
density ratios the agreement with the numerical solution is slightly worse.
However, for the largest density ratios the first-order analytical solutions
remain valid up to $t^{*}\approx 6$ (Figures 2 (g), (h) and (i)). For all the
density ratios one may observe the good agreement between the analytical and
numerical solutions out of the range of applicability of the analytical
solution. Also, this agreement seems to be dependent on $\bar{\rho}$. The
explanation lies in the magnitude of the coupling terms between translation
and rotation in equations 10 and 11. At first-order this coupling terms are
found to be proportional to $\mathcal{A}^{2}\mathcal{B}\mathcal{C}$ and
$\mathcal{A}\mathcal{B}^{2}\mathcal{C}$ in equations 14 and 15.

Figure 3: Leading order terms coupling translation and rotation. $-$ :
$\chi=2$, $--$ : $\chi=8$, $-\cdot-$ : $\chi=32$, $\cdot\cdot\cdot$ :
$\chi\rightarrow\infty$.

Figure 3 illustrates the variation of $\mathcal{A}^{2}\mathcal{B}\mathcal{C}$
and $\mathcal{A}\mathcal{B}^{2}\mathcal{C}$ as a function of $\chi$ and
$\bar{\rho}$. Both coupling terms have the same order of magnitude and are
smaller than approximatively $0.8$. The behaviour of
$\mathcal{A}^{2}\mathcal{B}\mathcal{C}$ and
$\mathcal{A}\mathcal{B}^{2}\mathcal{C}$ as function of $\chi$ and $\bar{\rho}$
disserve also some comments. Since both quantity evolves comparably, we will
focus on $\mathcal{A}^{2}\mathcal{B}\mathcal{C}$. This quantity increases as a
function of $\chi$ for moderate $\chi$. Hence one may expect a better match of
the theory for small cylinders. Then $\mathcal{A}^{2}\mathcal{B}\mathcal{C}$
becomes nearly independent of $\chi$ for $\chi\approx 20$. This is in line
with the results of appendix A in which the added mass coefficients become
almost independent of $\chi$ for this aspect ratio.
$\mathcal{A}^{2}\mathcal{B}\mathcal{C}$ behaves non-monotonously with
$\bar{\rho}$. It increases for $\bar{\rho}\leq 5$ and then decreases for
larger $\bar{\rho}$. The value of $\bar{\rho}$ for which we observed the
maxima is weakly dependent on $\chi$. All these trends can be easily obtained
by looking more closely at the behaviour of
$\mathcal{A}^{2}\mathcal{B}\mathcal{C}$ as a function of $\bar{\rho}$. First
one may note that for $\bar{\rho}\rightarrow 1$ we have
$\mathcal{A}^{2}\mathcal{B}\mathcal{C}\sim(\bar{\rho}-1)^{3}(A_{q}^{*}-A_{p}^{*})/[(1+A_{p}^{*})^{2}(1+A_{q}^{*})(J_{q}^{*}+D_{p}^{*})]$
which explains the strong increase of $\mathcal{A}^{2}\mathcal{B}\mathcal{C}$
for $\bar{\rho}\leq 5$. In the opposite limit $\bar{\rho}\gg 1$ one obtain
$\mathcal{A}^{2}\mathcal{B}\mathcal{C}\sim
1/\bar{\rho}\times(A_{q}^{*}-A_{p}^{*})/J_{q}^{*}$ which explains the decay of
$\mathcal{A}^{2}\mathcal{B}\mathcal{C}$ for large density ratio.

It is also interesting to discuss the range of applicability of the asymptotic
results in experiments. Our observations suggest that the asymptotic expansion
is valid up to $t\sim(L/g)^{1/2}$ and even larger value for large density
ratios. Thus one may expect this solution to be valid as long as
$(L/g)^{1/2}\ll D^{2}/\nu$ which can be written $Ar\gg(\bar{\rho}-1)\chi$
where $Ar=(\rho_{p}-\rho)\rho gD^{3}/\mu^{2}$ is the Archimedes number. This
implies that for sufficiently inertial regimes, the proposed analytical
solution is expected to be observed in experiments. We are not aware of any
experiments specifically dedicated to testing the validity of the theory.
However, the qualitative observations provided in Toupoint  _et al._ [25] may
provide a first comparison with the present results. Toupoint  _et al._ [25]
investigated the settling of cylinder in inertia dominated regime $Ar\approx
40000$. Their density ratio was fixed to $\bar{\rho}=1.16$ and the elongation
ratio of the cylinders ranged between $2$ and $20$. In particular, the time
related to viscous diffusion in their experiments with a cylinder of diameter
$2$mm and $\chi=5$ is $t_{\nu}\approx 4s$ while $T\approx 0.03s$. Hence the
present theory should apply to their experimental results. They observed that
the distance needed for the cylinders to rotate to their equilibrium position
was less than $10d$ even when the cylinder was released vertically. From our
analytical results, we obtain the time needed for the cylinder to change its
orientation from nearly vertical ($\phi^{(0)}\ll 1$) to horizontal
($\phi(t^{*})=\pi/2$)
$t\sim(L/g)^{1/2}(6\pi/(\phi^{(0)}\mathcal{A}\mathcal{B}\mathcal{C}))^{1/4}$.
With $\phi^{(0)}=\pi/10$ one get $t\approx 0.2s$. Injecting this estimate in
the velocity equation and integrating one obtains the normalised distance
$l/D$ for which the cylinder changes its orientation:
$l/D=\mathcal{A}gt^{2}/2\approx 11.2$. This estimate is very close to the
experimental prediction and gives confidence on the non-viscous flow origin on
which this model is grounded.

##  IV Quasi-steady models

The particle equations of motion (3 \- 5) are coupled to the Navier-Stokes
equations via the boundary conditions and the hydrodynamic loads. This make
the problem very hard to solve and simplifying assumptions are required to
make analytical progress as explained in the introduction. Since the
characteristic time needed for the body to change its orientation scales like
$\Omega^{-1}$, the unsteady term in the Navier-Stokes equation scales as $\rho
U\Omega$ where $\Omega$ and $U$ may be taken as the nominal scales for the
angular and settling velocities. Since the inertial term scales as $\rho
U^{2}/L$ the unsteady term is negligible in comparison to the inertial term as
long as $\Omega L/U\ll 1$. This condition is satisfied when the time needed
for the vorticity to diffuse from the body for a Reynolds number of unity is
much smaller than the rotation time scale [15]. In this limit, one may
consider the fluid unsteady term and as a consequence the history loads to be
negligible. To be fully consistent with the neglect of the unsteady term we
will also neglect the added mass loads. It remains to consider the inertial
correction to the loads proportional to $\Omega U$. One may expect the $\Omega
U$ corrections to be smaller than the inertial $U^{2}$ corrections as long as
$\Omega L/U\ll 1$. Hence under the assumption, we postulate that both the
history terms and the $\Omega U$ contribution are negligible with respect to
the other load contribution as long as $\Omega L/U\ll 1$. The choice of the
length of the particles ensures that the particle length is the relevant scale
in the limit of a slender fibre [12, 26]. The validity of the assumption
$\Omega L/U\ll 1$ will be evaluated through comparisons with experimental
results and direct numerical simulations in Section V.

###  IV.1 Quasi-steady models for an arbitrarily oriented cylinder

Under the assumption $\Omega L/U\ll 1$, the loads can be approximated as their
quasi-steady counterparts, disregarding the $\Omega U$ terms:
$F_{p}^{\omega}\approx F_{p}$, $F_{q}^{\omega}\approx F_{q}$ and
$-U_{p}U_{q}(A_{q}-A_{p})+T_{r}^{\omega}\approx T_{r}^{i}+T_{r}^{\Omega}$
where $F_{p}$ and $F_{q}$ are the quasi-steady forces in the longitudinal and
perpendicular directions to the cylinder. The inertial torque, $T_{r}^{i}$,
drives the cylinder towards its equilibrium position, while the hydrodynamic
torque, $T_{r}^{\Omega}$, resists rotation. We will distinguish two different
configurations depending on the cylinder aspect ratio: moderately long rods
$2\leq\chi\leq 30$ and long fibres $\chi>30$. For $\chi\leq 30$, expressions
derived through slender-body theory and direct numerical simulations by
Kharrouba  _et al._ [27], Pierson  _et al._ [18], Fintzi  _et al._ [19] will
be used, as they provide more accurate predictions compared to those by Khayat
and Cox [12]. Conversely, the expression by Khayat and Cox [12] will be used
for $\chi>30$. To the author’s knowledge, finite-inertia effects have not been
considered for a slender body rotating in a fluid at rest at infinity. Thus,
the most accurate expression for this case under the Stokes flow assumption by
Pierson  _et al._ [18] will be used. The expressions for $F_{p}$, $F_{q}$,
$T_{r}^{i}$, and $T_{r}^{\Omega}$ are outlined in Appendix C. In both
configurations, the linearized approximation is employed to describe the force
as a function of the particle velocity. This approximation, which is exact in
the Stokes flow regime, has been demonstrated to be accurate up to $Re\approx
1$ for $10\leq\chi\leq 30$ and up to $Re\approx 10$ for $\chi<10$ by Fintzi
_et al._ [19] where $Re=\rho DU/\mu$ is the Reynolds number based on the body
diameter. The experimental works of Lopez and Guazzelli [28] with
$11.5\leq\chi\leq 34.5$ and Roy  _et al._ [16] with $20\leq\chi\leq 100$ have
also validated its relevance for larger aspect ratios.

Equations 3 \- 5 are normalized by defining dimensionless (starred) quantities
as $U_{p}=UU_{p}^{*}$, $U_{q}=UU_{q}^{*}$, $\Omega_{r}=\Omega\Omega^{*}$ and
$t=\Omega^{-1}t^{*}$ where the characteristic angular and velocity scales are
a priori unknown. In the small inertia limit making use of the linearized
approximation, the forces can be expressed as $F_{p}=-\mu
UU_{p}^{*}LF_{p}^{*}(Re_{L}^{*},\chi)$ and $F_{q}=-\mu
UU_{q}^{*}LF_{q}^{*}(Re_{L}^{*},\chi)$ where the expressions for $F_{p}$ and
$F_{q}$ can be obtained from appendix C.
$Re_{L}^{*}=Re_{L}(U_{p}^{*2}+U_{q}^{*2})^{1/2}$ is the Reynolds number based
on the instantaneous settling velocity while $Re_{L}=\rho UL/(2\mu)$ is the
characteristic Reynolds number based on the body half length. The inertial
torque reads $-\rho
U^{2}U_{p}^{*}U_{q}^{*}L^{3}T_{i}^{*}(Re_{L}^{*},\chi,\theta)$ while the
torque resisting rotation reads
$-\mu\Omega\Omega_{r}^{*}L^{3}T_{\Omega}^{*}(Re_{\Omega}^{*},\chi)$ where
$Re_{\Omega}^{*}=Re_{\Omega}|\Omega_{r}^{*}|$ and $Re_{\Omega}=\rho\Omega
D^{2}/\mu$. Injecting all those scaling in equations 3, 4 and 5 one obtain

$\displaystyle\bar{\rho}Re_{\Omega}\left(\frac{dU_{p}^{*}}{dt^{*}}-\Omega_{r}^{*}U_{q}^{*}\right)$
$\displaystyle=-\frac{4}{\pi}F_{p}^{*}(Re_{L}^{*},\chi)U_{p}^{*}+\frac{(\rho_{p}-\rho)gD^{2}}{\mu
U}\cos\phi,$ (18)
$\displaystyle\bar{\rho}Re_{\Omega}\left(\frac{dU_{q}^{*}}{dt^{*}}+\Omega_{r}^{*}U_{p}^{*}\right)$
$\displaystyle=-\frac{4}{\pi}F_{q}^{*}(Re_{L}^{*},\chi)U_{q}^{*}-\frac{(\rho_{p}-\rho)gD^{2}}{\mu
U}\sin\phi,$ (19)
$\displaystyle\bar{\rho}\left(\frac{Re_{\Omega}}{Re}\right)^{2}J_{q}^{*}\frac{d\Omega_{r}^{*}}{dt^{*}}$
$\displaystyle=-\frac{4}{\pi}\left(T_{i}^{*}(Re_{L}^{*},\chi,\theta)U_{p}^{*}U_{q}^{*}+\frac{\mu\Omega}{\rho
U^{2}}T_{\Omega}^{*}(Re_{\Omega}^{*},\chi)\Omega_{r}^{*}\right).$ (20)

From the above set of equations and for moderate density ratios
$\bar{\rho}\sim 1$ it appears that all the unsteady terms vanish if
$Re_{\Omega}\ll 1$ and $Re_{\Omega}\ll Re$. In contrast for large density
ratios, one may take into account the particle inertia [15]. We have to recall
that to derive equations 18 \- 20 we have to assume $\Omega L/U\ll 1$, which
is equivalent to $Re_{\Omega}\ll\chi Re$. This a rather restrictive condition
for the applicability of the model but we shall see hereafter that it is valid
for a non-negligible range of dimensionless parameters. Since the last terms
in equation 18 \- 20 must be of order one the velocity and angular velocity
scales can be readily obtained as $U\sim D^{2}(\rho_{p}-\rho)g/\mu$ and
$\Omega=\rho D^{4}(\rho_{p}-\rho)^{2}g^{2}/\mu^{3}$. Inserting these scalings
yields

$\displaystyle\bar{\rho}Ar^{2}\left(\frac{dU_{p}^{*}}{dt^{*}}-\Omega_{r}^{*}U_{q}^{*}\right)$
$\displaystyle=-\frac{4}{\pi}F_{p}^{*}(Re_{L}^{*},\chi)U_{p}^{*}+\cos\phi,$
(21)
$\displaystyle\bar{\rho}Ar^{2}\left(\frac{dU_{q}^{*}}{dt^{*}}+\Omega_{r}^{*}U_{p}^{*}\right)$
$\displaystyle=-\frac{4}{\pi}F_{q}^{*}(Re_{L}^{*},\chi)U_{q}^{*}-\sin\phi,$
(22) $\displaystyle\bar{\rho}Ar^{2}J_{q}^{*}\frac{d\Omega_{r}^{*}}{dt^{*}}$
$\displaystyle=-\frac{4}{\pi}\left(T_{i}^{*}(Re_{L}^{*},\chi,\theta)U_{p}^{*}U_{q}^{*}+T_{\Omega}^{*}(Re_{\Omega}^{*},\chi)\Omega_{r}^{*}\right).$
(23)

The dimensionless nonlinear system of coupled differential equations,
represented by Equations 21 \- 23, describe the velocity and angular velocity
of a cylinder under an external force. Analytical solutions to this system can
only be obtained in the limit of low Archimedes number, as reported by Cox
[5]. Otherwise, numerical solutions must be sought. Note that
$\bar{\rho}Ar^{2}$ can be seen as a Stokes number ($St$) which is a
dimensionless measure of particle inertia. The system exhibits two distinct
regimes depending on $St$. For moderate particle inertia ($St\sim 1$ or
equivalently $\bar{\rho}\sim 1/Ar^{2}$), as in the case of small solid
particles settling in air, the left-hand sides (LHS) of equations 21 \- 23 are
non-negligible. The resulting system then consists of three ordinary
differential equations which are solved thanks to a Runge-Kutta 4 algorithm.
Conversely, for small particle inertia ($St\ll 1$) as for small particles with
densities close to that of the fluid the (LHS) can be neglected. The
simplified system is a non-linear function composed of three equations, which
can be solved using classical root finding methods. It reads

$\displaystyle 0$
$\displaystyle=-\frac{4}{\pi}F_{p}^{*}(Re_{L}^{*},\chi)U_{p}^{*}+\cos\phi,$
(24) $\displaystyle 0$
$\displaystyle=-\frac{4}{\pi}F_{q}^{*}(Re_{L}^{*},\chi)U_{q}^{*}-\sin\phi,$
(25) $\displaystyle 0$
$\displaystyle=T_{i}^{*}(Re_{L}^{*},\chi,\theta)U_{p}^{*}U_{q}^{*}+T_{\Omega}^{*}(Re_{\Omega}^{*},\chi)\Omega_{r}^{*}.$
(26)

The validity of the two sets of equations is restricted to the condition
$\Omega L/U\ll 1$. Since $\Omega L/U\sim Ar\chi$ we obtain the very
restrictive condition $Ar\ll 1/\chi$. The next section aim to demonstrate that
the present set of equations remain valid up to moderate Archimedes numbers,
despite its seemingly limited range of validity. It is also important to note
that these scalings are a priori limited to small fluid inertial effects as
the forces and resistive torque are viscous in nature. However, Fintzi  _et
al._ [19] and Pierson  _et al._ [18] have demonstrated that this viscous based
laws remain valid for moderate inertial effects ($Re_{L}\leq 1$). The behavior
for larger inertial effects is complex and non-linearly dependent on the
Reynolds number, as described in Appendix C.

###  IV.2 Damped oscillations around $\phi=\pi/2$

Upon reaching its equilibrium position $\phi=\pi/2$, which is perpendicular to
its direction of motion, an additional regime of interest can be expected for
the particle. We linearize the equation of motion around this position. By
writing $\phi(t)=\pi/2+\epsilon\phi^{*}(t)$, where $\epsilon\ll 1$, the
leading-order expressions for $\cos\phi$ and $\sin\phi$ are found to be
$-\epsilon\phi^{*}$ and $1$, respectively. In order for the velocity terms in
the momentum equation to be non-trivial, they must scale as $\epsilon U$ and
$U$ for the parallel and perpendicular velocity, respectively, such that
$U_{p}=\epsilon UU_{p}^{*}$ and $U_{q}=UU_{q}^{*}$. By balancing the viscous
torque with the inertial torque, we get
$\Omega_{r}=\epsilon\Omega\Omega_{r}^{*}$ with
$\Omega_{r}^{*}=d\phi^{*}/dt^{*}$, resulting in the following system of
equations

$\displaystyle
St\left(\frac{dU_{p}^{*}}{dt^{*}}-\frac{d\phi^{*}}{dt^{*}}U_{q}^{*}\right)$
$\displaystyle=-\frac{4}{\pi}F_{p}^{*}(Re_{L}^{*},\chi)U_{p}^{*}-\phi^{*},$
(27) $\displaystyle
St\left(\frac{dU_{q}^{*}}{dt^{*}}+\epsilon^{2}\Omega_{r}^{*}U_{p}^{*}\right)$
$\displaystyle=-\frac{4}{\pi}F_{q}^{*}(Re_{L}^{*},\chi)U_{q}^{*}-1,$ (28)
$\displaystyle StJ_{q}^{*}\frac{d^{2}\phi^{*}}{dt^{*2}}$
$\displaystyle=-\frac{4}{\pi}\left(T_{i}^{*}(Re_{L}^{*},\chi,\theta)U_{p}^{*}U_{q}^{*}+T_{\Omega}^{*}(Re_{\Omega}^{*},\chi)\frac{d\phi^{*}}{dt^{*}}\right).$
(29)

Despite assuming a small angle around the equilibrium position, the resulting
system comprises three coupled nonlinear equations, rendering analytical
solutions for arbitrary Reynolds numbers impractical. However, various
simplifications can be made. First, the second term on the left-hand side of
equation 28 becomes negligible in the limit of $\epsilon\ll 1$. Second, using
the velocity estimate mentioned above, $Re_{L}^{*}\sim Re_{L}|U_{q}^{*}|$ at
leading order, resulting in equation 28 being decoupled from the other two
equations. Third, in the long-time limit $t^{*}\gg St$, the left-hand side of
equation 28 becomes negligible, leading to
$U_{q}^{*}\sim-\pi/4/F_{q}^{*}(Re_{L}^{*},\chi)$. Fourth, we consider the
small inertia limit $Ar\ll 1$. Under this limit, the system reads

$\displaystyle
St\left(\frac{dU_{p}^{*}}{dt^{*}}+\frac{\pi}{4}\frac{1}{F_{q}^{*}}\frac{d\phi^{*}}{dt^{*}}\right)$
$\displaystyle=-\frac{4}{\pi}F_{p}^{*}U_{p}^{*}-\phi^{*},$ (30) $\displaystyle
StJ_{q}^{*}\frac{d^{2}\phi^{*}}{dt^{*2}}$
$\displaystyle=\frac{T_{i}^{*}}{F_{q}^{*}}U_{p}^{*}-\frac{4}{\pi}T_{\Omega}^{*}\frac{d\phi^{*}}{dt^{*}}.$
(31)

The resulting set of equations is made of two coupled linear ordinary
differential equations (ODE) since in the limit $Ar\ll 1$, $T_{i}^{*}$,
$F_{q}^{*}$, $T_{\Omega}^{*}$ are independent of $Re_{L}^{*}$,
$Re_{\Omega}^{*}$ and $\theta$. By substituting the expression for $U_{p}$
from Equation 31 into Equation 30, the resulting equation becomes a third-
order linear ODE with constant coefficients. Although an analytical solution
may be obtained by summing three exponential functions, the coefficients
inside these functions are complex and obey a third-order polynomial. Instead,
we assume that the left-hand side of Equation 30 is negligible compared to the
right-hand side. This assumption is not satisfactory from a pure asymptotic
perspective but allows for clearer disentanglement of the physical mechanism
underlying the particle dynamic. In this limit,
$U_{p}^{*}\approx-\pi/4/F_{p}^{*}(Re_{L}^{,}\chi)\phi^{*}$ and

$StJ_{q}^{*}\frac{d^{2}\phi^{*}}{dt^{*2}}+\frac{4}{\pi}T_{\Omega}^{*}\frac{d\phi^{*}}{dt^{*}}+\frac{\pi}{4}\frac{T_{i}^{*}}{F_{q}^{*}F_{p}^{*}}\phi^{*}=0.$
(32)

The above equation is a classical damped harmonic oscillator. This equation
shares many similarities with the one used by Gustavsson  _et al._ [2, 3] who
studied the effect of particle inertia on the orientation of spheroids in
turbulent flows. Although the derivation presented in Gustavsson  _et al._ [3]
is different from the present one they nicely complement each other. Two
natural solutions exist depending on the sign of discriminant related to the
second order polynomial of this linear ODE. If the discriminant is positive
the solution is overdamped and leads to a decay of the motion without
oscillation while for a negative discriminant the solution is underdamped and
oscillatory motion is observed. Hence to obtain underdamped solution the
particle must obey the following criteria

$St>\frac{16}{\pi^{3}}\frac{T_{\Omega}^{*2}F_{q}^{*}F_{p}^{*}}{J_{q}^{*}T_{i}^{*}},$
(33)

which means that the particle inertia must be sufficiently large to lead to an
oscillating regime. Under the condition $Ar\ll 1$, one may expect underdamped
oscillation only in the case where the density ratio is very high.
Specifically, oscillatory motion can only be expected for high-density
particles, such as steel, settling in air.

Figure 4: Regime map showing the overdamped versus underdamped configuration
in the $(St,\chi)$ plane. Solid line : analytical criterion 33.

Figure 4 presents a map of the two regimes. The critical Stokes number
separating the two regimes is a decreasing function of $\chi$. This may be
proved using scaling arguments for $\chi\gg 1$. Indeed under this assumption
$F_{p}^{*}\sim F_{q}^{*}\sim 1/\ln(\chi)$, $T_{\Omega}^{*}\sim 1/\ln(\chi)$
and $T_{i}^{*}\sim 1/\ln^{2}(\chi)$. Hence for $\chi\gg 1$ the critical Stokes
number decreases as $1/\ln^{2}(\chi)$. Hence for a given Stokes number, the
underdamped regime will be promoted by using longer particle. If the criterion
33 is satisfied the solution of equation 32 is given by

$\phi^{*}(t)=Ae^{-\delta
t^{*}}\cos(\zeta+\omega_{0}t^{*})\quad\text{where}\quad\delta=\frac{2}{\pi}\frac{1}{St}\frac{T_{\Omega}^{*}}{J_{q}^{*}}\quad\text{and}\quad\omega_{0}=\left(\frac{\pi}{4}\frac{1}{St}\frac{T_{i}^{*}}{F_{q}^{*}F_{p}^{*}J_{q}^{*}}-\delta^{2}\right)^{1/2}.$
(34)

$A$ and $\zeta$ are constants given by the initial boundary conditions while
$\omega_{0}$ is the dimensionless natural frequency of the oscillator and
$\delta$ is the dimensionless damping coefficient. In the limit $St\gg 1$
which is relevant limit to expect underdamped oscillations one get

$\omega_{0}=\left(\frac{\pi}{4}\frac{1}{St}\frac{T_{i}^{*}}{F_{q}^{*}F_{p}^{*}J_{q}^{*}}\right)^{1/2}.$
(35)

As the value of the Stokes number increases, both $\omega_{0}$ and $\delta$
decrease, albeit with a less pronounced decrease observed for $\omega_{0}$.
Figure 5 displays the behavior of both quantities, which approach a near-
constant value for $\chi\gg 1$. Specifically, as $\chi$ increases, $\delta\sim
1/(St\ln\chi)$ while $\omega_{0}\sim(1/St)^{1/2}$. The results for $\chi$
larger than 30 were not plotted; however, it is expected that the behavior for
larger $\chi$ will remain consistent as all expressions for the loads converge
to the slender body theory results for $\chi\gg 1$.

$(a)$ $(b)$

Figure 5: (a) : Evolution of $\delta$ (equation 34) as function of $\chi$.
$(b):$ Evolution of $\omega_{0}$ (equation 34) as function of $\chi$.
$-$:$St=2000$, $--$:$St=1000$, $-\cdot$:$St=500$, $+$ : equation 35 for
$St=2000$.

One should recall the many assumptions required to derive the model. In
particular, we have assumed that the left-hand side (LHS) of Equation 30 is
negligible. Upon substitution of our estimate for $U_{p}$, the LHS of Equation
30 simplifies to $\pi/4Std\phi^{*}/dt(1/F_{q}^{*}-1/F_{p}^{*})$. For $\chi\leq
10$, it has been reported that $F_{q}^{*}\approx F_{p}^{*}$ [19], leading to a
near-cancellation of the LHS. However for longer fiber this assumption is
discutable. Additionally, the model assumes that both the history effect and
the coupling effect between translation and rotation can be neglected. The
former is questionable as the current model exhibits a constant period of
oscillation while $U_{p}$ is a decreasing function of time. Hence, in the
long-time limit, the unsteady term in the Navier-Stokes equation may be much
larger than the advective term. Therefore, the proposed model is not designed
to provide an exact value for the damping coefficient and natural frequency
for arbitrary dimensionless parameters especially larger fluid inertia, but
rather serves as a useful framework for qualitatively understanding the origin
of the damped oscillations observed in the subsequent section. Newsom and
Bruce [15] investigated the settling of elongated graphite particles in air
($\bar{\rho}\approx 1442,Ar\approx 0.054$, $\chi\geq 100$) and did not observe
oscillations around their stable orientation. However, in their study
$St\approx 4.2$, which is probably too low to expect such a phenomenon (see
Figure 4). The recent study of Bhowmick  _et al._ [29] (not published yet) who
investigated the settling of ellipsoidal particle with larger inertia in air
indicate that this solution is relevant.

##  V Comparison of the quasi-steady models with experiments and direct
numerical simulations

###  V.1 Experiments for $Re\ll 1$

In this section, we conduct a comparison between the experimental results
presented in Cabrera  _et al._ [20] and Roy  _et al._ [16] with those provided
by the models presented in the previous section. Despite having previously
compared our quasi-steady results to those of Cabrera  _et al._ [20] in Fintzi
_et al._ [19], we present the comparison again as we aim to further discuss
the effect of the particle inertia.

$\chi=8$ $\chi=8$

$(a)$ $(b)$
$\chi=16$ $\chi=16$

$(c)$ $(d)$

Figure 6: Dimensionless sedimenting velocities. (a), (c) : Dimensionless
sedimenting velocity as function of $\phi$. (b), (d) : Dimensionless angular
velocity divided by the square of the particle velocity as function of
$\theta$ (Figure 1). $\bullet$ : Cabrera  _et al._ [20] experiments for
$Ar\approx 0.147$, $\bar{\rho}\approx 12.2$. $-$ : prediction from the
unsteady equations 21 \- 23 for $Ar\approx 0.147$, $-\cdot-$ : prediction from
the unsteady equations 21 -23 for $Ar\approx 0.147$, $\bar{\rho}=1000$, $--$ :
prediction from the steady equations 24 \- 26 for $Ar\approx 0.147$,
$\bar{\rho}\approx 12.2$. In the unsteady computations, we have not displayed
the transient behaviour of the particles starting from rest. Top panel
corresponds to $\chi=8$ and bottom panel to $\chi=16$.

The sedimentation velocity and angular velocity measured in the study of
Cabrera  _et al._ [20] for particles with aspect ratios of $\chi=8$ and
$\chi=16$ are presented in Figures 6 (a) and (c). For $\chi=8$ the Reynolds
number varies in the range $0.026\leq Re\leq 0.035$ and is slightly higher for
$\chi=16$ ($0.032\leq Re\leq 0.047$). The model predictions are consistent
with the experimental data, with good agreement for the angular velocity for
$\chi=16$ but less so for $\chi=8$ (Figures 6 (b) and (d)). The experimental
data for $\chi=8$ exhibits considerable scatter, which the authors possibly
attribute to the presence of particle defects, particularly mass
inhomogeneities. The gravitational torque caused by mass inhomogeneities along
the body axis, $T_{g}$, scales as $T_{g}\sim mgL$ [16]. The ratio of inertial
torque to gravitational torque, scales as
$T_{i}/T_{g}\sim(\bar{\rho}-1)Ar\chi$. Therefore, for low Archimedes number
and moderate aspect ratio, as in the present case, mass inhomogeneities may
impact the angular velocity. Additionally, the results obtained from both sets
of equations, (21 \- 23) and (24 \- 26), are indistinguishable, thus particle
inertia can be safely neglected in the experiments of Cabrera  _et al._ [20].
To further examine the impact of particle inertia, we conducted computations
of the same physical configuration reported by Cabrera  _et al._ [20] (with
$Ar\approx 0.147$ and $\chi=8$ or $\chi=16$), but with a significantly larger
density ratio $\bar{\rho}=1000$, which is representative of particles settling
in air (represented by the dashed-dotted line in Figures 6). Our results
indicated that particle inertia has minimal effect on the outcome. This can be
explained by the fact that the maximum value of $Re_{\Omega}^{*}$ was found to
be approximately $\text{max}(Re_{\Omega}^{*})\approx 9\cdot 10^{-5}$ implying
that the inertial term in the momentum equations which scales as
$\bar{\rho}Re_{\Omega}$ remains much smaller than unity even with
$\bar{\rho}=1000$.

$(a)$ $(b)$ $(c)$

$(d)$ $(e)$ $(f)$

Figure 7: Dimensionless sedimenting velocities. $(a),$ $(d)$ : vertical
sedimenting velocity, $(b)$, $(e)$ : drift velocity and $(c)$, $(f)$ : angular
velocity. $\bullet$ : Roy  _et al._ [16] experiments for $Ar\approx 0.76$ and
$\bar{\rho}\approx 1.16$. $-$ : prediction from equations 21 \- 23, $--$ :
prediction from equations 24 \- 26, $-\cdot$ : prediction from the unsteady
equations 21 \- 23 for $Ar\approx 0.76$, $\bar{\rho}=1000$. Top panel
corresponds to $\chi\approx 20.5$ and bottom panel to $\chi\approx 101$.

The sedimenting velocities and angular velocities as a function of $\phi$,
from the experiments of Roy  _et al._ [16], are shown in Figures 7. Although
the angular velocity is not directly given in their work, it can be deduced
from their plot of the inertial torque, as they utilize a quasi-steady balance
for angular momentum. The results obtained in this study perfectly match the
experimental data. Furthermore, the difference between the sets of equations
(21 \- 23) and (24 \- 26) is once again negligible, indicating that particle
inertia is negligible in the experiments of Roy  _et al._ [16]. However, when
considering more inertial particles ($\bar{\rho}=1000$), deviations from the
experimental data are observed. Specifically, the dimensionless velocities
remain largely unchanged in amplitude, but a noticable shift to larger $\phi$
is observed. The system of ODEs represented by Equations 21 \- 23 can be
visualized as a phase portrait, as shown in Figures 7 (c) and (f). Notably,
the stable fixed point located at $\phi=\pi/2$ is identified as a focus or
spiral.

Figure 8: Evolution of $|\phi-\pi/2|$ as function of time with $Ar\approx
0.76$, $\bar{\rho}=1000$ and $\chi\approx 20.5$. $-$ : prediction from
equations 21 \- 23. $--$ : prediction from equation 34 with $Re\approx 0.156$
and $Re_{\Omega}\ll 1$.

The oscillatory behaviour of the angle as a function of time is evidenced in
figure 8 for $\chi\approx 20.5$. Exponential decay of $\phi$ over time is
observable for both the simplified model, expressed in Equation 34, and the
complete numerical solution. However, one may observe discrepancies between
the two models both for the decrease slope and oscillation period. This is not
surprising due to the many assumptions required to derive the model. The
analytical solution only offers a qualitative depiction of the underdamped
behaviour in this inertial configuration.

###  V.2 Direct numerical simulations for $Re\sim 1$

In the preceding subsection, we have demonstrated the unexpected agreement of
the quasi-steady equations with various literature findings for $Re\ll 1$.
However, no experiments of such nature have been conducted for higher Reynolds
numbers, and to the best of our knowledge no available simulations. Therefore,
we conducted direct numerical simulations to validate the applicability of the
model at higher Reynolds numbers ($Re\sim 1$).

####  V.2.1 Numerical methodology

Computations are carried out with the PeliGRIFF code Wachs  _et al._ [30].
This code was used in the past to investigate the settling of spherical and
angular particles [31, 32]. In brief, the code solves the three-dimensional
Navier-Stokes equations using a finite-volume discretization on a staggered
grid. The time-stepping strategy for the fluid phase is done thanks to a
second-order time-accurate Crank-Nicolson and Adams-Bashforth schemes. In
order to enforce the rigid-body motion inside the solid region a Lagrange
Multiplier/Fictitious Domain (DLM/FD) is used. We make use of a uniform
distribution of the lagrangian points along the surface of the cylindrical
body as detailed in [33]. More details on the numerical methods can be found
in Wachs  _et al._ [30].

$\mathbf{g}$$\mathcal{L}$$2\mathcal{L}$$\mathcal{L}$$x$$y$$z$ Figure 9: Scheme
of the computational domain (not to scale).

We explore the effect of $\chi$, $Ar$ and $\bar{\rho}$ within the range
$\chi=\\{2,4\\}$, $Ar=\\{24,96\\}$ and $\bar{\rho}=\\{1.5,10\\}$. The cylinder
is initiated with an angle $\phi=5^{\circ}$ and it is released at rest. In
addition, a case with a larger density ratio was computed, starting at
$\phi=60^{\circ}$, to investigate the damped oscillation regime. The
computational domain for this problem is depicted in Figure 9. Its
characteristic size length depends on the aspect ratio: for $\chi=2$,
$\mathcal{L}=12D$ while for $\chi=4$ $\mathcal{L}=15D$. Note that the value
prescribed for $\chi=2$ is larger than the one used in [32]. We have checked
that the present domain was sufficiently large by increasing its size of
$25\%$ in all directions finding less than $2\%$ error on $\Omega L/U$ for the
more challenging configuration ($Ar=24$, $\chi=4$). Obviously, larger domains
might be used to avoid this small effect of the boundaries at the expense of a
substantially larger numerical cost. However, we have to recall that the
present paper is not aimed to provide domain size perfectly independent
simulations of the problem but rather to provide bound to the quasi-steady
assumptions. To this aim, we compare the numerical results to the model of
section IV which heavily relies on fit and also contains some error with
respect to the simulation they are extracted from. Regarding the boundary
conditions, the domain is biperiodic in the lateral directions while zero
velocity and outflow conditions are imposed on the upstream and downstream
boundaries, respectively. A uniform cell distribution is imposed in a
rectangular region depicted in Figure 9 as an orange volume. This region is
located $2D$ below the middle of the domain and extends up to $6D$ in the
direction of gravity. In this flow region, $25$ cells are distributed per body
diameter which is sufficient to accurately compute the dynamic for the range
of Reynolds number investigated [33]. The number of cells per mesh is 26
millions for $\chi=2$ and 46 millions for $\chi=4$. The time step is imposed
such that the CFL always fall below $0.25$. To prevent the cylinder to exit
the numerical domain the computational domain is moved in the vertical
direction so as to keep the particle at least at a distance of $10D$ of the
upstream boundary. More details on the domain translation technique can be
found in Rahmani and Wachs [31], Seyed-Ahmadi and Wachs [32] The magnitude of
this domain translation is chosen to be the grid size ($D/25$). The
simulations are run up until the cylindrical particle reaches its equilibrium
position.

####  V.2.2 Numerical results

To investigate the impact of moderate inertia on the rod motion we first
consider the situation $Ar=24$. In this scenario, the Reynolds number will
fluctuate between $1.7$ and $1.9$ when $\chi$ equals $2$, and between $2$ and
$2.7$ when $\chi$ is $4$.

$(a)$ $(b)$ $(c)$

$(d)$ $(e)$ $(f)$

Figure 10: Dimensionless sedimenting velocities for $Ar=24$. $(a),$ $(d)$ :
vertical sedimenting velocity, $(b)$, $(e)$ : drift velocity and $(c)$, $(f)$
: angular velocity. $-$ : direct numerical simulation results with
$\bar{\rho}=1.5$, $-\cdot-$: direct numerical simulation results with
$\bar{\rho}=10$, $--$ : prediction from equations 21 \- 23 with
$\bar{\rho}=1.5$, $\cdot\cdot$ : prediction from equations 21 \- 23 with
$\bar{\rho}=10$. Top panel corresponds to $\chi=2$ and bottom panel to
$\chi=4$.

For the shortest cylinder with aspect ratio $\chi=2$, there is a noticeable
agreement between the results obtained from equations 21 \- 23 (including
particle inertia) and those from direct numerical simulations, as illustrated
in Figure 10 (a) - (c). Varying the density ratio from 1.5 to 10 has a
negligible effect on the sedimentation velocities in this particular case. In
Figure 10 (d) - (f), we present the results for $\chi=4$. Although there is a
qualitative agreement between the model and the numerical results, the model
overestimates the drift velocity and the angular velocity. Moreover, for
$\chi=4$, the effect of particle inertia on the sedimentation velocities is
more pronounced. In particular, this case exhibits underdamped motion. To
investigate this oscillating behavior further, we have carried out one
simulation of the settling of a rod with the same dimensionless parameters
($Ar=24$, $\chi=4$), but with a much larger density ratio ($\bar{\rho}=50$).
In order to maintain reasonable computational time, the simulation is
initiated with $\phi=60^{\circ}$.

$(a)$ $(b)$

Figure 11: Comparison between the simulation results and the quasi-steady
model with $Ar=24$, $\chi=4$ and $\bar{\rho}=50$. (a) : Phase space diagram of
the oscillator. (b) : Time evolution of the inclination angle. $-$ : direct
numerical simulation results, $--$ : prediction from equations 21 \- 23.

There is a good agreement between the simulation and the model even if the
model overpredicts the angular velocity (Figure 11). We recall that we neglect
the history loads in equations 21 \- 23 which may affect the initial transient
since the rod starts from rest. One may also observe a good agreement between
the oscillating period, the decrease in amplitude given by the numerical
results and the model made of equations 21 \- 23 (Figure 11 (b)). This is
somehow surprising since there is a priori no reason for neglecting the
history loads in the underdamped regime. They appear to have a negligible
effect on this regime.

$(a)$ $(b)$ $(c)$

$(d)$ $(e)$ $(f)$

Figure 12: Dimensionless sedimenting velocities for $Ar=96$. $(a),$ $(d)$ :
vertical sedimenting velocity, $(b)$, $(e)$ : drift velocity and $(c)$, $(f)$
: angular velocity. $-$ : direct numerical simulation results with
$\bar{\rho}=1.5$, $-\cdot-$: direct numerical simulation results with
$\bar{\rho}=10$, $--$ : prediction from equations 21 -23 with
$\bar{\rho}=1.5$, $\cdot\cdot$ : prediction from equations 21, 22 and 23 with
$\bar{\rho}=10$. Top panel corresponds to $\chi=2$ and bottom panel to
$\chi=4$.

We now consider much larger inertia effects ($Ar=96$). Within this regime, the
Reynolds number is found to vary between $5.3$ and $6.2$ for $\chi=2$ and
between $5.8$ and $8.6$ for $\chi=4$. Surprisingly, the quasi-steady model
exhibits good agreement with numerical results for $\chi=2$ (as depicted in
Figure 12 (a)-(c)), despite the significant role played by inertia effects.
However, for $\chi=4$, the model overestimates both the drift velocity and
angular velocity (as shown in Figure 12 (d)-(f)). We observe that all
configurations result in underdamped regimes. This regime is much more
pronounced for $\chi=4$ and $\bar{\rho}=10$ than for $\chi=2$ and
$\bar{\rho}=1.5$ for which the amplitude of oscillation is small and nearly
indistinguishable. Finally, we compare the numerical and model-based results
(derived from Equations 21 \- 23) for $Ar=96$, $\chi=4$, and $\bar{\rho}=10$,
by examining the oscillation period and the decay of amplitude. The results
are displayed on Figure 13. The model slightly underpredicts the amplitude
decay as well as the oscillation period.

Figure 13: Time evolution of the inclination angle ($Ar=96$,
$\chi=4,\bar{\rho}=10$). $-$ : direct numerical simulation results, $--$ :
prediction from equations 21 \- 23.

##  VI Discussion and conclusion

Based on the results presented above, two main conclusions can be drawn.
Firstly, it can be inferred that the quasi-steady assumption holds true across
a broad range of dimensionless parameters, as it is supported by both the
experiments conducted by Cabrera  _et al._ [20] and Roy  _et al._ [16] and by
the direct numerical simulations. Although these experiments and simulations
do not strictly adhere to the condition $Ar\ll 1/\chi$, the quasi-steady
assumption is still valid for a significantly larger range of values than
initially predicted. This can be attributed to the magnitude of the angular
velocity $\Omega$, which has the correct scaling but is at least three order
of magnitude smaller than the anticipated value. As a result, $\Omega L/U\ll
1$ for all configurations except the highest inertial simulations (see Table
1). Aditionnally, these findings support the notion that the quasi-steady
assumption remains valid as long as $\Omega L/U\ll 1$, thereby strengthening
our analysis.

| Cabrera  _et al._ [20] | Roy  _et al._ [16] | Direct numerical simulations
---|---|---|---
| $Ar\approx 0.147$ | $Ar\approx 0.76$ | $Ar=24$ | $Ar=96$
| $\chi=8$ | $\chi=16$ | $\chi=20$ | $\chi=100$ | $\chi=2$ | $\chi=4$ | $\chi=2$ | $\chi=4$
max$(\Omega_{r}L/|\mathbf{U}|)$ | 0.024 | 0.035 | 0.12 | 0.11 | 0.13 | 0.25 | 0.2 | 0.33
Table 1: Values of max$(\Omega_{r}L/|\mathbf{U}|)$ in the experiments of [20],
[16] and in the simulation results.

Secondly, it can be concluded based on the results of the previous section
that particle inertia plays no significant role on the magnitude of the
sedimenting velocities in the system under investigation. Hence, a correct
estimate of $\Omega L/U$ can be obtained by disregarging the particle inertia
in the equations of motion and considering the system made of equations 24 \-
26. We have solved this system for $\chi\in[2;500]$ and $Ar\in[0.001;150]$ by
using the semi-empirical expression for the loads proposed by Fintzi  _et al._
[19] for $\chi\leq 30$ and Khayat and Cox [12] theory for $\chi>30$. In each
of the run we have computed the maximum value of $\Omega_{r}L/|\mathbf{U}|$
(Figure 14).

$(a)$ $(b)$

Figure 14: Evolution of max($\Omega_{r}L/|\mathbf{U}|)$ with $\chi Ar$.
Equations 24 \- 26 are solved using a root-finding algorithm for 20 values of
$\phi$ ranging between $0$ and $\pi/2$.

For $Ar\chi\leq 10$ except for the smallest aspect ratio $\Omega L/U$ scales
as $Ar\chi$ in agreement with our scaling analysis. For larger $\chi Ar$ the
rate of increase of $\Omega L/U$ decreases strongly, especially for $\chi\geq
12$. Indeeed it can be noted that qualitatively different behavior may be
observed for moderate inertial effects ($Re_{L}\sim 1$) and very large aspect
ratios ($\chi\gg 1$). Fintzi  _et al._ [19] observed that for $\chi=30$, the
inertial torque scales as $T_{i}/(\mu UL^{2})\sim Re_{L}^{1/3}$, yielding
$T_{i}\sim\rho^{1/3}\mu^{2/3}U^{4/3}L^{7/3}$. By balancing this inertial
torque with the resistive torque, we obtain
$\Omega\sim\rho^{1/3}\mu^{-1/3}U^{4/3}L^{-2/3}$, and $\Omega L/U\sim
Re_{L}^{1/3}$. Using the same scaling as in section IV for the velocity, we
obtain $\Omega L/U\sim Ar^{1/3}\chi^{1/3}$. Therefore, the rate of increase of
$\Omega L/U$ is slower as the Reynolds number increases for elongated
particle. This trend is even more pronounced for very long fibers (Figure 14
(b)). In this regime, Khayat and Cox [12] demonstrated that the inertial
torque decreases with the Reynolds number for $Re_{L}\geq 4$. However, the
validity of their theory for such high Reynolds numbers may be questionable,
as their asymptotic solution requires $Re_{L}\ll\ln(\chi)$. Nonetheless,
simulations performed by Khair and Chisholm [34] and Shin  _et al._ [17]
indicate that the solutions proposed by Khayat and Cox [12] remain valid for
$Re_{L}\approx 10$ and $Re_{L}\approx 5$, respectively, for the longitudinal
force on a long spheroid aligned with the flow direction and the torque on a
long fiber, with $\chi=100$. Moreover, Khayat and Cox [12] models are in very
good agreement with Roy  _et al._ [16] experimental results for $\chi=100$ and
$Re_{L}\approx 7.6$. Consequently, the solution provided by Khayat and Cox
[12] is considered to provide quantitatively accurate results up to
$Re_{L}\approx 10$ as long as the fiber is adequately elongated ($\chi\geq
100$). If we consider that the quasi-steady assumption fails for values of
$\Omega L/U$ larger than $0.2$, then we can expect the quasi-steady models to
be accurate for $\chi Ar\approx 200$ if $\chi=2$, and for $\chi Ar\approx 40$
if $2<\chi\leq 30$. However, for more elongated fibers, the quasi-steady
assumption should remain valid as long as the underlying assumptions made in
the derivation of Khayat and Cox [12] models are fullfilled, particularly if
$Re\ll 1$.

This paper also focuses on the underdamped regime and its characteristics,
such as oscillation period and amplitude decrease rate, in comparison to
simplified models. It was found that the most simple model i. e. the
underdamped oscillator qualitatively reproduces numerical results but not
quantitatively. This may have significant implications for atmospheric flows
where particle orientation is driven by this solution[2, 3]. To improve the
model, particle inertia could be included in the parallel velocity equation,
resulting in a third-order linear ordinary differential equation that can be
easily solved. A natural perspective to this model might be to study rods wake
instability, like fluttering motion [25]. However, in such high inertial flow,
there is no reason for the quasi-steady assumption to remain valid Fabre  _et
al._ [35]. Moreover, even under the quasi-steady limit, the loads on such
bodies at high Reynolds numbers are only known for very few configurations
[33, 27].

##  VII Acknowledgements

ANR MUSCATs financial support is greatly appreciated. We thank Bernhard Mehlig
for stimulating discussions and for pointing out the former studies on the
underdamped oscillator. The author is indebted to Greg Voth for providing the
experimental data of Roy  _et al._ [16]. We also thank Jacques Magnaudet for
fruitful discussions on the added mass torque.

##  Appendix A Added mass coefficients

Figure 15: Dimensionless added mass coefficients as a function of $\chi$.
($\blacksquare$) : Loewenberg [9] results, $\bullet$ : present results based
on JADIM code [22].

In this appendix, we compute the coefficient $A_{p}$, $A_{q}$ and $D_{q}$ by
using the JADIM code for $1\leq\chi\leq 15$. The numerical details concerning
the code as well as the mesh properties can be found elsewhere [22, 27, 18].
At the time $t=0$ we impose a constant linear or angular acceleration to a
cylinder initially at rest. Since in the very short time limit viscous and
rotational contributions are negligible in comparison to potential flow
contribution [21] one can easily recover the added mass coefficients by
computing the loads on the body. Figures 15 display the dimensionless added
mass coefficient $A_{p}^{*}=A_{p}/(\rho V),A_{q}^{*}=A_{q}/(\rho V)$ and
$D_{q}^{*}=D_{q}/(\rho VL^{2})$ as a function of the aspect ratio. A good
agreement is observed between the present results and Loewenberg [9] results
obtained using potential flow calculations. In the limit $\chi\gg 1$, $\chi
A_{p}^{*}$ and $A_{q}^{*}$ tend toward a constant value. For the particular
case $A_{q}^{*}$ this constant is simply $1$, i.e. the added mass coefficient
on an infinitely long cylinder perpendicular to the flow direction. As a
result, we propose the following correlation

$A_{p}^{*}=\frac{1}{\chi}\left(0.655-\frac{0.141}{1+\chi^{1.17}}\right),$ (36)
$A_{q}^{*}=1-\frac{0.828}{1+\chi^{1.12}}.$ (37)

To the best of our knowledge, the functional dependency of $D_{q}^{*}$ with
respect to $\chi$ has not been published yet in the literature at least for
$\chi\geq 1$. Figure 15 shows that in the limit $\chi\gg 1$ $D_{q}^{*}$ tends
toward a constant value which is a priori unknown. However, in the case of a
Rankine ovoid, for $\chi\gg 1$, $D_{q}^{*}\approx 1/12$ Howe [23]. Since in
the limit of large aspect ratio one may assume that the rounded ends of the
Rankine ovoid have little effect on the added mass coefficient, we propose the
following empirical correlation

$D_{q}^{*}=\frac{1}{12}-\frac{0.11}{(1+\chi^{0.8})},$ (38)

which matches closely the numerical results.

##  Appendix B Short-time asymptotic expansion

We calculate the leading-order terms of the short-time asymptotic expansion.
Since
$\cos(\phi^{(0)}+\epsilon\phi^{(1)})\sim\cos\phi^{(0)}-\epsilon\phi^{(1)}\sin\phi^{(0)}$
and
$\sin(\phi^{(0)}+\epsilon\phi^{(1)})\sim\sin\phi^{(0)}+\epsilon\phi^{(1)}\cos\phi^{(0)}$
at zero-th order equations 10, 11, 12 and 13 simplifies to

$\displaystyle\frac{dU_{p}^{*(0)}}{dt^{*}}$
$\displaystyle=\mathcal{A}\cos\phi^{(0)},$ (39)
$\displaystyle\frac{dU_{q}^{*(0)}}{dt^{*}}$
$\displaystyle=-\mathcal{B}\sin\phi^{(0)},$ (40)
$\displaystyle\frac{d\Omega_{r}^{*(0)}}{dt^{*}}$
$\displaystyle=-\mathcal{C}U_{p}^{*(0)}U_{q}^{*(0)},$ (41)
$\displaystyle\frac{d\phi^{(0)}}{dt^{*}}$ $\displaystyle=0.$ (42)

The zero-th order solution is easily obtained and reads

$\displaystyle U_{p}^{*(0)}$ $\displaystyle=\mathcal{A}t^{*}\cos\phi^{(0)},$
(43) $\displaystyle U_{q}^{*(0)}$
$\displaystyle=-\mathcal{B}t^{*}\sin\phi^{(0)},$ (44)
$\displaystyle\Omega_{r}^{*(0)}$
$\displaystyle=\frac{\mathcal{A}\mathcal{B}\mathcal{C}}{3}t^{*3}\cos\phi^{(0)}\sin\phi^{(0)},$
(45) $\displaystyle\phi^{(0)}$ $\displaystyle=\phi(t^{*}=0),$ (46)

where we have assumed that the cylinder starts from rest. At first order
equations 10, 11, 12 and 13 give

$\displaystyle\frac{dU_{p}^{*(1)}}{dt^{*}}$
$\displaystyle=\frac{\mathcal{A}}{\mathcal{B}}\Omega_{r}^{*(0)}U_{q}^{*(0)}-\mathcal{A}\phi^{(1)}\sin\phi^{(0)},$
(47) $\displaystyle\frac{dU_{q}^{*(1)}}{dt^{*}}$
$\displaystyle=-\frac{\mathcal{B}}{\mathcal{A}}\Omega_{r}^{*(0)}U_{p}^{*(0)}-\mathcal{B}\phi^{(1)}\cos\phi^{(0)},$
(48) $\displaystyle\frac{d\Omega_{r}^{*(1)}}{dt^{*}}$ $\displaystyle=0,$ (49)
$\displaystyle\frac{d\phi^{(1)}}{dt^{*}}$ $\displaystyle=\Omega_{r}^{*(0)},$
(50)

which leads to

$\displaystyle U_{p}^{*(1)}$
$\displaystyle=-\frac{1}{12}\mathcal{A}^{2}\mathcal{B}\mathcal{C}t^{*5}\cos\phi^{(0)}\sin^{2}\phi^{(0)},$
(51) $\displaystyle U_{q}^{*(1)}$
$\displaystyle=-\frac{1}{12}\mathcal{A}\mathcal{B}^{2}\mathcal{C}t^{*5}\cos^{2}\phi^{(0)}\sin\phi^{(0)},$
(52) $\displaystyle\Omega_{r}^{*(1)}$ $\displaystyle=0,$ (53)
$\displaystyle\phi^{(1)}$
$\displaystyle=\frac{\mathcal{A}\mathcal{B}\mathcal{C}}{12}t^{*4}\cos\phi^{(0)}\sin\phi^{(0)}.$
(54)

##  Appendix C Quasi-steady loads

The expression of the quasi-steady loads used in this paper can be found
below.

###  C.1 Moderately long rods : $\chi\leq 30$

##### Expression of $F_{p}$ :

The expression of $F_{p}$ reads [19]

$\displaystyle F_{p}(Re_{L}^{*},\chi,\theta)=-2\pi\mu|\mathbf{U}|L\cos\theta$
$\displaystyle\left(\frac{A_{Re=0}^{(1)}+A^{(1)}(Re_{L}^{*})}{\ln(2\chi)}+\frac{A_{Re=0}^{(2)}+A^{(2)}(Re_{L}^{*})}{\ln^{2}(2\chi)}+\frac{A_{Re=0}^{(3)}+A^{(3)}(Re_{L}^{*})}{\ln^{3}(2\chi)}\right.$
$\displaystyle\left.+\frac{A_{Re=0}^{(4)}+A^{(4)}(Re_{L}^{*})}{\ln^{4}(2\chi)}+\frac{2.34}{\chi^{2/3}(\chi-\frac{1}{2})^{1.75}}\right),$
(55)

where $A_{Re=0}^{(1)}=1$, $A_{Re=0}^{(2)}\approx 0.807$,
$A_{Re=0}^{(3)}\approx 0.829$, $A_{Re=0}^{(4)}\approx 1.45$ [27]. The first
order inertial correction is null, $A^{(1)}(Re_{L}^{*})=0$ , while the second,
third and fourth order inertial functions read

$\displaystyle A^{(2)}(Re_{L}^{*})$
$\displaystyle=\frac{1}{2}\left(\frac{E_{1}(2Re_{L}^{*})+\ln(2Re_{L}^{*})-e^{-2Re_{L}^{*}}+\gamma+1}{2Re_{L}^{*}}+E_{1}(2Re_{L}^{*})+\ln(2Re_{L}^{*})+\gamma-2\right),$
(56) $\displaystyle A^{(3)}(Re_{L}^{*})$
$\displaystyle=A_{A}^{(3)}(Re_{L}^{*})+A_{B}^{(3)}(Re_{L}^{*})+2A^{(2)}(Re_{L}^{*})\ln{(2)},$
(57) $\displaystyle A^{(4)}(Re_{L}^{*})$
$\displaystyle=3\ln{(2)}\left(A_{A}^{(3)}(Re_{L}^{*})+A_{B}^{(3)}(Re_{L}^{*})\right)+3A^{(2)}(Re_{L}^{*})\ln{(2)}^{2}-0.636Re_{L}^{*0.762},$
(58)

where $\gamma$ is the Euler constant,
$E_{1}(x)=\int_{x}^{\infty}\frac{e^{-t}}{t}dt$ the exponential integral
function and here $Re_{L}^{*}=\rho|\mathbf{U}|L/(2\mu)$. One may observe that
$|\mathbf{U}|\cos\theta=U_{p}$.

##### Expression of $F_{q}$ :

The expression of $F_{q}$ can be found in [19] and reads

$\displaystyle F_{q}(Re_{L}^{*},\chi,\theta)=4\pi\mu|\mathbf{U}|L\sin\theta$
$\displaystyle\left(\frac{B_{Re=0}^{(1)}+B^{(1)}(Re_{L}^{*})}{\ln(2\chi)}+\frac{B_{Re=0}^{(2)}+B^{(2)}(Re_{L}^{*})}{\ln^{2}(2\chi)}+\frac{B_{Re=0}^{(3)}+B^{(3)}(Re_{L}^{*})}{\ln^{3}(2\chi)}\right.$
(59)
$\displaystyle+\left.\frac{B_{Re=0}^{(4)}+B^{(4)}(Re_{L}^{*})}{\ln^{4}(2\chi)}-\frac{0.568}{\chi^{2/3}(\chi-\frac{1}{2})^{1.75}}\right).$
(60)

where $B_{Re=0}^{(1)}=1$, $B_{Re=0}^{(2)}\approx-0.193$,
$B_{Re=0}^{(3)}\approx 0.214$, $B_{Re=0}^{(4)}\approx 0.387$ Kharrouba  _et
al._ [27]. $B^{(1)}(Re_{L}^{*})=0$ and

$\displaystyle B^{(2)}(Re_{L}^{*})$
$\displaystyle=E_{1}\left(Re_{L}^{*}\right)+\ln{\left(Re_{L}^{*}\right)}-\frac{e^{-Re_{L}^{*}}-1}{Re_{L}^{*}}+\gamma-1$
(61) $\displaystyle,B^{(3)}(Re_{L}^{*})$
$\displaystyle=2\ln{(2)}B^{(2)}(Re_{L}^{*})+B_{e}^{(3)}(Re_{L}^{*}),$ (62)
$\displaystyle B^{(4)}(Re_{L}^{*})$
$\displaystyle=3\ln{(2)}^{2}B^{(2)}(Re_{L}^{*})+3\ln{(2)}B_{e}^{(3)}(Re_{L}^{*})+B_{e}^{(4)}(Re_{L}^{*}).$
(63)

We have $|\mathbf{U}|\sin\theta=-U_{q}$.

##### Expression of $T_{r}^{i}$ :

The expression of $T_{r}^{i}$ can be found in [19]

$\displaystyle
T_{r}^{i}(Re_{L}^{*},\chi,\theta)=\rho|\mathbf{U}|^{2}L^{3}\sin{(2\theta)}\frac{5\pi}{48(1+Re_{L}^{*1.991})^{0.331}}$
$\displaystyle\left(\frac{1}{\ln^{2}(3\chi)}+\frac{2.244-1.813Re_{L}^{*0.543}}{\ln^{3}(3\chi)}\right.$
(64)
$\displaystyle\left.-\frac{3.603+8.854Re_{L}^{*0.538}}{\ln^{4}(3\chi)}-\frac{14.301(Re_{L}^{*}/\chi)^{0.448}}{\ln^{5}(3\chi)}\right).$
(65)

One may note that $|\mathbf{U}|^{2}\sin{(2\theta)}=-2U_{p}U_{q}$.

##### Expression of $T_{r}^{\Omega}$ :

An expression for the resisting torque due to the particle rotation can be
found in Pierson  _et al._ [18] and may be expressed as

$\displaystyle T_{r}^{\Omega}=$
$\displaystyle-\frac{-\pi\mu\Omega_{r}L^{3}}{3}\left[\frac{1}{\ln(2\chi)}+\frac{1}{\ln^{2}(2\chi)}\left(\frac{11}{6}-\ln
2+f(\chi,,Re_{\Omega}^{*})\right)+\frac{1}{\ln^{3}(2\chi)}\left(\frac{161}{36}-\frac{\pi^{2}}{12}-\frac{11}{3}\ln
2+(\ln 2)^{2}\right)\right.$
$\displaystyle\left.+\frac{1}{\ln^{4}(2\chi)}\left(1-\frac{1}{(2\chi)^{1.2}}\right)^{5}\left(-\frac{5}{4}\zeta(3)+\frac{1033}{72}-\ln^{3}(2)+\frac{11}{2}\ln^{2}(2)-\frac{161}{12}\ln
2-\pi^{2}\left(\frac{11}{24}-\frac{1}{4}\ln 2\right)\right)\right],$ (66)

with $f(\chi,Re_{\Omega}^{*})=0.018\chi^{2.3}Re_{\Omega}^{*0.9}$ and
$Re_{\Omega}^{*}=\rho|\boldsymbol{\Omega}|D^{2}/\mu$.

###  C.2 Long rods : $\chi>30$

For sufficiently long rods Khayat and Cox [12] expressions for the loads are
accurate. In the following we present the Khayat and Cox [12] expression still
making use of the linearization proposed by [28] for the forces. Also we make
use of the $1/\ln\chi$ original expansion proposed by Khayat and Cox [12]
rather than the $1/\ln(2\chi)$ expansion.

##### Expression of $F_{p}$ :

$F_{p}(Re_{L}^{*},\chi,\theta)=\frac{-2\pi\mu|\mathbf{U}|L\cos\theta}{\ln\chi}\left(1-\frac{A^{(2)}(Re_{L}^{*})-4\ln
2+3}{\ln\chi}\right)^{-1}.$ (67)

##### Expression of $F_{q}$ :

$F_{q}(Re_{L}^{*},\chi,\theta)=\frac{4\pi\mu|\mathbf{U}|L\sin\theta}{\ln\chi}\left(1-\frac{B^{(2)}(Re_{L}^{*})+1/2-\ln
4}{\ln\chi}\right)^{-1}.$ (68)

##### Expression of $T_{r}^{i}$ :

$\displaystyle T_{r}^{i}(Re_{L}^{*},\chi,\theta)=-\mu
UL^{2}\frac{\pi}{2}\left(\frac{1}{\ln\chi}\right)^{2}\left[\cos\theta\left(P(X)-Q(X)+P(Y)-Q(Y)\right)\right.\left.+P(Y)-P(X)\right]\sin\theta,$
(69)

with $Q(x)=\frac{E_{1}(x)+\ln(x)+\gamma}{x}$,
$P(x)=\frac{2}{x}\left(1+\frac{e^{-x}-1}{x}\right)$,
$X=Re_{L}^{*}\left(1-\cos\theta\right)$ and
$Y=Re_{L}^{*}\left(1+\cos\theta\right)$.

##### Expression of $T_{r}^{\Omega}$ :

Since to the best of our knowledge there is no expression for this torque for
finite inertia effect and very long fibers we make use of expression 66
neglecting the inertial correction $f$.

## References

  * Poulain  _et al._ [2018] M. Poulain, M. J. Mercier, L. Brach, M. Martignac, C. Routaboul, E. Perez, M. C. Desjean, and A. Ter Halle, Small microplastics as a main contributor to plastic mass balance in the north atlantic subtropical gyre, Environmental science & technology 53, 1157 (2018).
  * Gustavsson  _et al._ [2019] K. Gustavsson, M. Sheikh, D. Lopez, A. Naso, A. Pumir, and B. Mehlig, Effect of fluid inertia on the orientation of a small prolate spheroid settling in turbulence, New Journal of Physics 21, 083008 (2019).
  * Gustavsson  _et al._ [2021] K. Gustavsson, M. Sheikh, A. Naso, A. Pumir, and B. Mehlig, Effect of particle inertia on the alignment of small ice crystals in turbulent clouds, Journal of the Atmospheric Sciences 78, 2573 (2021).
  * Batchelor [1970] G. K. Batchelor, Slender-body theory for particles of arbitrary cross-section in stokes flow, Journal of Fluid Mechanics 44, 419–440 (1970).
  * Cox [1965] R. Cox, The steady motion of a particle of arbitrary shape at small reynolds numbers, Journal of Fluid Mechanics 23, 625 (1965).
  * Howe [1995] M. Howe, On the force and moment on a body in an incompressible fluid, with application to rigid bodies and bubbles at high and low reynolds numbers, The Quarterly Journal of Mechanics and Applied Mathematics 48, 401 (1995).
  * Ern  _et al._ [2012] P. Ern, F. Risso, D. Fabre, and J. Magnaudet, Wake-induced oscillatory paths of bodies freely rising or falling in fluids, Annual Review of Fluid Mechanics 44, 97 (2012).
  * Kim and Karrila [2013] S. Kim and S. J. Karrila,  _Microhydrodynamics: principles and selected applications_ (Courier Corporation, 2013).
  * Loewenberg [1993] M. Loewenberg, Stokes resistance, added mass, and basset force for arbitrarily oriented, finite-length cylinders, Physics of Fluids A: Fluid Dynamics 5, 765 (1993).
  * Kabarowski and Khair [2020] J. K. Kabarowski and A. S. Khair, The force on a slender particle under oscillatory translational motion in unsteady stokes flow, Journal of Fluid Mechanics 884 (2020).
  * Lawrence and Weinbaum [1988] C. t. Lawrence and S. Weinbaum, The unsteady force on a body at low reynolds number; the axisymmetric motion of a spheroid, Journal of Fluid Mechanics 189, 463 (1988).
  * Khayat and Cox [1989] R. Khayat and R. Cox, Inertia effects on the motion of long slender bodies, Journal of Fluid Mechanics 209, 435 (1989).
  * Rubinow and Keller [1961] S. I. Rubinow and J. B. Keller, The transverse force on a spinning sphere moving in a viscous fluid, Journal of Fluid Mechanics 11, 447 (1961).
  * Lovalenti and Brady [1993] P. M. Lovalenti and J. F. Brady, The hydrodynamic force on a rigid particle undergoing arbitrary time-dependent motion at small reynolds number, Journal of Fluid Mechanics 256, 561 (1993).
  * Newsom and Bruce [1994] R. Newsom and C. Bruce, The dynamics of fibrous aerosols in a quiescent atmosphere, Physics of Fluids 6, 521 (1994).
  * Roy  _et al._ [2019] A. Roy, R. J. Hamati, L. Tierney, D. L. Koch, and G. A. Voth, Inertial torques and a symmetry breaking orientational transition in the sedimentation of slender fibres,  875, 576 (2019).
  * Shin  _et al._ [2006] M. Shin, D. L. Koch, and G. Subramanian, A pseudospectral method to evaluate the fluid velocity produced by an array of translating slender fibers, Physics of Fluids 18, 063301 (2006).
  * Pierson  _et al._ [2021] J.-L. Pierson, M. Kharrouba, and J. Magnaudet, Hydrodynamic torque on a slender cylinder rotating perpendicularly to its symmetry axis, Physical Review Fluids 6, 094303 (2021).
  * Fintzi  _et al._ [2023] N. Fintzi, L. Gamet, and J.-L. Pierson, Inertial loads on a finite-length cylinder embedded in a steady uniform flow, Physical Review Fluids (submitted) (2023).
  * Cabrera  _et al._ [2022] F. Cabrera, M. Z. Sheikh, B. Mehlig, N. Plihon, M. Bourgoin, A. Pumir, and A. Naso, Experimental validation of fluid inertia models for a cylinder settling in a quiescent flow, Phys. Rev. Fluids 7, 024301 (2022).
  * Mougin and Magnaudet [2002] G. Mougin and J. Magnaudet, The generalized kirchhoff equations and their application to the interaction between a rigid body and an arbitrary time-dependent viscous flow, International journal of multiphase flow 28, 1837 (2002).
  * Kharrouba [2020] M. Kharrouba,  _Etude analytique et numérique des efforts s’éxerçant sur des particules cylindriques allongées._ , Ph.D. thesis (2020).
  * Howe [2006] M. S. Howe,  _Hydrodynamics and sound_ (Cambridge University Press, 2006).
  * Lamb [1953] H. Lamb,  _Hydrodynamics_ (At The University Press, 1953).
  * Toupoint  _et al._ [2019] C. Toupoint, P. Ern, and V. Roig, Kinematics and wake of freely falling cylinders at moderate reynolds numbers, Journal of Fluid Mechanics 866, 82 (2019).
  * Dabade  _et al._ [2015] V. Dabade, N. K. Marath, and G. Subramanian, Effects of inertia and viscoelasticity on sedimenting anisotropic particles, Journal of Fluid Mechanics 778, 133 (2015).
  * Kharrouba  _et al._ [2021] M. Kharrouba, J.-L. Pierson, and J. Magnaudet, Flow structure and loads over inclined cylindrical rodlike particles and fibers, Physical Review Fluids 6, 044308 (2021).
  * Lopez and Guazzelli [2017] D. Lopez and E. Guazzelli, Inertial effects on fibers settling in a vortical flow, Physical Review Fluids 2, 024306 (2017).
  * Bhowmick  _et al._ [2023] T. Bhowmick, J. Seesing, K. Gustavsson, J. Guettler, A. Pumir, B. Mehlig, Y. Wang, and G. Bagheri, Inertial angular dynamics of non-spherical atmospheric particles, arXiv preprint arXiv:2303.04299 (2023).
  * Wachs  _et al._ [2015] A. Wachs, A. Hammouti, G. Vinay, and M. Rahmani, Accuracy of finite volume/staggered grid distributed lagrange multiplier/fictitious domain simulations of particulate flows, Computers & Fluids 115, 154 (2015).
  * Rahmani and Wachs [2014] M. Rahmani and A. Wachs, Free falling and rising of spherical and angular particles, Physics of Fluids 26, 083301 (2014).
  * Seyed-Ahmadi and Wachs [2019] A. Seyed-Ahmadi and A. Wachs, Dynamics and wakes of freely settling and rising cubes, Physical Review Fluids 4, 074304 (2019).
  * Pierson  _et al._ [2019] J.-L. Pierson, F. Auguste, A. Hammouti, and A. Wachs, Inertial flow past a finite-length axisymmetric cylinder of aspect ratio 3: Effect of the yaw angle, Physical Review Fluids 4, 044802 (2019).
  * Khair and Chisholm [2018] A. S. Khair and N. G. Chisholm, A higher-order slender-body theory for axisymmetric flow past a particle at moderate reynolds number, Journal of Fluid Mechanics 855, 421 (2018).
  * Fabre  _et al._ [2011] D. Fabre, P. Assemat, and J. Magnaudet, A quasi-static approach to the stability of the path of heavy bodies falling within a viscous fluid, Journal of fluids and structures 27, 758 (2011).