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

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 time scale, 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 (MOI) and centre-of-mass (COM) offset systematically and beyond what is feasible experimentally. For certain buoyant density ratios, the particle dynamics exhibits 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 MOI. Remarkably, depending on the input parameters, resonance can be observed for COM offsets as small as a few per cent of the particle diameter, showing that the particle dynamics can be highly sensitive to this parameter.

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, Chong & Lim 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 that is generally 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 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, e.g. the shedding frequency, path amplitude, and most practically relevant, the drag. These changes to particle dynamics and kinematics are important for our understanding of, e.g. the mixing behaviour in chemical reactors and in waste/resource extraction by flotation and sedimentation (Alméras et al. 2015; Chan, Ng & Krug 2023) or in natural processes such as sedimentation in rivers or oceans (Meiburg & Kneller 2010), particle-turbulent interaction of falling snowflakes in the atmosphere (Nemes et al. 2017) or, as previously mentioned, the spread of micro plastics in our oceans (Sutherland et al. 2023).

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 (2-D) cylinders by introducing a centre-of-mass (COM) offset. This approach is motivated by the 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. Stronger rotational dynamics is also observed for anisotropic particles. However, in this case inertial (normal) forces on the particle also contribute a net torque, such that the driving in this case is not exclusively via the rotational–translational coupling introduced by the COM. The observed phenomena, for the case of spheres, could be explained via the analogy to a simple driven harmonic oscillator and expressing the results in terms of a time-scale 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 $\varGamma$, 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 2-D cylinders with COM offset, to show that the underlying physics are similar between the 2-D and three-dimensional (3-D) 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, the effects of forced rotations were examined in recent work by Bourguet (2023) for elastically mounted cylinders at subcritical Reynolds numbers (${\textit {Re}} \leq 30$), uncovering significant alterations in the flow structure and amplitude of oscillation depending on ${\textit {Re}}$ and the rotational magnitude and frequency. For the case of freely rising and settling 2-D cylinders, a critical mass density ratio ($\varGamma$), the ratio of particle-to-fluid mass density, was encountered, marking the threshold between reduced coupling at high $\varGamma$ 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 MOI was also investigated as a relevant parameter, governing the dynamics of rising and settling 2-D cylinders in the numerical work of 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 $\varGamma$ 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 2-D assumption in this work effectively corresponds to the limiting case of particle motion for very long cylinders settling or rising in a 3-D environment. Beyond a certain Reynolds number, the flow will become inherently three dimensional, even for a cylinder of infinite length; for a fixed cylinder, this is found to occur around ${\textit {Re}} \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 three dimensional, showing both horizontal (azimuthal) cylinder rotation (around the vertical axis) and fluttering motion (around a horizontal axis) (Toupoint, Ern & Roig 2019). Therefore, we would like to preface this work by stating that the principle goal is not to perfectly predict dynamics of 3-D cylinders settling or rising in a quiescent fluid but rather to better understand the physics underlying the behaviour resulting from different combinations of the four control parameters for both cylindrical and, more relevantly, spherical geometries.

1.1. Problem definition and equations of motion

In this work we concern ourselves with the dynamics and kinematics of freely rising and settling 2-D 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 the schematic in figure 1a) 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-\varGamma |gD}$.

Figure 1. (a) Schematic of the cylinder with the COM ($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 ($\kern0.08em 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 = O(1)$, the value on the left $y$ axis is indicative of body rotation.

The dynamics of freely rising and settling buoyancy-driven particles is governed by the linear and angular momentum balances. In the present work, all particle dynamics is 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

(1.1)\begin{equation} \varGamma \dfrac{\textrm{d}\boldsymbol{v}_G}{\textrm{d}t} = \dfrac{\boldsymbol{F}_f}{m_f} + (1-\varGamma)g\boldsymbol{e}_y, \end{equation}

with $\varGamma =\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.1) is independent of the COM offset. The angular momentum balance for a 2-D cylinder with COM offset reads

(1.2)\begin{equation} I_C \dfrac{\textrm{d}^2 \theta}{\textrm{d} t^2} = T_f - \dfrac{1}{2} m_p \gamma D (g\boldsymbol{e}_y +\boldsymbol{a}_C) \times \boldsymbol{p}, \end{equation}

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 $(g\boldsymbol {e}_y+\boldsymbol {a}_C) \times \boldsymbol {p}$ on the right-hand side results from the COM offset. In the current work we treat $I_C$, the MOI of the particle around point $C$, as the governing parameter characterising the body's rotational inertia. With $I_G$ following from the parallel axis theorem: $I_G = I_C - m_p \ell ^2$. Furthermore, $\gamma$ denotes the ratio $\gamma = 2\ell /D$. The rotational dynamics described in (1.2) resembles that of a forced pendulum (Will & Krug 2021b), for which, when linearised, the natural frequency is given by

(1.3)\begin{equation} f_p = \tau_p^{-1}\ = \dfrac{1}{\rm \pi}\sqrt{\dfrac{\gamma g}{D I^*}}, \end{equation}

with $I^* = I_C/I_\varGamma$ and $I_\varGamma = m_p D^2/8$ a reference value of a homogeneous cylinder with identical mass $m_p$. Therefore, the physically relevant range is $0 \leq I^* \leq 2$, i.e. $0 \leq I_C \leq 1/4 m_p D^2$ (all mass in the centre, all mass on the edge, respectively). In our work, however, we consider cases beyond $I^* = 2$ to unravel the role of this parameter. We can further write (1.2) in dimensionless form by introducing a dimensionless time scale $\tilde {t} = t/\tau _v$, where $\tau _v = D/V_b$ is the vortex shedding time scale and, therefore, represents the typical time scale of the forcing in (1.2). Since the geometry is cylindrical and only shear forces contribute to the torque around $C$, we use viscous scaling to non-dimensionalise the viscous torque term as $T_f = \mu D L V_b T^*_f$ (Jordan & Fromm 1972; Bouchet, Mebarek & Dus̆ek 2006). Here, $L$ is the length of the cylinder that is set to unity in the present work. On the other hand, since the contribution of the body acceleration is related to lateral, pressure induced, forcing, we use inertial scaling in order to non-dimensionalise this term writing the particle acceleration as $\boldsymbol {a}_C \sim \boldsymbol {F}_f/m_p$, and using inertial scaling $\|\boldsymbol {F}_f\| \sim \rho _f L D V_b^2$ to obtain $\boldsymbol {a}_C = V_b^2/(D \varGamma ) \boldsymbol {a}_C^*$. Applying this to (1.2) we obtain

(1.4)\begin{equation} \dfrac{\textrm{d}^2 \theta}{\textrm{d} \tilde{t}^2} = \dfrac{32}{\rm \pi}\dfrac{1}{Ga \varGamma I^*} T^*_f - \dfrac{4\gamma}{I^*}\left( \dfrac{1}{|1-\varGamma|} \boldsymbol{e}_y + \dfrac{1}{\varGamma} \boldsymbol{a}_C^* \right) \times \boldsymbol{p}. \end{equation}

Here, the Galileo number is defined as ${\textit {Ga}} = \sqrt {|1-\varGamma |gD^3}/\nu$. The dimensionless prefactor $\gamma /(|1-\varGamma |I^*)$ in front of the pendulum term is proportional to the square of the ratio $\mathcal {T}$ of the vortex shedding to the pendulum time scale as defined in Will & Krug (2021b). For a 2-D cylinder, this time-scale ratio is equal to

(1.5)\begin{equation} \mathcal{T} = \dfrac{\tau_v}{\tau_p} = \dfrac{1}{\rm \pi}\sqrt{\dfrac{\gamma}{|1-\varGamma| I^*}}. \end{equation}

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 first outline the numerical approach used to obtain the results as well as the data processing applied to the dataset. 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)

(2.1a)\begin{gather} \frac{\partial\boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} [\boldsymbol{u}(\boldsymbol{u}-\boldsymbol{v}_C)] =-\boldsymbol{\nabla} p + \dfrac{1}{{\textit{Ga}}} \nabla^2 \boldsymbol{u} +\boldsymbol{f}, \end{gather}

(2.1b)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}

with $\boldsymbol {u}$ the fluid velocity vector, $p$ the kinematic pressure and $\boldsymbol {f}$ the boundary forcing from the immersed boundary method 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 2-D 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 transition from the fine constant spacing surrounding the cylinder to the linearly expanding mesh, which is required to avoid numerical artefacts, 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.

Table 1. Overview of the grids. The first column denotes the Galileo number ${\textit {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 (IBPM) originally developed by Taira & Colonius (2007) and Lcis, Taira & Bagheri (2016) of which a brief overview is presented. We solve (2.1a) and (2.1b) on a staggered grid, where the spatial gradients are computed using a conservative second-order central finite difference scheme. The nonlinear term of (2.1a) is advanced in time via the explicit second-order Adams–Bashforth scheme and the viscous terms via the second-order implicit Crank–Nicolson scheme, i.e.

(2.2a)$$\begin{gather} \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{\boldsymbol{\mathsf{G}}}\boldsymbol{\varphi}^{n+1/2} \nonumber\\ -\hat{\boldsymbol{\mathsf{G}}}\boldsymbol{p}^n+ \frac{1}{2{\textit{Ga}}}\hat{\boldsymbol{\mathsf{L}}} (\boldsymbol{u}^{n+1} +\boldsymbol{u}^n) + \hat{\boldsymbol{\mathsf{H}}} \boldsymbol{f}^{n+1/2}, \end{gather}$$

where

(2.2b)\begin{equation} \hat{\boldsymbol{\mathsf{D}}}\boldsymbol{u}^{n+1}=0 \quad \text{and} \quad \hat{\boldsymbol{\mathsf{E}}} \boldsymbol{u}^{n+1}=\boldsymbol{v}_C^{n+1}+\boldsymbol{\omega}^{n+1}\times \boldsymbol{\mathcal{L}}. \end{equation}

Here, $\hat {N}(\boldsymbol {u},\boldsymbol {v}_C)$ denotes the nonlinear operator, $\hat {\boldsymbol{\mathsf{G}}}$ the gradient operator, $\hat {\boldsymbol{\mathsf{L}}}$ the Laplace operator, $\hat {\boldsymbol{\mathsf{D}}}$ the divergence operator, $\hat {\boldsymbol{\mathsf{H}}}$ the spreading (regularisation) operator, $\hat {\boldsymbol{\mathsf{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 (2.1a). The interpolation and regularisation matrices $\hat {\boldsymbol{\mathsf{E}}}$ and $\hat {\boldsymbol{\mathsf{H}}}$ make use of a discrete three-point $\delta$ function (Roma, Peskin & Berger 1999).

For the particle, the non-dimensional Newton–Euler equations are advanced in time via

(2.3)\begin{equation} \dfrac{1}{\Delta t} (\boldsymbol{v}_C^{n+1}-\boldsymbol{v}_C^n) =\boldsymbol{\mathsf{C}}_B(\boldsymbol{\mathsf{N}}_B\tilde{\boldsymbol{f}}^{n+1/2} + \Delta\boldsymbol{q}_B) + \boldsymbol{g}^n - \boldsymbol{\zeta}^n, \end{equation}

with

(2.4a,b)\begin{align} \boldsymbol{\mathsf{C}}_B = \dfrac{4}{{\rm \pi}\varGamma} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\0 & 0 & 8/I^*\end{bmatrix},\quad \boldsymbol{\mathsf{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}, \end{align}

$\Delta \boldsymbol {q}_B\equiv \boldsymbol{\mathsf{Q}} (\boldsymbol {u}^n -\boldsymbol {u}^{n-1})/\Delta t$ and $\boldsymbol{\mathsf{Q}}$ the matrix that interpolates the velocity inside the cylinder (cf. Kempe & Fröhlich 2012). The time step is limited with a Courant–Friedrichs–Lewy number of 0.4.

Vector $\boldsymbol {g}^n$ contains the buoyancy force and the torque induced by the particle weight. The Newton equation in (1.1) is solved with respect to the geometric centre and we make use of $\boldsymbol {v}_C= \boldsymbol {v}_G + \gamma /2\,\boldsymbol {\omega }\times \boldsymbol {p}$ for the transformation of (1.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$,

(2.5a,b)\begin{equation} \boldsymbol{g}^n \equiv \begin{Bmatrix} 0\\ \textrm{sgn}(1-\varGamma)/\varGamma\\ - 4\gamma/I^* \sin\theta^n \end{Bmatrix},\quad \boldsymbol{\zeta}^n \equiv \dfrac{\gamma }{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) / I^* \end{Bmatrix}, \end{equation}

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 (2.2a) together with constraints (2.2b) and (2.3) can be rewritten as

(2.6) \begin{equation} \begin{bmatrix} \boldsymbol{\mathsf{A}} & \boldsymbol{\mathsf{0}} & \boldsymbol{\mathsf{G}} & \boldsymbol{\mathsf{E}}^T \\ \boldsymbol{\mathsf{0}} & \boldsymbol{\mathsf{I}}_B & \boldsymbol{\mathsf{0}} & \boldsymbol{\mathsf{N}}_B \\ \boldsymbol{\mathsf{G}}^T & \boldsymbol{\mathsf{0}} & \boldsymbol{\mathsf{0}} & \boldsymbol{\mathsf{0}} \\ \boldsymbol{\mathsf{E}} & \boldsymbol{\mathsf{N}}_B^T & \boldsymbol{\mathsf{0}} & \boldsymbol{\mathsf{0}} \end{bmatrix} \begin{Bmatrix} \boldsymbol{q}^{n+1} \\ \boldsymbol{v}_C^{n+1} \\ \boldsymbol{\varphi}^{n+1/2} \\ \tilde{\boldsymbol{f}}^{n+1/2} \end{Bmatrix} = \begin{Bmatrix} \boldsymbol{r}^n \\ \boldsymbol{r}_B^n \\\boldsymbol{0}\\ \boldsymbol{0} \end{Bmatrix}, \end{equation}

with $\boldsymbol {q}^{n+1}=\boldsymbol{\mathsf{R}}\boldsymbol {u}^{n+1}$, $\boldsymbol{\mathsf{E}} = \hat {\boldsymbol{\mathsf{E}}}\boldsymbol{\mathsf{R}}^{-1}$, diagonal matrix $\boldsymbol{\mathsf{R}}\equiv [\Delta y_j, 0; 0, \Delta x_i]$, $\boldsymbol{\mathsf{A}}=\hat {\boldsymbol{\mathsf{M}}}[\boldsymbol{\mathsf{I}}/\Delta t -\hat {\boldsymbol{\mathsf{L}}}/(2{\textit {Ga}})]$, and $\boldsymbol {r}^n$, $\boldsymbol {r}_B^n$ containing the explicit terms. We approximate the inverse of $\boldsymbol{\mathsf{A}}$ as $\boldsymbol{\mathsf{A}}^{-1}\approx \boldsymbol{\mathsf{B}} = \Delta t \boldsymbol{\mathsf{M}}^{-1}$ (Lcis et al. 2016, cf. § 3), with $\boldsymbol{\mathsf{M}}=\hat {\boldsymbol{\mathsf{M}}}\boldsymbol{\mathsf{R}}^{-1}$, $\hat {\boldsymbol{\mathsf{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 (2.6) is performed via a block-LU decomposition following a three step procedure

(2.7a)\begin{gather} \begin{bmatrix} \boldsymbol{\mathsf{A}} & \boldsymbol{\mathsf{0}}\\ \boldsymbol{\mathsf{0}} & \boldsymbol{\mathsf{I}}_B \end{bmatrix} \begin{Bmatrix} \boldsymbol{q}^* \\ \boldsymbol{v}_C^* \end{Bmatrix} = \begin{Bmatrix} \boldsymbol{r}^n \\ \boldsymbol{r}_B^n \end{Bmatrix}, \end{gather}

(2.7b)\begin{gather} \begin{bmatrix} \boldsymbol{\mathsf{G}}^T\boldsymbol{\mathsf{BG}} & \boldsymbol{\mathsf{G}}^T\boldsymbol{\mathsf{B}} \boldsymbol{\mathsf{E}}^T \\ \boldsymbol{\mathsf{EBG}} & \boldsymbol{\mathsf{EBE}}^T + \boldsymbol{\mathsf{N}}_B^T \boldsymbol{\mathsf{I}}_B^{-1}\boldsymbol{\mathsf{N}}_B \end{bmatrix} \begin{Bmatrix} \boldsymbol{\varphi}^{n+1/2} \\ \tilde{\boldsymbol{f}}^{n+1/2} \end{Bmatrix} = \begin{Bmatrix}\boldsymbol{\mathsf{G}}^T\boldsymbol{q}^* \\ \boldsymbol{\mathsf{E}} \boldsymbol{q}^* + \boldsymbol{\mathsf{N}}_B^T\boldsymbol{v}_C^* \end{Bmatrix}, \end{gather}

(2.7c)\begin{gather} \begin{Bmatrix} \boldsymbol{q}^{n+1} \\ \boldsymbol{v}_C^{n+1} \end{Bmatrix} = \begin{Bmatrix} \boldsymbol{q}^* \\ \boldsymbol{v}_C^* \end{Bmatrix} -\begin{Bmatrix}\boldsymbol{\mathsf{BG}}\boldsymbol{\varphi}^{n+1/2}+ \boldsymbol{\mathsf{BE}}^T \tilde{\boldsymbol{f}}^{n+1/2} \\ \boldsymbol{\mathsf{I}}_B^{-1} \boldsymbol{\mathsf{N}}_B\kern0.06em \tilde{\boldsymbol{f}}^{n+1/2} \end{Bmatrix}. \end{gather}

Here, $\boldsymbol {q}^*$ in (2.7a) 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 to solve (2.7b). This combination of solvers was found to be robust and converge within 12–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) to

(2.8)\begin{equation} \boldsymbol{p}^{n+1} = \boldsymbol{p}^{n} + \boldsymbol{\varphi}^{n+1/2} - \dfrac{\Delta t}{2{\textit{Ga}}}\hat{\boldsymbol{\mathsf{L}}}\boldsymbol{\varphi}^{n+1/2}. \end{equation}

The additional solving routines for (2.7b) and (2.7c) were tested to ensure that they yield machine precision solutions satisfying the divergence free and no-slip condition (defined in (2.2b)). The overall solution procedure was found to provide a first-order convergence rate in the $L_2$ norm for the velocity field and first order in time (owing to the approximation of $\boldsymbol{\mathsf{A}}^{-1}\approx \Delta t \boldsymbol{\mathsf{M}}^{-1}$). Multiple validations for fixed and freely rising cylinders showed good agreement with previous numerical and experimental studies (see Appendix A). The method was found to be stable, even for density ratios as low as $\varGamma =0.001$. This stability is inherent to the coupling between the particle and the flow via the matrix definitions defined in (2.7c). This means that we may use explicit finite differences for the predictor velocity $\boldsymbol {v}_C^*$ in $\boldsymbol {r}^n_B$, since the corrector velocity $\boldsymbol {v}_C^{n+1}$ is solved for simultaneously with $\boldsymbol {q}^{n+1}$.

2.3. Dataset and processing

In this work, a total of 938 cases have been simulated for different combinations of the four control parameters: $\mathcal {T}$, ${\textit {Ga}}$, $\varGamma$ 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, ${\textit {Ga}}$ in the range between 50 up to 2000, $\varGamma$ between 0.001 and 5, and $I^*$ from 0.5 to 16. Compiled input and output parameters for all cases in our dataset are included in the supplementary materials available at https://doi.org/10.1017/jfm.2024.30.

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 with 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 $Ga =200$, this typically is the case after a transient time of $60 D/V_b$, which is short compared with 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.

Due to the intrinsic unsteady and non-regular motion of these bodies, 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,\varGamma,\mathcal {T},I^*)$, or one full oscillation time. Thus, we obtain a ‘centreline’ ($\boldsymbol {x}_{cl}(t)$, with horizontal mean drift velocity $v_d = \langle | \textrm {d}\kern0.06em \boldsymbol {x}_{cl, x} /\textrm {d}t |\rangle$, documented in the supplementary data) that 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^\circ$ to $180^\circ$. 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

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$, $\varGamma = 0.5$ and $I^* = 1$.

Figure 2. Snapshots of particle trajectories and wake structures of rising cylinders with Galileo number ${\textit {Ga}}=200$ and density ratio $\varGamma =0.5$ for six different COM offsets $\gamma$ (and $\mathcal {T}$) (a–f). The offset increases from left to right as indicated by the listed parameters at 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 _zD/V_b$).

The cylinder with zero COM offset in figure 2(a) is seen to rise almost straight, with regular vortex shedding occurring at double the frequency of 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, e.g. $\mathcal {T} = 0.201$ shown in figure 2(c), 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 movie 1. For this case, it is observed that the wake structure intermittently switches between several modes: (i) path oscillations with no significant shedding events (denoted by $0$ in figure 2), (ii) one single vortex pair per oscillation as in the 2S regime or (iii) four vortices per oscillation cycle; in two pairs of two shed when the body is changing direction; 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 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 figure 2(d), we observe that the large amplitude path oscillations persist albeit with a reduced wavelength. Furthermore, 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 movie 1 for $\mathcal {T} >0.3$ cases). The results shown here are representative of the $\varGamma$ range where the resonance phenomenon is present. In the following, we evaluate how this resonance behaviour depends on all of the other governing parameters.

3.2. On the importance of fluid inertia

In order to investigate the effect of fluid inertia, we first consider the mean rotational amplitude ($\hat {\theta }$) for a constant density ratio ($\varGamma = 0.6$) as a function of the time-scale 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}$.

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

When comparing results across different ${\textit {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 ${\textit {Ga}}$ decreases. Such a variation with ${\textit {Ga}}$ is not surprising since the definition of $\mathcal {T}$ does not incorporate any viscous effects. However, for low values of ${\textit {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 time scale. 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 time-scale ratio of the system become

(3.1a)\begin{equation} \tilde{f}_p = \dfrac{1}{\rm \pi}\sqrt{\dfrac{\gamma g}{D(I^* + I^*_a/\varGamma)}} \end{equation}

and

(3.1b)\begin{equation} \tilde{\mathcal{T}} = \dfrac{1}{\rm \pi}\sqrt{\dfrac{\gamma}{|1-\varGamma| \left( I^* + I^*_a/\varGamma \right) }}, \end{equation}

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\varGamma$. 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

(3.2)\begin{equation} I^*_a({\textit{Ga}}) = \dfrac{8c_1}{\sqrt{{\textit{Ga}}}} + \dfrac{24c_1^2}{{\textit{Ga}}} + \dfrac{32 c_1^3}{{\textit{Ga}}^{3/2}} + \dfrac{16 c_1^4}{{\textit{Ga}}^2}, \end{equation}

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 ${\textit {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 ${\textit {Ga}}$ are provided in figure 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 two times that of the displaced fluid. Beyond ${\textit {Ga}} = O(10^3)$, the value of $I^*_a$ changes much more slowly, explaining the weak $Ga$ dependence at higher ${\textit {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 B.

3.3. Who's driving?

When considering the right-hand side of (1.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 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 time scale $\tau _v$ and the pendulum time scale $\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 into 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 (1.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$, $\varGamma = 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 returns to that 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.

Figure 4. (a,b) Results for $Ga = 200$ and $\varGamma = 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 ($\textrm {rad}$) versus dimensionless time. During these runs at $t =0$ the $\boldsymbol {a}_C \times \boldsymbol {p}$ term for (1.2) is turned off, showing that in the absence of this coupling term the dynamics of particles with offset almost completely reverts back to that of particles without offset. (c) Phase lag $\Delta \psi$ between the rotational particle acceleration ($\alpha$) and the viscous torque ($T_f$).

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 $\varGamma$ 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^\circ$, such that the viscous torque predominantly acts as damping in these cases. In essence, these trends also hold for higher ${\textit {Ga}}$. However, the dynamics becomes somewhat more chaotic at higher ${\textit {Ga}}$, as will be shown in § 6, resulting in slightly lower values of $\Delta \psi$ on the order of $120\unicode{x2013}150^\circ$.

4. The effects of density ratio on COM offset

4.1. Frequency of oscillation

In the following sections we 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 is forced. In figure 5(a) we plot the data in the form of the Strouhal number

(4.1)\begin{equation} {\textit{Str}} = \dfrac{f\kern0.08em D}{V_b} \end{equation}

as a function of $\mathcal {T}$ and $\varGamma$. The marker colour in the figure indicates the exact $Str$ obtained from the simulations as can be read from the legend, the isocontours 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 COM, 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 = \text {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}$ and $\varGamma$. In (a) the marker colour indicates the dimensionless frequency ($Str = f\kern0.08em D/V_b$) according to the continuous variation in the top part of 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, and the colour of the regions between the isocontours corresponds to the discrete increments at the bottom of the colour bar. Dashed white lines represent isocontours of $\tilde {\mathcal {T}}$, the time-scale 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 $\varGamma$ 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$ versus the time-scale ratio $\tilde {\mathcal {T}}$. Here the marker colour indicates $\varGamma$ 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 $\varGamma$ effects in terms of this parameter. The grey shaded region in this figure indicates the frequency lock-in regime. 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 ${\textit {Str}}$ varies quite significantly from ${\textit {Str}} = 0.195$ at high density ratios ($\varGamma \geq 0.2$) to ${\textit {Str}} = 0.127$ for $\varGamma \leq 0.1$. This transition appears to be quite sudden, suggesting the existence of a critical density ratio as previously observed for rising and settling blunt bodies (Namkoong, Yoo & Choi 2008; Horowitz & Williamson 2010b; Mathai et al. 2017; Auguste & Magnaudet 2018; Will & Krug 2021a). The change in the path frequency at the lowest $\varGamma$ 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 $\varGamma$ (i.e. moving vertically in the figure). Depending on $\varGamma$, three distinct effects of increased offset can be observed. First of all, for $\varGamma \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 (1.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-\varGamma |^{-1}$ and, therefore, goes to infinity for $\varGamma \to 1$. This is not the case for the driving term that scales according to $\varGamma ^{-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 $\varGamma$ 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 switches 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 \varGamma \leq 0.8$), ${\textit {Str}}$ stays constant at approximately 0.195 (upper branch). However, as the offset increases, there is a sharp jump to the lower branch of ${\textit {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 $\varGamma = 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 $\varGamma \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 $\varGamma$, there is no effect of increasing offset, at intermediate density ratios we see a large difference between different $\mathcal {T}$, and at low $\varGamma$ 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 ${\textit {Str}}$ in this region approximately follow the lines of constant $\tilde {\mathcal {T}}$ (white dashed lines) and, in particular, the minimum of ${\textit {Str}}$ coincides roughly with $\tilde {\mathcal {T}} = 0.08$. The correlation between ${\textit {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 $\varGamma$ 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, ${\textit {Str}}$ appears to return to the trend at large $\varGamma$ where higher density ratios have a slightly higher ${\textit {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 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 $\varGamma$. The case with the lowest density of $\varGamma = 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/\varGamma$, it will likely dominate the pendulum torque, which does not diverge for small $\varGamma$.

As mentioned at the beginning of this section, the frequency of the path oscillations is important for the driving of the rotational dynamics through (1.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 time-scale ratio $\tilde {\mathcal {T}}$. Curves of constant ${\textit {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 ${\textit {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 3-D flow environment.

4.2. On the transition to resonance

In this section we discuss the transition to resonance in the intermediate $\varGamma$ regime, i.e. for $0.2 \leq \varGamma \leq 0.8$ in more detail. In the range $0.3 \leq \varGamma \leq 0.7$, the transition from the high ${\textit {Str}}$ number mode to the low ${\textit {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 $\varGamma =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 ${\textit {Str}} \approx 0.1$ and ${\textit {Str}} \approx 0.2$, respectively. The former peak also dominates for the intermediate case $\mathcal {T} = 0.195$, however, a weaker secondary peak at ${\textit {Str}} \approx 0.21$ is also seen to emerge at this offset value. Similar trends can be observed across the range $0.3 \leq \varGamma \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.

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}\ {\rm s}^{-1}$) versus time for $\varGamma = 0.8$ and $\mathcal {T} = 0.327$. We see that the dynamics exhibits a cyclical behaviour on a time scale 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 $\varGamma = 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 § 4.1 we mentioned that $\varGamma =0.2$ and $0.8$ were on the edges of the $\varGamma$ range for which a sharp transition to the resonance regime was encountered. We investigate these cases in more detail here. For $\varGamma = 0.2$, the range of $\mathcal {T}$ where multiple modes are observed widens significantly as compare to $0.3 \leq \varGamma \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 $\varGamma$ close to 0 and is similar to the cases of $\varGamma \leq 0.1$. However, looking at figure 5(c), we still observe that the cases at $\varGamma = 0.2$ and $\gamma \approx 0$ follow the upper branch in terms of $f/f_p$ and Str, making the behaviour transitional between the $\varGamma$ regimes.

At $\varGamma = 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 time scales much larger than that of the path oscillations. This behaviour stands in stark contrast to the rest of the cases that exhibit very regular periodic motion. Specifically, the parameter combinations that show this characteristic behaviour are: $\varGamma = 0.8$ with $\mathcal {T} = 0.327$ and $\mathcal {T} = 0.365$, and $\varGamma = 0.7$ with $\mathcal {T} = 0.411$. A video showing this behaviour along with the vorticity in the wake can be found in supplementary movie 2.

To quantify the frequency variations, we plot the instantaneous ${\textit {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 ${\textit {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 buildup attached wake, quickly reducing its length. This state is similar to the dynamics observed at higher $\varGamma$. 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 $\varGamma = 0.8$ and $\mathcal {T} = 0.327$ at this Galileo number and is indicative of the density regime transitions, where the dynamics exhibits signs of both regimes. These observations highlight that in the transition range, multiple states can coexist.

4.3. Drag coefficient and Magnus force

In this section we 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 2-D cylinder, this results in

(4.2)\begin{equation} C_d = \dfrac{{\rm \pi} |1-\varGamma|gD}{2\langle v_y \rangle^2}. \end{equation}

Note that in this definition $C_d$ solely reflects variations in the vertical velocity of the particle. When plotted as a function of $\varGamma$ and $\mathcal {T}$ (figure 7a), 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 $\varGamma = 0.2$. Moving to the right in the figure, i.e. towards increasing $\gamma$, the previously (§ 4.1) defined $\varGamma$ regimes can again be noted. For $\varGamma \geq 0.8$, no increase in $C_d$ is observed as a consequence of the COM offset. However, for $\varGamma \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 $\varGamma \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 $\varGamma$. 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. $\varGamma \geq 0.8$), the drag decreases slightly. It appears that the magnitude of the decrease is inversely correlated with the mass density, resulting in a larger reduction for lighter particles (figure 7c). 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).

Figure 7. (a) Particle drag coefficient as a function of the particle-to-fluid density ratio $\varGamma$ and the time-scale ratio $\mathcal {T}$ for rising particles at $Ga = 200$ and $I^* = 1$. (b) Drag coefficient plotted explicitly versus $\tilde {\mathcal {T}}$, the time-scale ratio including fluid inertia effects. Capturing the maximum drag trend reasonably well. (c) Zoomed in version of (b) showing the slight reduction of drag beyond the resonance peak. (d) Phase lag $\Delta \phi$ between the horizontal Magnus force and the horizontal component of instantaneous acceleration versus $\tilde {\mathcal {T}}$ alongside the experimental results for spheres (Will & Krug 2021b). The inset shows the correlation between $C_d$ and $\Delta \phi$.

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 dataset in the inset of figure 7(d). In the main panel of figure 7(d), the phase lag is plotted explicitly versus $\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 two dimensions is indeed comparable to that in three dimensions.

4.4. Drag correlations

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 3-D 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'^2 \rangle _{{\textit {rms}}}} /V_b$, where $\boldsymbol {v}_C' = \boldsymbol {v}_C - \langle \boldsymbol {v}_C \rangle$, presented in figure 8(a), and dimensionless root-mean-squared rotation rate $\omega ^* = \langle \omega \rangle _{{{\textit {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 ($\varGamma <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 C). 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 8c). 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.0036 v^{*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.

Figure 8. (a) Normalised velocity fluctuations $v^* = \sqrt {\langle \boldsymbol {v}_C'^2 \rangle _{{\textit {rms}}}} /V_b$ and (b) normalised rotational fluctuations $\omega ^* = \langle \omega \rangle _{{\textit {rms}}} D/V_b$ (in deg.) as a function of $\varGamma$ and $\mathcal {T}$. Correlations between the particle drag coefficient and $v^{*2}$ (c) and $\omega ^{*2}$ (d). The density ratio ($\varGamma$) in represented by the line colour. The marker edge colour, white or grey, represents $\Delta \phi < 0$ or $\Delta \phi \geq 0$, respectively.

4.5. Settling particles

Up until this point, the focus was exclusively on light (rising) 2-D 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 $\varGamma >1$. To explore this, the density ratio range from $\varGamma = 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. Results for settling ($\varGamma >1$) 2-D 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}}$.

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 with that observed for rising bodies (from around 1.1 up to 4). Furthermore, the drag increase is more pronounced at larger $\varGamma$ 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 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 9c) although this effect is much smaller here as compared with the resonance mode encountered for $\varGamma <1$. Consistent with the trends established for light cylinders at zero offset, the lateral (figure 9d) and rotational (figure 9e) amplitudes are also elevated in this parameter region. This behaviour can also be observed in supplementary movie 3 showing six values of $\mathcal {T}$ for $\varGamma = 2.5$. While the magnitude of all the path oscillations remains significantly lower than those encountered in the resonance regime for $\varGamma <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 primarily results 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 $\varGamma = 1$. This behaviour is consistent with the previous results for $\varGamma <1$ and the explanation provided in § 4.1, namely the divergence of the pendulum term.

5. Effects of varying MOI

In this section we explore the effects of 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 $\varGamma = 0.4$, $Ga=200$ and $\mathcal {T}\in [0,0.5]$. The corresponding results for $C_d$, $\Delta \phi$, ${\textit {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.

Figure 10. Investigation on the effect of the dimensionless MOI $I^*$ in combination with the time-scale ratio ($\mathcal {T}$) for $Ga = 200$ and $\varGamma = 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, it is plotted against the modified time-scale ratio $\tilde {\mathcal {T}}$ to include effects of the rotational added mass due to the Stokes layer.

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 are plotted directly versus $\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 dataset. 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 are 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}} = \text {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 C at a matching Galileo number ($Ga =500$) and overlapping range of $I^*$ and $\varGamma$. 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 the 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

In this section we revisit the Galileo number, previously discussed in § 3.2; however, here we take a broader scope and 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 $\varGamma = 0.6$ and $I^* = 1$. These results are presented in figure 11. Furthermore, for ${\textit {Ga}}$ = 500 and 700, we also varied $\varGamma$ 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.

Figure 11. Exploration of the combined effects of COM offset and Galileo number for fixed density ratio $\varGamma = 0.6$ and dimensionless MOI $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 ${\textit {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 figure 11(a) the drag coefficient is shown as a function of ${\textit {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 ${\textit {Ga}}$ is related to the viscous dominance in this regime and, for higher ${\textit {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 at near constant $\tilde {\mathcal {T}}$. Additionally, the range of offsets where $C_d$ is affected extends to larger $\tilde {\mathcal {T}}$ for increasing ${\textit {Ga}}$ in a similar way as decreasing $\varGamma$ 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$; see 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 more than 80 % reduction in $I^*_a$ for $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}$, $\varGamma$, $I^*$ and ${\textit {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 ${\textit {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 ${\textit {Ga}}$ and small values of $\mathcal {T}$. We can see that as ${\textit {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 11c). This behaviour is akin to that observed for spheres as reported in Auguste & Magnaudet (2018) and Will & Krug (2021a), where a $\varGamma$ threshold is encountered demarcating the transition between a vertical and chaotic rise mode, notably in the absence of strong rotation. The $\varGamma$ value of this threshold was shown to increase with ${\textit {Ga}}$, which is what we find here as well since, for ${\textit {Ga}} =200$, this transition was encountered for $0.1 < \varGamma < 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.

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}'/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}'/D \leq 0.1$ even for the transitional cases discussed in § 4.2. This stresses that behaviour at higher ${\textit {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}'/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 supplementary movies 4 and 5 showing results for ${\textit {Ga}} = 500$ and $700$, respectively.

Figure 12. (a) Amplitude of the path oscillations for ${\textit {Ga}} = 700$ and $I^* = 1$ as a function $\mathcal {T}$ and $\varGamma$. (b) Plot showing the standard deviation of the peak path amplitude for the same parameter space as in (a). (d,e) Trajectories and wake structure for $Ga= 50$, $\varGamma = 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.

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 COM offset effects for a freely rising or settling cylinder in a quiescent fluid via direct numerical simulations employing the IBPM. The non-dimensional parameter characterising the COM offset is given by the time-scale ratio (1.5) $\mathcal {T}\equiv \tau _v/\tau _p$, with $\tau _v$ defining the vortex shedding frequency time scale (set by the buoyancy velocity and particle diameter), and $\tau _p$ a time scale 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 terms of this time-scale ratio, which depends on both $\varGamma$ 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 two dimensions, i.e. cylinders, allowed us to examine the four-dimensional control parameter space that is not numerically feasible for three dimensions. The thus studied parameter space ranges from $0 \leq \mathcal {T} \leq 0.6$ for COM offset for Galileo numbers ranging from 50 up to 2000, $0.001 \leq \varGamma \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 $\varGamma$, constant (or large) ${\textit {Ga}}$ and constant $I^*$ appeared to be well described by $\mathcal {T}$ (as was the case for Will & Krug 2021b). However, for lower $\varGamma$ or ${\textit {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 time-scale ratio resulting in (3.1b), that was previously not explored. For rising cylinders, for which the resonance phenomenon is present, this modified time-scale 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 per cent 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 $\varGamma$ 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 (3.1b). 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 $\varGamma > 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 $\varGamma > 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 $\varGamma$ is around unity. This can be explained by the fact that the ratio of the pendulum torque over the driving torque ($\varGamma /|1-\varGamma |$), (3.1b), goes to infinity for $\varGamma \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.