# Rising and settling 2D cylinders with
centre-of-mass offset

Martin P. A. Assen1Equally contributing authors  Jelle B. Will1,2  Chong Shen
Ng1  Detlef Lohse1  Roberto Verzicco1,3,4  Dominik Krug1<EMAIL_ADDRESS>1Physics of Fluids Group, Max Planck UT Center for Complex Fluid Dynamics,
Faculty of Science and Technology, and J.M. Burgers Centre for Fluid Dynamics,
University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands 2
Department of Mechanical Engineering, University of Melbourne, Victoria 3010,
Australia 3 Dipartimento di Ingegneria Industriale, University of Rome ‘Tor
Vergata’,
Via del Politecnico 1, Roma 00133, Italy 4 Gran Sasso Science Institute -
Viale F. Crispi, 7 67100 L’Aquila, Italy

###### Abstract

Rotational effects are commonly neglected when considering the dynamics of
freely rising or settling isotropic particles. Here, we demonstrate that
particle rotations play an important role for rising as well as for settling
cylinders in situations when mass eccentricity, and thereby a new pendulum
timescale, is introduced to the system. We employ two-dimensional simulations
to study the motion of a single cylinder in a quiescent unbounded
incompressible Newtonian fluid. This allows us to vary the Galileo number,
density ratio, relative moment of inertia, and Centre-Of-Mass offset (COM)
systematically and beyond what is feasible experimentally. For certain buoyant
density ratios, the particle dynamics exhibit a resonance mode, during which
the coupling via the Magnus lift force causes a positive feedback between
translational and rotational motions. This mode results in vastly different
trajectories with significantly larger rotational and translational amplitudes
and an increase of the drag coefficient easily exceeding a factor two. We
propose a simple model that captures how the occurrence of the COM offset
induced resonance regime varies, depending on the other input parameters,
specifically the density ratio, the Galileo number, and the relative moment of
inertia. Remarkably, depending on the input parameters, resonance can be
observed for centre-of-mass offsets as small as a few percent of the particle
diameter, showing that the particle dynamics can be highly sensitive to this
parameter.

###### keywords:

flow–structure interactions, vortex-shedding

##  1 Introduction

One of the striking characteristics of freely rising and settling particles at
sufficiently large Reynolds numbers is the existence of non-vertical paths.
This type of motion is a common characteristic for the dynamics of blunt
bodies in general (Ern et al., 2012) and is related to the presence of a
(periodically) oscillating wake with vortex shedding, akin to that forming
behind a fixed geometry (Gerrard, 1966; Perry et al., 1982; Williamson, 1996).
The phase regimes of freely rising or settling bodies are more complicated
than their fixed counterparts, however, due to the inherent coupling between
the motion of the body-fluid interface and the surrounding flow morphology.
This results in a complex, often only quasi-periodic motion which generally is
difficult to predict a priori. Therefore, a complete solution or model of
these systems remains elusive in many cases such that new insights rely on
empirical results and parameter studies based on experiments or numerical
simulations. Understanding and modelling of particle dynamics is important in
many fields, for instance to predict the the spread of (plastic) pollutants in
the ocean (Sutherland et al., 2023).

Properties of the paths and dynamics observed for buoyancy/gravity driven
motion are determined by the strength of the coupling between particle and
fluid. When this coupling is weak (Horowitz & Williamson, 2006) or when the
degrees of freedom of motion are limited (Williamson & Govardhan, 2004), then
the fluid response will be similar to that of a fixed geometry. On the
contrary, regimes exist where particle kinematics are strongly affected by the
fluid motion, and vice versa leading to alterations in flow morphology and
particle trajectory and dynamics, the shedding frequency, but also the drag,
and possibly larger path amplitudes.

Previous studies have often focused solely on the translational dynamics in
relation to the particle-to-fluid density ratio, often disregarding body
rotation, as noted by Mathai et al. (2017). In the present work, we primarily
focus on the effects of rotational coupling. To this end, we vary the internal
mass distribution of freely rising and settling two-dimensional (2D) cylinders
by introducing a Centre-Of-Mass (COM) offset. This approach is motivated by
recent work of Will & Krug (2021b), where the COM of freely rising and
settling spheres was varied experimentally. This study found that the
rotational dynamics of the spheres are strongly affected by the internal mass
distribution, which in turn strongly affects all other aspects of particle
motion. The observed phenomena could be explained via the analogy to a simple
driven harmonic oscillator and expressing the results in terms of a timescale
ratio between the ‘pendulum’ frequency, induced by the offset, and the driving
frequency, which is determined by the vortex-shedding. This model captured
several key features of the particle dynamics when COM offset was present. It
further predicted additional dependencies on the particle-to-fluid density
ratio $\Gamma$, the dimensionless Moment Of Inertia (MOI) of the particle
$I^{*}$, and implicitly on the Galileo number Ga, which is similar to the
Reynolds number $Re$ as it is a measure of the ratio between inertial to
viscous forces but uses an a priori defined buoyancy velocity scale instead of
a dynamical one, and governs the onset and transitions in between the various
wake-structure topologies. However, due to experimental and physical
constraints, the parameter space available in Will & Krug (2021b) is
insufficient to test these conclusively. Therefore, we aim to systematically
investigate these dependencies by means of numerical simulations of 2D
cylinders with COM offset, to show that the underlying physics are similar
between the 2D and the three-dimensional (3D) case, and that the results
presented here can shed light on the remaining open questions.

The problem of freely rising or settling cylinders is an extension of the
classical case of vortex-induced vibrations (Bearman, 1984; Parkinson, 1989;
Govardhan & Williamson, 2000; Williamson & Govardhan, 2004), where a cylinder
is placed in a free stream with only limited degrees of freedom. The applied
restrictions indicated that the degrees of freedom and, therefore, the amount
of particle motion (tuned by body inertia, spring stiffness, and damping of
the supports) strongly influence the wake structure and coupled dynamics.
Williamson & Roshko (1988) presented qualitatively similar results for a
cylinder that was forced periodically in a free stream and classified the
resulting wake patterns as a function of the driving amplitude and frequency.
Building on this, the work by Jeon & Gharib (2004) showed that the type of
vortex-shedding depends on transverse and streamwise oscillations as well as
on their relative phase. Similarly, for elastically mounted cylinders at
subcritical Reynolds numbers ($\Rey\leq 30$) the effects of forced rotations
were examined in recent work by Bourguet (2023), uncovering significant
alterations in flow structure and amplitude of oscillation depending on $\Rey$
and rotational magnitude and frequency. For the case of freely rising and
settling 2D cylinders, a critical mass density ratio ($\Gamma$), the ratio of
particle to fluid mass density, was encountered, marking the threshold between
reduced coupling at high $\Gamma$ where the particle dynamics and its wake
barely influence each other. On the contrary, below this threshold, particles
show large path amplitudes and substantial alterations in the wake vortex
shedding frequency (Horowitz & Williamson, 2006, 2010b). Similar density ratio
related transitions in the regime of motion have also been observed for
spheres (Horowitz & Williamson, 2010a; Auguste & Magnaudet, 2018; Will & Krug,
2021a). Following the same train of thought, the rotational moment of inertia
was also investigated as a relevant parameter, governing the dynamics of
rising and settling 2D cylinders in the numerical work by (Mathai et al.,
2017), as well as experimentally for spheres (Mathai et al., 2018; Will &
Krug, 2021a). Due to these previous observations, we investigate the effects
of $\Gamma$ and MOI separately and in combination with effects induced by a
COM offset.

Before proceeding with the problem definition, caution is warranted when
interpreting the results as the 2D assumption in this work effectively
corresponds to the limiting case of particle motion for very long cylinders
settling or rising in a three-dimensional (3D) environment. Beyond a certain
Reynolds number, the flow will become inherently 3D, even for a cylinder of
infinite length; for a fixed cylinder, this is found to occur around
$\Rey\approx 190$ (Henderson, 1997; Williamson & Brown, 1998; Aleksyuk & Heil,
2023). Moreover, the cylinder length and associated end-effects play an
important role (Inoue & Sakuragi, 2008), such that the motion of these
cylinders becomes inherently 3D, showing both horizontal (azimuthal) cylinder
rotation (around the vertical axis) and fluttering motion (around a horizontal
axis) (Toupoint et al., 2019).

###  1.1 Problem definition and equations of motion

Figure 1: (a) Schematic of the cylinder with the centre-of-mass ($G$)
displaced by a distance $\ell$ from the volumetric centre ($C$). The pointing
vector $\boldsymbol{p}$ is a unit vector in the direction from $G$ to $C$ and
$\theta$ is the angle between $\boldsymbol{p}$ and the vertical
($y$-direction). The forces acting on the body are buoyancy
($\boldsymbol{F}_{b}$) and the remaining fluid forces $\boldsymbol{F}_{f}$ (in
$C$) and gravity $\boldsymbol{F}_{g}$ (in $G$). (b) Schematic depicting the
direction of the Magnus lift force the horizontal component of which is used
together with the horizontal particle acceleration to calculate the phase lag
$\Delta\phi$. (c) Time signals of the horizontal component of the Magnus force
($\boldsymbol{F}_{m}\sim-\omega_{z}v_{y}$, blue) and horizontal particle
acceleration ($a_{x}/g$, red) for three different offset cases, showing
different phase lags. Note that, since $\langle v_{y}\rangle=\mathcal{O}(1)$,
the value on the left $y$-axis is indicative of body rotation.

In this work, we concern ourselves with the dynamics and kinematics of freely
rising and settling 2D cylinders in an otherwise quiescent, infinite fluid.
The fluid phase has a constant mass density $\rho_{f}$ and a kinematic
viscosity $\nu$. The motion of the cylinder is confined to move within the
$xy$-plane. Here, the $y$-axis is anti-parallel to the gravity vector (which
has magnitude $g$). Subscripts $x$ and $y$ are assigned to vector components
in this plane. We denote particle position, velocity, and acceleration by
$\boldsymbol{x}_{p}$, $\boldsymbol{v}$, and $\boldsymbol{a}$, respectively.

The particle (see schematic in figure 1 (a)) has a circular diameter $D$, and
an effective mass $m_{p}$, mass density $\rho_{p}$, and a volume
$\mathcal{V}$, per unit length. The COM of the cylinder (designated with point
$G$) is displaced by a distance $\ell$ from the geometric centre $C$.
Subscripts $C$ and $G$ throughout will refer to these points. A unit pointing
vector $\boldsymbol{p}$ is defined between these points from $G$ to $C$. From
here, we define the angles $\theta$ between $\boldsymbol{p}$ and
$\boldsymbol{e}_{y}$ (the vertical unit vector), and $\theta_{v}$ between
$\boldsymbol{p}$ and $\boldsymbol{v}_{C}$ (the instantaneous particle velocity
at the centre). The buoyancy force acts upwards through point $C$, while the
gravitational force acts downwards through point $G$. The relevant velocity
scale characterising this system is the buoyancy velocity, i.e.
$V_{b}=\sqrt{|1-\Gamma|gD}$.

The dynamics of freely rising and settling buoyancy-driven particles are
governed by the linear and angular momentum balances. In the present work, all
particle dynamics are for unconstrained motion, implying that the only two
forces acting on the geometry are a body force due to gravity
($\boldsymbol{F}_{g}=-\rho_{p}\mathcal{V}g\boldsymbol{e}_{y}$) and the force
exerted by the fluid on the particle $\boldsymbol{F}_{F}$. For convenience we
split this total fluid force into a contribution due to buoyancy and a time-
varying part $\boldsymbol{F}_{F}=\boldsymbol{F}_{b}+\boldsymbol{F}_{f}$, where
$\boldsymbol{F}_{b}=\rho_{f}\mathcal{V}g\boldsymbol{e}_{y}$. Therefore, the
conservation of linear momentum for the cylinder is given by

$\Gamma\dfrac{\textrm{d}\boldsymbol{v}_{G}}{\textrm{d}t}=\dfrac{\boldsymbol{F}_{f}}{m_{f}}+(1-\Gamma)g\boldsymbol{e}_{y},$
(1)

with $\Gamma=\rho_{p}/\rho_{f}$ the mass density ratio, and $m_{f}$ the mass
of the fluid displaced by the particle. It is important to note that (1) is
independent of the COM offset. The angular momentum balance for a 2D cylinder
with COM offset reads

$I_{C}\dfrac{\textrm{d}^{2}\theta}{\textrm{d}t^{2}}=T_{f}-\dfrac{1}{2}m_{p}\gamma
D(\boldsymbol{a}_{C}+g\boldsymbol{e}_{y})\times\boldsymbol{p},$ (2)

where the balance is constructed with respect to the _geometric centre_ of the
particle. Here, $T_{f}$ is the fluid torque due to viscous stresses on the
particle and the additional terms
$(\boldsymbol{a}_{C}+g\boldsymbol{e}_{y})\times\boldsymbol{p}$ on the right-
hand side results from the COM offset. The moment of inertia of the particle
around point $C$ ($I_{C}$) is related to that around point $G$ by the parallel
axis theorem $I_{C}=I_{G}+m_{p}\ell^{2}$. Further, $\gamma$ denotes the ratio
$\gamma=2\ell/D$. The rotational dynamics described in (2) resemble those of a
forced pendulum (Will & Krug, 2021b), for which, when linearised, the natural
frequency is given by

$f_{p}=\tau_{p}^{-1}\ =\dfrac{1}{\pi}\sqrt{\dfrac{\gamma g}{DI^{*}}},$ (3)

with $I^{*}=I_{C}/I_{\Gamma}$, and $I_{\Gamma}=m_{p}D^{2}/8$ a reference value
of a homogeneous cylinder with identical mass $m_{p}$. We can further write
(2) in dimensionless form by introducing a dimensionless timescale
$\tilde{t}=t/\tau_{v}$, where $\tau_{v}=D/V_{b}$ is the vortex shedding
timescale and therefore represents the typical timescale of the forcing in
(2). We further non-dimensionalise the viscous torque term as $T_{f}=\mu
DLV_{b}\,T^{*}_{f}$ (Jordan & Fromm, 1972; Bouchet et al., 2006) and write the
particle acceleration as $\boldsymbol{a}_{C}\sim\boldsymbol{F}_{f}/m_{p}$.
Finally, using $||\boldsymbol{F}_{f}||\sim\rho_{f}LDV_{b}^{2}$ to obtain
$\boldsymbol{a}_{C}=V_{b}^{2}/(D\Gamma)\,\boldsymbol{a}_{C}^{*}$ yields

$\dfrac{\textrm{d}^{2}\theta}{\textrm{d}\tilde{t}^{2}}=\dfrac{1}{Ga\Gamma
I^{*}}\,T^{*}_{f}-\dfrac{\gamma}{I^{*}}\left(\dfrac{1}{|1-\Gamma|}\boldsymbol{e}_{y}+\dfrac{1}{\Gamma}\boldsymbol{a}_{C}^{*}\right)\times\boldsymbol{p},$
(4)

The dimensionless prefactor $\gamma/(|1-\Gamma|I^{*})$ in front of the
pendulum term is proportional to the square of the ratio $\mathcal{T}$ of the
vortex shedding to the pendulum timescale as defined in Will & Krug (2021b).
For a 2D cylinder, this timescale ratio is equal to

$\mathcal{T}=\dfrac{\tau_{v}}{\tau_{p}}=\dfrac{1}{\pi}\sqrt{\dfrac{\gamma}{|1-\Gamma|I^{*}}}.$
(5)

Note that the control parameter $\mathcal{T}$ is solely dependent on the
prescribed particle and fluid properties. It was shown by Will & Krug (2021b)
that this control parameter governs the rotation dynamics of spheres with a
COM offset. We will validate and expand upon this finding in the present work.

In § 2 we will first outline the numerical approach used to obtain the results
as well as the data processing applied to the data set. Then, in § 4 – § 6,
our results are presented and discussed. This is split into four sections,
where § 3 contains a general discussion on effects of COM offset, the
resonance mechanism, and the effect of fluid inertia for freely rising
cylinders. Next, in § 4 the role of the particle-to-fluid density ratio and
how it affects COM-induced phenomena is discussed along with the differences
between rising and settling, followed by a discussion on MOI effects in § 5
and Galileo number effects in § 6. Finally, in § 7 the primary results and
findings are summarised.

##  2 Numerical framework

###  2.1 Fluid phase

The incompressible Navier–Stokes equations describe the flow of an unbounded
Newtonian fluid around the particle and satisfy the boundary conditions at the
body surface and infinity. An approximate computational strategy to model this
configuration is achieved by solving the Navier–Stokes equations on a finite
size domain in a moving reference frame attached to the body. For a perfectly
circular particle, this reference frame does not need to rotate due to the
body’s inherent symmetry. A co-moving frame also allows for a configuration
where the gravity vector is directed towards the outlet such that the wake can
gently leave the domain without disturbing the particle dynamics. The
incompressible Navier–Stokes equations are non-dimensionalised with $V_{b}$
and $D$. For the translating frame, the momentum and continuity equations are
given by (see e.g. Mougin & Magnaudet, 2002; Jenny & Dušek, 2004)

$\frac{d\boldsymbol{u}}{dt}+\nabla\cdot[\boldsymbol{u}(\boldsymbol{u}-\boldsymbol{v}_{C})]=-\nabla
p+\dfrac{1}{\textit{Ga}}\nabla^{2}\cdot\boldsymbol{u}+\boldsymbol{f},$ (6a)
$\nabla\cdot\boldsymbol{u}=0,$ (6b)

with $\boldsymbol{u}$ the fluid velocity vector, $p$ the kinematic pressure
and $\boldsymbol{f}$ the boundary forcing from the immersed boundary method
(IBM) that enforces the no-slip condition (described in detail § 2.2). Note
that the hydrostatic component is absent from $p$ as it is explicitly added to
the forces acting on the cylinder.

The velocity at the inflow is set to zero to simulate an asymptotically
quiescent fluid. At the outflow, a convective boundary is imposed (see e.g.
Kim & Choi, 2006). The side walls of the two-dimensional domain are periodic.
The domain size is set to $60D$ in the gravity direction and $20D$ in the
transverse direction, which was found to be sufficiently large to avoid box
size effects. The grid spacing was kept constant within a square of size $2D$
adjacent to the cylinder. Outside this domain, the mesh spacing expands
linearly towards the edges of the domain. A heat equation is solved to
smoothen the mesh, and the final grids employed for the various cases are
presented in table 1. These grids were chosen such that the particle boundary
layer was resolved by three to four grid points. Tests confirmed that halving
the respective grid spacing did not alter the overall statistics.

Ga | $D/\Delta x$ | $N_{x}\times N_{y}$
---|---|---
50-200 | 50 | $420\times 960$
500 | 84 | $520\times 1104$
700 | 95 | $560\times 1152$
2000 | 161 | $700\times 1440$
Table 1: Overview of the grids. The first column denotes the Galileo number
Ga. The second column represents the number of grid points per diameter of the
cylinder. The third column is the grid resolution for the fluid phase.

###  2.2 Numerical method

The numerical scheme closely follows the immersed boundary projection method
as formulated by Lǎcis et al. (2016) of which a brief overview is presented.
We solve (6a) and (6b) on a staggered grid, where the spatial gradients are
computed using a conservative second-order central finite difference scheme.
The non-linear term of (6a) is advanced in time via the explicit second-order
Adams–Bashforth scheme and the viscous terms via the second-order implicit
Crank–Nicolson scheme:

$\displaystyle\begin{split}\frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{n}}{\Delta
t}&+\frac{3}{2}\hat{N}(\boldsymbol{u}^{n},\boldsymbol{v}_{C}^{n})-\frac{1}{2}\hat{N}(\boldsymbol{u}^{n-1},\boldsymbol{v}_{C}^{n-1})=-\hat{G}\boldsymbol{\varphi}^{n+1/2}\\\
&-\hat{G}\boldsymbol{p}^{n}+\frac{1}{2\textit{Ga}}\hat{L}(\boldsymbol{u}^{n+1}+\boldsymbol{u}^{n})+\hat{H}\boldsymbol{f}^{n+1/2},\end{split}$
(7a)
$\textrm{where}\quad\hat{D}\boldsymbol{u}^{n+1}=0\quad\text{and}\quad\hat{E}\boldsymbol{v}^{n+1}=\boldsymbol{v}^{n+1}+\boldsymbol{\omega}^{n+1}\times\boldsymbol{\mathcal{L}}.$
(7b)

Here, $\hat{N}(\boldsymbol{u},\boldsymbol{v}_{C})$ denotes the non-linear
operator, $\hat{G}$ the gradient operator, $\hat{L}$ the Laplace operator,
$\hat{D}$ the divergence operator, $\hat{H}$ the spreading (regularisation)
operator, $\hat{E}$ the interpolation operator, $\boldsymbol{\varphi}^{n+1/2}$
the discrete incremental pressure, $\boldsymbol{p}^{n}$ the discrete pressure,
$\boldsymbol{\mathcal{L}}$ the Lagrangian marker coordinates with respect to
geometric centre and $\boldsymbol{f}^{n+1/2}$ the discrete analogue of
$\boldsymbol{f}$ in (6a). The interpolation and regularisation matrices
$\hat{E}$ and $\hat{H}$, respectively make use of a discrete three-point
$\delta$-function (Roma et al., 1999).

The Newton–Euler equations are advanced in time via

$I_{B}(\boldsymbol{v}_{C}^{n+1}-\boldsymbol{v}_{C}^{n})=N_{B}\tilde{\boldsymbol{f}}^{n+1/2}+\Delta\boldsymbol{q}_{B}+\boldsymbol{g}^{n}-\boldsymbol{\zeta}^{n},$
(8)

with

$\quad I_{B}=\dfrac{1}{\Delta t}\begin{bmatrix}m_{p}&0&0\\\ 0&m_{p}&0\\\
0&0&I_{C}\end{bmatrix},\quad N_{B}=-\begin{bmatrix}1&\ldots&1&0&\ldots&0\\\
0&\ldots&0&1&\ldots&1\\\
-\mathcal{L}_{y_{1}}&\ldots&-\mathcal{L}_{y_{n}}&\mathcal{L}_{x_{1}}&\ldots&\mathcal{L}_{x_{n}}\end{bmatrix},$
(9)

$\Delta\boldsymbol{q}_{B}\equiv
Q(\boldsymbol{v}^{n}-\boldsymbol{v}^{n-1})/\Delta t$ and $Q$ the matrix that
interpolates the velocity inside the cylinder (cf. Kempe & Fröhlich, 2012).
The time step is limited with $\textrm{CFL}=0.4$.

Vector $\boldsymbol{g}^{n}$ contains the buoyancy force and the torque induced
by the particle weight. The Newton equation in (1) is solved with respect to
the geometric centre and we make use of
$\boldsymbol{v}_{C}=\boldsymbol{v}_{G}+\boldsymbol{\omega}\times\gamma/2D\boldsymbol{p}$
for the transformation of (1). Additionally, the equation of angular
conservation is solved with respect to the geometric centre (see § 1) yielding
an additional $\boldsymbol{a}_{C}\times\boldsymbol{p}$ term. The terms from
the latter contributions are collected in $\boldsymbol{\zeta}^{n}$. Finally,
we have for $\boldsymbol{g}^{n}$ and $\boldsymbol{\zeta}^{n}$:

$\boldsymbol{g}^{n}\equiv\begin{Bmatrix}0\\\ (\Gamma-1)g\\\
-m_{p}gD\gamma/2\sin\theta^{n}\end{Bmatrix},\quad\boldsymbol{\zeta}^{n}\equiv\dfrac{\gamma
m_{p}D}{2\Delta
t}\begin{Bmatrix}\omega^{n}\cos\theta^{n}-\omega^{n-1}\cos\theta^{n-1}\\\
\omega^{n}\sin\theta^{n}-\omega^{n-1}\sin\theta^{n-1}\\\
a_{x}^{n}\cos\theta^{n}+a_{y}^{n}\sin\theta^{n}\end{Bmatrix},$ (10)

with $a_{x}^{n}\equiv v_{x}^{n}-v_{x}^{n-1}$ and $a_{y}^{n}\equiv
v_{y}^{n}-v_{y}^{n-1}$, discretised components relating to
$\boldsymbol{a}_{C}$, respectively.

Equation (7a) together with constraints (7b), and (8) can be rewritten as

$\begin{bmatrix}A&0&G&E^{T}\\\ 0&I_{B}&0&N_{B}\\\ G^{T}&0&0&0\\\
E&N_{B}^{T}&0&0\end{bmatrix}\begin{Bmatrix}\boldsymbol{q}^{n+1}\\\
\boldsymbol{v}_{C}^{n+1}\\\ \boldsymbol{\varphi}^{n+1}\\\
\tilde{\boldsymbol{f}}^{n+1/2}\\\
\end{Bmatrix}=\begin{Bmatrix}\boldsymbol{r}^{n}\\\ \boldsymbol{r}_{B}^{n}\\\
\boldsymbol{0}\\\ \boldsymbol{0}\end{Bmatrix},$ (11)

with $\boldsymbol{q}^{n+1}=R\boldsymbol{v}^{n+1}$, $E=\hat{E}R^{-1}$, diagonal
matrix $R\equiv[\Delta y_{j},\,0;0,\,\Delta x_{i}]$, $A=\hat{M}\left[I/\Delta
t-\hat{L}/(2\textit{Ga})\right]$, and $\boldsymbol{r}^{n}$,
$\boldsymbol{r}_{B}^{n}$ containing the explicit terms. We approximate the
inverse of $A$ as $A^{-1}\approx B=\Delta tM^{-1}$ (Lǎcis et al., 2016, cf. §
3), with $M=\hat{M}R^{-1}$, $\hat{M}\equiv[(\Delta x_{i}+\Delta
x_{i-1})/2,\,0;0,\,(\Delta y_{j}+\Delta y_{j-1})/2)]$ a diagonal matrix. The
solution procedure of (11) is performed via a block-LU decomposition following
a three step procedure

$\begin{bmatrix}A&0\\\
0&I_{B}\end{bmatrix}\begin{Bmatrix}\boldsymbol{q}^{*}\\\
\boldsymbol{v}_{C}^{*}\end{Bmatrix}=\begin{Bmatrix}\boldsymbol{r}^{n}\\\
\boldsymbol{r}_{B}^{n}\end{Bmatrix},$ (12a)
$\begin{bmatrix}G^{T}BG&G^{T}BE^{T}\\\
EBG&EBE^{T}+N_{B}^{T}I_{B}^{-1}N_{B}\end{bmatrix}\begin{Bmatrix}\boldsymbol{\varphi}^{n+1/2}\\\
\tilde{\boldsymbol{f}}^{n+1/2}\end{Bmatrix}=\begin{Bmatrix}G^{T}\boldsymbol{q}^{*}\\\
E\boldsymbol{q}^{*}+N_{B}^{T}\boldsymbol{v}_{C}^{*}\end{Bmatrix},$ (12b)
$\begin{Bmatrix}\boldsymbol{q}^{n+1}\\\
\boldsymbol{v}_{C}^{n+1}\end{Bmatrix}=\begin{Bmatrix}\boldsymbol{q}^{*}\\\
\boldsymbol{v}_{C}^{*}\end{Bmatrix}-\begin{Bmatrix}BG\boldsymbol{\varphi}^{1/2}+BE^{T}\tilde{\boldsymbol{f}}^{n+1/2}\\\
I_{B}^{-1}N_{B}\tilde{\boldsymbol{f}}^{n+1/2}\end{Bmatrix}.$ (12c)

Here, $\boldsymbol{q}^{*}$ in (12a) is solved for via a well tested
factorisation procedure (see e.g. Verzicco & Orlandi, 1996). The solution of
$\boldsymbol{\varphi}^{n+1/2}$ and $\tilde{\boldsymbol{f}}^{n+1/2}$ are
obtained via the PETSc library (Balay et al., 1997, 2019) using the algebraic
multigrid method BoomerAMG as the preconditioner and the general minimum
residual method (GMRES) to solve (12b). This combination of solvers was found
to be robust and converge within 12 to 17 iterations depending on the selected
grid size and time step. A relative tolerance was set to $10^{-13}$, to obtain
solutions that satisfy the divergence free condition and no slip condition up
to the limit of double-precision calculations. Once the solution vector is
found, we update the pressure field (cf. Verzicco & Orlandi, 1996)

$\boldsymbol{p}^{n+1}=\boldsymbol{p}^{n}+\boldsymbol{\varphi}^{n+1/2}-\dfrac{\Delta
t}{2\textit{Ga}}\hat{L}\boldsymbol{\varphi}^{n+1/2}.$ (13)

The additional solving routines for equations (12b) and (12c) were tested to
ensure that they yield machine precision solutions satisfying the divergence
free, and no-slip condition (defined in 7b). The overall solution procedure
was found to provide first-order convergence rate in the $L_{2}$ norm for the
velocity field and first-order in time (owing to the approximation of
$A^{-1}\approx\Delta tM^{-1}$). Multiple validations for fixed, and freely
rising cylinders showed good agreement with previous numerical and
experimental studies. We refer the reader for further details on the immersed
boundary projection approach to Taira & Colonius (2007); Lǎcis et al. (2016).

###  2.3 Data set and processing

In this work, a total of 938 cases have been simulated for different
combinations of the four control parameters: $\mathcal{T}$, Ga, $\Gamma$ and
$I^{*}$. The main goal is to investigate the effect of the COM offset in
combination with the other parameters. For this, we varied $\mathcal{T}$
between 0 and a maximum of 0.6, Ga in the range between 50 up to 2000,
$\Gamma$ between 0.001 and 5, and $I^{*}$ from 0.5 to 16. Compiled input and
output parameters for all cases in our data set are included in the
supplementary materials.

All results presented in this study are obtained after a statistically steady
state has been reached. To ensure this, first, a moving average is computed of
the time trace of the vertical velocity with an averaging window much larger
than a single period of the typical fluctuations. This processed signal is
compared to the terminal velocity of that case (determined by the average of
the last 10% of the time signal). The initial transient is considered to have
ended once the filtered time signal deviates less than 5% from this terminal
velocity. For $\textit{Ga}=200$, this typically is the case after a transient
time of $60D/V_{b}$, which is short compared to the total average simulation
time of $(1.9\times 10^{3}\,D/V_{b})$.

A number of different properties are derived from the simulations to
characterise kinematics and dynamics of the particle path and the surrounding
flow field. In the following, we describe the procedures used to extract these
in detail.

The frequency $f$ of the horizontal path oscillations is determined by the
peak of the power spectrum of $v_{x}(t)/V_{b}$, for which we applied local
peak fitting in order to increase the accuracy of the estimated $f$. In the
case of multiple peaks, the most prominent one is used in subsequent analysis
and data visualisation. Some specific cases featuring multiple peaks are
discussed in § 4.2. The obtained values of $f$ were cross-checked with an
autocorrelation analysis of $v_{x}(t)/V_{b}$, which was found to yield almost
identical results in all cases.

Some additional processing is required to obtain oscillation amplitudes of the
particle rotation and translation due to drift present in the time signals of
$\boldsymbol{x}_{p}(t)$ and $\theta$. The reference $\theta=0$ is either
defined by the direction of the offset or by the initial orientation in the
case of zero offset without loss of generality. To correct for the slow drift
present for some of the cases, we employ a moving averaging filter on the
signal with a window size of approximately $1/f(Ga,\Gamma,\mathcal{T},I^{*})$,
or one full oscillation time. Thus, we obtain a ‘centre-line’
($\boldsymbol{x}_{cl}(t)$, with horizontal mean drift velocity
$v_{d}=\langle|\textrm{d}\boldsymbol{x}_{cl,\,x}/\textrm{d}t|\rangle$,
documented in the supplementary data) which is subtracted from the actual
position and orientation time signal to remove any low frequency effects. The
absolute value of the signal processed this way is used to determine a list of
the individual peak amplitudes ($A$) for the path and ($\theta$) for the
rotational oscillations, the mean of which is denoted by $\hat{A}$ and
$\hat{\theta}$, respectively. Note that, as a consequence, this can mean that
the particle does not exhibit rotational oscillations around $\theta=0$ (where
$\boldsymbol{p}$ is pointing upwards). Instead, especially for small offsets
we observed a behaviour where $\theta$ might drift away from $\theta=0$
followed by a large rotation back to the reference state when the rotational
amplitude becomes large.

The phase lag $\Delta\phi$ between the horizontal component of the Magnus lift
force $\boldsymbol{F}_{m}$ and the horizontal body acceleration
$\boldsymbol{a}_{x}$ is calculated via cross-correlation of these quantities.
Similarly, $\Delta\psi$ denotes the phase lag between the angular acceleration
$\alpha$ and the fluid torque $T_{f}$. The lag obtained from the cross-
correlation is divided by the length of an oscillation period $1/f$ and then
expressed as a phase angle ranging from -180∘ to 180∘. The respective
components that define $\Delta\phi$ are illustrated in figure 1(b). Figure
1(c) provides three examples of signals with varying $\mathcal{T}$ showing a
negative, zero and positive value of $\Delta\phi$.

##  3 General effect of the COM offset on dynamics and kinematics

###  3.1 Particle kinematics and wake structures

Figure 2: Snapshots of particle trajectories and wake structures of rising
cylinders with Galileo number $\textit{Ga}=200$ and density ratio $\Gamma=0.5$
for six different centre-of-mass offsets $\gamma$ (and $\mathcal{T}$) (a–f).
The offset increases from left to right as indicated by the listed parameters
in the top left of each subfigure. Particle trajectories are indicated by the
black lines, the grid spacing has dimensions of the particle diameter $D$.
Coloured contours represent the normalised vorticity field
($\omega_{z}D/V_{b}$).

We present figure 2 to give an impression of how the wake patterns and
particle kinematics change in the presence of COM offset. These snapshots
display the non-dimensionalised fluid vorticity field
($\omega_{z}=\partial_{y}u_{x}-\partial_{x}u_{y}$) along with particle tracks
(black lines) for six different COM offsets in the range $\gamma\in[0,\,1.23]$
(increasing from left to right). All cases here are for $\textit{Ga}=200$,
$\Gamma=0.5$, and $I^{*}=1$.

The cylinder with zero COM offset in figure 2a is seen to rise almost
straight, with regular vortex shedding occurring at the same frequency as the
path oscillations. This vortex pattern, where two single vortices of opposite
vorticity are shed during a single oscillation cycle, is the so-called “2S”
mode (Williamson & Roshko, 1988). No visible effect of the COM offset is
observed for cases up to $\mathcal{T}=0.174$ (figure 2b), but beyond this
value, $\mathcal{T}=0.201$ shown in panel figure 2c, remarkably different
kinematics are encountered. Both amplitude and wavelength of the path
oscillations are significantly larger in this case, and the wake now exhibits
an irregular vortex shedding pattern as can be seen in supplementary video 1.
For this case it is observed that the wake structure intermittently switches
between several modes: (i) path oscillations with no significant shedding
events, (ii) one single vortex pair as in the 2S-regime per oscillation, or
(iii) even two vortex pairs as typically found in the so-called “2P”-mode.
Note that these different modes appear to alternate without any noticeable
long time scale pattern. This chaotic shedding pattern occurs for cases close
to what we will call “resonance”, where the rotational forcing induced by the
path oscillations occurs at the same frequency as the inherent pendulum time
scale. For even higher values of $\mathcal{T}$ beyond resonance, represented
by $\mathcal{T}=0.285$ in panel figure 2d, we observe that the large amplitude
path oscillations persist albeit with a reduced wavelength. Further, the
vortex shedding returns again to an unperturbed 2S mode now with staggered
vortex cores due to the strong path oscillations. Finally, with even larger
offsets, figure 2(e, f), the amplitude of the path oscillations begins to
gradually reduce, returning to a state very much like that for the zero offset
case (see supplementary video 1 for $\mathcal{T}>0.3$ cases). The results
shown here are representative of the $\Gamma$-range where the resonance
phenomenon is present. In the following, we will evaluate how this resonance
behaviour depends on all of the other governing parameters.

###  3.2 On the importance of fluid inertia

Figure 3:  (a) Mean rotational amplitude $\hat{\theta}$ as a function of
Galileo number versus the timescale ratio $\mathcal{T}$, here $\Gamma=0.6$ and
$I^{*}=1$. (b) Schematic showing the parameters of the fluid inertia model.
(c) $\hat{\theta}$ for the same cases as in (a) plotted against the modified
timescale ratio $\tilde{\mathcal{T}}$, which includes the effects of a Galileo
number dependent added fluid inertia as per (14b). (d) Thickness of the fluid
inertia layer $\delta$ and (e) added inertia as a function of Ga, based on
empirical collapse of the data.

In order to investigate the effect of fluid inertia, we first consider the
mean rotational amplitude ($\hat{\theta}$) for a constant density ratio
($\Gamma=0.6$) as a function of the timescale ratio $\mathcal{T}$ as shown in
figure 3(a). Focusing initially on the case $\textit{Ga}=200$ (corresponding
to figure 2), we see that at $\mathcal{T}=0$ the rotational amplitude is small
$\hat{\theta}=0.4^{\circ}$. Introducing a small amount of offset results in a
marginal increase in this amplitude up to $\hat{\theta}=1.4^{\circ}$ at
$\mathcal{T}=0.16$. However, around $\mathcal{T}=0.2$, there is a strong
increase reaching a maximum amplitude of more than $\hat{\theta}=35^{\circ}$
at $\mathcal{T}=0.225$. This rapid increase is associated with the resonance
phenomenon that was already visible in figure 2(c). Beyond this point, the
amplitude decreases gradually with increasing $\mathcal{T}$.

When comparing results across different Ga numbers, the same characteristic
behaviour is observed for all cases in figure 3(a). However, the value of
$\mathcal{T}$ at which resonance appears (marked by the steep increase in
$\hat{\theta}$) is consistently shifted towards higher values as Ga decreases.
Such a variation with Ga is not surprising since the definition of
$\mathcal{T}$ does not incorporate any viscous effects. However, for low
values of Ga one would expect the Stokes layer surrounding the particle to
contribute significantly to the total rotational inertia of the body and
thereby also to modify the particle pendulum timescale. We can account for
this effect by additionally including the rotational inertia $I_{a}$ resulting
from the Stokes layer with thickness $\delta$ as illustrated in figure 3(b) in
our analysis. As a result, the modified pendulum frequency and timescale ratio
of the system become

$\tilde{f}_{p}=\dfrac{1}{\pi}\sqrt{\dfrac{\gamma
g}{D(I^{*}+I^{*}_{a}/\Gamma)}},$ (14a) and
$\tilde{\mathcal{T}}=\dfrac{1}{\pi}\sqrt{\dfrac{\gamma}{|1-\Gamma|\left(I^{*}+I^{*}_{a}/\Gamma\right)}},$
(14b)

respectively. Here $I^{*}_{a}$ is the dimensionless fluid inertia defined as
$I^{*}_{a}\equiv 8I_{a}/(m_{f}D^{2})$, the ratio of the Stokes layer’s
rotational inertia to that of the displaced fluid. The total rotational
inertia is thus given by $I^{*}+I^{*}_{a}\Gamma$. We assume that the thickness
of this Stokes layer scales as $\delta\sim 1/\sqrt{\textit{Ga}}$ (Williamson &
Brown, 1998; Schlichting & Gersten, 2003; Mathai et al., 2018), which for a
cylinder leads to

$I^{*}_{a}(\textit{Ga})=\dfrac{8c_{1}}{\sqrt{\textit{Ga}}}+\dfrac{24c_{1}^{2}}{\textit{Ga}}+\dfrac{32c_{1}^{3}}{\textit{Ga}^{3/2}}+\dfrac{16c_{1}^{4}}{\textit{Ga}^{2}},$
(15)

with $c_{1}$ as the only free parameter. We find that choosing $c_{1}=2.3$
results in a reasonable collapse of the resonance regime for different Ga when
plotting $\hat{\theta}$ against $\tilde{\mathcal{T}}$ as shown in figure 3(c).
The corresponding thickness of the Stokes layer and magnitude of the added
fluid inertia as a function of Ga are provided in figures 3 (d, e),
respectively. For $\textit{Ga}=200$, the thickness of the fluid layer is
approximately 0.33 particle radii and the rotational inertia amounts to about
2 times that of the displaced fluid. Beyond $\textit{Ga}=\mathcal{O}(10^{3})$,
the value of $I^{*}_{a}$ changes much more slowly, explaining the weak $Ga$
dependence at higher Ga observed in figure 3(a) as well as in previous work on
spheres (Will & Krug, 2021b). Note, however, that $I_{a}^{*}$ is still 0.72 of
the displaced fluid mass at $\textit{Ga}=1000$ for cylinders and therefore by
no means negligible. We performed an estimate of the history torque to confirm
that the obtained values for $I_{a}^{*}$ are realistic. A complete discussion
of this for both cylinders and spheres is provided in Appendix A.

###  3.3 Who’s driving?

Figure 4:  (a, b) Results for $Ga=200$ and $\Gamma=0.5$ for two cases; one
without offset (grey line) and one with offset (red line). (a) Dimensionless
horizontal velocity of the cylinder ($v_{x}$) and (b) dimensionless rotation
rate (rad.) versus dimensionless time. During these runs at $t=0$ the
$\boldsymbol{a}_{C}\times\boldsymbol{p}$ term for (2) is turned off, showing
that in absence of this coupling term the dynamics of particles with offset
almost completely revert back to those of particles without offset. (c) Phase
lag $\Delta\psi$ between the rotational particle acceleration ($\alpha$) and
the viscous torque ($T_{f}$).

When considering the right hand side of (2), there are two potential drivers
of the rotational motion, the viscous torque $T_{f}$ and the translational-
rotational coupling term $\boldsymbol{a}_{C}\times\boldsymbol{p}$, the latter
being a consequence of the COM offset. Here, we will investigate their
respective role with respect to the resonance behaviour. We know from the
analysis in §3.2 that the maximum in rotational amplitude is related to
resonance between the vortex shedding timescale $\tau_{v}$ and the pendulum
timescale $\tau_{p}$. However, both the viscous and translational driving will
occur at the vortex shedding frequency, such that a distinction of their
effects is not possible on this basis only. Answering the question of ‘who’s
driving’ also provides insight in the effectiveness of COM offset in specific
regimes of motion.

In order to untangle the effects of both contributions, simulations were
performed where, after a statistically steady state had been reached, the
($\boldsymbol{a}_{C}\times\boldsymbol{p}$)-term was turned off in the
integration of (2). In figure 4 (a), the horizontal component of the particle
velocity $v_{x}$ is shown as a function of time for these runs at
$\textit{Ga}=200$, $\Gamma=0.5$, and $\mathcal{T}=0$ (grey line) and
$\mathcal{T}=0.285$ (red line). At $t=0$, the coupling term
$\boldsymbol{a}_{C}\times\boldsymbol{p}$ is turned off for the case with
offset. Figure 4 (b) displays the rotation rate $\omega$ for the same
simulations. These results clearly indicate that as the coupling-term is
turned off, the particle dynamics return to those of the case without offset.
Note here that the pendulum term ($~{}\boldsymbol{e}_{y}\times\boldsymbol{p}$)
is still present for $t\geq 0$, but evidently it has no effect without
translational driving of the rotational dynamics. Therefore, we conclude that
the rotational resonance phenomenon is linked to the translational coupling.
As a consequence, we expect COM offset to have no impact on particle dynamics
for cases where no horizontal path oscillations (i.e. no horizontal
accelerations) are present, e.g. at low $Ga$. This also suggests that the
resonance behaviour might also be triggered by outside periodic forcing, as
would be present in a turbulent flow environment. It would be interesting to
study how the settling/rising velocities of low Galileo number bodies with COM
offset are affected in turbulence via this mechanism.

On the role of $T_{f}$, it is further instructive to consider the phase lag
$\Delta\psi$ between $T_{f}$ and the rotational acceleration $\alpha$, which
is shown in figure 4 (c) for the full range for $\Gamma$ and $\mathcal{T}$ at
$\textit{Ga}=200$. For zero or very small offsets, $T_{f}$ is driving the
(weak) rotational dynamics as evidenced by $\alpha$ and $T_{f}$ being close to
in phase. However, as the offset increases towards resonance and beyond,
$\Delta\psi$ switches swiftly to values close to 180∘, such that the viscous
torque predominantly acts as damping in these cases. In essence, these trends
also hold for higher Ga. However, the dynamics become somewhat more chaotic at
higher Ga, as will be shown in §6, resulting in slightly lower values of
$\Delta\psi$ on the order of $120$ to $150^{\circ}$.

##  4 The effects of density ratio on COM offset

###  4.1 Frequency of oscillation

In the following sections, we will discuss how the effect of the COM offset
varies with the density ratio. In doing so, we focus on the representative
case of $\textit{Ga}=200$ and $I^{*}=1$. We first consider the frequency of
the path oscillations ($f$), as this parameter also corresponds to the
frequency at which the rotational dynamics are forced. In figure 5(a), we plot
the data in the form of the Strouhal number

$\mbox{{Str}}=\dfrac{fD}{V_{b}}.$ (16)

as a function of $\mathcal{T}$ and $\Gamma$. The marker colour in the figure
indicates the exact $Str$ obtained from the simulations as can be read from
the legend, the iso-contours and background colours represent a linear
interpolation of this data. The transparent white area bordered by the black
dashed line indicates the region where $\gamma>\sqrt{0.5}$. This corresponds
to the theoretical state where $I_{G}$, the MOI of the particle around the
centre-of-mass, is zero in accordance with the parallel axis theorem
($I_{C}=I_{G}+m_{p}\ell^{2}$) as a consequence of keeping $I_{C}\equiv
m_{p}D^{2}/8=const.$ (i.e. $I^{*}=1$). Furthermore, we also add a line where
$\gamma=1$, i.e. when the point $G$ lies on the particle edge. Results within
this marked region are therefore not physically viable, yet still satisfy the
governing equations.

Figure 5: Results for the path oscillation frequency ($f$) of rising particles
at $Ga=200$ as a function of $\mathcal{T}$) $\Gamma$. In (a), the marker
colour indicates the dimensionless frequency ($Str=fD/V_{b}$) according to the
colour bar provided below. The marker type indicates the different regimes in
terms of the resonance behaviour discussed in the following. The isocontours
are based on a linear interpolation of the data. Dashed white lines represent
isocontours of $\tilde{\mathcal{T}}$, the timescale ratio including effects of
fluid inertia. (b) Horizontal particle position over the cylinder diameter
($x/D$) as a function of dimensionless time grouped in three values of
$\mathcal{T}$ for three values of the $\Gamma$ as indicated by the line
colours showing characteristic behaviour for each. (c) Ratio of the frequency
($f$) of the path oscillations over the pendulum frequency $f_{p}$ vs. the
timescale ratio $\tilde{\mathcal{T}}$. Here the marker colour indicates
$\Gamma$ as listed in the legend below the figure. The two dashed black lines
show a constant value of $Str$. Both of these show the collapse of COM and
$\Gamma$ effects in terms of this parameter. We further see that the results
also collapse with the results from spheres with COM offset (Will & Krug,
2021b) shown as black symbols. The inset of the figure shows the same data as
(a) plotted as $Str$ vs. $\tilde{\mathcal{T}}$.

Considering the zero offset ($\mathcal{T}=0$) cases in figure 5(a) first, we
find that Str varies quite significantly from $\mbox{{Str}}=0.195$ at high
density ratios ($\Gamma\geq 0.2$) to $\mbox{{Str}}=0.127$ for $\Gamma\leq
0.1$. This transition appears to be quite sudden, suggesting the existence of
critical density ratio as previously observed for rising and settling blunt
bodies (Namkoong et al., 2008; Horowitz & Williamson, 2010b; Mathai et al.,
2017; Auguste & Magnaudet, 2018; Will & Krug, 2021a). The change in the path
frequency at the lowest $\Gamma$ is also associated with an increase in the
path amplitude as is evident from the trajectories at $\mathcal{T}=0.15$
(which resemble those at $\mathcal{T}=0$) in figure 5(b).

Now let us consider the effects of COM offset for varying $\Gamma$ (i.e.
moving vertically in the figure). Depending on $\Gamma$, three distinct
effects of increased offset can be observed. First of all, for $\Gamma\geq
0.9$ (marked by square symbols throughout), we find that increasing COM offset
has almost no effect on the oscillation frequency. There is only a slight
decrease for $Str$ at extreme offset as can best be seen in the inset of
figure 5(c). The general lack of response to COM offset in this regime can be
explained by considering the rotational equation of motion as presented in
(4). Remember here that the system is similar to a driven damped harmonic
oscillator where the pendulum term is analogous to the spring stiffness, the
viscous torque is the damping term, and the accelerated reference frame
($\boldsymbol{a}_{C}\times\boldsymbol{p}$) provides the driving. The restoring
torque is proportional to $|1-\Gamma|^{-1}$ and therefore goes to infinity for
$\Gamma\to 1$. This is not the case for the driving term which scales
according to $\Gamma^{-1}$. Therefore, when the body becomes close to
neutrally buoyant, the pendulum torque goes to infinity and as a result the
forcing can not rotate the body significantly enough to induce any
circulation. Therefore, there will be no Magnus force and no rotational-
translational feedback loop leading to resonance. Thus, for the cases where
$\Gamma$ is close to unity, the oscillation frequency (as well as other output
parameters) of the body remain unaffected by the offset.

The second regime is characterised by a sharp transition in particle dynamics
where in a narrow range of $\mathcal{T}$ the dynamics switch between the base
state (near identical to $\gamma=0$) and the resonance state. This is best
shown in the inset of figure 5(c) where we see that at low values of
$\tilde{\mathcal{T}}$ for intermediate density ratios ($0.2\leq\Gamma\leq
0.8$) Str stays constant at approximately 0.195 (upper branch). However, as
the offset increases there is a sharp jump to the lower branch of Str. The
upper branch corresponds to a system state with minimal body rotation and
translation, and in the lower branch the vortex shedding latches on to body
motion and is affected by the pendulum frequency. The cases $\Gamma=0.2$ and
$0.8$ are edge cases and show characteristics of their neighbouring regimes.

Finally, the third regime is characterised by a gradual transition to the
resonance state and is encountered for $\Gamma\leq 0.1$. Here we find that
even at zero offset they are already following the lower branch in figure 5(c.
Since the particle is already exhibiting path oscillations and minute
rotational oscillations even at zero offset, no critical threshold of offset
needs to be exceeded for the coupling to begin occurring. For these cases,
even at $\tilde{\mathcal{T}}$ below resonance, we already see offset affecting
the particle dynamics. The footprint of these three regimes is also evident in
the amplitude and spatial path frequency as shown in figure 5(b). For high
$\Gamma$, there is no effect of increasing offset, at intermediate density
ratios we see a large difference between different $\mathcal{T}$, and at low
$\Gamma$, we observe large path oscillations even at small/zero offset.

Beside the jump at the onset of resonance, $Str$ also varies significantly
beyond the resonance state. The isocontours of Str in this region
approximately follow the lines of constant $\tilde{\mathcal{T}}$ (white dashed
lines) and in particular the minimum of Str coincides roughly with
$\tilde{\mathcal{T}}=0.08$. The correlation between Str and
$\tilde{\mathcal{T}}$ is explicitly shown in the inset of figure 5(c). This
plot also highlights the existence of two branches of the system state and the
fact the COM offset can trigger a transition between the two, indicated for
each $\Gamma$ by the coloured dashed section of the lines. The collapse of the
data on these two curves is not trivial and underlines the validity of the
Stokes layer argument at the core of the definition of $\tilde{\mathcal{T}}$.
As $\tilde{\mathcal{T}}$ becomes very large, Str appears to return to the
trend at large $\Gamma$ where higher density ratios have a slightly higher
Str. It is further clear that while $\tilde{\mathcal{T}}$ is the relevant
parameter to describe the behaviour after the transition from the low $Str$ to
the high $Str$ state, the transition itself does not coincide with
$\tilde{\mathcal{T}}=\textrm{const.}$, but occurs at lower values of
$\tilde{\mathcal{T}}$ for lower $\Gamma$. The case with the lowest density of
$\Gamma=0.001$ spans only a tiny range in terms of $\tilde{\mathcal{T}}$ even
for the largest offsets. This could explain why there is no noticeable
variation in the particle behaviour for this density ratio even at large
offsets. However, since the driving term is proportional to $1/\Gamma$, it
will likely dominate the pendulum torque, which does not diverge for small
$\Gamma$.

As mentioned at the beginning of this section, the frequency of the path-
oscillations is important for the driving of the rotational dynamics through
equation (2). As with any harmonic oscillator the parameter of prime
importance is the ratio of the driving to natural frequency of the system
$f/\tilde{f}_{p}$, which we show in the main panel of figure 5(c) as a
function of the timescale ratio $\tilde{\mathcal{T}}$. Curves of constant Str
corresponding to the two different states are indicated by the black dashed
lines. Importantly, we find that $f/\tilde{f}_{p}=1$ occurs around
$\tilde{\mathcal{T}}=0.11$, which corresponds to the bold white dashed line in
figure 5(a). We further see in figure 5(c) that the path oscillation frequency
of the body appears to be drawn towards $f_{p}$ as it begins to deviate from
$\mbox{{Str}}=0.127$ to meet $f/\tilde{f}_{p}=1$, consistent with the so
called lock-in phenomenon (Bishop & Hassan, 1964; Bearman & Obasaju, 1982).
The region of (approximate) frequency lock-in, ranging from
$0.09\leq\tilde{\mathcal{T}}\leq 0.12$ is indicated by a grey shaded area in
the figure background throughout this work. Finally, we included the results
for rising and settling spheres with COM offset from the work by Will & Krug
(2021b) as black circles in figure 5(c). The good agreement with the present
results suggests that the underlying physics of the resonance mechanism are
indeed the same and that results and trends presented here are also relevant
for spherical bodies in a 3D flow environment.

###  4.2 On the transition to resonance

Figure 6:  (a) Single sided amplitude spectrum ($\mathcal{F}$) based on the
particle horizontal particle velocity $(v_{x}/V_{b})$ normalised by the
maximum amplitude. (b) Dimensionless horizontal velocity and rotation rate
($\textrm{rad.}/s$) versus time for $\Gamma=0.8$ and $\mathcal{T}=0.327$. We
see that the dynamics exhibit a cyclical behaviour on a timescale much greater
than that of the vortex-shedding dynamics. This behaviour is split into three
parts as indicated by the colours in the background of the figure. (c)
Instantaneous Strouhal number as a function of time for $\Gamma=0.8$ and
$\mathcal{T}=0.327$, calculated based on the peak-to-peak times of
$|v_{x}/V_{b}|$. (d–f) Vortex shedding and path oscillations correlating to
the modes in (b, c). (d) Very low frequency oscillations of minimal amplitude,
attached wake is very large. (e) The buildup vorticity is rapidly shed in the
wake at a high frequency, resulting in small amplitude high frequency path
oscillations. (f) Slower periodic vortex shedding with larger amplitude path
oscillations, the attached vorticity slowly grows throughout this phase until
the cycle begins anew.

In this section, we will discuss the transition to resonance in the
intermediate $\Gamma$ regime, i.e. for $0.2\leq\Gamma\leq 0.8$ in more detail.
In the range $0.3\leq\Gamma\leq 0.7$, the transition from the high Str number
mode to the low Str one occurs within a narrow band of $\mathcal{T}$. This is
also visible from figure 6(a), where the power spectra normalised by the
maximum amplitude ($\mathcal{F}$) of $v_{x}/V_{b}$ are shown for
$\mathcal{T}=0.159$, $0.195$, and $0.225$ in the vicinity of the transition
point at $\Gamma=0.6$ (yellow diamonds, pink pentagons and purple hexagons
respectively). For both, $\mathcal{T}=0.159$ and $\mathcal{T}=0.225$, the
spectra feature singular peaks only at $\mbox{{Str}}\approx 0.1$ and
$\mbox{{Str}}\approx 0.2$, respectively. The former peak also dominates for
the intermediate case $\mathcal{T}=0.195$, however, a weaker secondary peak at
$\mbox{{Str}}\approx 0.21$ is also seen to emerge at this offset value.
Similar trends can be observed across the range $0.3\leq\Gamma\leq 0.7$ with
varying ratios of relative peak height, suggesting that the transition between
modes happens in a narrow band of $\mathcal{T}$, but is not entirely discrete.

In §4.1, we mentioned that $\Gamma=0.2$ and $0.8$ were on the edges of the
$\Gamma$-range for which a sharp transition to the resonance regime was
encountered. We will investigate these cases in more detail here. For
$\Gamma=0.2$, the range of $\mathcal{T}$ where multiple modes are observed
widens significantly as compare to $0.3\leq\Gamma\leq 0.7$. We observe
multiple peaks in the spectra for $0.138\leq\mathcal{T}\leq 0.225$ as
exemplified by the case shown in figure 6(a) (green circles). This extended
range is most likely due to the intrinsic rotational and transitional
oscillations present at $\gamma=0$ for $\Gamma$ close to 0 and is similar to
the cases of $\Gamma\leq 0.1$. However, looking at figure 5(c), we still
observe that the cases at $\Gamma=0.2$ and $\gamma\approx 0$ follow the upper
branch in terms of $f/f_{p}$ and Str, making the behaviour transitional
between the $\Gamma$-regimes.

At $\Gamma=0.8$, we find only a single case ($\mathcal{T}=0.327$) for which
the system state jumps to the lower branch in figure 5(c). In figure 6(a)
(blue squares), we observe that this jump is accompanied by a wide range of
frequency peaks. The occurrence of multiple peaks originates from a very
unique frequency and amplitude-modulation cycle in the fluid-structure
interaction, the signature of which is shown in terms of the time-evolution of
$v_{x}/V_{b}$ and $\omega D/V_{b}$ in figure 6(b). For both quantities, a
modulation of the amplitude but also of the frequency is evident at timescales
much larger than that of the path oscillations. This behaviour stands in stark
contrast to the rest of the cases which exhibit very regular periodic motion.
Specifically, the parameter combinations that show this characteristic
behaviour are: $\Gamma=0.8$ with $\mathcal{T}=0.327$ and $\mathcal{T}=0.365$,
and $\Gamma=0.7$ with $\mathcal{T}=0.411$. A video showing this behaviour
along with the vorticity in the wake can be found in supplementary video 2.

To quantify the frequency variations, we plot the instantaneous Str in figure
6(c) based on determining the distance between the maxima and minima in the
signal $v_{x}/V_{b}$. Instances where the frequency is relatively low are
indicated by the blue regions in figure 6(b,c). The corresponding particle
kinematics, showing a marked reduction in the lateral amplitude, and the
associated vortical structures for the low Str mode are illustrated in figure
6(d). In this case, the cylinder rises almost vertically and the length of the
attached wake is at its maximum extent. After this, in the transitional period
marked in yellow, the path amplitude remains low but the frequency of the
oscillations increases. Figure 6(e) shows that this is linked to the rapid
shedding of the build-up attached wake, quickly reducing its length. This
state is similar to the dynamics observed at higher $\Gamma$. Finally, in the
period indicated by red shading, the lateral amplitude is relatively large and
the oscillation frequency is intermediate. In the corresponding figure 6(f),
it is seen that over the course of a number oscillation periods the attached
wake slowly grows again until this cycle repeats. Due to the large amplitude
and longest duration of this phase, the red region manifests as the strongest
peak in the Fourier spectrum and thus the result for $Str$. This behaviour is
characteristic for the cases near $\Gamma=0.8$ and $\mathcal{T}=0.327$ at this
Galileo number and is indicative of the density regime transitions, where the
dynamics exhibit signs of both regimes. These observations highlight that in
the transition range, multiple states can coexist.

###  4.3 Drag coefficient and Magnus force

Figure 7: (a) Particle drag coefficient as a function of the particle-to-fluid
density ratio $\Gamma$ and the timescale ratio $\mathcal{T}$ for rising
particles at $Ga=200$ and $I^{*}=1$. (b) Correlation between $C_{d}$ and the
phase lag $\Delta\phi$ between the horizontal Magnus force and the horizontal
component of instantaneous acceleration. (c) Non-dimensional mean magnitude of
the fluctuating velocity component
$v^{*}=\sqrt{\langle\boldsymbol{v}_{C}^{\prime 2}\rangle}/V_{b}$ as a function
of the timescale ratio $\mathcal{T}$.

In this section, we will concern ourselves with the mean vertical velocity,
i.e. the terminal rising/settling velocity, which is of particular practical
relevance. We define the drag coefficient $C_{d}$ obtained from the time-
averaged vertical force balance between the buoyancy and the drag force, given
that the particle has reached terminal velocity. For a two-dimensional
cylinder, this results in

$C_{d}=\dfrac{\pi|1-\Gamma|gD}{2\langle v_{y}\rangle^{2}}.$ (17)

Note that in this definition $C_{d}$ is essentially a dimensionless rising or
settling velocity. When plotted as a function of $\Gamma$ and $\mathcal{T}$
(figure 7(a)), the drag coefficient exhibits considerable variations almost up
to a factor of 4 across the parameter space that are predominantly induced by
COM effects. At $\mathcal{T}=0$, $C_{d}$ is found to be lowest for
$\Gamma=0.2$. Moving to the right in the figure, i.e. towards increasing
$\gamma$, the previously (§4.1) defined $\Gamma$ regimes can again be noted.
For $\Gamma\geq 0.8$, no increase in $C_{d}$ is observed as a consequence of
the COM offset. However, for $\Gamma\leq 0.7$ the resonance behaviour
manifests itself in a strong increase in $C_{d}$. These trends are explicitly
plotted in terms of $\tilde{\mathcal{T}}$ in figure 7 (b, c). For $\Gamma\leq
0.7$, the resonance behaviour reaches a maximum for
$\tilde{\mathcal{T}}\approx 0.11$. This is indicated in figure 7(a) by the
bold white dashed line, and even more evident from the location of the peak in
$C_{d}$ in figure 7(b). The value of $\tilde{\mathcal{T}}=0.11$ corresponds to
$f/f^{*}_{p}=1$ in figure 5(c), i.e. where the driving frequency $f$ and the
pendulum frequency are identical. The magnitude of the peak drag monotonically
increases with $\Gamma$. Beyond $\tilde{\mathcal{T}}=0.11$, $C_{d}$ gradually
decreases again in all resonance cases. Finally, figure 7(c) also shows that
for all cases at large offsets, even those that did not exhibit resonance
(i.e. $\Gamma\geq 0.8$), the drag decreases slightly. It appears that the
magnitude of the decrease is inversely correlated with the mass density,
resulting in larger reduction for lighter particles (figure 7(c)). This
phenomenon does not appear to occur at fixed values of $\mathcal{T}$ or
$\tilde{\mathcal{T}}$. It is noteworthy, though, that a similar drag reduction
was encountered around similar values of $\tilde{\mathcal{T}}$ for settling
and rising spheres with COM offset (Will & Krug, 2021b).

In previous work by Will & Krug (2021b), the connection was made between the
drag increase and the maximum in the enhancement of horizontal particle
acceleration ($a_{x}$) through the rotation induced Magnus lift force
($F_{m,\,x}\sim-\omega_{z}v_{y}$). Besides $F_{m,\,x}$, lateral accelerations
can also be driven by pressure fluctuations induced by the vortex shedding. To
study the enhancement of the path oscillations via the Magnus force, we
consider the phase lag ($\Delta\phi$) between $F_{m,\,x}$ and $a_{x}$.
Examples of time series with different phase lags for three values of
$\mathcal{T}$ are shown in figure 1(c). When $F_{m,\,x}$ and $a_{x}$ are in
phase ($\Delta\phi=0$), the enhancement of the horizontal particle motion by
the Magnus force is maximum. The connection between $\Delta\phi$ and $C_{d}$
is established for the current data set in the inset of figure 7(d). In the
main panel of figure 7(d), the phase lag is plotted explicitly vs.
$\tilde{\mathcal{T}}$. Again there is excellent collapse of the data onto two
branches, representing the oscillating and non-oscillating states, identical
to those encountered for $Str$ in §4.1. We further see that
$\Delta\phi=0^{\circ}$ around $\tilde{\mathcal{T}}=0.11$ for all cases where
resonance is present, whereas acceleration and Magnus force are significantly
out of phase ($\Delta\phi\approx-60^{\circ}$) in the same range on the lower
branch. This point is emphasised by the inset of figure 7(d), where the peaks
in $C_{d}$ are seen to align with $\Delta\phi=0$. The good agreement with the
experimental sphere data of Will & Krug (2021b), included in figure 7(d),
again underlining the fact that the resonance phenomenon in 2D is indeed
comparable to that in 3D.

###  4.4 Drag correlations

Figure 8:  (a) Normalised velocity fluctuations
$v^{*}=\sqrt{\langle\boldsymbol{v}_{C}^{\prime
2}\rangle_{\textrm{rms}}}/V_{b}$ and (b) normalised rotational fluctuations
$\omega^{*}=\langle\omega\rangle_{\textrm{rms}}D/V_{b}$ (in degrees) as a
function of $\Gamma$ and $\mathcal{T}$. Correlations between the particle drag
coefficient and $v^{*2}$ (c) and $\omega^{*2}$ (d). The density ratio
($\Gamma$) in represented by the line colour. The marker infill, white of grey
represents $\Delta\phi<0$ or $\Delta\phi\geq 0$, respectively.

The question to what extent the drag of freely rising and settling bodies is
correlated to path-oscillations and/or particle rotations is subject to an
ongoing discussion. It is indisputable that the presence of horizontal
oscillations affects the overall drag coefficient (Horowitz & Williamson,
2010a). However, the presence of rotations was also clearly found to play a
prominent role (Namkoong et al., 2008; Auguste & Magnaudet, 2018; Mathai et
al., 2018). In the work on spheres with COM offset by Will & Krug (2021b), it
was shown that the drag correlated better with the mean rotation rate than
with the amplitude of the path oscillations or the horizontal velocity
fluctuations for cases at or beyond resonance. On the other hand, for zero
offset, the drag appeared to correlate equally well with both in the three-
dimensional (3D) chaotic regime when varying the MOI of spheres (Will & Krug,
2021a) and both of these quantities did not result in an adequate prediction
of drag for the spiralling regime.

To add to this, we investigate how variations in $C_{d}$ relate to the
presence and strength of rotational and translational fluctuations in the
present data. To this end, we define the dimensionless fluctuating velocity
$v^{*}=\sqrt{\langle\boldsymbol{v}_{C}^{\prime
2}\rangle_{\textrm{rms}}}/V_{b}$, where
$\boldsymbol{v}_{C}^{\prime}=\boldsymbol{v}_{C}-\langle\boldsymbol{v}_{C}\rangle$,
presented in figure 8(a), and dimensionless root-mean-squared rotation rate
$\omega^{*}=\langle\omega\rangle_{\textrm{rms}}D/V_{b}$, shown in figure 8(b),
for all rising cases at $\textit{Ga}=200$. Since particle velocity
fluctuations are dominated by their horizontal component, they are
qualitatively similar to the amplitude of the path oscillations $\hat{A}/D$. A
striking difference between the distributions of $v^{*}$ and $\omega^{*}$
concerns the region at low density ratios ($\Gamma<0.2$) and small offsets
($\mathcal{T}<0.2$), where significant velocity fluctuations, and hence path-
oscillations, are present in the almost complete absence of body rotation (a
feature that is different from the findings of Mathai et al. (2017) as
discussed in Appendix B). In figure 8(c, d), we plot $C_{d}$ vs. $v^{*}$ and
$\omega^{*}$, respectively. In all panels of figure 8, the marker edge colour
is white if $\Delta\phi<0$ and black if $\Delta\phi\geq 0$. For a subset of
the markers labelled with a white border, an approximately linear scaling of
$C_{d}$ with $v^{*2}$ can be observed in figure 8(c)). This linear range
corresponds to the region of small offsets and low density ratios featuring
path oscillations but importantly no rotations. These cases adhere to the
scaling $C_{d}(v^{*})=2.0036v^{*2}+1.1$ as indicate by the dashed black line
in figure 8(c). However, once rotation begins to become significant, we find
that it dominates the drag behaviour. This is demonstrated in figure 8(d),
from which it is clear that for the markers with black borders $C_{d}$ is
approximately proportional to $\omega^{*2}$. For cases beyond resonance,
results are reasonably well represented by the fit
$C_{d}(\omega^{*2})=0.0024\omega^{*2}+1.15$. We find similar quadratic
relationships for the higher Galileo number cases examined in this work, but
not in the case of settling particles.

###  4.5 Settling particles

Figure 9: Results for settling ($\Gamma>1$) 2D cylinders at $\textit{Ga}=200$
and $I^{*}=1$. With (a) showing the drag coefficient, (b) the Strouhal number,
(c) the Phase angle between Magnus force and particle horizontal acceleration,
(d) the path amplitude and (e) the rotational amplitude. White lines in (a)
represent isocontours of $\tilde{\mathcal{T}}$.

Up until this point, the focus was exclusively on light (rising) 2D cylinders.
For heavy particles, it was already demonstrated in Will & Krug (2021b) that
the feedback between the Magnus lift force and particle acceleration becomes
negative, effectively suppressing the resonance mode. This implies that the
strong coupling between rotation and translation and the associated drag
increase are absent, but not necessarily that the COM offset has no effect at
$\Gamma>1$. To explore this, the density ratio range from $\Gamma=1.1$ up to 5
was studied with $\mathcal{T}$ ranging from 0 to the contour $I_{G}=0$ at
$\textit{Ga}=200$, 500 and 700 and $I^{*}=1$. We present only the results for
$\textit{Ga}=200$ in figure 9(a-e) since the trends for the higher Galileo
numbers are similar.

Figure 9(a) confirms that the drag coefficient does vary as a function of
$\mathcal{T}$ also for settling particles. Yet, the magnitude of the increase
in $C_{d}$ (from around 1.2 to 1.8) is much smaller compared to that observed
for rising bodies (from around 1.1 up to 4). Furthermore, the drag increase is
more pronounced at larger $\Gamma$ and contrary to rising cylinders, the
contours of constant $C_{d}$ do not align well with isocontours of either
$\mathcal{T}$ or $\tilde{\mathcal{T}}$ (dashed white lines), suggesting that
the mechanism of drag increase here is not resonance related. This is further
evidenced by figure 9(b), where for the the phase lag $\Delta\phi$ is shown
for the same cases as in figure 7(d). Unlike for rising particles, there is no
monotonic increasing trend between offset and $\Delta\phi$ and no regime where
$\Delta\phi=0$ can be identified. In fact, the phase lag is strongly positive
(between $90^{\circ}$ and $135^{\circ}$) in the regions of elevated $C_{d}$,
which implies that $F_{m}$ (at least in part) counteracts $a_{x}$. The drag
behaviour for settling particles rather seems correlated with a reduction of
the rotational inertia around point $G$, as the latter tends to zero (black
dashed line) for increasing offset due to the fact we maintain $I^{*}=1$. In
the inset of figure 9(a) we show this explicitly by rescaling the horizontal
axis according to $\mathcal{T}/\mathcal{T}|_{I_{G}=0}$. Doing so reveals that
the drag is maximum for low, but non-zero, rotational inertia (around
$\mathcal{T}/\mathcal{T}|_{I_{G}=0}=0.9$). Similarly to rising cylinders, the
increase in drag coincides with a decrease in $Str$ (see figure 9(c)) although
this effect is much smaller here as compared to the resonance mode encountered
for $\Gamma<1$. Consistent with the trends established for light cylinders at
zero offset, the lateral (figure 9(d)) and rotational (figure 9(e)) amplitudes
are also elevated in this parameter region. This behaviour can also be
observed in supplementary video 3 showing six values of $\mathcal{T}$ for
$\Gamma=2.5$. While the magnitude of all the path-oscillations remains
significantly lower as those encountered in the resonance regime for
$\Gamma<1$, surprisingly the rotational amplitudes are on a similar level.

The relation between drag and path/rotational oscillations for settling
cylinders, shown in figure 8(c,d), is qualitatively similar to that discussed
in §4.4 for rising bodies. However, the exact scaling of $C_{d}$ with
$\omega^{*2}$ is not exactly identical. Furthermore, due to the absence of
rotational-translational coupling, the observed increase in body rotation is
not reflected in a similar increase in the path-oscillations. We suspect that
for settling the increase in drag is primarily resulting from the rotational
motion given the minute increase in the translational dynamics in this case.
Finally, we observe that the magnitudes of $C_{d}$, $\hat{A}/D$, and
$\hat{\theta}$ all decrease towards $\Gamma=1$. This behaviour is consistent
with the previous results for $\Gamma<1$ and the explanation provided in §4.1,
namely the divergence of the pendulum term.

##  5 Effects of varying moment of inertia

Figure 10: Investigation on the effect of the dimensionless moment of inertia
$I^{*}$ in combination with the timescale ratio ($\mathcal{T}$) for $Ga=200$
and $\Gamma=0.4$ on the drag coefficient (a), Strouhal number (b), phase lag
(c), translational amplitude (d), and the rotational amplitude (e). In these
figures the solid and dashed black lines, respectively, represent contours
along which $\gamma=1$ and $I_{G}=0$. In (a) the inset shows the same data as
the main panel, however is plotted against the modified timescale ratio
$\tilde{\mathcal{T}}$ to include effects of the rotational added mass due to
the Stokes layer.

In this section, we will explore the effects variations in the particle MOI
around $G$, which thus far has been kept constant. Since particle rotation
proved to be a critical aspect in the preceding analysis, it is anticipated
that the variations of the MOI will also affect the overall dynamics. We
investigate this by varying the dimensionless MOI around the geometric centre
($I^{*}$) in the range $I^{*}\in[0.5,16]$ for the cases of $\Gamma=0.4$,
$Ga=200$ and $\mathcal{T}\in[0,0.5]$. The corresponding results for $C_{d}$,
$\Delta\phi$, Str, $\hat{A}/D$, and $\hat{\theta}$ are presented in figure 10.
In these figures, physically feasible boundaries are represented by the solid
black line, which indicates the line where $\gamma=1$ and the dashed black
line indicates where $I_{G}=0$. The grey shaded region marks parameter
combinations beyond both these two criteria. This region is probed less
extensively and therefore no linear interpolation of the data is provided
there.

The results for $C_{d}$ as a function of $I^{*}$ and $\mathcal{T}$ in figure
10(a) clearly underline the need to include the fluid inertia in the analysis
of the problem. This is obvious from the fact that the $\mathcal{T}$ values at
the maximum in $C_{d}$ show significant variation as a function of $I^{*}$,
while inclusion of the fluid inertia in the definition of
$\tilde{\mathcal{T}}$ resolves this dependence. The latter can be seen from
the white dashed lines in figure 10(a), but is even more evident from the
inset, where the same data is plotted directly vs. $\tilde{\mathcal{T}}$; the
maxima in $C_{d}$ collapse onto $\tilde{\mathcal{T}}=0.11$. Besides the
dependence on $\tilde{\mathcal{T}}$, our data further show that higher values
of $I^{*}$ lead to an increased peak drag coefficient at resonance.

In figure 10(b), $\Delta\phi$ is shown for the same data set. Identically to
the results described in §4.3, we find that peak drag occurs for
$\Delta\phi=0$ when the system is in resonance. The phase data is the best way
to asses the validity of the inclusion of fluid inertia as shown in figure
10(b), the isocontours of $\tilde{\mathcal{T}}$ almost exactly match the
interpolated $\Delta\phi$ data proving the efficacy of this model. Note that
this match is obtained with no additional fitting such that this validates
also our choice for the value of $I_{a}^{*}$, that was obtained based on $Ga$
trends in §3.2.

In figure 10 c–e, we present corresponding results for the Strouhal number and
for the translational and rotational amplitudes, respectively. Consistent with
the observations for $C_{d}$ and $\Delta\phi$, isocontours of all these
quantities also line up with lines for which $\tilde{\mathcal{T}}=const$. Also
in line with the $C_{d}$ results, the resonance induced changes become
stronger with increasing $I^{*}$ in all quantities considered. This behaviour
can be understood by considering the systems as a driven and damped harmonic
oscillator noting that when the rotational inertia of the system increases,
the damping ratio decreases resulting in larger rotational amplitudes
$\hat{\theta}$. This, in turn, then affects the other parameters ($C_{d}$,
$Str$ and $\hat{A}/D$), for which the known effects of rotation become
enhanced.

As an aside, we would like to remark on the zero offset case,
$\gamma=\mathcal{T}=0$, which was studied in detail in Mathai et al. (2017).
Based on their simulations, these authors identified a transition in particle
dynamics and vortex shedding mode as a function of $I^{*}$. We were unable to
reproduce such a transition in our simulations. In order to clarify this
difference, a direct comparison of these two contradictory results is provided
in Appendix B at matching Galileo number ($Ga=500$) and overlapping range of
$I^{*}$ and $\Gamma$. Finally, we would like to point out that the conclusion
based on the present data, i.e. that rotation plays a very marginal role in
affecting regime transition in absence of a COM offset, is also consistent
with the experimental study on the effect of varying MOI for rising spheres by
Will & Krug (2021a).

##  6 Galileo number effects

Figure 11: Exploration of the combined effects of COM offset and Galileo
number for fixed density ratio $\Gamma=0.6$ and dimensionless moment of
inertia $I^{*}=1$. With (a) showing the drag coefficient $C_{d}$ and the inset
highlighting the scaling in terms of $\tilde{\mathcal{T}}$, (b) the phase lag
$\Delta\phi$, (c) the mean rotational amplitude $\hat{\theta}$, (d) the
Strouhal number Str, and (e) the mean path-amplitude $\hat{A}/D$ . The white
dashed lines indicate isocontours in $\tilde{\mathcal{T}}$, the black dashed
line indicates $\gamma=\sqrt{0.5}$.

In this section we revisit the Galileo number, previously discussed in §3.2,
however here we will take a broader scope and will look beyond the effects of
fluid inertia. In this work seven Galileo numbers, ranging from 50 up to 2000,
were examined for varying COM offset at fixed $\Gamma=0.6$, and $I^{*}=1$.
These results are presented in figure 11. Furthermore, for Ga = 500 and 700 we
also varied $\Gamma$ from 0.001 up to 5 identical to what was previously
presented for $\textit{Ga}=200$. The only difference in the results for higher
Galileo number settling particles is that for $0.16<\mathcal{T}<0.19$ a
significant horizontal drift is encountered $v_{d}/V_{b}>0.1$, resulting in
oblique trajectories. Results for the drift velocity for all cases can be
found in the supplementary data.

In figure 11 (a) the drag coefficient is shown as a function of Ga and
$\mathcal{T}$. For all cases depicted here, an increase in drag associated
with COM offset is observed. The magnitude of this increase in drag is found
to become larger with increasing Galileo number. At $\gamma=0$, a minimum drag
($C_{d}=1.2$) occurs for $\textit{Ga}=200$. The increase in $C_{d}$ towards
lower Ga is related to the viscous dominance in this regime and for higher Ga
the increase in $C_{d}$ is associated with increasing path and rotational
oscillations. Furthermore, the isocontours of $\tilde{\mathcal{T}}$, including
fluid inertia effects, capture the essential features of Galileo number
dependence of the $C_{d}$ variation reasonably well. This is highlighted in
the inset of figure 11 (a), where we rescale the horizontal axis to
$\tilde{\mathcal{T}}$. This also reveals that the onset of resonance occurs
near at constant $\tilde{\mathcal{T}}$. Additionally, the range of offsets
where $C_{d}$ is affected extends to larger $\tilde{\mathcal{T}}$ for
increasing Ga in a similar way as decreasing $\Gamma$ or $I^{*}$ would.

Figure 11 (b) shows $\Delta\phi$ for the same parameter range. Here, the
$\tilde{\mathcal{T}}=\textrm{const.}$ contours match with constant
$\Delta\phi$ predominantly over the range $-45^{\circ}<\Delta\phi<-15^{\circ}$
($0.08<\tilde{\mathcal{T}}<0.11$). This is a consequence of our choice to base
the value of $I^{*}_{a}$ and more specifically the fitting coefficient $c_{1}$
on collapsing the rising edge of $\langle\hat{\theta}\rangle$, figure 3 (c).
An alternate choice would have been to fit $c_{1}$ to align $\Delta\phi=0$
across all $Ga$. Doing so results in a value $c_{1}$ of approximately $0.5$
(note that this yields a more than 80% reduction in $I^{*}_{a}$ for
$\textit{Ga}<2000$), resulting in $\tilde{\mathcal{T}}$ contours overlapping
the $\Delta\phi=0$ isocontour in figure 11 (b). This does not alter the
conclusion that added rotational mass is responsible for the observed
behaviour. Previously discussed results collapse in a similar way for both
$c_{1}=0.5$ or $c_{1}=2.3$, but the differences highlight the fact that the
actual value depends on the particle dynamics and kinematics for a given
parameter combination of $\mathcal{T}$, $\Gamma$, $I^{*}$, and Ga and is, in
fact, not constant in time altogether.

For all Galileo numbers a reduction in the Strouhal is encountered around
$\tilde{\mathcal{T}}=0.11$, see figure 11 (d). This drop is contingent on the
presence of rotational oscillations: when rotation is absent the path-
oscillation frequency is high and when it is present it is much lower
identical to the behaviour observed in §4.1. Based on the results here
rotational amplitudes of approximately $5^{\circ}$ at high Galileo numbers are
sufficient to cause a drastic drop in Str.

In figure 11(e) the path-amplitude response is depicted. The behaviour here is
similar to that of $C_{d}$ and $\hat{\theta}$ except for the large increase in
path-amplitude at high Ga and small values of $\mathcal{T}$. We can see that
as Ga increases, the path amplitude at $\gamma=0$ also grows for
$\textit{Ga}\geq 500$, whereas the rotational amplitude remains low in the
same parameter region (see figure 11 (c). This behaviour is akin to that
observed for spheres as reported in Auguste & Magnaudet (2018); Will & Krug
(2021a), where a $\Gamma$ threshold is encountered demarcating the transition
between a vertical and chaotic rise mode, notably in the absence of strong
rotation. The $\Gamma$ value of this threshold was shown to increase with Ga
which is what we find here as well since for $\textit{Ga}=200$ this transition
was encountered for $0.1<\Gamma<0.2$, see §4.1. And indeed, for
$\textit{Ga}\geq 500$ the behaviour of the 2D cylinders is chaotic with large
fluctuations in both path and rotational amplitudes without period-to-period
regularity.

Figure 12: (a) Amplitude of the path oscillations for $\textit{Ga}=700$ and
$I^{*}=1$ as a function $\mathcal{T}$ and $\Gamma$. (b) For the same parameter
space shows the standard deviation of the peak path-amplitude. (d–e)
Trajectories and wake structure for $Ga$ 50, $\Gamma=0.6$, $I^{*}=1$ and
$\mathcal{T}=0$, 0.3, and 0.6, respectively. The colour gradient in the wake
indicates non-dimensional fluid vorticity ($\omega_{f}D/V_{b}$) as indicated
by the colour bar.

To highlight this period-to-period irregularity and to demonstrate the effect
of COM offset on these dynamics, we show the path-amplitude ($\hat{A}/D$) for
$\textit{Ga}=700$ in figure 12 (a) as well as the standard deviation of
$\hat{A}^{\prime}/D$ of this quantity in figure 12 (b). For these higher
Galileo number cases the magnitude of the fluctuations in the amplitude is
comparable to the path-amplitude itself. The irregularity is also much higher
than at $\textit{Ga}=200$, for which $\hat{A}^{\prime}/D\leq 0.1$ even for the
transitional cases discussed in §4.2. This stresses that behaviour at higher
Ga is in fact chaotic and it is therefore quite remarkable that the mean
quantities are reasonably well behaved and in-line with results for lower
Galileo numbers. This can in part be explained by a second observation
pertaining to figure 12 (b), namely that when the resonance threshold is
exceeded, i.e. $\tilde{\mathcal{T}}>0.11$, the value of $\hat{A}^{\prime}/D$
drops drastically. The periodicity imposed by the pendulum frequency becomes
dominant and stabilises the chaotic motion resulting from the body-fluid
interaction. Qualitatively, this chaotic behaviour at low offset and the
stabilising effect of large offsets is visible in the supplementary video 4
and supplementary video 5 showing results for $\textit{Ga}=500$ and $700$,
respectively.

Finally, we focus on the lowest Galileo number ($\textit{Ga}=50$), where the
zero offset case exhibited a vertical (non-oblique) rise mode with
superimposed path oscillations. For these cases, the discrete vortex shedding
that is characteristic of high $Re$ flow around a blunt body is no longer
observed. Instead, an oscillating wake is encountered as depicted for the
$\gamma=0$ case in figure 12 (c), where the particle path and wake pattern (in
terms of fluid vorticity) are visualised. Importantly, however, we find that
even the small pressure asymmetry induced by the oscillating wake and the
associated small path-oscillations ($\hat{A}/D=0.06$ at $\gamma=0$) combined
with COM offset are enough to trigger resonance behaviour similar to that
observed for the $\textit{Ga}=200$ case. In figure 12 (d), a snapshot of the
simulations for the resonance case is shown, featuring larger amplitude path
oscillations and a wider wake. Similar to the $\textit{Ga}=200$ case,
increasing the offset beyond resonance again leads to a reduction of the
amplitude of the oscillations as shown in figure 12 (e). Based on this result,
we conclude that COM offset will affect the dynamics as long as the base state
at $\gamma=0$ exhibits path-oscillations.

##  7 Summary and conclusions

In this work, we have systematically studied Centre-Of-Mass (COM) offset
effects for a freely rising or settling cylinder in a quiescent fluid via
direct numerical simulations employing the Immersed Boundary Projection
Method. The non-dimensional parameter characterising the COM offset is given
by the timescale ratio (5) $\mathcal{T}\equiv\tau_{v}/\tau_{p}$, with
$\tau_{v}$ defining the vortex shedding frequency timescale (set by the
buoyancy velocity and particle diameter), and $\tau_{p}$ a timescale
originating from the “pendulum”-like restoring torque resulting from the
offset between the centres of mass and buoyancy (cf. Will & Krug, 2021b). The
main goal of this work has been to confirm that the behaviour of COM offset
can be predicted in term of this timescale ratio, which depends on both
$\Gamma$ and $I^{*}$. Additionally, a dependence on the Galileo number is
expected but this is not reflected in the definition of $\mathcal{T}$ by Will
& Krug (2021b). These dependencies could not be adequately confirmed
experimentally in previous work due to physical constraints, therefore, a
numerical study was desirable since then one can exactly set all of the
control parameters. Simplification to 2D, i.e. cylinders, allowed us to
examine the 4-dimensional control parameter space which is not feasible for
3D. The thus studied parameter space ranges $0\leq\mathcal{T}\leq 0.6$ for COM
offset for Galileo numbers ranging from 50 up to 2000, $0.001\leq\Gamma\leq
5$, and $0.5\leq I^{*}\leq 16$.

First of all, we found that for rising particles the dynamics and response to
the offset was qualitatively similar to that of spheres; a resonance mode was
encountered at a particular offset for which the particle rotation and drag
are significantly enhanced. For increasing offset, this effect slowly reduces
towards a case where no rotation is present. This behaviour at larger
$\Gamma$, constant (or large) Ga, and constant $I^{*}$ appeared to be well
described by $\mathcal{T}$ (as was the case for Will & Krug (2021b)). However,
for lower $\Gamma$ or Ga we found that an additional effect was playing a
role. This was identified as the contribution of rotational fluid inertia
($I_{a}$). We modelled this contribution as an annulus surrounding the
cylinder, the thickness of which scales according to a boundary layer
($1/\sqrt{\textit{Ga}}$, which is identical to $1/\sqrt{\textit{Ga}}$ when
$C_{d}$ is constant). This hence introduces a Galileo (or Reynolds) number
dependence in the definition of the timescale ratio resulting in (14b), that
was previously not explored. For rising cylinders, for which the resonance
phenomenon is present, this modified timescale ratio ($\tilde{\mathcal{T}}$)
was seen to capture the trends in the observed resonance behaviour with
resonance occurring at $\tilde{\mathcal{T}}=0.11$.

Secondly, we find that body rotation is of crucial importance when considering
the dynamics and kinematics of a body moving through a fluid. When altering
COM offset, only the rotational equation of motion of the body is affected,
and indeed we note a substantial increase in rotation around
$\tilde{\mathcal{T}}=0.11$. But more importantly, this also affects the
frequency of oscillation, path-amplitude, and drag coefficient (terminal
velocity). Altering the COM a couple of percent can induce a more than three-
fold increase in $C_{d}$. This increase in drag can be attributed to an
increase in both rotational and path oscillations. However, the effect of
rotation typically is way more significant and instances exhibiting
translational oscillations without rotation featured only relatively small
drag increases. While $\tilde{\mathcal{T}}$ describes the behaviour of the
output parameters relative to the offset, the magnitude of the variation (such
as the drag increase) still depends on the full parameter combination; for
instance the magnitude of the drag increase in resonance still depends on both
$\Gamma$ and $I^{*}$.

Thirdly, we determined that the driving of the rotational dynamics for bodies
with COM offset originates from the torque generated by horizontal path-
oscillations, i.e. the $\boldsymbol{a}_{C}\times\boldsymbol{p}$-term in (14b).
When switching only this coupling term off in the equations of motion while
leaving the pendulum term unchanged, there was almost no difference in the
particle behaviour with or without offset for both rising and settling
particles. Conversely, the (viscous) fluid torque almost entirely acts as a
damping term in the presence of an offset limiting the rotational
oscillations. Generally, the magnitude of the viscous torque was found to be
too low to drive significant rotations. This applies also at zero offset where
the present data did not reproduce the regime transition for varying
rotational inertia reported in Mathai et al. (2017).

Finally, we confirmed that for $\Gamma>1$, i.e. settling cylinders, the
resonance phenomenon is no longer present. As outlined by Will & Krug (2021b),
the feedback between the rotation induced Magnus lift force and the direction
of horizontal acceleration becomes negative for heavy particles. Nevertheless,
some effects of the offset also exists for $\Gamma>1$, however, they occur at
larger offsets than expected based on the light counterparts and the effect on
$C_{d}$ is significantly smaller (and importantly does not scale according to
$\tilde{\mathcal{T}}$). Instead of the resonance mechanism the explanation for
this behaviour appears to be related to the reduction of $I_{G}$ resulting
from us enforcing $I^{*}=1$. Surprisingly, we also uncover that for both
rising and settling no effect of COM occurs when $\Gamma$ is around unity.
This can be explained by the fact that the ratio of the pendulum torque over
the driving torque ($\Gamma/|1-\Gamma|$), in (14b) goes to infinity for
$\Gamma\to 1$. Thus the driving can not overcome the restoring force of the
pendulum and rotations remain too small to engage the feedback mechanism.

To summarise, we have given a complete overview of how COM offset depends on
the four control parameters governing the system for both rising and settling
cylinders in a quiescent fluid. The dynamics and kinematics uncovered here in
terms of $\tilde{\mathcal{T}}$ qualitatively match those uncovered by Will &
Krug (2021b) for rising and settling spheres. This suggests that the present
findings largely transfer to the behaviour of spherical particles. Especially,
this concerns the behaviour at low Galileo numbers, where also for spheres
fluid inertia will become increasingly important which can be accounted for
analogously via an added mass term.

## Acknowledgements

We acknowledge PRACE for awarding us access to MareNostrum at Barcelona
Supercomputing Center (BSC), Spain (Project 2020225335). This project has
further received funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No. 950111, BU-PACT), as well as from the Netherlands Organisation
for Scientific Research (NWO) under VIDI Grant No. 13477 and through the
research programme of the Foundation for Fundamental Research on Matter with
project number 16DDS001, which is financially supported by the NWO.

Declaration of interests. The authors report no conflicts of interest.

##  Appendix A On the value of $I^{*}_{a}$

In this appendix we address how the added inertia $I_{a}^{*}$ may be
estimated, with $I_{a}^{*}$ the added inertia due to the surrounding Stokes
layer. The torque induced by this layer is estimated for a body that starts
spinning in a quiescent viscous fluid such that squares of the velocity field
may be neglected. For a sphere it may be shown that this layer yields an
effective torque that can be analytically expressed. For a cylinder the
analytical solution takes the form of an infinite series expansion (Basset,
1888). For our analysis we shall assume the torque of the cylinder to take the
same form as that of a sphere. The analysis further applies under the
condition that the relevant timescale is small enough such that effectively
only the history torque contributes (see Feuillebois & Lasek, 1978; Auguste &
Magnaudet, 2018). Here, we introduce a fitting parameter which we estimate
through a series of numerical experiments to match the history torque to that
of the cylinder.

###  A.1 History torque of a sphere

Here, we will present the main results for the history torque of a sphere,
which is the dominating contribution of the Stokes layer (see Auguste &
Magnaudet (2018)), and then find an approximate analogue history torque for a
cylinder. For small time $t$ the leading contribution of the history torque
for a sphere is (see e.g. Feuillebois & Lasek, 1978)

$T_{h}\approx-\dfrac{1}{6}\sqrt{\pi}\mu
D^{4}\nu^{-1/2}\int_{0}^{t}\dfrac{d\omega/d\tau}{\sqrt{t-\tau}}d\tau,$ (18)

with the rotational velocity $\omega$ assumed to be continuously
differentiable. Discretising the integral in (18) for small $t$ one finds

$T_{h}\approx-\dfrac{1}{3}\rho_{f}\sqrt{\pi\nu}D^{4}\Delta\omega/\sqrt{\Delta
t},$ (19)

with $\Delta\omega\equiv\omega^{n+1}-\omega^{n}$, and the difference between
time levels $n$ and $n+1$ being equal to $\Delta t$. If one then introduces

$\delta=\sqrt{\nu\Delta t/(\pi d^{2})}$ (20)

one can write for the rotational equation (cf. Auguste & Magnaudet, 2018)

$\dfrac{1}{10}\Gamma\dfrac{\Delta\omega}{\Delta
t}\sim-2\delta\dfrac{\Delta\omega}{\Delta t}.$ (21)

###  A.2 History torque of a cylinder

Our goal is to find the history torque $T_{h}$ acting on a cylinder. Here, we
assume that $T_{h}$ approximately takes a similar form as that of a sphere
described in § A.1. To this point we consider a rotating cylinder in a viscous
fluid that is at rest initially. The flow is assumed to be axisymmetric and
squares of the velocity field are neglected, yielding the expression for the
azimuthal velocity component $\hat{u}_{\theta}$ (non-dimensionalised with
length scale $D$ and velocity scale $\omega D$)

$\Rey_{\theta}\partial_{\hat{t}}\hat{u}_{\theta}=\partial_{\hat{r}}^{2}\hat{u}_{\theta}+\hat{r}^{-1}\partial_{\hat{r}}\hat{u}_{\theta}-\hat{u}_{\theta}/\hat{r}^{2},$
(22)

with $\Rey_{\theta}\equiv 0.25\omega D^{2}/\nu$. Here, we solve (22)
numerically, by using the fourth order Runge–Kutta scheme for the time-
discretisation and a second order central difference scheme for spatial
gradients. In addition, we point out that $\omega$ was given a fixed value
when solving (22). We selected multiple $\Rey_{\theta}\approx O(1)$ and
integrated up to times such that $t\ll D^{2}/\nu$ (the used time interval
ranged from $10^{-4}D^{2}/\nu$ down to $10^{-8}D^{2}\nu$). In this time-
interval we assume the history torque

$T_{h}\approx-0.25\pi D^{3}\mu c_{2}\,\dfrac{\omega}{\sqrt{\pi\nu t}},$ (23)

to be dominating the torque on the cylinder. The torque on the cylinder in
(23) has the fitting parameter $c_{2}$, which we find by calculating the
actual torque on the cylinder from our numerical experiment via
$T_{h}=0.25D^{2}\mu\omega\int_{0}^{2\pi}[\partial_{\hat{r}}\hat{u}_{\theta}-\hat{u}_{\theta}/{\hat{r}}]_{\hat{r}=D/2}\,d\theta$.
With this procedure we then find the fitting parameter to be equal to

$c_{2}\approx 1.00.$ (24)

Now that we have an approximate form of the torque on the cylinder, we can
find the analogue form of (21) for a cylinder. We then plug the obtained
expression for the torque in the angular momentum balance and find

$\Gamma\dfrac{\pi}{32}D^{4}\ddot{\theta}\sim-0.25\sqrt{\dfrac{\pi\nu}{t}}\omega
D^{3}.$ (25)

To convert the expression for continuously differentiable $\omega$ we apply
Duhamel’s principle and find (by approximating
$\Delta\omega=\omega^{n+1}-\omega^{n}$)

$\Gamma\dfrac{1}{8}\dfrac{\Delta\omega}{\Delta
t}\sim-2.00\delta\dfrac{\Delta\omega}{\Delta t}.$ (26)

The expression in (26) teaches us that the added inertia $I_{a}^{*}$ takes the
form of

$I_{a}^{*}=16\delta.$ (27)

###  A.3 Dataset fit and comparison

The moment of inertia for an annulus is obtained via

$I_{a}=\dfrac{1}{2}\rho\pi(r_{2}^{4}-r_{1}^{4}),$ (28)

with $r_{1}=0.5D$ and $r_{2}=(0.5+c_{1}/\sqrt{\textit{Ga}})D$, and $c_{1}$ a
constant. Plugging in the latter radii and calculating $I^{*}_{a}\equiv
I_{a}/I_{f}$ yields

$I^{*}_{a}=\dfrac{8c_{1}}{\sqrt{\textit{Ga}}}+\dfrac{24c_{1}^{2}}{\textit{Ga}}+\dfrac{32c_{1}^{3}}{\textit{Ga}^{3/2}}+\dfrac{16c_{1}^{4}}{\textit{Ga}^{2}}.$
(29)

We fitted $c_{1}$ such that the rotational data depicted in figure 3(a)
collapses with respect to $T^{*}$ (presented in 3b). For this fit we found
$c_{1}\approx 2.3$ to yield good results for cylinders. To compare with the
history torque value from the analysis in § A.2, we follow Mathai et al.
(2018) and set the timescale as half the oscillation time yielding $\Delta
t=0.5\mbox{{Str}}D/V_{b}$. From this it follows that
$\delta=(2\pi\mbox{{Str}}\textit{Ga})^{-1/2}$. By plugging this relationship
into (27) and assuming $\mbox{{Str}}\approx 0.11$ (the value for resonance
where rotation is dominant), we find a value $c_{1}=2.4$, which does match the
leading term of our fit ($c_{1}=2.3$) in (29) surprisingly well. In previous
work by Mathai et al. (2018) only the leading order term was taken into
account, however the inclusion of the additional higher order terms in (29)
results in a large mismatch in the actual value of $I^{*}_{a}$; more than a
factor 2 at a Galileo number of 50. However, we found that inclusion of the
higher order terms resulted in a better collapse of the data presented in the
current work, their inclusion is therefore recommended.

##  Appendix B Comparison to previous work

Figure 13: (a) Path and (b) rotational amplitude for Galileo 500 cases without
COM offset and $I^{*}=1$ compared with results from Mathai et al. (2017)
(extracted from their figures).

For the cases at $\textit{Ga}=500$ with zero offset, we examined a parameter
space spanning $\Gamma$ and $I^{*}$ that was already explored extensively in
the work by Mathai et al. (2017). Here, we take a closer look into the
differences between that work and the present one.

We compare our result for $Ga=500$ with $\gamma=0$ and $\Gamma$ ranging from
0.001 to 0.99 to results of Mathai et al. (2017) at identical parameters. A
note on the difference in convention: in the work by Mathai et al. (2017) the
parameter $I^{*}$ is equal to $I^{*}\Gamma$ in the present work. We extracted
the results from Mathai et al. (2017, figures 2 a, b), where, due to the
difference in the definitions of the nondimensional rotational inertia, our
results lie on the diagonal $m^{*}=I^{*}$, i.e. the line from (A) to (D). One
of the main findings in Mathai et al. (2017) was change in the vortex shedding
mode from a 2S mode at high $\Gamma$ and $I^{*}$ to a 2P mode at low values of
these two parameters. We found no such transition as all cases in the present
work exhibited a 2S mode. Furthermore, Mathai et al. (2017) reported that this
transition was accompanied by large increases in the path and rotational
amplitudes, $\hat{A}/D$ and $\hat{\theta}$, respectively. A direct comparison
between the results for these two parameters is presented in figure 13 (a, b).
While the general trends of a gradual increase for decreasing $\Gamma$ (i.e.
low $I^{*}\Gamma$) for both amplitudes is consistent between the works, there
are large deviations in the magnitudes of both translational and rotational
amplitudes for identical cases, especially towards lower density ratios.

These differences for identical parameters combinations raise the question if
their employed virtual mass approach (Schwarz et al., 2015) could explain
these variations. The use of the virtual mass approach (VMA) was required to
stabilise the explicit scheme in Mathai et al. (2017). We tested this
hypothesis by modifying equations (1), (2) to include a virtual mass
contribution on both sides scaled by coefficient $C_{v}$ (for which we
discretised the added time derivatives on the right hand side with a forward
Euler scheme)

$(\Gamma+C_{v})\dfrac{\textrm{d}\boldsymbol{v}_{G}}{\textrm{d}t}=\dfrac{\boldsymbol{F}_{f}}{m_{f}}+(1-\Gamma)g\boldsymbol{e}_{z}+C_{v}\dfrac{\textrm{d}\boldsymbol{v}_{G}}{\textrm{d}t},$
(30a)
$\dfrac{1}{8}(\Gamma+C_{v})D^{2}I^{*}\dfrac{\textrm{d}^{2}\theta}{\textrm{d}t^{2}}=\frac{1}{m_{f}}T_{f}+\dfrac{1}{8}C_{v}D^{2}I^{*}\dfrac{\textrm{d}^{2}\theta}{\textrm{d}t^{2}}.$
(30b)

A value of $C_{v}=0$ corresponds to the present approach, while typical values
to stabilise explicit schemes are on the order of the added mass, i.e.
$C_{v}=1$ for a cylinder (Schwarz et al., 2015). We did not observe
appreciable changes in the particle dynamics when varying $C_{v}$ in the range
$C_{v}\in[0,5]$ and it therefore appears that the discrepancies between our
work and Mathai et al. (2017) are not related to the VMA and may be caused by
other unknown factors.

## References

  * Aleksyuk & Heil (2023) Aleksyuk, A. I. & Heil, M. 2023 On the onset of long-wavelength three-dimensional instability in the cylinder wake. J. Fluid Mech. .
  * Auguste & Magnaudet (2018) Auguste, F. & Magnaudet, J. 2018 Path oscillations and enhanced drag of light rising spheres. J. Fluid Mech. 841, 228–266.
  * Balay et al. (2019) Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W. D., Karpeyev, D., Kaushik, D., Knepley, M. G., May, D. A., McInnes, L. Curfman, Mills, R. Tran, Munson, T., Rupp, K., Sanan, P., Smith, B. F., Zampini, S., Zhang, H. & Zhang, H. 2019 PETSc users manual. Tech. Rep. ANL-95/11 - Revision 3.11. Argonne National Laboratory.
  * Balay et al. (1997) Balay, S., Gropp, W. D., McInnes, L. Curfman & Smith, B. F. 1997 Efficient management of parallelism in object oriented numerical software libraries. In Modern Software Tools in Scientific Computing (ed. E. Arge, A. M. Bruaset & H. P. Langtangen), pp. 163–202. Birkhäuser Press.
  * Basset (1888) Basset, A. B. 1888 A treatise on hydrodynamics: with numerous examples, , vol. 2. Deighton, Bell and Company.
  * Bearman (1984) Bearman, P. W. 1984 Vortex shedding from oscillating bluff bodies. Annu. Rev. Fluid Mech. 16 (1), 195–222.
  * Bearman & Obasaju (1982) Bearman, P. W. & Obasaju, E. D. 1982 An experimental study of pressure fluctuations on fixed and oscillating square-section cylinders. J. Fluid Mech. 119, 297–321.
  * Bishop & Hassan (1964) Bishop, R. E. D. & Hassan, A. Y. 1964 The lift and drag forces on a circular cylinder oscillating in a flowing fluid. Proc. R. Soc. Lond. Ser. A 277 (1368), 51–75.
  * Bouchet et al. (2006) Bouchet, G., Mebarek, M. & Dus̆ek, J. 2006 Hydrodynamic forces acting on a rigid fixed sphere in early transitional regimes. Eur. J. Mech. B/Fluids 25 (3), 321–336.
  * Bourguet (2023) Bourguet, Rémi 2023 Forced rotation enhances cylinder flow-induced vibrations at subcritical reynolds number. J. Fluid Mech. 955, R3.
  * Ern et al. (2012) Ern, P., Risso, F., Fabre, D. & Magnaudet, J. 2012 Wake-induced oscillatory paths of bodies freely rising or falling in fluids. Annu. Rev. Fluid Mech. 44 (1), 97–121.
  * Feuillebois & Lasek (1978) Feuillebois, F. & Lasek, A. 1978 On the rotational historic term in non-stationary stokes flow. Q. J. Mech. Appl. Math. 31 (4), 435–443.
  * Gerrard (1966) Gerrard, J. H. 1966 The mechanics of the formation region of vortices behind bluff bodies. J. Fluid Mech. 25 (2), 401–413.
  * Govardhan & Williamson (2000) Govardhan, R. & Williamson, C. H. K. 2000 Modes of vortex formation and frequency response of a freely vibrating cylinder. J. Fluid Mech. 420, 85–130.
  * Henderson (1997) Henderson, R. D. 1997 Nonlinear dynamics and pattern formation in turbulent wake transition. J Fluid Mech. 352, 65–112.
  * Horowitz & Williamson (2006) Horowitz, M. & Williamson, C. H. K. 2006 Dynamics of a rising and falling cylinder. J. fluids Struct. 22 (6-7), 837–843.
  * Horowitz & Williamson (2010a) Horowitz, M. & Williamson, C. H. K. 2010a The effect of Reynolds number on the dynamics and wakes of freely rising and falling spheres. J. Fluid Mech. 651, 251–294.
  * Horowitz & Williamson (2010b) Horowitz, M. & Williamson, C. H. K. 2010b Vortex-induced vibration of a rising and falling cylinder. J. Fluid Mech. 662, 352–383.
  * Inoue & Sakuragi (2008) Inoue, O. & Sakuragi, A. 2008 Vortex shedding from a circular cylinder of finite length at low Reynolds numbers. Phys. Fluids 20 (3), 033601.
  * Jenny & Dušek (2004) Jenny, M. & Dušek, J. 2004 Efficient numerical method for the direct numerical simulation of the flow past a single light moving spherical body in transitional regimes. J. Comp. Phys. 194 (1), 215–232.
  * Jeon & Gharib (2004) Jeon, D. & Gharib, M. 2004 On the relationship between the vortex formation process and cylinder wake vortex patterns. J. Fluid Mech. 519, 161–181.
  * Jordan & Fromm (1972) Jordan, S. K. & Fromm, J. E. 1972 Oscillatory drag, lift, and torque on a circular cylinder in a uniform flow. Phys. Fluids 15 (3), 371–376.
  * Kempe & Fröhlich (2012) Kempe, T. & Fröhlich, J. 2012 An improved immersed boundary method with direct forcing for the simulation of particle laden flows. J. Comput. Phys. 231 (9), 3663–3684.
  * Kim & Choi (2006) Kim, D. & Choi, H. 2006 Immersed boundary method for flow around an arbitrarily moving body. J. Comput. Phys. 212 (2), 662 – 680.
  * Lǎcis et al. (2016) Lǎcis, U., Taira, K. & Bagheri, S. 2016 A stable fluid–structure-interaction solver for low-density rigid bodies using the immersed boundary projection method. J. Comp. Phys. 305, 300–318.
  * Mathai et al. (2017) Mathai, V., Zhu, X., Sun, C. & Lohse, D. 2017 Mass and moment of inertia govern the transition in the dynamics and wakes of freely rising and falling cylinders. Phys. Rev. Lett. 119 (5), 054501.
  * Mathai et al. (2018) Mathai, V., Zhu, X., Sun, C. & Lohse, D. 2018 Flutter to tumble transition of buoyant spheres triggered by rotational inertia changes. Nat. Comm. 9 (1), 1792.
  * Mougin & Magnaudet (2002) Mougin, G. & Magnaudet, J. 2002 The generalized kirchhoff equations and their application to the interaction between a rigid body and an arbitrary time-dependent viscous flow. Int. J. Multiphase Flow 28 (11), 1837–1851.
  * Namkoong et al. (2008) Namkoong, K., Yoo, J. Y. & Choi, H. G. 2008 Numerical analysis of two-dimensional motion of a freely falling circular cylinder in an infinite fluid. J. Fluid Mech. 604, 33–53.
  * Parkinson (1989) Parkinson, G. V. 1989 Phenomena and modelling of flow-induced vibrations of bluff bodies. Prog. Aerosp. Sci. 26 (2), 169–224.
  * Perry et al. (1982) Perry, A. E., Chong, M. S. & Lim, T. T. 1982 The vortex-shedding process behind two-dimensional bluff bodies. J. Fluid Mech. 116, 77–90.
  * Roma et al. (1999) Roma, A. M., Peskin, C. S. & Berger, M. J 1999 An adaptive version of the immersed boundary method. J. Comput. Phys. 153 (2), 509–534.
  * Schlichting & Gersten (2003) Schlichting, Herrmann & Gersten, Klaus 2003 Boundary-layer theory. Springer Science & Business Media.
  * Schwarz et al. (2015) Schwarz, S., Kempe, T. & Fröhlich, J. 2015 A temporal discretization scheme to compute the motion of light particles in viscous flows by an immersed boundary method. J. Comput. Phys. 281, 591–613.
  * Sutherland et al. (2023) Sutherland, B. R., DiBenedetto, M., Kaminski, A. & van den Bremer, T 2023 Fluid dynamics challenges in predicting plastic pollution transport in the ocean: A perspective. Phys. Rev. Fluids 8 (7), 070701.
  * Taira & Colonius (2007) Taira, K. & Colonius, T. 2007 The immersed boundary method: A projection approach. J. Comput. Phys. 225 (2), 2118–2137.
  * Toupoint et al. (2019) Toupoint, C., Ern, P. & Roig, V. 2019 Kinematics and wake of freely falling cylinders at moderate Reynolds numbers. J. Fluid Mech. 866, 82–111.
  * Verzicco & Orlandi (1996) Verzicco, R. & Orlandi, P. 1996 A finite-difference scheme for three-dimensional incompressible flows in cylindrical coordinates. J. Comput. Phys. 123 (2), 402–414.
  * Will & Krug (2021a) Will, J. B. & Krug, D. 2021a Dynamics of freely rising spheres: the effect of moment of inertia. J. Fluid Mech. 927, A7.
  * Will & Krug (2021b) Will, J. B. & Krug, D. 2021b Rising and sinking in resonance: Mass distribution critically affects buoyancy-driven spheres via rotational dynamics. Phys. Rev. Lett. 126, 174502.
  * Williamson (1996) Williamson, C. H. K. 1996 Vortex dynamics in the cylinder wake. Annu. Rev. Fluid Mech. 28 (1), 477–539.
  * Williamson & Brown (1998) Williamson, C. H. K. & Brown, G. L. 1998 A series in 1/$\sqrt{Re}$ to represent the strouhal–reynolds number relationship of the cylinder wake. J. Fluids Struct. 12 (8), 1073–1085.
  * Williamson & Govardhan (2004) Williamson, C. H. K. & Govardhan, R. 2004 Vortex-induced vibrations. Annu. Rev. Fluid Mech. 36, 413–455.
  * Williamson & Roshko (1988) Williamson, C. H. K. & Roshko, A. 1988 Vortex formation in the wake of an oscillating cylinder. J. Fluids Struct. 2 (4), 355–381.