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 Reference Mougin and Magnaudet2002; Jenny & Dušek Reference Jenny and Dušek2004)

(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 Reference Kim and Choi2006). 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 (Reference Taira and Colonius2007) and Lcis, Taira & Bagheri (Reference Lcis, Taira and Bagheri2016) 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 Reference Roma, Peskin and Berger1999).

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 Reference Kempe and Fröhlich2012). 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. Reference Lcis, Taira and Bagheri2016, 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 Reference Verzicco and Orlandi1996). The solution of $\boldsymbol {\varphi }^{n+1/2}$ and $\tilde {\boldsymbol {f}}^{n+1/2}$ are obtained via the PETSc library (Balay et al. Reference Balay, Gropp, Curfman McInnes and Smith1997, Reference Balay2019) 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 Reference Verzicco and Orlandi1996) 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.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 Reference Williamson and Roshko1988). 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 Reference Williamson and Brown1998; Schlichting & Gersten Reference Schlichting and Gersten2003; Mathai et al. Reference Mathai, Zhu, Sun and Lohse2018), 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 Reference Will and Krug2021b). 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.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 Reference Will and Krug2021b) 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 Reference Namkoong, Yoo and Choi2008; Horowitz & Williamson Reference Horowitz and Williamson2010b; Mathai et al. Reference Mathai, Zhu, Sun and Lohse2017; Auguste & Magnaudet Reference Auguste and Magnaudet2018; Will & Krug Reference Will and Krug2021a). 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 Reference Bishop and Hassan1964; Bearman & Obasaju Reference Bearman and Obasaju1982). 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 (Reference Will and Krug2021b) 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 Reference Will and Krug2021b).

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 Reference Will and Krug2021b). The inset shows the correlation between $C_d$ and $\Delta \phi$.

In previous work by Will & Krug (Reference Will and Krug2021b), 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 (Reference Will and Krug2021b), 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 Reference Horowitz and Williamson2010a). However, the presence of rotations was also clearly found to play a prominent role (Namkoong et al. Reference Namkoong, Yoo and Choi2008; Auguste & Magnaudet Reference Auguste and Magnaudet2018; Mathai et al. Reference Mathai, Zhu, Sun and Lohse2018). In the work on spheres with COM offset by Will & Krug (Reference Will and Krug2021b), 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 Reference Will and Krug2021a) 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. (Reference Mathai, Zhu, Sun and Lohse2017) 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 (Reference Will and Krug2021b) 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.